第二章 弹性力学平面问题有限元法

第二章弹性力学平面问题有限元法,2.1有限元法的基本思想及优越性,1.在应用有限元法时，我们首先将一个连续的弹性体看作由许多尺寸有限的小单元-有限元组成。这就是所谓区域划分，在数学上称为“离散化”。2.根据计算对象的简化模型，单元的形状，取成平面三角形或四边形，四面体或六面体等。单元与单元之间，通过若干个称为“节点”的点铰接相连，由此组合成整体。3.以一个个小单元为计算单位，首先进行单元分析，然后把它们组装起来，进行整体分析，最后求出结构的近似解。,这种把复杂结构看成有限个单元组成的整体，就是有限元法的基本思想。,有限元法是从基于能量的变分法发展而来的。如应用最小势能原理的雷利-里兹法，当按位移求解时，它首先要寻找一个满足整个弹性体几何边界条件的位移函数，这对工程实际问题往往有困难。而用有限元法时，将结构进行离散，从一个个单元入手，只要假设单元上的分片插值函数，然后综合起来，代替整个域上的位移函数，这就使问题大为简便和灵活。,因此，有限元法是以变分原理和分片插值为基础的。得到的是近似解。,（1）有限元法直接在力学模型上进行离散化（剖分），物理概念清晰，明白易懂。（2）有限元法有较好的适应性。对于简单问题和复杂问题基本上同等处理。（3）有限元法的各个计算步骤，如单元分析，总体分析和方程解算等都较易标准化和程式化,有一套比较固定的分析顺序，目前已发展成各种通用程序，便于掌握和使用。,有限元法应用于应力分析，按所选取的未知量不同可分为三类：(1)位移法-取节点位移作为基本未知量；(2)力法-取节点力，作为基本未知量；(3)混合法-取一部分节点位移和一部分节点力作为基本未知量。在推导有限元方程时，主要有两种方法：直接法（如直接刚度法）；变分法（如固体力学中的最小势能原理和最小余能原理）把问题归结为求泛函的极值问题。作为初步介绍，我们将以直接刚度法来讨论弹性力学平面问题中的有限元法概念。,有限元模型是真实系统理想化的数学抽象。,节点和单元,节点：空间中的坐标位置，具有一定自由度和存在相互物理作用。,单元:一组节点自由度间相互作用的数值矩阵描述(称为刚度或系数矩阵)单元有线、面或实体以及二维或三维的单元等种类。,载荷,有限元模型由一些简单形状的单元组成，单元之间通过节点连接，并承受一定载荷。,1.每个单元的特性是通过一些线性方程式来描述的。2.作为一个整体，单元形成了整体结构的数学模型。3.尽管梯子的有限元模型低于100个方程（即“自由度”），然而在今天一个小的ANSYS分析就可能有5000个未知量，矩阵可能有25,000，000个刚度系数。,节点自由度是随连接该节点单元类型变化的,J,I,I,J,J,K,L,I,L,K,I,P,O,M,N,K,J,I,L,三维杆单元(铰接),UX,UY,UZ,三维梁单元,二维或轴对称实体单元,UX,UY,三维四边形壳单元,UX,UY,UZ,三维实体热单元,TEMP,J,P,O,M,N,K,J,I,L,三维实体结构单元,ROTX,ROTY,ROTZ,ROTX,ROTY,ROTZ,UX,UY,UZ,UX,UY,UZ,用有限元法对弹性力学平面问题进行应力分析，不仅具有实际意义，而且带有一定的典型性。通过它可以看到:(1)一般情况下此处理问题的方法(2)有限元法的特点(3)使用中的应注意的问题可为今后进一步的深入研究打下基础。,当取节点位移为基本未知量时，有限元法的解题步骤归纳如下：下面我们就按上述顺序介绍。,2.2弹性体的剖分,作为用有限元法解决弹性力学问题的第一步，必须先对弹性体区域进行剖分。对于平面问题来说，最简单的方法是用直线将弹性体区域剖分为有限个三角形或四边形单元。本章将只讨论三个节点的三角形单元。,图2-1,剖分要一直进行到弹性区域的边界上当边界是直线段时，就取其为三角形单元的一条边；当边界是曲线时，则在每小段上用相应的直线近似地代替曲线而作为三角形单元的一边，如图2-1,单元分得越小，计算结构越精确。因此，应当在计算机容量的允许的范围内，尽可能地提高工程上的精确要求，适当地确定单元的大小和数目。,单元的大小和数目要根据精度的要求和计算机容量来确定。,1）任意一个三角形单元的顶点，必须同时也是其相邻的三角形的单元的顶点（如图2-2a)，而不能是其相邻三角形的内点（如图2-2b)。,具体进行剖分时，一般应注意以下几点：,2）尽可能使同一个三角形单元各边的长度相差不太大。此外在三角形单元中最好不要出现钝角。因此，在图2-3a、b两种剖分式中，虽然都涉及到同样的四个顶点，但我们通常都采用a，而不采用b。,3）在事先估计应力较为集中、应力变化较大的地方，例如孔洞附近以及形状突变的角点等处，单元应分得小一些；在应力变化比较平缓的地方，如离开孔洞一定的距离处，单元可以分得比较大一点，如图2-4。,图2-4,有时应力情况事先无法估计，可先采用比较均匀的剖分法进行一次初算，然后经初算的结果重新合理剖分，再进行第二次计算，或用光弹性的方法事先对应力场作一个大概的了解，再在此基础上作合理的剖分和计算，这也是一种常用的方法。,4）在厚度或材料常数有突变的地方，除了应把这些部位的单元分得较小，较密一些以外，还必须把突变线作为单元的分界线。也就是说，在一个单元内部，只能包含一个厚度和一种材料常数。,5）当整个弹性体区域在几何上具有对称轴，而载荷又对称于该轴或反对称于该轴时，则其位移和应力也必然具有这种对称性质。为了减少计算量，只需取其一部分作为求解区域进行单元剖分和计算即可。图2-4实际上只取了整个弹性体区域的四分之一。,作了这样的剖分之后，再以三角形单元的顶点作为节点（注意，如果边界上有集中力，则一般将其作用点选定为节点），然后对单元和节点分别进行编号。编号的顺序不影响计算结果，原则上是可以任意的。但用直接法求解有限元法的基本方程时，从压缩计算机存储量的角度来看，在对节点编号时应注意，单元的两个相邻节点编号之差应尽可能地小。因为这个差值就反映在方程组的系数矩阵（总刚度矩阵）的带宽上，它直接决定了系数矩阵元素的存储数量。有关问题，以后还要作详细的说明。,2.3单元分析,在进行了弹性体的剖分后，可任取一单元作为研究对象。设某三角形单元e的节点编号为i,j,m,（为了在以后的计算中使三角形的面积不致为负值，规定i,j,m的次序为逆时针方向。）并设三个节点i,j,m在右手坐标系的坐标值分别为（xi,yi），（xj，yj），（xm，ym），如图2-5所示。,图2-5,对于平面问题，三个节点的位移分别为：,单元的节点位移列阵为：,所谓单元分析，就是建立节点位移（基本未知量）和单元内任意一点的：位移f，单元应变，单元应力单元节点力Fe之间的关系，使f，Fe等都用节点位移e来表示。如此，则基本未知量e一经求得，其它各量皆可随之而定。,1）节点位移e和单元内任意一点位移f关系首先，我们要确定三角形单元内各点的位移变化规律。即当节点位移确定时，单元内各点的位移应如何插值？设单元内任一点的位移是该点坐标（x,y)的线性函数。对于采用三角形单元的平面问题来说，当单元取得足够小时，取线性位移插值函数是合理的。,即：式中是待定常数，它们可以由单元的边界条件，即节点的位移值来确定。为此，只要将节点的坐标值代入式（a)，就得到节点的位移值：,用克莱姆法则求解线性方程组（b）,（c），得:,式中：而:是三角形单元的面积,将式(d)、(e)代入式(a)，即得单元的位移插值函数：,进行整理后得：若令：为单元的形函数。由式(f)可知，单元内任意点的位移与单元的节点位移是通过形函数来联系的，而形函数则是点的坐标的线性函数。,引入式(2-5)后，式(f)可以表示为：,写成矩阵形式就是：式中为二阶单位矩阵,若则单元上的位移插值函数可表示为：f=Nee应当指出，任意两个相邻的三角形单元，如图中的i,j,m及p,j,i,它们在i和j点具有相同的位移。,我们已假定了位移分量在每一单元中是坐标的线性函数，则在公共边ij上，位移也必然是按同样的线性变化的。因此，在上述两个单元中，公共边ij上各点也都具有相同的位移。这就保证了相邻单元在公共边界上位移的连续性，也即弹性体在受力变形后，各单元的边界线上的材料不致产生空隙或重叠的现象。,2)节点位移e和单元应变分量的关系从上一节可知，在取定单元的位移插值函数以后，只要求得各个节点的位移值，则每个单元内各点的位移（因而也是整个弹性体内各点的位移）即可确定。这一节，我们将联系到平面问题的几何方程和物理方程，来进一步确定单元的应变和应力。,由(2-4)式，即可得到用节点位移表示单元任一点的应变表达式：,写成矩阵形式：,简写为：=Be其中：,称为单元的应变矩阵。,由于单元的面积以及各几何参数bi,ci,cm的值，都可以由节点坐标直接确定，而且均为常数。因此在每一个单元中，应变分量x,y,xy都是常量。故线性位移插值函数的单元又称为常应变单元。,3）节点位移e和单元应力分量的关系由弹性力学已知，平面应力情况下，应力与应变之间的关系可表达为：,写成矩阵形式：,或缩写成：=D(2-12)式中：D弹性矩阵,在平面应变情况下，弹性矩阵为：,将（2-10）式代入（2-12）式，即可得到用节点位移表示的单元内任意一点应力的表达式：=DBe或写成:=Se式中:S=DB称为应力矩阵。,4）节点位移e和单元节点力Fe的关系-单元刚度矩阵k所谓有限元法的“位移法”，就是把节点位移作为未知数。而把单元内各点的位移、应力和应变都表达成节点位移的函数。这样，就把全部问题都归纳为求取节点位移的问题了。在弹性力学中，以应力形式表示的物体内部各部分之间的相互作用力，在有限元法中就相应地以节点力的形式来代替。,所谓节点力就是单元周围部分对我们所考虑的那个单元的作用。设三角形单元i,j,m三个节点上的节点力分量分别为：,则单元的节点力列阵为：,与应力-应变成线性关系相似，就弹性系统而言，节点位移和节点力之间，也同样保持线性关系：,用矩阵形式表示：,或简写成：Fe=kee(2-16b)这就是单元节点位移与节点力的表达式。式中：,称为单元刚度矩阵。,关于它的性质，可以讨论如下：单元刚度矩阵每一列的意义：（为了讨论方便，设i=1,j=2,m=3）令u1=1,而v1=u2=v2=u3=v3=0，由（2-16a）式可得，,由此可见,单元刚度矩阵的第一列元素表示:当节点1在x方向有单位位移（u1=1）而其它位移为零（v1=u2=v2=u3=v3=0）时，各节点上产生的节点力。如k61，就表示当节点1在x方向有单位位移时，在节点3,y方向产生的节点力。因单元在这些力作用下处于平衡，所以x方向和y方向的节点力之和分别为零，其总和亦为零，从而有：k11+k21+k31+k41+k51+k61=0,对于其它各列也有类似的性质，元素kij下标j表示产生单位位移的节点序号和方向，i表示产生节点力的节点序号和方向。,单元刚度矩阵主对角线上的元素为正例如k11,表示节点1在x方向产生单位位移而其它位移均为零时，必须在节点1,x方向施加的力。该力显然应和位移方向一致。因而k11应为正值。,单元刚度矩阵是对称矩阵它的对称性是由弹性结构的反力互等定理（第j个单位位移分量引起的第i个节点力等于第i个单位位移分量引起的第j个节点力）得到的：即：,通过单元刚度矩阵ke,建立了单元节点力列阵Fe与单元节点位移列阵e之间的关系式（2-16)。在有限元法中，必须逐一求出每个单元的单元刚度矩阵。这是计算中必不可少的重要一环。为了求取单元刚度矩阵，我们应用虚功原理推演它的具体表达式。任取某一单元，虚功原理可简述如下：当处于平衡状态的单元体发生约束条件所允许的微小虚位移时，则节点力在虚位移上所作的虚功等于整个单元体的虚变形能。,如以*表示节点虚位移，*表示与虚位移相应的虚应变，即:,虚应变*与节点虚位移*之间的关系仍然如公式（2-10）所示，即:*=B*节点力Fe在节点虚位移*上所作的功为：整个单元体的虚应变能:(其中t为薄板厚度）,按虚功原理，（b）和（c）相等，即:将（a）式和（2-14）代入上式的右端:,根据矩阵相乘的逆序法则：（B*)T=*TBT便有：由于单元节点虚位移*是任意给的微小量，故*T可以提到积分符号外面去。,于是得：与(2-16b)式:比较可得：,由于弹性矩阵D完全决定于弹性常数E和，而在选取线性位移插值函数的条件下，应变矩阵B中的元素又都是常量，因此可把（h）式中的常量都提到积分号外面来，并注意到:是三角形ijm的面积，于是:ke=BTDBt,用矩阵形式表示为：,算例：设某三角形单元ijm如下图(2-8)所示。已知：弹性常数E=1.5,v=0.25;厚度t=0.4试计算:应变矩阵B应力矩阵S单元刚度矩阵ke1.常数计算和面积计算面积=1/2（bicjbjci)=0.5,2.应变矩阵B和应力矩阵S:,3.单元刚度矩阵ke=BTDBt=BTSt,根据上述条件，试计算图(2-9)所示三角形ijm的单元刚度矩阵：,答：,由公式(2-2):可以看出，bi,bj,bm和ci，cj，cm等只与单元的方位有关，不随坐标的平移而变化。可以利用这一性质适当选择原点的位置以简化单元刚度矩阵的计算。,即对于右图。下述两种坐标的选取，所得ke一样。,当把单元刚度矩阵写成2x2的分块子矩阵时，ke可写成下列形式：,式中:,如为平面应变问题，将(2-19)式中的:E换成E/(1-2)换成/(1-)于是得到：,按(2-20a）式或（2-20b）式先依次计算出分块矩阵中的四个元素可同样得到单元刚度矩阵中的36个元素。,在实际计算中，每个单元的节点都有一个编号，称为总体编码（或实际编码）。如图（2-10）。如果按实际编码进行运算，就会给编制程序带来困难；因为每个单元的实际编码并无一定规律。但是如赋予每个单元以局部编码（1）、（2）、（3）这样就可使用同一个过程体（子程序）。使之能十分方便地对局部编码的（1）、（2）、（3）三角形单元进行运算。,总体编码123如单元局部编码总体编码245单元局部编码其余类推,5）等效节点载荷作用在弹性体上的载荷不外乎是：*体力（自重、惯性力）*面力*集中力三种。用有限元法解题时，既然全部问题都归结到节点来处理，那么，当单元上作用有外载荷时，也应把它们移置到节点上来，成为节点载荷。这种移置必须按照静力等效原则进行,所谓“静力等效”，系指原载荷与移置后的节点载荷，在弹性体的任何虚位移过程中的虚功相等。当插值函数已经确定时，这种移置的结果是唯一的。在取线性位移插值时，符合刚体静力等效原则，即：载荷与节点载荷在任一轴上投影之和相等，对任一轴的力矩之和也相等，也就是说，原载荷与节点载荷将具有相同的主矢和主矩。,通常，我们总将集中力的作用点取为节点，不需要移置。因此，下面只讨论：体力和面力的移置,(1)体力的移置,以重力为例来说明这个问题。设有匀质、等厚度、编号为e的三角形单元，三个节点为i,j,m，重力作用在形心c上，如图2-11所示。由初等几何知，,首先，我们求移置到i节点上的垂直节点载荷。为了便于计算虚功，假想该单元在节点i处沿y方向产生一个单位虚位移，而其它两点不动；这相当于在j点及m点安置了铰支座，在i点安置了水平连杆支座，如图2-11。由于我们在三角形单元中采用的位移插值函数是线性的，所以任一条直线上各点位移都呈线性变化。,现在m点及j点的位移都等于0，所以在边上各点位移都等于0；线上各点的垂直位移也按线性变化，在b点等于0，在i点为1。因此c点的垂直位移将为1/3。按静力等效原则，体力载荷W的虚功应等于的虚功，即有：或Yei=-W/3负号表示Yei的方向与图上所画的方向相反。,用同样的方法可以得到Yej=-W/3,Yem=-W/3,下面来求移置到节点i上的水平节点载荷Xei。与前面一样，在形心c处有W作用，假设节点i只沿x方向产生一个单位的虚位移，而其它两点不动。由图2-12可知，c点的垂直位移等于0，水平位移等于1/3。按静力等效原则，有：故:,用同样的方法可以得到：Xej=0,Xem=0,写成分量形式:,由此可以得出如下结论：对于匀质、等厚度的三角形单元，当考虑自重时，只需把1/3的重量移置到每个节点上，就完成了重力载荷的移置，而不必再去列出虚功相等的条件。这也完全符合对刚体的静力等效原则。但必须指出：上述结果是由于我们采用线性位移插值函数造成的。如果位移插值函数是非线性的，例如是坐标的二次函数，那就不满足1/3的关系，也就不能按简单的刚体静力等效原则来处理，而必须用虚功方程来建立普遍的表达式。,(2)面力的移置1）设等厚度的三角形单元e的三个节点为i,j,m其边界ij上受有垂直均匀分布的面力载荷。载荷集度（单位长度上的力）为q，如图2-13所示。,仍然采用线性位移插值函数方法。则根据静力等效原则，将此均匀分布的面力载荷移置到两侧节点i,j上时，等效节点载荷为：式中l为的边长，Rei,Rej仍与原载荷平行，,故此时单元的节点载荷列阵为：,或:,2)若ij边上受有三角形分布的面力载荷，在i点上的载荷集度为,j点上为0，则其单元的节点载荷列阵是：或,3）若ij边上受有垂直的梯形分布的面力载荷，在i,j点上载荷集度分别为和，如图2-15所示，则其等效节点载荷为：,或写成分量形式：,式中分别是在x,y方向上的分量，其方向与x,y轴正向一致为正，反之为负。,2.3整体分析,1.基本方程总刚度矩阵K的形成(1)节点力的组合以上我们分析了一个单元的情况，现在进而研究单元的组合。为了说明问题，今选用一个包含9个节点8个单元的平面问题来分析。如图2-16所示，除1、3、7、9四个节点外，其余五个节点均联接着四个单元。,对于联结着n个单元的节点，一个节点的位移当然涉及到n个单元，与节点位移相应的节点力将是n个单元的综合效应。如节点5的位移涉及到(2)、(3)、(6)、(7)单元。与节点5相应的节点力将与上述四个单元有关。首先，我们按(2-16a)式逐个建立单元节点力和节点位移的关系。下面，选写其中(2)、(3)、(6)、(7)四个单元。,对于单元(2)，节点编号2、5、4,对于单元(3)，节点编号2，6，5,对于单元(6),节点编号4，5，8,对于单元(7),节点编号5，6，8,同样可写出(1)、(4)、(5)、(8)单元的节点力和节点位移的关系。(学习者可自己完成）。对于第5个节点，其X方向的节点力即:U5=U5(2)+U5(3)+U5(6)+U5(7),将上式代入(a),式中:,所以式中的,上式可简写成：F=K式中K称为总刚度矩阵。它是一个对称正定阵(其对角线上各元素为正值，且有元素Kij=Kji)。总刚度矩阵由单元刚度矩阵累加而成，每个单元刚度矩阵对总刚度矩阵都有一定贡献。,K=ke,假定结构离散化后有n个节点，每个节点有两个方程。因此总刚度矩阵为(2nx2n)的矩阵。可将单元刚度矩阵用补零的方法由6x6扩大到(2nx2n)的方阵(图中虚点上的元素均为0）如图所示，则单元刚度矩阵中各元素在总刚度矩阵中的位置即可确定。,例如将第3单元刚度矩阵中的元素填入总刚度矩阵（亦即将该单元刚度矩阵用补零的方法扩大成总刚度矩阵）,在计算程序编制中，我们的做法是先都按局部编码，使用同一过程体算出各单元的刚度矩阵；然后转换成总体编码，最后将相同编码的元素合并成总刚度矩阵中的元素。转换过程示意如下：,(2)节点载荷组合当进行单元组合时，除了考虑节点力的组合外，同时还应进行节点载荷的组合。设结构上的载荷（如体力、面力、集中力等）均已移置到节点上，则单元节点载荷列阵为：,Yme,Xme,Yie,Xie,Xje,Yje,图2-17,x,y,对于联结着两个以上单元的节点，把相同方向上的载荷迭加起来，显然有：,这就是该节点载荷在x和y方向的两个分量。（表示对环绕节点i的单元求和）,若各节点上的Xi和Yi（i=1，2，3n）均已经求出，并按节点编码的顺序排列起来，就得到弹性体总的节点载荷列阵：,(3)平衡方程在求得了节点外力矩阵以后，我们就可以写出位移法中位移分量必须满足的平衡条件。在有限元法中，也就是节点位移必须满足的节点平衡条件。根据公式（2-21）和（2-22），表示所有节点内力与外力平衡的数学表达式为：,通常写成：,等式右端为总节点载荷列阵，当弹性体上的外载荷确定时，它是已知的；总刚度矩阵K由单元刚度矩阵ke集合而成。因此，式（2-23a）是一个以为未知量，以K为系数的线性代数方程组，这是有限元法的基本方程。,其矩阵展开式为：,顺便提一下，为了编程方便，对各矩阵元素的下标，均按其所在位置标定：,即为节点1，x方向的位移u1，为节点1，y方向位移v1,.余此类推。R1为节点1，x方向的节点载荷x1，R2为节点1，y方向的节点载荷y1,.余此类推。,2.关于总刚度矩阵的性质和常用的存贮方法：,在有限元法中，结构总刚度矩阵的性质和常用的存贮方法：在有限元法中，结构总刚度矩阵的阶数为节点总数和自由度数的乘积。如平面问题，自由度为2，若有n个节点，则总刚度矩阵为:(2nx2n)阶。,若一个具有200个节点的小型平面问题，其总刚度矩阵的元素就有400 x400=160000个，为一般中小型电子计算机的内存容量所不允许的。因此，如何缩小总刚度矩阵所需的存贮单元，是有限元法程序编制中一个需要考虑的突出问题。我们可以根据总刚度矩阵的某些特性，寻求节省存贮量的途径。,1)总刚度矩阵是对称阵。利用对称性，我们只需要存贮总刚度矩阵的上三角形部分（或下三角形部分）中的元素。2)当合理编排节点编码时，总刚度矩阵可呈带状的稀疏阵。其中有不少元素为零，而非零元素对称地分布于主对角线的两旁，形成一带状阵。,下面介绍两种常见的压缩存贮方法：,(1)半宽带存放：仍以图2-16为例，分析与该图相应的总刚度矩阵，其元素的分布表示于图2-18中,当把对角线向右平移到第十个元素时，剩下的就全部是零元素了。这是因为节点编排时，就一个单元而言，节点编码最大相差五个号码，而每个节点又有两个方向的位移所致。我们把自对角线开始向右侧平移至最后一根带有非零元素的斜线为止，两线之间一行内包含的元素个数（包括两线之间的零元素）称为最大半带宽。,而刚度矩阵中每一行的半带宽均取决于一个单元中节点编号差和节点的自由度（即有几个方向的位移）。如以NBD表示最大半带宽，以d来表示各单元中最大节点差值。（对于三角形单元，也就是相邻节点编号的最大差值），则对具有两个自由度的平面问题，最大半带宽的计算公式为：,NBD=（d+1）x2,对于以上所述的这样一个对称的，具有NBD半带宽的总刚度矩阵，由于它在半带宽以外有许多零元素，而对角线一侧的元素又和另一侧元素对称相同，因此为了节省计算机内存，可以采用长方形压缩存储法只将半带宽内的元素存放起来，如图（2-18b）所示，这种存贮总刚度矩阵的方法也称为“半带宽存贮”或“等带宽存贮”。,图2-18,如节点总数为n，而每一节点有两个方程，则方程总数为：,NEQ=2n,我们可以用一个两维数组SK1：NEQ;1：NBD即图（2-18b）来存放总刚度矩阵中必须保存的元素。,总刚度矩阵K和长方形压缩存贮的数组SK之间有如下的对应关系：,矩阵K,数组SK,对角线,第一列,r行,r行,r行s列元素,r行(s-r+1)列元素,由K形成SK时，元素的行号相同。新的列号等于原先的列号减去行号加1，即:,新列号=原列号-行号+1,为了节约内存，我们力求减小半带宽NBD，这就是要求在编排节点号码时，应使相邻节点的节点差值d尽可能小。SK右下角三角块中的元素不和K中的任何元素对应，那些元素应是零。,考虑到总刚度矩阵K中各行的半带宽并不相同，有时，由于结构几何形状等原因，某些行的半带宽特别大，而其它行又较小，这种情况下如以NBD为半带宽（即等带宽）存贮，就可能把许多零元素也存贮起来，这对节省计算机存储量来说是不利的。,(2)一维压缩存贮法,一维压缩存贮法是将总刚度矩阵K的下三角形中每一行从第一个非零元素开始，逐行存放入一维数组K1：S中，S是元素的总个数）举例如下，设有一对称正定阵：,在一维数组K1：13中依次存放为：5，3，4，6，0，3，7，2，8，9，0，0，8共13个数。但是，仅使用这样一个一维数组并不能将元素在K中的位置确定下来。,为此，还须将主对角线上的元素在一维压缩存贮中的序号用另一个一维数组N1：2n0存放起来（n0是节点总数）。如对于上述矩阵数组N1：6中存放的是：,1，3，6，8，9，13,即指出，在一维数组K1：13中，第1，第3，第6.个元素是对角元，显然，这两个数组完全确定了各元素在K中的位置。,2.4方程求解,1.引入位移约束条件上一节，我们建立了有限元法的基本方程式：有了这个方程还不能立即求解节点位移，因为到现在为止，我们还没有考虑到弹性体的几何边界条件，即边界位移的约束条件.,很明显，如弹性体的边界没有位移约束，则在外载荷的作用下，它将有产生刚体运动的可能性，反映在基本方程上，其系数矩阵K将是一个奇异阵（对应的行列式的值等于零）。逆矩阵不存在，方程将具有不定解。在进行应力分析时，为了使解具有唯一性（即排除刚体运动），必须根据弹性体具体的边界位移约束条件，对基本方程加以处理，方能求解。,1）引入约束条件的原因（1）结构实际上可能存在若干约束条件。它们应该加以考虑，否则计算的结果将与实际不符。如一端固定，一端铰支的静不定梁，如图2-19所示。,图2-19,节点1，2，3及15既不允许有x方向的位移也不允许有y方向的位移。,(2)有些约束是由于考虑到结构和载荷的对称性（或反对称性）可取其中一部分作为计算对象而附加的。如图2-20所示一对角受压的方形薄板。,图2-20,由于结构和载荷的均对称，可以只计算薄板的四分之一。图2-20b因为变形对称于对角线，水平对角线nn上不可能有y方向的位移，所以附加的约束条件应是：（节点4，5，6为垂直的可动铰支座）；,垂直对角线mm不可能有x方向的位移，所以附加的约束条件应是:（节点1，2，4为水平的可动铰支座），薄板的中心点O，任何方向的位移均为零，故节点4处为固定铰支座。,2）引入约束条件的处理：（1）零位移约束条件的处理举例说明，对于在剖分后有n个节点的弹性体，假定第n个节点处有约束，其位移为零，即：un=0vn=0此外，还已知第（n-1）个节点沿y轴方向有约束，其y方向位移为零，即：vn-1=0（这对平面问题来说，是不产生刚体运动的最低限度的边界位移约束条件）。,则应对基本方程作这样的处理：把刚度矩阵K的最后三行和最后三列划去，得到一个（2n-3）阶的方阵，称为总刚度矩阵的缩聚或降阶。研究表明，降阶以后的总刚度矩阵K是一个非奇异的对称正定矩阵。相应地，分别将节点位移列阵和节点载荷列阵R的最后三行划去，得到（2n-3）维列阵和。这样，经过零位移边界条件处理后的基本方程就成为:,这样得到的结果，不但满足平衡条件和相容条件，而且满足全部边界条件，按上面的假定及处理方法，方程的具体形式是：,式中Kij是第i行，j列的元素。,如果零位移的节点编号不在最后而在中间，也可以用同样的方法划去想相应的行和列，而得到降阶后的式（2-25）：,但是在计算机的计算程序中，我们采用的方法是使总刚度矩阵K中要划去的行和列除对角线元素充成1以外，其余的元素均充成零。如此得到的矩阵不降阶，仍为2n阶方阵，记为，并且也是一个非奇异的对称正定矩阵。与此同时，将节点载荷列阵中相应的元素也充成零，如此得到的2n维列阵记为,因而基本方程就变为：,（2-26）,式中为包含所有节点位移的列阵。显然，在求解式（2-26）时，原来节点位移为零的仍保持为零。式（2-26）与式（2-25）实际上是等价的。按前面同样的假定，式（2-26）的具体形式应该是：,可以明显看出，如果将上式乘开，则后三行的结果是：也即：vn-1=0,un=0,vn=0,表示了边界位移约束条件，由此可见，式（b）与式（a）等价。,有几个约束条件，就应重复几次上面的修改过程。处理边界节点零位移条件，还可以采用其它方法，这里不一一列举。,（2）位移不为零的约束处理,如果某些节点的位移值是已知的，如：我们可将基本方程作如下修改：将第i行的对角元改为1第i行的其它元素改为零将右端项改为,式（2-27a）,i行,i列,这样修改后，第i个方程变为：,为了保持总刚度矩阵的对称性，以便于以后求解，我们把K的第i列也改为零，这时，节点的载荷列阵要作如下修改：,i行,i列,式（2-27b）,当将（2-27a）式和（2-27b）式乘开时，两组代数方程组完全相同。下面举一简例说明:,已知3=u0，由（2-27a）式：（节点位移和节点载荷按列阵中的位置标号）。,由(2-27b)式:,方程（b）等同于方程(a)。(2-27b)式保持系数矩阵对称，目的是使求解方便。,2.基本方程求解-高斯消去法在把位移约束条件引入基本方程（2-23）式后，问题就归结为线性代数方程的求解了。一般来说，求解的方法可分成两大类。一类是直接法，本节要介绍的高斯消去法就是其中的一种；另一类是迭代法，它是按照一定的迭代程序，由假设的初始值，通过多次迭代去逐步逼近方程的精确解。,直接法和迭代法相比，具有速度较快的优点，但所需的存贮量比迭代法大，在计算机存贮量许可的条件下，通常都用直接法求解。下面，我们通过一个简单的算例，说明消去法的基本思想和步骤。,求解线性代数方程组：(a)写成矩阵形式为：,首先，可将（a）式第一式中X1的系数化为1，为此，用2去除第一式的两边,得到：(b),其次，我们把（b）式作为主导方程，使（a）的第二、三式X1的系数化为0,为此，只需用（+1）乘式（b），被第二式减，用（+3）乘式（b）被第三式减，如此，线性代数方程组可化为:,写成矩阵形式为：,同理，可将（c）中第二式X2的系数化为1，得到:(d),以（d）式为主导方程，将（c）式第三式X2的系数消去，得：,最后得到下列形式的线性代数方程组：,写成矩阵形式：,其中系数矩阵的特点是:对角线系数为1，对角线下侧的元素全部为零，这样的方阵称为上三角阵。,由（e）式第三式得：,代入（e）第二式得：,将x3及x2代入（e）第一式得：,至此，可把消去法求解线性代数方程组归纳成两个主要步骤：,第一步是通过“化1消零”把线性代数方程组的系数矩阵化为上三角阵，这一过程称为消元。当这一步骤完成时，方程组中最后一个未知数已求得；第二步称为反代，从最末第二个方程开始，将已求出的未知数代入，逐个求出上一个未知数。,前面已指出：总刚度矩阵K是一个对称正定阵。其主对角线上所有元素均为正值。且有:,对于对称正定阵，其逆矩阵一定存在。因此，以对称正定阵K为系数的线性代数方程组：,必有唯一解。,（1）消元和反代过程：,用矩阵表示：,当用高斯消元时，可将上式第一方程两边各除以K11，使第一方程成为：,（g）为主导方程。,上式两边：乘以K21后与（f）式第二方程相减，乘以K31后与（f）式第三方程相减，乘以K41后与（f）式第四方程相减，便将方程化作如下形式：,此时方程组（h）的系数矩阵的第一行第一列元素为1外，第一列的元素全部为零。,用矩阵形式表示：,可以证明，第一循环后，自第二行和第二列开始的方程仍是一对称正定阵。,因此，按同样的方法，可使第二方程X2前的系数为1，而第三方程及其以下各方程在X2项的系数为零，依次类推，消元过程中的系数矩阵可表示如下：,用矩阵形式表示消完元后的方程：,1.，最后一个右端项的值就是最后一个未知量的解；,2.从倒数第二个方程开始回代，每一个方程中都只有一个未知量，依次可求得全部节点位移值；,（2）递推公式：,上面的消元反代作法同样适用于系数矩阵为n阶的线性代数方程组，由上过程可知，消第一元时，系数矩阵从f1变到f2，也就是说，每消元一次，要作四件事（建立四个递推公式）。,（1）修改非主导方程的系数项（2）修改主导方程的系数项（3）修改非主导方程的右端项（4）修改主导方程的右端项,其中,相应于(2)(3)(4)项的递推公式容易建立,唯有非主导方程系数项的处理涉及到有关元素如何从刚度矩阵K中过渡到SK中的问题。,递推公式：,值得注意的是：,1.每消元一次，需修改的方程只涉及到NBD个；主导方程1非主导方程2，3，4，.NBD,2.随着被修改方程的行数增加，方程中需修正的系数逐个减少（呈三角形，见下图）。3.在总刚K中的对角元，在SK中变为第一列。,注意:（1）最后一个右端项就是最后一个未知量的解。（2）反代从第NEQ-1个方程开始，倒推。反代结束。R阵是节点位移。（3）每反代一个方程涉及该方程的NBD-1列元素和NBD个右端项。,2.5应力计算,全部位移求得后，各单元的节点位移确定后，单元内的应力可由公式（2-14）求得：,其中应力矩阵S=DB,在计算程序中,计算某一单元的应力时，首先要确定该单元的三个节点的六个位移分量是已求得的全部节点位移中的那几个。相对应的编号关系建立后，就可以把这些位移分量代入应力分量关系式（2-14）求得x，y，xy。,如仍以图（2-16）为例，当消元的反代过程结束后，节点位移分量全部求得：,如欲确定单元（4）的应力，按（2-14）式:,需从全部位移分量中找出单元（4）相应的位移分量。,