结构有限元法的数学基础.doc

结构有限元法的数学基础（航空学院 航空结构工程系 王生楠）一、弹性力学的基本方程在弹性力学问题里，通常是已知物体的边界、物体的弹性常数、物体所受的面力以及位移边界条件，求解其应力分量、应变分量和位移分量。为了由已知量求出未知量，必须建立这些已知量与未知量之间的关系，以及未知量之间的关系，从而导出求解的方程。设在卡氏直角坐标系中，体积为的弹性体中任意一点的坐标参数为，该点的位移参数为、应力分量为以及应变分量为。由线弹性力学理论，可以得到如下的用于描述一个弹性静力学小位移变形问题的基本方程式。（1）力的平衡方程用张量形式描述： （在内） （1-1）式中表示体力，表示应力分量对坐标分量的偏导数（以下相同）。（2）应变位移关系式（几何关系）， ， ， ， ， 用张量形式描述： （在内） （1-2）注意：应变张量与工程切应变之间的关系为。变形协调方程（从应变位移关系式中消去位移分量）： ，或 （3）应力应变关系式、本构方程（物理关系） （1-3） （1-4）式中为弹性模量系数，为劲度系数，和都具有对称性。（4）在弹性体的边界上，表面S可划分为两部分：外力已知的边界及位移为已知的边界，即，前者称为面力边界，后者称为位移边界。在面力边界上， （1-5）式中为已知边界力，为的边界外法线向量与坐标轴夹角的方向余弦。在位移边界上， （1-6）式中为已知边界位移。（1-5）式和（1-6）式统称为“边界条件”。上述的诸方程共有15个，即3个平衡方程，6个应变位移关系方程，6个物理关系方程。而未知变量也共计15个：6个应力分量，6个应变分量和3个位移分量。因此基本方程是可以求解的。这里没有考虑变形协调方程，原因是位移已经作为基本未知量。对于任意的单值连续的位移函数，如果设其有三阶的连续导数，则变形协调方程仅仅是几何方程微分的结果，自然地满足，所以位移作为基本未知量时，则不需要考虑变形协调方程。要使基本方程有确定的解，还要有对应的面力或位移边界条件。几何方程物理方程平衡微分方程位移分量ui应变分量ij应力分量ij体力 F变形协调方程位移边界条件面力边界条件混合边界条件已知位移 ui已知面力Ti二、弹性力学问题求解的一般性原理弹性力学问题就是：在给定的边界条件下，求解基本方程，求出未知量。当然，具体求解弹性力学问题时，并不需要同时求解15个基本未知量，可以而且必须做出必要的简化。根据几何方程和本构方程可见，位移、应力和应变分量之间不是相互独立的。假如已知位移分量，通过几何方程可以得到应变分量，然后通过物理方程可以得到应力分量。反之，如果已知应力分量，也可通过物理方程得到应变分量，再由几何方程的积分求出位移分量，不过这时的应变分量必须满足一组补充方程，即变形协调方程。基于上述的理由，为简化求解的难度，选取部分未知量作为基本未知量。位移解法。若以位移函数作为基本未知量求解，称为位移解法。以位移分量作为基本未知函数（未知量），从基本方程和边界条件中消去应力分量和形变分量（应变分量），导出只含位移分量的方程和相应的边界条件，并由此解出位移分量，然后再由几何方程求出形变分量，由物理方程求出应力分量。此类解法类似于结构力学中的位移法。通过几何方程将位移函数表达为应变分量，再通过物理方程将其表达为应力分量，代入平衡微分方程即可得到位移解法的基本方程。应力解法。若以应力函数作为基本未知量，称为应力解法。以应力分量作为基本未知函数（未知量），从基本方程和边界条件中消去位移分量和形变分量（应变分量），导出只含应力分量的方程和相应的边界条件，并由此解出应力分量，然后再由物理方程求出形变分量，由几何方程求出位移分量。此类解法类似于结构力学中的力法。以应力分量作为基本未知函数求解弹性力学问题时，应力分量必须满足平衡微分方程和面力边界条件。但是仅此还不够，仅仅满足平衡微分方程和面力边界条件的应力分量并不是真正的应力，因为这组应力分量求出的应变分量代入几何方程，将可能得到一组矛盾方程，不可能求出单值连续的位移分量。要使这组方程不矛盾，则要求应力分量不仅满足平衡微分方程和面力边界条件，而且要求应力分量对应的应变分量必须满足变形协调方程。这个问题也可以从物理上解释，应力分量满足平衡微分方程和面力边界条件，只能保证物体的平衡，但是不能保证物体的连续。只有这组应力分量求出的应变分量满足变形协调方程时，才能保证变形后的物体是连续的。当位移分量作为基本未知函数求解时，变形协调方程是自然满足的。只有当应力作为基本未知函数求解时，变形协调方程作为一组补充方程是必须的。因此，对于应力解法，应力分量必须满足平衡微分方程和变形协调方程。混合解法。以部分位移分量和部分应力分量作为基本未知量，称为混合解法。弹性力学问题归结为，在给定的边界条件下，求解偏微分方程组，数学上称为偏微分方程的边值问题。按照不同的边界条件，弹性力学有三类边值问题：第一类边值问题： 已知弹性体内的体力和其表面的面力，求平衡状态的弹性体内各点的应力分量和位移分量，这时的边界条件为面力边界条件。第二类边值问题：已知弹性体内的体力以及表面的位移分量，求平衡状态的弹性体内各点的应力分量和位移分量，这时的边界条件为位移边界条件。第三类边值问题：已知弹性体内的体力分量，以及物体表面的部分位移分量和部分面力分量，求平衡状态的弹性体内各点的应力分量和位移分量。这时的边界条件在面力已知的部分用面力边界条件，位移已知的部分用位移边界条件，称为混合边值问题。以上三类边值问题，代表了一些简化的实际工程问题。若不考虑物体的刚体位移，则三类边值问题的解是唯一的。三、弹性力学问题的变分解法弹性力学问题只有通过数学运算才能精确求解，在这里，可用的数学方法通常使这种求解限于过分简单的问题。为了克服实际的连续体问题的不易处理性，数学家和工程师们不时地提出了各种离散化方法，这些方法都包含着这样的一种近似：“当离散变量的数目增加时，它如所希望的那样逼近于真实的连续解”。在实现连续体问题的离散化上，数学家和工程师采用了不同的方法。数学家从连续体问题的微分方程出发，建立了可直接应用于这些方程的一般方法，如有限差分法，各种加权残值法，以及求适当定义的泛函的极值的近似方法。而工程师则是通过建立实际离散单元与连续区域的有限部分之间的模拟这种更直观的方法来处理连续体问题。“有限单元”一词的产生正是来自于这种工程的“直接模拟”的观点。变分原理就是一种最有成效的近似解法，就其本质而言，是把弹性力学的基本方程的定解问题，转换为求解泛函的极值或者驻值问题，这样就可以将基本方程由偏微分方程的边值问题转换为线性代数方程组。变分原理不仅是弹性力学近似解法的基础，而且也是数值计算方法例如有限元方法等的理论基础。3.1 等价于微分方程的积分表达形式变形体的虚功原理如果弹性体处于静力平衡状态，即平衡方程（1-1）在体积中的每个点处都成立，并且面力边界条件（1-5）在Ss上也同时满足，则对任意的一组函数ui，下面的积分表达式成立：利用分部公式，代入前面的积分表达式中，有如果我们事先要求ui满足几何方程(1-2)和位移边界条件(1-6)，即要求 （在内） （在上）则有上式表明：如果弹性体处于静力平衡状态，对于任意一组满足变形连续条件的虚位移（ui）而言，外力在虚位移上所做的虚功，等于真实应力分量在对应的虚应变上所做的虚功（即虚应变能）。这就是变形体的虚功原理。上式就是虚功原理的虚功方程。注意到在位移边界Su上，虚位移是恒等于零的，所以在上述面积分中仅需要在面力边界S上完成。反之，如果虚功方程成立，即有利用几何方程(1-2)利用分部公式，上式等号右边积分可变换为上式在推导中应用了在位移边界Su上，ui0的条件。现在将上式回代到虚功方程，整理可得因为虚位移ui是任意的，因此上式的成立，要求在弹性体内而在面力已知的边界S上,要求满足即，如果虚功方程成立，则弹性体处于静力平衡状态。显然，虚功原理等价于平衡微分方程和面力边界条件，它满足了静力平衡的要求。应该指出：虚功原理的推导并没有涉及任何材料性质，因此适用于任何材料。当然，由于使用了小变形假设，即线性的几何方程，因此虚功原理必须是在小变形条件下适用于任何材料。将虚功方程展开后，有3.2 位移变分方程-最小位能原理下面根据虚功方程，推导出弹性体的最小位能（或势能）原理。设弹性体的应变能密度函数U0是位移分量（应变分量）的函数，应用公式，则应变能密度函数的一阶变分为将上式代入虚功方程，则上式表示外力虚功等于弹性体应变能U的一阶变分（虚应变能）。定义为外力位能。注意到虚位移与真实的应力无关，因为在发生虚位移过程中外力保持不变，即变分与外力无关。而且积分和变分两种运算次序可以交换的，所以外力位能的一阶变分可以写作回代可得p称为弹性体的总位能，它等于弹性体的应变能U与外力位能V的和，是应变分量的泛函。由于应变分量通过几何方程可以用位移分量表示，所以总位能又是位移分量的泛函。以上公式表明，在所有几何可能的位移中，真实位移将使弹性体总位能p的一阶变分为零，因此真实位移使总位能取驻值。以下将证明：对于弹性体的稳定平衡状态，总位能p将取最小值。将几何可能位移和与之对应的应变代入总位能表达式，可以得到与几何可能位移对应的总位能将上式减去与真实位移分量及其应变分量对应的总位能，得将U0(ij+ij) 按泰勒级数展开，并略去二阶以上的小量，有回代可得由于总位能的一阶变分为零，因此而总位能的二阶变分为由于由于应变能密度函数U0为正定函数，即只有在所有的应变分量全部为零时其才可能为零，否则总是大于零的。因此，有所以以上证明了在所有的可能位移场中，真实位移场的总位能取最小值。这一原理称为最小位能原理。数学描述即总位能的一阶变分为零，其二阶变分是正定的（大于零）。必须强调指出的是，真实位移与其他的可能位移之间的差别在于是否满足静力平衡条件（即V内的应力平衡微分方程和给定面力边界上的平衡方程），所以说最小位能原理是用变分形式表达的平衡条件。通过总位能的一阶变分为零，可以推导出平衡微分方程和面力边界条件，这和虚功原理是相同的，即最小位能原理也等价于平衡微分方程和面力边界条件。最小位能原理完整的叙述为：在所有几何可能位移中，真实位移使总位能取最小值。最小位能原理是以位移函数作为基本未知量求解弹性力学问题的。当然，选择的位移函数必须是连续的，且在位移已知的边界上满足位移边界条件。例3.1 图示的杆承受均匀的轴向分布力作用，杆的1端为固持，2端为自由，杆长为，截面拉伸刚度为常数，为弹性模量，为断面积。试写出该杆的总位能泛函，并导得其平衡微分方程。解：设轴向位移为u(x)，满足位移边界条件。该系统的总位能为其一阶变分为由最小位能原理p=0，考虑到位移边界条件，得到平衡微分方程为和面力边界条件为（即自由端的轴力为零）例3.2 图示直梁，分布载荷q(x)作用在轴线所在的铅垂平面内，设梁的抗弯模量为EI。用最小位能原理推导问题的平衡微分方程和面力边界条件。解：设w(x)表示梁的挠度，则梁的应变能为梁的外力位能对应于横向分布载荷的位能则梁的总位能为对上式作一阶变分并且令其为零，有其中，代入上式，整理可得因此有， ， ， 上述关系式的第1式即问题的平衡微分方程，第2、3和4式为梁的边界条件，第2式代表在边界x=0和x=l处挠度等于零，第3式代表在边界x=0处挠度转角等于零，第4式代表在边界x=l处弯矩等于零。3.3 最小位能原理的应用-弹性力学问题的变分解法最小位能原理的主要用途并非推导平衡微分方程和面力边界条件，它是弹性力学问题近似解法的基础。如果要使得某个原理要应用于实际问题，必须有对应的求解方法。本节介绍基于最小位能原理的两种近似解法：瑞利-里茨（Rayleigh- Ritz）法和伽辽金（）法。根据最小位能原理，如果能够列出所有的几何可能位移，那么使总位能p取最小值的那一组位移就是真实位移。问题是列出所有几何可能的位移是非常困难的，甚至是不可能的。因此，对于实际问题的计算，只能凭借经验和直觉缩小寻找范围，在这个范围内的一族几何可能的位移中，找到一组位移使得总位能p最小。虽然这一组位移一般地说并不是真实的，但是可以肯定，它是在这个缩小的给定范围内部，与真实位移最为接近的一组位移，由此解答可以作为近似解。从上述思想出发，在一般情况下，可以将位移分量选择为如下的形式其中，Am、Bm和Cm均为任意的常数。函数u0、v0和w0以及um、vm和wm都是坐标的已知函数，并且在位移边界Su上，有， ， 和， ， 这样构造的位移试函数，不论系数Am、Bm和Cm取何值，总是满足位移边界条件的。而且对于连续函数，必然满足几何方程。因此这样构造的位移试函数满足几何可能位移的条件。现在的问题是将要如何选择待定系数Am、Bm和Cm，使得总位能p在位移表达式表示的这一族位移中取最小值。为此，将位移表达式代入几何方程求得应变分量，然后代入总位能p的表达式，注意到应变能密度函数是应变分量的齐二次函数，因此总位能p表达式的第一个积分成为待定系数Am、Bm和Cm的齐二次函数，而第二和第三个积分为Am、Bm和Cm的一次函数。于是，总位能p原本是位移函数的泛函，现在变成为待定系数Am、Bm和Cm的二次函数。这样就把求解泛函的极值问题，转化成为求解函数的极值问题。总位能p取极值的条件为，总位能p取极值的条件又可以写作上述公式是一组以Am、Bm和Cm（m=1，2，3）为未知数的线性非齐次代数方程组，求解方程可得待定系数，回代就可以得到近似位移解答。这一方法称为瑞利-里茨法。注意到，由虚功方程，可以得到如果选择的位移试函数不仅在位移边界上满足位移边界条件，而且在面力边界上满足面力边界条件，即位移试函数满足全部的边界条件，则上式可以进一步简化为上式展开可以写作将位移试函数表达式代入几何方程求得应变分量，再根据物理方程求出应力分量代入上式，并且注意到， ， 将上述结果代入虚功方程，可得由于Am, Bm和Cm 彼此独立而且是完全任意的，所以上式成立的条件为由于应力分量为Am、Bm和Cm的线性函数，所以上述公式为Am、Bm和Cm的线性非齐次代数方程组。 解出待定系数代入公式就得到位移函数的近似解答，这种方法称为伽辽金法。例3.3 对于图示右端受集中力F作用下的悬臂梁，如取挠度曲线为试用最小位能原理（瑞利-里茨法）求，的值。解：该梁的总位能的表达式为将代入，积分后得由最小位能原理，可知，：解得， 例3.4 两端简支的等截面梁，受均匀分布载荷q作用如图所示。试求解梁的挠度w(x)。解：首先使用瑞利-里茨法求解。为了满足梁的位移边界条件，即简支梁两端的约束条件：在x=0和l 处，w=0，取位移试函数，即挠度曲线方程为本问题的总位能为即根据，有并由此求得回代到位移公式，可得挠度曲线表达式是无穷级数，它给出了本问题的精确解答。这个级数收敛很快，只要取少数几项就可以得到足够的精度。最大挠度在梁的中点即x = l/2处，因此如果取一项，有。这一结果与精确值十分接近。由于上述位移试函数表示的挠度曲线方程在求二阶导数后仍为正弦函数，所以二阶导数在x=0 和x=l处仍旧为零。下面应用伽辽金法来求解本问题。本问题的静力边界条件是梁在两个铰支处的弯矩为零，所以该挠度表达式也满足面力边界条件，因此这一试函数也可以应用于伽辽金法求解。注意到将位移试函数公式代入上式并且积分，可以得到与瑞利-里茨法相同的结果。例3.5 求图示一端固支，另一端简支的梁在横向均布载荷和简支端轴向力作用挠度的近似值。解：假设挠度为 （a）因为有的作用，至少应该是的四次多项式。（a）式中的、为待求常数。下面考察位移边界条件。显然，（a）式满足处的位移边界条件：由处的位移边界条件：可得于是 （b）现在，（b）所表示的近似挠度满足全部的位移边界条件。但是，我们还可以进一步要求它满足力的边界条件 （c）则有， 以提高解的精确度。于是有 （d）并由（d）式可得梁的弯曲应变能为梁的外力位能则由轴向力的位能以及横向分布载荷的位能组成式中是梁由于弯曲而引起的轴向缩短。系统的总位能为将上面获得的各式代入上式，可得由于变分可取任意值，上式等价于即由此解得 （e）从（e）式可以看出，当时，挠度将达到无限大，即使此时。此值即为等截面杆在一端固持、另一端铰支情况下的临界压载荷的近似值，而临界压载荷的精确值为近似值的相对误差6.388。此外，（e）式给出的近似挠度函数的最大值（）为而精确值为相对误差为0.199。最后，如果取近似挠度函数为（b）式，而不是（d）式，即不同时要求满足力的边界条件，则最大挠度的相对误差为7.7，而临界压载荷的相对误差则高达41.85。例3.6 图示矩形薄板，四边固定，受有平行于板面的体力作用。设坐标轴如图所示，试用瑞利-里茨法求解。解：设位移试函数为 式中m和n为正整数，在边界x=0，a上，和y=0，b上，u=v=0，所以试函数满足位移边界条件。由于问题属于平面应力问题，所以因此将位移试函数代入上述公式，求导数后再积分，并且注意到方程则由此可见，只要体力的分布是已知的，通过积分即可以求得待定系数Amn和Bmn，从而位移分量可以求解，根据几何方程可以得到应变分量，再由物理方程求出应力分量。例3.7 图示矩形薄板，三边固定，而另外一条边的位移给定为受有平行于板面的体力作用。设坐标轴如图所示，试用伽辽金法求解。解：设位移试函数为位移试函数满足位移边界条件。由于问题没有面力边界条件，因此我们可以认为位移试函数满足面力边界条件，即可以采用伽辽金方法求解。由于问题属于平面应力问题，有将位移试函数代入上式，积分后可得积分后，求解关于Amn和Bmn的线性方程组则问题可解。如果，则问题与例3.6完全相同。本问题当然可以采用瑞利-里茨法求解。但是，一般地讲，使用伽辽金法求解相对的工作量要小一些。四、有限元法的数学原理基础知识传统变分原理的试函数采用整体插值，而在有限元方法中，是把整个弹性体物理分割为有限个分区（元素）的集合体，试函数采用的是分区插值，即试函数的选取不是整体的，而是在分区内完成的，因此试函数形式简单统一。有限元法是目前工程上应用最为广泛的结构数值分析方法，它的理论基础仍然是弹性力学的变分原理。有限元法将单元内部位移用节点位移表示，这可以使用插值函数构造单元内位移函数。并且通过单元内位移描述单元的应力和应变分量。通过最小位能原理建立单元节点位移与单元节点力的关系，即构造单元平衡方程。对于由单元集合得到的弹性体整体，应用最小位能原理构造整体平衡方程。这个方程是一个线性方程组，求解可以得到弹性体的位移，以及单元的应力和应变分量。4.1 有限元法与变分原理的关系弹性力学问题的本质是求解偏微分方程的边值问题。由于偏微分方程边值问题的复杂性，只能采取各种近似方法或者渐近方法求解。变分原理就是将弹性力学的基本方程偏微分方程的边值问题转换为代数方程求解的一种方法。有限元法的理论基础仍然是弹性力学的变分原理。那么，为什么变分原理在工程上的应用有限，而有限元原理却应用广泛。有限元原理与一般的变分原理求解方法有什么不同呢。问题在于变分原理用于弹性体分析时，不论是瑞利-里茨法还是伽辽金法，采用整体建立位移试函数或者应力试函数的方法。由于试函数要满足一定的条件，导致对于实际工程问题求解仍然困难重重。有限元法选取的试函数不是整体的，而是在弹性体内的每个分区（元素）上完成的，因此试函数形式简单统一。当然，这使得转换后的代数方程的阶次比较高。但是，面对强大的计算机处理能力，线性方程组的求解不再有任何困难。因此，有限元原理成为目前工程结构分析的重要工具。4.2 有限元法的数学基础变分原理实际是把求解偏微分方程边值问题转换为求解某一泛函的最小值问题。例如对于最小位能原理，变分方程除了满足给定的位移边界条件之外，等价于平衡微分方程和面力边界条件。当然这个转换过程的实质是在能量基础上对于边界条件的放松。因此，最小位能原理和偏微分方程边值问题仅仅是形式的不同，实质是相同的。本节将从位移变分方程引出有限元方法的基本概念。对于最小位能原理，物体的总位能为用矢量形式表达如果将物体分解为若干个有限尺寸的单元，则物体总位能为所有单元体总位能的和。有其中e为单元序号，m为单元总数。而任意一个单元体的总位能为这里Ve、(S)e分别表示第e个单元的体积和面力边界。显然如果选取的位移试函数是连续的，物体的总位能可以用单元总位能的和表示。有限元法中，物体的位移是由单元位移确定的，因此位移连续需要单元选取的位移试函数保证单元的边界位移与所有相邻单元位移相同。有限元法就是采用单元选取位移函数，而选取的位移函数在弹性体内又是连续的这一基本思想的变分方法。可见，有限元法的主要内容包括单元分析和整体分析两个部分。 （1）单元分析。将物体划分为有限个较小单元的集合体，即进行所谓的离散化。单元分析的任务是建立单元内位移模式，通过几何方程建立单元内应变位移之间的关系，通过物理方程建立单元内应力与位移之间的关系，分析单元的内力与位移之间的关系，建立单元的平衡方程，得到单元的刚度矩阵。在有限元素法中，单元的内位移模式是以元素节点上的位移自由度来表示。设元素的内位移场为 （A）其中，是元素内位移场的待定的广义位移参数，只与坐标有关。元素的内位移场在元素节点上应等于节点位移，由此可确定出广义位移参数，再代会到（A）式，可以得到元素的内位移与节点位移之间的关系式 （B）式中为联系节点位移与元素内位移的矩阵，称为元素形函数矩阵。由几何方程，元素的广义应变与节点位移之间的关系为 （C）式中：为元素广义应变列阵；为元素几何矩阵。由物理方程，元素应力列阵为 （D）式中为元素弹性矩阵。元素的总位能泛函表示为 （E）将（B）、（C）、（D）三式代入（E）式，可得 （F）其中，为作用在元素节点上的力的列阵， （G） （H）为作用在元素表面上的已知面力的节点等效节点力列阵，为联结了元素节点位移与元素节点力的矩阵，称之为元素的刚度矩阵。现在对（F）式表示的元素总位能取驻值，即，得 （I）（I）式就是元素的平衡方程，它是由元素节点位移表示的，称之为元素刚度方程。可见，最小位能原理导致“位移法”。从虚功原理，也可以得出（I）式。假设元素的节点虚位移为，相应的元素的虚内位移和虚应变为则由元素的虚功方程为考虑到的任意性，可得上式即为（I）式。 （2）整体分析。即将各单元又集合成原来的结构，要求各单元满足原结构的变形协调条件和平衡条件，从而建立整个结构的刚度方程，以求解原结构的位移和内力。设一个连续的弹性体被划分为若干(如M)个分区（元素），如右图所示。任一分区的体积为、体积力为、表面为，一般由三部分组成：其中，为中包含给定位移的边界面，为中包含给定表面力的边界面，为与相邻分区的交接面属于分区的边界。在有限元法的变分原理中，每个分区均为位能区，每个分区中独立变分的量只有一类变量（）。在每个分区的边界面上要求满足给定的位移边界条件，在分区的交接面上要求满足交接面上位移的位移相容条件。现在用和表示两个任意的相邻元素，用表示和的交接面，如图所示。另外引用两个符号和来区别交接面属于的还是属于的（这里表示的整个边界）。设每个元素的内位移表示为，； 如果选择的每一个元素的位移函数满足下列要求：（）在元素内，是连续的和单值的；（）在元素的交接面上，满足位移相容条件，即在上，（）如若元素的边界包含有，则该元素的位移函数应满足位移边界条件：（在上）。则这些位移函数的集合可以作为最小位能原理泛函的容许位移函数，那么，集合体的最小位能原理的泛函由下式给出式中经受变分的独立量是，是由（F）式确定的元素e的总位能泛函，则设整体被划分为N个节点，这N个节点的节点位移列阵为()，与之对应的节点外载荷列阵为P，则整体的总位能可以写成利用总位能原理，即，则有，或 称之为全结构平衡方程或总刚度方程，K 称之为全结构刚度矩阵，简称为总刚度矩阵，它是关于节点位移的一组代数方程组。接下来的问题，就是如何求解总刚度方程，求出节点位移。五、平面杆元素与平面桁架有限元法求解下面通过最简单的平面杆元素和平面桁架，逐步介绍有限元法的基本原理和计算过程。桁架本身就是离散结构，因此，不存在离散化问题。 对于图示平面桁架，编号为1、2、3、4、5、6的铰节点称为节点，两节点之间的链杆称为杆元素，如杆元素12、杆元素23等。在图示坐标系中，由于每一杆元素的方位不尽相同，为具普遍性，任取其中一杆元素，首先来研究杆元素的平衡关系。5.1 单元分析如图所示的杆元素ij，建立元素局部坐标系，轴沿杆元素的轴线由i节点指向j节点，杆长为Li j 。节点上的节点位移分别记为和，与节点位移相对应的节点力分别记为和，节点位移和节点力一律以顺坐标系的正方向为正。定义为杆元素在局部坐标系下的节点位移列阵和节点力列阵。问题的关键：将杆元素的其他物理量（如元素内位移、应变、应力、节点力，应变能等）用节点位移表示。（1）元素的内位移模式元素内各点的位移叫做内位移。对杆元素来说，取内位移模式为其中位移参数a、b为待定系数。内位移满足杆元素的边界（即节点）位移条件：， 则有由此求出位移参数a、b为， 代会到内位移公式中，得到用节点位移表示的内位移场式中，Ni (x)、Nj (x) 称为位移形状函数；N(x) 称为元素位移形状函数矩阵。对于杆元素，其位移形状函数具有如图所示的形状：（2）变形协调条件与几何矩阵利用变形协调条件，求元素应变，并用节点位移表示： 式中，B称为元素几何矩阵。（3）物理关系与应力矩阵利用物理关系，对于杆元素，D = E，则有式中，S称为元素应力矩阵。（4）元素平衡条件与元素刚度矩阵设杆的横截面积为A。建立杆元素的总位能为利用总位能原理，即有由此得到，将B，D代入上式，假设杆的横截面积A为常数，积分后得，或 上式就是杆元素的平衡方程，也称为元素刚度方程，它给出了元素节点力与节点位移之间的关系式。其中称为杆元素在局部坐标系中的刚度矩阵。元素刚度矩阵将元素的节点位移列阵和节点力列阵联系了起来。kii、kij、kji 和 kjj 称为刚度矩阵系数，简称为刚度系数。将上述平衡方程展开后，得到刚度系数的物理意义为 元素刚度方程也可以通过虚功原理导出。假设节点位移发生虚位移，对应的虚应变为，则有：外力虚功：虚应变能：由虚功原理：，由于虚位移是任意的，故上式成立，只有， 或记为其中： 称为元素在局部坐标系中的刚度矩阵。5.2 元素刚度矩阵的坐标变换由于结构是由许多不同元素组成的，而各个元素的局部坐标系又是不全相同的，用位移法求解节点位移时，必须规定统一的坐标系，各节点位移的矢量必须按统一的坐标系来定义，便于建立全结构的平衡方程。因此，由各个元素局部坐标系定义的元素节点位移和元素刚度矩阵必须向一个统一的坐标系转换，统一的坐标系称之为“总体坐标系”或“结构坐标系”。现考察图示平面杆元素的情况，将x、y坐标系定义为总体坐标系，总体坐标系与局部坐标系之间的夹角为（以逆时针方向为正）。将元素在局部坐标下的节点位移列阵、节点力列阵升级为：， 与之对应地，元素在局部坐标下的刚度矩阵扩阶后，变为：局部坐标系中的节点位移与总体坐标系中的节点位移，有以下转换关系：将其写成矩阵形式：记：，则有：T或 T 称为坐标变换矩阵。元素在总体坐标系下的节点位移列阵和节点力列阵，记为则对节点位移，有相应地，对节点力，有下面给出元素刚度矩阵的转换。元素在局部坐标中的作功与元素在总体坐标中的作功是相等的，据此，有：元素在总体坐标系中的刚度方程为：则元素在总体坐标系下的刚度矩阵的变换式为： 需要指出：由T矩阵中仅仅包含坐标的倾角，当坐标平移时，对刚度矩阵没有影响。因而如果仅平行移动坐标轴，刚度矩阵中元素值不变。 平面杆元素在总体坐标下的刚度矩阵的展开式为：式中，5.3 整体分析，全结构平衡方程及结构总刚度矩阵利用节点的力系平衡，建立节点平衡方程，从而得到全部节点的平衡方程。例如图示平面桁架，考虑节点2的平衡。作用在节点2上的沿总体坐标系轴正向的外载荷为X 2、Y 2，与节点2相连的杆元素的节点力如图所示，杆元素给节点2提供的力与杆元素的节点力大小相等，但方向相反，在节点2处形成平衡的共点力系。节点2上的平衡方程为或写成，利用杆元素的平衡方程，可以得到杆元素的节点力。例如杆元素1-2的平衡方程为：由此得到杆元素1-2在节点2处的节点力为：结构在节点2处的平衡方程，就可写作：一般地，对于结构中的节点 i 的平衡方程，可写作： 如果全结构有N个节点，将N个节点的平衡方程联立起来，组成全结构的平衡方程为：简记为：， 或 5.4 总刚度矩阵K 的组集总刚度矩阵K实际上是由各个元素的刚度矩阵组集而成的。元素的刚度矩阵只需按其节点位移在总刚度矩阵中的位置“按号就座”。设全结构有 N 个节点，这 N 个节点的节点位移列阵为，按节点位移列阵排列节点外载荷列阵P，即， 对于元素ij的刚度矩阵，遵照“按号就座”的原则，放置在总刚度矩阵K 中的 i 和 j所在行和列的交叉点处。对于元素im的刚度矩阵，放置在总刚度矩阵K 中的 i 和 m所在行和列的交叉点处。刚度系数在此迭加。要注意的是：刚度系数和在 i 和 i所在行和列的交叉点处迭加。5.5 总刚度方程的求解总刚度方程在没有引入位移约束条件之前，是不能求解的。这是因为，结构在没有引入位移约束条件之前是自由体，整体上存在3个(空间体6个)静力平衡条件，节点力列阵P中的项并不是完全独立的，或刚度矩阵K中的相关行和列有线性相关性，导致刚度矩阵的逆阵不存在，是奇异矩阵。在结构有限元分析中，采用引入边界条件的办法删去相关的行(列)。如果一个结构的刚体位移被限制，则与限制位移相对应的支反力自行保持结构的平衡，也即这些支反力必然是其他外力的线性组合，而它所对应的行(列)，恰恰是线性相关的行列。位移边界条件处理之后，线性相关的行(列)就删去，处理之后的总刚度矩阵K变成可逆的，此时，总刚度方程才能求解。对零位移边界条件，从总刚度矩阵中删去零位移所对应的行和列，这种方法称为删行删列法。5.6 平面桁架的有限元法求解例5.1 如图所示的长度为的杆元素1-2，设节点1处的截面积为，节点2处的截面积为，截面积沿轴线性变化，节点1和节点2处的位移和节点力分别为、和、。试求元素的刚度矩阵。解：元素刚度矩阵为例5.2 图所示的三节点杆单元，已知单元长度为，单元的横截面面积为，材料弹性模量为。单元节点位移分别为、，单元节点力分别为、，假设单元的内位移函数为，试导出：（1） 单元的形状函数矩阵N（即内位移与节点位移之间的关系矩阵）；（2） 单元的几何矩阵B（即内应变与节点位移之间的关系矩阵）；（3） 应力矩阵S；（即内应力与节点位移之间的关系矩阵）；（4） 单元的刚度矩阵K（即节点力与节点位移之间的关系矩阵）。解：（1）单元的形状函数矩阵由节点位移条件，：；：；：，即， 求得，代入位移函数表达式，得 （2）单元的几何矩阵（3）应力矩阵（4）单元的刚度矩阵由，可以导出单元的刚度矩阵的显式表达式（略）。例5.3 试用有限元法求图示桁架中3点的位移以及元素内力和支反力。设全部杆件的拉伸刚度均为EA。解：以此例说明位移法的解题步骤。1、节点编号、元素编号。2、建立总体坐标系xoy 。 3、建立元素在总体坐标系xoy 中的刚度矩阵。元素：节点编号为1-2，长度为L，=270o=-90o，元素：节点编号为1-3，长度为，=315o=-45o，元素：节点编号为2-3，长度为L，=0o，4、组集总刚度矩阵，建立总刚度方程。按编号顺序列出节点位移和节点力列阵：5、引入位移约束条件：，消除刚体位移。简单做法：从总刚度方程中，删去等式两端零位移所在的行和列。6、解引入位移约束条件后的总刚度方程，求出节点位移。7、求支反力。可以从总刚度方程中，求出支反力。8、求元素内力，即杆轴力。方法1：将已求出的节点位移代入到元素刚度方程中，可得到元素的节点力，然后合成得到杆的轴力。例如：对于元素i-j ：注意：这里需要根据节点力的正负，来判断杆子是受拉、还是受压。方法2：将总体坐标系中的节点位移，转化到元素局部坐标系中，在元素坐标系中计算杆的轴力。例如：对于元素i-j ：并且：Nij 0，杆受拉； Nij 0，杆受压。例5.4 用有限元法求图示桁架中节点1、3的位移，以及元素内力，设全部杆件的拉伸刚度均为EA。 解：1、节点编号、元素编号。2、建立总体坐标系xoy ，考虑到对称性，取一半结构来分析。3、建立元素刚度矩阵。元素2-1，元素2-3，元素1-3，4、组集总刚度矩阵，建立总刚度方程。按编号顺序列出节点位移和节点力列阵：考虑到位移约束条件：。因为要从总刚度方程中，删去等式两端零位移所在的行和列。因此，只需要针对 v1 和 v3 建立总刚度方程就可以了。5、总刚度方程，求处节点位移。6、求元素内力。六、平面梁元素与平面刚架有限元法求解6.1 一维梁元素的刚度方程及刚度矩阵如图所示的一维梁元素，同时受横向分布载荷作用，利用总位能原理求该梁元素的刚度方程及刚度矩阵。设节点位移分别为，下标1代表节点1处的量，下标2代表节点2处的量，节点位移列阵为 （6-1）与之对应的节点力列阵为 （6-2）因为节点位移有四个，我们以三次多项式表达梁的挠度， （6-3）或 （6-4）式中：，显然，式（6-3）给出的位移阶次不满足梁的平衡微分方程式。利用节点位移（6-1）式，可得 （6-5）则（6-4）式化为 （6-6）式（6-6）中的矩阵为位移插值函数，其物理涵意表述如下。下面由式（6-6）式导出几何矩阵，梁的弯曲应变为 （6-7）（6-7）式中的阵为 （6-8）将（6-6）式中的代入（6-8）式，可求出几何矩阵为 （6-9）梁的总位能泛函为 （6-10）式中为梁的抗弯模量，为梁横截面关于轴的惯性矩。由泛函p的驻值条件，即，可得 （6-11）式中 （6-