弹性力学平面问题的有限元法.doc

第二章 平面问题的有限元法 实际工程问题，在进行适当简化后，可看成平面问题。平面问题分平面应力问题和平面应变问题。 平面问题比较简单，容易理解。因此，本章以平面问题为研究对象，阐明有限元法的基本概念、理论和求解的一般步骤。第一节 两种平面问题 任何一个实际结构或构件都是空间物体，外力也都是空间力系。因此，严格地说，任何实际问题都是空间问题，应该考虑所有的应力、应变和位移分量。但是，如果所研究的结构或构件具有特殊的几何形状，承受某种特殊的外力和几何约束，就可以对其进行简化，如简化为平面问题。在平面问题中，忽略一些应力或应变分量，其余的则只是坐标x、y的函数。这样，可以使分析计算工作量大为减少。 一、平面应力问题 属于平面应力问题的弹性体，其几何形状是等厚度薄板，受到平行于板面且沿板厚t均匀分布的外力作用，前、后板面为自由表面。研究这种薄板时，坐标面xoy总是取在平分板厚的中面内，z轴垂直于板面，如图2-1a所示。 （a） (b) 图2-1由于板面是自由表面，故在z=t/2的前、后表面处，有。又由于板很薄，外力不沿厚度t变化，因此可认为板内所有各点处都有。在板内垂直于x和y轴的平面上，只剩下平行于xoy平面的三个应力分量， ，板内的点处于图2-1b所示的平面应力状态。而且，这三个应力分量沿板厚t不变，即与点的z坐标无关，只是坐标x、y的函数。平面应力问题的应力列阵简化为： （2-1） 在式（1-9）中，由可知：，但它不是独立的，取决于。因此，在平面应力问题中所需要考虑的应变分量只有，它的应变列阵简化为： （2-2） 从式（1-9）中的第一、二、四式，可解得用应变分量表达应力分量的物理方程为： （2-3）式中，平面应力问题的弹性矩阵D为： （2-4）几何方程（1-8）简化为： （2-5）在工程实际中，许多机械零件都可作为平面应力问题处理，例如，发动机连杆、直齿圆柱齿轮、平面凸轮等等。二、平面应变问题设想有无限长的截面柱体，受到平行于横截面且不沿轴线变化的外力作用，如图2-2a所示。显然，任一横截面的情况都相同，都是一个对称面，所以柱内的应力、应变、位移分量都不沿z轴变化，只是x、y的函数；而且不存在垂直于对称面的位移，w=0，每个横截面内只有沿x和y向的位移分量u、v存在，因而这种问题称为平面位移问题。 （a） （b） 图2-2根据上述位移的特点，由几何方程（1-8）可知：。因而在6个应变分量中，只剩下平行于xoy平面的3个应变分量，故称为平面应变问题。显然，平面应变问题的几何方程仍可简化为式（2-5）。在物理方程（1-9）中，由，可知 由，可知，但不是独立的。柱体内任一点的应力状态如图2-2b所示，已不是平面应力状态了。但仍然只需考虑平行于xoy平面的三个应力分量和三个应变分量。物理方程的展开式可化简为： （2-6）其中，平面应变问题的弹性矩阵为：比较式（2-4）和式（2-7）可见，只要把平面应力问题弹性矩阵中的E换成，把换成，就可得出平面应变问题的弹性矩阵。 机械工程中的长滚柱、某些轴类零件、长花键轴，以及承受均布内压或外压的长厚壁筒等，都可简化为平面应变问题来计算。实践证明，虽然它们不是无限长，但对于离开两端稍远的中间部分，按平面应变问题进行计算，所得结果与实际情况很接近，因而是可行的。 平面应力和平面应变问题虽各有特点，但有许多共同之处，即它们的应力、应变及位移分量都只是坐标x、y的函数，与z坐标无关，几何方程完全相同，物理方程的形式也相同。因此，它们的解法相同，故统称为平面问题。第二节 有限元法的概念和平面问题的离散化一、有限元法的概念在弹性力学的解析方法中，把弹性体划分成无限多个微单元体的组合体来研究。通过对其中任一微单元体的平衡、几何及物理关系的分析，建立起对弹性体内任一点都适用的一系列基本微分方程，然后，设法找出满足所有基本方程和边界条件的解析解。这些解是点的坐标的表达式，可给出弹性体内任一点处所要求的未知量。从实质上来说，这种分析方法可以称为“无限小单元法”。然而，在大多数工程实际问题中，由于结构几何形状不规则和载荷情况的复杂性，要求得完全满足边界条件的解析解，在数学上非常困难。与弹性力学解析法不同，有限元法是把连续的弹性体划分成有限多个彼此只在有限个点相联接的、有限大小的单元组合体来研究的，也就是用一个离散结构来代替原结构，作为真实结构的近似力学模型，以后的数值计算就在这个离散结构上进行。把实际连续体划分为离散结构的过程，叫做有限单元离散化。有限大小的单元，称为有限单元，常简称为单元。各单元之间相连接的点，称为节点。在平面问题中，把节点看成为光滑铰链。 目前，最广泛应用的有限元法实际上是有限元位移法。它取节点位移作为基本未知量，把原来具有无限多自由度的连续弹性体简化为有限个单元组成的离散结构，从而避免了解微分方程的麻烦。从单元分析入手，找出单元内的位移、应变、应力以及节点作用力与单元节点位移的关系，建立每个单元的刚度方程。然后，进行结构的整体分析，即组集联系整个结构的节点位移和节点载荷的总刚度方程。总刚度方程是包含有限个未知节点位移分量的线性代数方程组。最后，根据所求得的各单元的节点位移，利用单元分析得到的关系，就可求出各单元内的应力。有限元分析的基本步骤：连续体的有限单元离散化、单元分析、整体分析。二、平面问题的离散化对实际工程构件，首先要简化它的几何形状，按比例画出其平面简图，然后划分单元。图2-3a是一个受集中载荷悬臂板，其厚度不大，可作为弹性力学平面应力问题处理，采用矩形单元。在固定边上各点的位移为零，所以在该边上各节点处设置固定铰链支座。这样，就得到了悬臂板的有限元计算简图，如图2-3b所示。 （a） （b）图2-3总之，把连续体离散化成有限元计算简图时，应使原构件或结构尽可能比较准确地得到模拟。这步工作做得是否恰当，关系到最后计算精度的高低。在划分单元时，应注意下列一些原则： 1) 单元小、网格密，则计算精度高，但要求计算机的容量大，计算时间长。若单元大、网格稀，则计算精度低，但要求计算机的容量小，计算时间短。因此，应在计算机容量的范围内，根据合理的计算时间，并考虑工程上对精度的要求，合理地决定单元大小。2) 对结构的不同部位可采取不同大小的单元。对边界曲折的部位、应力或位移变化剧烈的重要部位，单元可划分得小一些(如在凹槽、孔洞等应力集中处)。对边界平直的部位、应力或位移变化平缓的次要部位，单元可划分得大些。 3) 对边界平直的构件，可用矩形单元，如图2-3，对有曲线边界的构件，宜用三角形单元，且应注意三角形单元的三个内角的大小不要相差太大，不要有钝角，最好不要小于30，因为应力和位移的误差都与单元最小内角的正弦成反比。 4) 应尽量利用结构和载荷的对称性，以减小计算工作量。 5) 在结构厚度或弹性性质的突变处，应把突变线作为单元的分界线，不要让突变线穿过任一单元。 6) 在分布载荷集度变化处和集中力作用处，应布置节点，以反映因此引起的应力变化，如图2-5c所示。 第三节 单元位移模式和解答的收敛性 图2-6设图2-6是从有限元计算简图中任取的一个单元。三个顶点取为三个节点，并规定按逆时针顺序给节点编局部号i、j、m。这三个节点的坐标分别为(xi，yi)、(xj、yj)，(xm，ym)。 平面问题中的每个节点有两个位移分量，也就是有两个自由度，可用列阵表示为：上式后的(i，j，m)表示式中下标i可以轮换为j、m，也即表示有三个列阵。把它们集合在一起即为如下的单元节点位移列阵： （2-8）右上标e表示单元。式中6个分量的正方向如图2-6所示。一、单元位移模式虽然有限元法把连续体离散成为有限个单元的集合体，但仍认为每一个单元是均匀、连续、各向同性的完全弹性体，弹性力学基本方程在各单元中都是成立的。如果知道单元内的位移分布规律，也就是知道位移函数u(x，y)、v(x，y), 那么由几何方程、物理方程就可确定单元内的应变和应力。如果只知道单元的节点位移，则不能直接求得单元内的应变和应力。有限元法克服这一困难的绝妙办法是：在一个单元的范围内，假定包含若干待定常数的近似位移函数u(x、y，)，v(x，y)，称为单元位移模式。在单元内，位移u,v应是x、y的连续函数，在各节点处，u,v的值应等于各节点的位移。只要u,v中所包含的待定常数的个数等于单元的自由度数，就可以利用单元节点位移分量确定出所有的待定常数，进而得到用节点位移表达单元内任一点位移的位移插值函数。因此，单元形状和节点数目选取不同时，所设位移模式也就应该不同。通常，把单元位移模式设为多项式形式。对于三节点三角形单元，因有6个自由度，所以，在其位移模式中应包含6个待定常数a1、a6。于是，应假设如下位移模式： （2-9）把单元的节点坐标分别代入式(2-9)，其左端的值应等于相应的节点位移，即有： (2-10),(2-11)解方程组(2-10)，可求得待定常数a1、a2、a3： 式中 (2-12),(2-13)显然，为了使式(2-13)中三角形单元面积A不为负值，单元节点局部编号i、j、m的次序必须逆时针转向。将式(2-12)中的ui、uj、um。换为vi、vj、vm，即是方程组(2-11)的解a4、a5、a6。把这样得到的常数a1、a6，代回式(2-9)中，加以整理，就得到由节点位移表达单元内任一点位移的插值公式，即位移模式的另一形式： (2-14)式中 (2-15) (2-16) 由式（2-14）可见，函数Ni（x，y）是节点i处ui=1或vi=1而另两节点位移为零时，单元内部各点处的位移u(x，y)和v（x，y）的分布形态；同理，可用来理解Nj、Nm，故称Ni、Nj、Nm为单元的形态函数，简称为形函数。通常把式（2-14）写成矩阵形式： (2-17)式中，N称为形函数矩阵： (2-18) 二、形函数的性质由式（2-13）和式（2-15），根据行列式的性质，可证明形函数有两个基本性质：1）形函数Ni在节点i处的值为1，在其余节点处之值为零，即 (2-19)注意，式（2-19）代表了9个等式。2）在单元内任一点的三个形函数之和等于1，即 由于形函数是x、y的线性函数，在单元内各点处的变化规律用几何图形直观表示时，各为一斜平面，在各边上为直线，如图2-7所示。由此图容易看出，三角形单元的形函数还有以下性质：3）在单元某一边上的形函数与第三个顶点的坐标无关。例如在ij边上 4）形函数在单元面积A上的二重积分之值，等于高为1、底为A的三角锥的体积，即：图2-7 5）形函数沿单元某一边的定积分之值，等于高为1、以该边为底的直角三角形的面积。例如，边的积分为： 三、解答的收敛性 为了使得当单元划分得越来越小时，所得解答能收敛于精确解，所设单元位移模式应满足一定的条件。 1)必须包含单元的常量应变。在一般情况下，单元中的应变包含着常量应变和变量应变。常量应变与单元中各点的位置无关，变量应变与单元中各点的位置有关。当单元的尺寸被分割得越来越小时，单元中各点的应变就会趋近于常量。可见，常量应变是单元应变的主要部分。为了正确反映单元的变形状态，所设单元位移模式必须包含单元的常量应变。把式(2-9)代入几何方程(2-5)中，得： 这就证明线性位移模式(2-9)中包含了单元的常量应变。 2)必须包含单元的刚体位移。在一般情况下，单元的位移包含着由本单元变形引起的变形位移和与本单元变形无关的刚体位移。例如，大家所熟悉的悬臂梁，在靠近自由端附近处，弯矩小、应力小、变形也很小，但位移却比弯矩大的固定端附近的位移大得多，如图2-8所示。显然，这是由于其它单元变形连带引起刚体位移的结果。由此可见，刚体位移往往是单元的主要位移。因此，所设位移模式必须包含这种位移。式(2-9)中所包含的刚体位移，可通过单元不变形时的特殊情况来研究。此时式(2-23)可知，入得单元的刚体位移分量为： 在以上两式中，、不随点的坐标变化，它们分别代表了单元沿x和y方向的刚体平动位移。当时，点的合成位移为： u合的方向如图2-9所示。它与y向的夹角为，则由此可见, 的方向与矢径R垂直,且u合的大小等于a5乘R，这表明a5，代表单元绕z轴的刚体转动角位移。 图2-9 在一般情况下，都不为零时，单元既有刚体平动位移，又有刚体转动位移。 3)要求位移在单元内部连续，在相邻单元公共边界上协调。后一要求意味着变形后相邻单元(见图2-10a)在公共边界处不会裂开或重叠(见图2-10b、c)。图2-10 由于形函数是坐标x、y的连续函数，这自然保证了位移在单元内的连续性。同时还能证明在边界上也是连续的。 在有限元法中，把满足条件1)和2)的单元称为完备单元；把满足条件3)的单元称为协调单元或保续单元。理论和实践证明，条件1)、2)是有限元解收敛于正确解的必要条件，再加上条件3)就是充分条件。线性位移模式(2-9)同时满足以上三个条件，因而三节点三角形单元是完备的协调单元，能保证当单元逐渐划小时，解答收敛于精确解。第四节 单元内的应变和应力一、单元应变与单元节点位移的关系把式(2-14)代入几何方程式(2-5)中，有： 再把式（2-15）代入上式，则单元应变可表示为： (2-24)上式中的矩阵B可写成分块形式：（2-26）式（2-24）表达了单元内应变和单元节点位移的关系，矩阵B就是用来表达的转换矩阵，故可称为应变转换矩阵。它的元素都是常量，因而单元内各点的应变分量也都是常量，所以，三节点三角形单元常称为常应变单元。二、单元应力与单元节点位移的关系把式（2-24）代入物理方程（2-3）中，即得到用单元节点位移表示的单元应力： （2-27）矩阵S是用表达单元应力的应力转换矩阵，它也可写成分块形式： （2-28）对于平面应力问题，其子矩阵为： （2-29）对于平面应变问题，只需把式（2-29）中的E换为，换为，即可得到相应的子矩阵。由式（2-29）可见，矩阵S中的各元素皆为常量，故单元中的应力也是常量。但即使各单元的S相同，由于不同，各单元的也就不同。这使相邻单元公共边界两侧的应力有突变，当然不符合实际情况。不过，随着单元的逐步取小，这种应力突变将会急剧减小。第五节 单元刚度矩阵把作用在每个单元上的截荷都静力等效地移置到各节点上之后，各单元所受到的力就只有通过节点传来的节点力。在单元节点i处，沿x和y方向有两个节点力分量，可组成列阵：每个单元共有6个节点力分量，其正方向如图2-11所示，它们组成如下单元节点力列阵： 对整个离散结构而言，单元节点力是结构内部节点和单元之间相互作用的内力。但是，对每个单元而言，节点作用于单元的节点力就是作用在单元上的外力。单元应力就是在这些节点力作用下产生的。利用虚功方程可推导出单元节点力和单元节点位移间的关系。设想单元发生了虚位移，相应的节点虚位移如图2-12所示，其6个分量可组成单元节点虚位移列阵：图2-11图2-12在单元内引起相应的虚应变为：由虚功方程（1-17）可知，单元节点力在节点虚位移上所作的虚功，等于单元内的应力在相应虚应变上所作的虚功之总和，即： （2-31）式中，A为单元面积，t为单元厚度。设单元的虚位移也具有线性位移模式，于是，根据式（2-17）和式（2-24）知，。并注意，将它们一起代入式（2-31）之后，把和提到积分号外，则式（2-31）变为：由于虚位移是任意的，因此，上式两边与相乘的矩阵应当相等，故有：可简记为： （2-32）上式就是联系单元节点力和单元节点位移的单元刚度方程，称为单元刚度矩阵。在三节点三角形单元中，矩阵B、D中的元素皆为常量；若材料是均质的，单元厚度t也为常量时，可以一起从积分号内提出，所以式（2-33）可变为： 把矩阵B的分块形式代入上式，可见由9个分块组成： （2-34）其中，每个子块是22的子矩阵： （2-35）式中，是在节点s处位移分量时，在节点r处沿x、y方向引起的节点力；是在节点s处位移分量时，在节点r处沿x、y方向引起的节点力。对平面应力问题，的展开式如下：对平面应变问题，只要把上式中的，即可得到相应的子矩阵。把式（2-35）代入式（2-34）中，可得单元刚度方程（2-32）的展开式 （2-37）单元刚度矩阵有如下一些性质：1）单元刚度矩阵中任一个元素都是一个刚度系数。2）单元刚度矩阵是对称矩阵。把式（2-36）中的下标r、s对调位置后，可见，就证明有。3）单元刚度矩阵是奇异矩阵。这可利用单元不受力作用时仍可作刚体移动来证明。此时，。由式（2-37）乘出为： 由于刚体移动时，、可以是任意值，故上式两括号内的值应等于零，于是的第一行元素之和为零。同理可证明的每行元素之和皆为零。据此，由行列式的性质可证明的行列式为零，所以，是一个奇异矩阵。而奇异矩阵的逆矩阵不存在。由此表明，若给定单元的节点位移，则可由式（2-37）确定出相应的；但是若给定，则不能由式（2-37）确定出。这是因为，单元的节点位移不仅由单元的变形引起，刚体运动也要引起节点位移。4）单元刚度矩阵只与单元的几何形状、大小、方位及材料的性质有关，而与单元的位置无关，不随单元或坐标轴的平移而改变。因为的各子块计算公式（2-36）中所包含的bi,ci等系数，只与单元节点的坐标差值有关，而与节点坐标的绝对值大小无关。因此，凡是形状、大小、相应节点坐标差值相同的单元，若材料、厚度也相同时，它们的就相同。第六节 截荷向节点的移置在有限元法中，为了便于建立节点的平衡方程，要把每个单元所受载荷按静力等效原则向节点移置，成为等效节点载荷。根据圣维南原理，移置引起的应力误差将局限于该单元的范围内，不会影响整个结构其它部分的应力。并且，对该单元应力分布的影响还会随着单元的细分而减小。静力等效可通过原载荷与移置后的节点载荷在任意虚位移上的虚功相等来实现。在一定的位移模式下，这样移置的结果是唯一的，符合通常理解的静力等效原则，即移置后的节点载荷在任一轴上投影之代数和与原载荷在同一轴上的投影相等，移置后的节点载荷对任一点的力矩之代数和与原载荷对同一点的力矩相等。 图2-13在图2-13所示单元上，M点作用着集中力；沿ij边作用着分布面力，在单位面积上的集度为分布体力在每单位体积内的集度为。把这些力移置到3个节点上，将有6个等效节点载荷分量，其正方向如图2-13所示，可集合成如下单元等效节点载荷列阵： (2-38) 设想单元发生了虚位移，单元的节点虚位移列阵为： 设单元内的虚位移仍按线性位移模式分布，则单元内的虚位移与的关系仍有： (2-39)根据移置前、后虚功相等，应有以下关系：式中，为集中力作用点M处的虚位移。把式(2-39)代入，并注意矩阵相乘的转置规则，有：由于虚位移的任意性，可得出单元等效节点载荷的计算公式： (2-40)可见，单元等效节点载荷与单元位移模式有关。若单元上只有集中力作用，则式(2-40)可简化为： (2-41)例如，在单元的ij边M点处，作用有一集中力如图2-14所示。由形函数的性质，可求得M点的形函数代入式(2-41)中，得： (2-42)若单元上只有ij边有分布面力作用，则式（2-40）可简化为： （2-43） 图2-15如果是均布面力，如图2-15所示，边长ij=l，把形函数矩阵代入式（2-43），并注意到：于是 (2-44)式中，边上总面力在x、y方向的分量。若单元上只有分布体力作用，则式（2-40）简化为：(2-45)如果为常量，注意到：由式(2-45)得：式中，为单元的总体力在x、y方向的分量。 第七节 总刚度矩阵的形成及其特点一、结构的总刚度方程设平面问题的离散结构有n个节点，其中，任一节点的位移可表示为：按节点整体编号由小到大的顺序排列，可写出结构的节点位移列阵为： （2-47）若已知，则各单元的节点位移就知道了。那么，根据前面单元分析的结果，就可求得各单元内的位移、应变和应力。只有通过对结构的整体分析，在整个离散结构中建立以节点位移为基本未知量的节点平衡方程组，求解这个方程组，才能求得各节点的位移。作用在节点上的力：一是各单元作用下节点的节点力；另一种外力，其中有各单元移置来的等效节点载荷（包括直接作用节点上的集中载荷）以及支座反力。在单元分析中，由单元刚度方程得到的是节点作用于单元的节点力，其分块形式如下：其中，任一节点i处的节点力为：上式表明，单元在节点i处的节点力与该单元所有的节点位移都有关，等于i、j、m处的节点位移在i处引起的节点力的代数和。同时，该单元也以大小相等、方向相反的节点力作用于节点i。通常，称同处于一个单元的节点为相关节点，环绕同一个节点的那些单元为相关单元。每个相关单元对所环绕的节点都作用有节点力。因此，作用在任一节点i上的节点力就等于它的所有相关单元作用于该节点的节点的代数，即： （2-48）式中，表示对环绕节点i的所有单元求和。再分析外力。作用在任一节点i上的外力，等于它的所有相关单元移置到节点i上的等效节点载荷的代数和。对受约束的节点，还应加上支反力。因此，任一节点i处的外力： （2-49）图2-16为简明起见，现以图2-16a中的简单离散结构为例。作用在节点2上的力有：相关单元、作用于节点2的节点力F2x,F2y,F2x,F2y；单元、移置到节点2上的等效节点载荷FP2x,FP2y,FP2x,FP2y。其受力图见图2-16d，因此，沿x和y方向的平衡方程是： 对任一节点i都可以建立这样两个平衡方程，可用求和符号表示成： 两式可合并写成矩阵形式： （2-50）或简写成 图2-16c、d中的1、2节点，其平衡方程为： 把式(2-48)、(2-49)代入式(2-50)，得任一节点的矩阵平衡方程； （2-51）在集合结构的所有节点平衡方程时，式(2-51)中各单元节点的局部编号i、j、m应换成对应的整体编号，并按整体编号从小到大的次序重新排列。例如，对图2-16a中单元、节点局部编号与整体编号的对应关系，可列出表2-1，、可改写成： 由式(2-51)可见，在每个节点平衡方程中，只包含与该节点相关的节点位移，而不包含与其不相关的节点位移。但是，为了写成整齐划一的矩阵方程，在每个节点平衡方程中应将不相关的节点位移也补充写出来，用22阶的零子块0乘以那些不相关的节点位移，是不会破坏各节点平衡的。对有n个节点的离散结构，则可列出n个这样的矩阵方程，依次集合在一起，即是结构接体的节点平衡方程组，可简记为： (2-52） 称为结构的总刚度方程。其中，结构的节点位移列阵是未知待求的，称为结构的节点载荷列阵，它们都是2n1阶的。K称为结构的总刚度矩阵，是2n2n阶的。由式(2-51)和上述1、2节点的平衡方程可以看出，K是由各ke扩大到2n2n阶后叠加形成的。结构整体分析的主要内容就是形成K。总刚度方程的分块形式为: 在平面问题中，每个节点有两个自由度，所以，它的K中任一子块都是22阶的。以子块为例:其中的每一个元素仍然是一个刚度系数。分别表示节点2的u2=1，其余节点位移分量全为零时，在节点1处沿x、y方向引起的节点力；分别表示v2=1，其余节点位移分量全为零时，在节点1处沿x、y方向引起的节点力。二、总刚度矩阵的组集 实际上，总刚度矩阵K是用刚度集成法或称直接刚度法组集而成的。组集办法是： 1）计算各单元刚度矩阵ke的各子块。 2）根据各单元节点的局部编号与整体编号的对应关系，把各ke子块的下标换成整体编号表示。这样，各子块的第一个下标就是该子块应放入K中的子块行号，第二个下标就是该子块应放入K中的子块列号。 3）“对号入座”地把各ke的子块送入K中相应的位置上。凡是下标相同的子块，应送入K中同一位置相叠加。把离散结构中所有K)e的子块全部“对号入座”之后，K就组集完毕了。 应注意，在K中，凡是相关节点号所对应的子块是非零子块，凡是不相关的节点号所对应的子块一定是零子块。 图2-16a中离散结构组集K为 从表2-2可见，图2-16中的二单元四节点结构的K，是由k和K分别扩大到88阶的矩阵后，叠加形成的。三、总刚度矩阵的特点1. 对称矩阵观察K矩阵的组装表2-2可见，由于各单元刚度矩阵是对称的，所以，总刚度矩阵也是对称的，即在主对角线两侧对称位置上的元素相等。利用K的对称性，可在计算机中只存放其上三角部分或下三角部分的元素，以节省计算机容量。 2. 稀疏，带状分布在一般的离散结构中，和一个节点相关的节点数最多也不会超过9个。若整个结构有200个节点，那么，在K矩阵每行中的非零子块数和该行总子块数之比不大于9200，小于5%，所以，K是一个稀疏矩阵。网格划分得越细，节点数越多，稀疏性就越突出。 另外，K中各子块的行、列号都是按节点的整体编号次序排列的，所以，其主对角线上的元素是每个节点本身发生单位位移时所引起的节点力。若节点的整体编号使相关节点之间的编号差值尽可能小，那么，非零子块就必然集中在主对角线附近，形成非零元素带状分布，通常，把K矩阵中每一行的第一个非零元素到该行主对角线上的元素的个数，叫做该行的半带宽，用d表示K的各行中的最大半带宽，如图2-17所示。在平面问题中， d(相关节点的最大编号差+1)2 (2-53)最大半带宽的大小除了与相关节点的个数有关之外，还与相关节点编号的最大差值有关。因此，在同一有限元网格中，不同的节点编号方式，会得到不同的d值。例如，在图2-18a、b中，相关节点的最大差值是5和3，因而图2-18a的d12，图2-18b的a8。这说明，为了减小最大半带宽，应尽量顺着网格的短边依次编号，而不要顺着长边编号。图 2-17 图2-l8 利用K矩阵疏稀、带状分布的特点，可设法在计算机中只存放K矩阵上三角或下三角部分中的非零元素，这样，可节省大量存贮量和运算工作量。 3. 奇异矩阵 可以用证明单元刚度矩阵是奇异矩阵的方法，来证明总刚度矩阵也是奇异矩阵。K的奇异性表明，已知节点位移可由总刚度方程(2-52)确定出各节点所受到的外力，反之若已知节点载荷，就不能由式(2-52)确定出节点位移，因为在节点位移中还可以包含任意的刚体位移。为了求出由于结构变形引起的节点位移，必须根据支座约束条件来修正总刚度方程，以消除K的奇异性。 第八节 边界约束条件的处理 根据不同的支承情况，约束条件可分为零位移约束和非零位移约束两类。 一、零位移约束 结构所含单元的多少与约束处理的方法无关。为了突出约束处理的方法，我们假设结构只有一个单元，如图2-19所示。由图中节点2、3的约束条件可知：，需要求解的未知量是：，已知节点载荷=0，其总刚度方程的展开式为： 由于K矩阵中与零位移对应的行和列上的元素，在求时是不起作用的，因而可把它们划去。同时，也应划去上式右边列阵中相应的未知反力。这样修正后，就降低了方程的阶数，使成为：显然，修正后的总刚度矩阵的行列式不等于零，其逆矩阵存在，即可由已知载荷解出节点位移。这种处理方法称为划行划列法或降阶法。特别适合于单元很少，可进行手算练习时采用。 为了编制计算程序更方便，常希望修改后的K矩阵保持原有阶数和排列顺序。为此，可采用在K矩阵中找出与零位移对应的行和列，在其主对角元素处置l，其余元素处置零；在载荷列阵中，在零位移对应的行元素处也置零，在节点位移列阵中，列出全部节点位移分量。修改后的总刚度方程变成： 这个方程完全等价于前面采用降阶法得到的方程和已知的边界约束条件，而且也消除了总刚度矩阵的奇异性。因此，由它可解出未知节点位移分量。 二、非零位移约束 若在图2-19中，支座2、3有升或降，如设。对这类有非零位移约束的情况，最方便和最常用的处理方法是主对角元素乘大数法或放大主对角元素法。 乘大数法的做法是：在K中，用一个很大的数，例如用1030(此数视计算机表示人数的能力而定)乘上已知位移对应的主对角元素，并用同样大的数1030乘以该已知位移和对应的主对角元素，以此代替载荷列阵FP中对应的行元素。修正后的总刚度方程为；把将方程中的第四行乘出来：在上式中，除包含大数1030的两项之外，其余各项与之相比甚小，可略去不计。于是，上式相当于给定位移条件。同理，乘出第六行有。在列阵的第三个元素位置上用零代替，表示位移约束条件。这说明乘大数法对处理零位移约束也是适用的。第九节 计算成果的整理单元的划分、单元及节点的编号、有限元计算模型的建立等是有限元分析的前处理工作，计算成果的整理是有限元分析的后处理工作，它主要包括位移和应力的整理。对于位移的整理，是根据计算机算出的节点位移，由节点的新位置画出结构变形后的网格图形（简称变形图）。对于应力计算结果的整理，工作量较大。由计算出来的是各单元内的应力分量，而不是单元内的主应力。在平面应力状态下，如图2-21a，应利用如下公式进一步计算单元内的主应力： 图2-21主应力和x轴间的夹角，可由图2-21b所示应力圆上的几何关系得出： （2-56）由式（2-55）求出的应力值，常被认为是三角形单元形心处的两个主应力。主应力的方向由式（2-56）确定。在每个单元形心处，沿两个主应力方向画出一定比例的线段，并规定带箭头的线段表示拉应力，带平头的线段表示压应力，如图2-22所示。这样就可画出结构的直观应力分布图。第十节 平面矩形单元当结构的外形都是由一些横线和竖线构成的规则边界时，采用平面艇形单元比较合理，因为它的位移模式次数比常应变三角形单元高，更接近于实际结构的位移和应力分布。图2-27所示平面矩形单元，边长为2a、 2b。以其对移轴作为局部坐标轴x、y。取4个角点为节点，其局部编号按图示位置依次为i、j、m、p。节点的坐标（）单元节点位移有8个分量，其列阵为： （2-59）单元节点力也有8个分量，其列阵为： （2-60）图2-27一、单元位移模式由于矩形单元有8个自由度，所以单元位移模式只能取包含8个待定常数a1a8的双线性形式：（2-61）把4个节点的坐标分别代入式（2-61），应等于相应的节点位移分量，它们形成两个四元线性联立方程组。由此，可解出a1a8。把上式中的代回式（2-61），整理后得：（2-62）可改写成矩阵形式： （2-63）式中，形函数矩阵为：（2-64）其中，形函数为： （2-65a）可统一写成：（2-65b）这些形函数具有如下特性：在双线性位移模式（2-61）中，单元边界上（和）的位移u,v分别是y和x的线性函数，它们可由边界两端的节点位移完全确定，因而，在相邻单元公共边界上的位移是连续的。式（2-61）中的常数项和x、y的一次项，反映了单元位移模式中包含着单元的刚体位移和常量应变，表明双线性位移模式完全满足解答的收敛性条件。二、单元应变与单元节点位移的关系把式（2-62）代入几何方程（2-5）中，得单元应变为：上式又简记为： （2-66）把式（2-65b）代入其子矩阵中，有： （2-67）式（2-67）表明，应变转换矩阵B中的元素包含有变量x、y，而不是常数，这与三角形常应变单元不同。三、单元应力与单元节点位移的关系把式（2-66）代入物理方程中，得单元应力： （2-68）应力转换矩阵S也可写成分块形式： （2-69）对于平面应力问题，其子矩阵为： （2-70）对于平面应变问题，只需把式（2-70）中的E换为把换为即可。 式（2-70）表明，应力转换矩阵S中的元素包含了变量x、y，因而矩形单元内各点的应力不是常量，较好地反映了实际应力的变化。四、单元刚度矩阵推导矩形单元刚度矩阵和推导三角形单元刚度矩阵的方法完全相同，仍是利用虚功方程式（2-31）导出式（2-33）的。不过，在矩形单元面积上进行二重积分时，积分上、下限可简单地写出来，即：（2-71）把B矩阵的分块形式代入式（2-71），可得单元风度矩阵的分块形式为： （2-72）对平面应力问题，其子矩阵为： （2-73）对平面应变问题，只要把式（2-73）中的E、作相应代换即可。五、载荷向节点的移置矩形单元的等效节点载荷列阵中有8个分量： （2-74）根据静力等效原则，由移置前后的载荷在任意虚位移上的虚功相等，仍可导出式（2-40）。单元受均布体力 作用时，仍可用式（2-45）计算只是矩形单元的积分区间与三角形单元不同：把矩形单元的形函数矩阵式（2-64）和形函数式（2-65a）代入上式，得： （2-75）式中，是单元总体力在x、y方向的分量，对于仅受自重力作用