有限元分析基础平面问题的有限元

第三章平面问题的有限元§3.1三角形单元§3.2矩形单元§3.3平面等参元§3.4平面单元应用比较分析第三章平面问题的有限元3.1三角形单元§3.1.1位移函数§3.1.2应变与应力矩阵§3.1.3形函数的性质§3.1.4单元刚度矩阵及其特点§3.1.5等效节点载荷§3.1.6整体刚度矩阵及其特点§3.1.7边界条件§3.1.8算例§3.1.9收敛准则§3.1.10从能量原理推导刚度矩阵本节主要研究三节点三角形单元，首先建立以单元节点位移表示单元内各点位移的关系式。设单元e的节点编号为i、j、m，如图所示。由弹性力学平面问题可知，每个节点在其单元平面内的位移有两个分量，所以整个三角形单元将有六个节点位移分量，即六个自由度。用列阵可表示为式中的子向量为

（i、j、m轮换）式中：——节点i在x轴方向的位移分量；

——节点i在y轴方向的位移分量。3.1.1位移函数在有限单元法中，虽然是用离散化模型来代替原来的连续体，但每一个单元体仍是一个弹性体，所以在其内部依然是符合弹性力学基本假设的，弹性力学的基本方程在每个单元内部同样适用。从弹性力学平面问题的解析解法中可知，如果弹性体内的位移分量函数已知，则应变分量和应力分量也就确定了。但是，如果只知道弹性体中某几个点的位移分量的值，那么就不能直接求得应变分量和应力分量。因此，在进行有限元分析时，必须先假定一个位移模式。由于在弹性体内，各点的位移变化情况非常复杂，很难在整个弹性体内选取一个恰当的位移函数来表示位移的复杂变化，但是如果将整个区域分割成许多小单元，那么在每个单元的局部范围内就可以采用比较简单的函数来近似地表示单元的真实位移，将各单元的位移模式连接起来，便可近似地表示整个区域的真实位移函数。这种化繁为简、联合局部逼近整体的思想，正是有限单元法的绝妙之处。3.1.1位移函数基于上述思想，可以选择一个单元位移模式，单元内各点的位移可按此位移模式由单元节点位移通过插值而获得。假设三节点三角形单元内任意点的位移是x、y的线性函数，如下是待定常数，因三角形单元共有六个自由度，且位移函数在三个节点处的数值应该等于这些点处的位移分量的数值。假设节点i、j、m的坐标分别为（式中～、）、（）、（），代入位移函数，得3.1.1位移函数从上式中可以解得式中：——三角形i、j、m的面积。为保证求得的面积为正值，节点i、j、m的编排次序必须是逆时针方向，如图所示。的计算公式，如下3.1.1位移函数将式所求得得系数代入位移函数，得式中：（i、j、m轮换）3.1.1位移函数同理可以求出单元内任意点的y方向位移，如下令

（i、j、m轮换）则单元内任意点的位移可以表示成如下形式式中：——单元内位移矢量3.1.1位移函数位移函数可以进一步写成如下形式式中：

（i、j、m轮换）是坐标的函数，它们反映了单元的位移状态，所以一般称之为形状函数，简称形函数。显然，形函数决定了单元内的“位移模式”，反映了i节点位移对单元内任意点位移的贡献。3.1.1位移函数有了单元的位移模式，就可以利用平面问题的几何方程求出单元内任意点的应变，式中：——单元应变列阵；

——单元应变转换矩阵，或称单元几何矩阵。3.1.2应变与应力矩阵单元应变转换矩阵，为3×6矩阵，其子矩阵为（i、j、m轮换）可以看出，表示i节点位移对单元应变的贡献率。当单元确定后，也就确定了，此时单元内的应变仅依赖于节点的位移。对于三节点三角形单元，由于面积和

等都是常量，所以单元应变转换矩阵中的诸元素都是常量，因而单元中各点的应变分量也都是常量，即单元内各点的应变相同，通常称这种单元为常应变单元。、3.1.2应变与应力矩阵求得应变之后，再将应变代入物理方程，即公式（2.19），便可推导出以节点位移表示的应力，如下式中：——单元应力矩阵。令则式中：——单元应力转换矩阵。3.1.2应变与应力矩阵的子矩阵为单元应力转换矩阵可以看出，表示i节点位移对单元应力的贡献率。当单元确定后，也就确定了，此时单元内的应力仅依赖于节点的位移。对于三节点三角形单元，由于面积和

等都是常量，所以单元应力转换矩阵中的诸元素都是常量，因而单元中各点的应力分量也都是常量，即单元内各点的应力相同，通常称这种单元为常应力单元。、3.1.2应变与应力矩阵可见，对于三节点三角形单元，由于所选取的位移模式是线性的，因而导致单元内各点的应变、应力相同，但相邻单元将具有不同的应力和应变，即在单元的公共边界上应力和应变的值将会有突变，但位移却是连续的。3.1.2应变与应力矩阵由前文可知，形函数决定了单元内的“位移模式”，反映了i节点位移对单元内任意点位移的贡献，形函数具有如下性质：1.形函数在各单元节点上的值，具有“本点是1、它点为零”的性质即

（i、j、m轮换）证明如下：3.1.3形函数的性质在节点j、m上类似有在节点i上3.1.3形函数的性质2.在单元内任意点上，三个形函数之和等于1即证明：将单元内任意点的坐标代入上式，得，因此容易得到此结论。由此可见，三个形函数中只有二个是独立的。可以证明，，3.1.3形函数的性质例如：假设单元各节点的位移均相同并等于，则根据此性质可得单元内任意点的位移为，显然，该性质反映了单元刚体位移，即当单元做刚体运动时，单元内任意点的位移均等于刚体位移。3.1.3形函数的性质3.三角形单元任意一条边上的形函数，仅与该边的两端节点坐标有关。例如在ij边上，有事实上，因ij

边的直线方程方程为3.1.3形函数的性质代入及，得则3.1.3形函数的性质另外，由可以求得利用形函数的这一性质可以证明，相邻单元的位移分别进行线性插值之后，在其公共边上将是连续的。例如，对图3.2所示的单元具有ij公共边，可知，在ij边上有这样，不论按哪个单元来计算，公共边ij上的位移均由下式表示3.1.3形函数的性质由此可见，在公共边上的位移u、v将完全由公共边上的两个节点i、j的位移所确定，因而相邻单元的位移是连续的。例3.1求如图3.3所示等边三角形的形函数。三角形的面积将以上数据代入式（3.10），得3.1.3形函数的性质①将各节点坐标代入以上三式，有即满足“本点是1、它点为零”的性质。3.1.3形函数的性质即满足“在单元内任意点上，三个形函数之和等于1”的性质。②将三个形函数相加，得，即在ij边上任意位置的位移，将完全由节点i、j的位移所确定。取其它边同样可以得出其结论，即满足“三角形单元任意一条边上的形函数，仅与该边的两端节点坐标有关”的性质。③以ij边为例，在ij边上，则3.1.3形函数的性质单元在节点处受力，单元会发生变形，也就是说单元在节点处受到的力与单元节点位移之间有必然的联系。单元间正是通过与节点间的相互作用力连接起来成为整体，而一个单元仍是一个弹性体，如果将每一个单元的受力－位移关系找到，则整体的受力－位移关系也容易清楚。为了推导单元的节点力和节点位移之间的关系，可应用虚位移原理对单元进行分析。单元在节点处受到的力称为单元节点力，单元在节点力的作用下处于平衡，节点力可以表示为

式中：——节点i沿x方向的节点力分量；

——节点i沿y方向的节点力分量。3.1.4单元刚度矩阵及其特点单元节点的虚位移为相应的单元内的虚位移场为

单元内的虚应变为

设单元仅在节点上受力，则单元节点力在节点虚位移上的虚功为

3.1.4单元刚度矩阵及其特点单元内的应力在虚应变上所做的功为式中：——在体积上积分。令单元刚度矩阵为由虚功原理可知，式（3.27）和（3.28）相等，再由虚位移的任意性，可得3.1.4单元刚度矩阵及其特点上式为单元刚度矩阵的定义式，其推导过程与形函数的具体表达形式、节点个数均无关，这说明该表达式具有普遍意义，则上式是单元节点力与单元节点位移之间的关系式，即单元节点平衡方程，或称单元刚度方程，具有普遍意义。对于三节点三角形单元，由于、中的元素都是常量，这里假设单元厚度均匀，则式中：——单元面积；

——单元厚度。3.1.4单元刚度矩阵及其特点上式进一步可以表示为将几何矩阵式（3.16）和弹性矩阵式（2.19）代入，可得每个子块的矩阵表达式，如下3.1.4单元刚度矩阵及其特点单元刚度矩阵具有如下特点：1.中的每一个元素都是一个刚度系数假设单元的节点位移为，由，得到节点力为

3.1.4单元刚度矩阵及其特点由此可得：表示i节点在水平方向产生单位位移时，在节点i的水平方向上需要施加的节点力；表示i节点在水平方向产生单位位移时，在节点i的垂直方向上需要施加的节点力。选择不同的单元节点位移，可以得到单元刚度矩阵中每个元素的物理含义，如下：表示s节点在水平方向产生单位位移时，在节点r的水平方向上需要施加的节点。表示s节点在水平方向产生单位位移时，在节点r的垂直方向上需要施加的节点力；表示s节点在垂直方向产生单位位移时，在节点r的水平方向上需要施加的节点力；表示s节点在水平方向产生单位位移时，在节点r的垂直方向上需要施加的节点力；因此，中的每一个元素都是一个刚度系数，表示单位节点位移分量所引起的节点力分量，中的每一个子块表示单元某个节点位移矢量对单元某个节点力矢量的贡献率。3.1.4单元刚度矩阵及其特点2.是对称矩阵由式（3.32）可知，是对称矩阵。根据对称性，可以减少计算量和存储量。

3.是奇异矩阵，即根据这一性质，当已知单元节点位移时，可以从式（3.31）中求出单元节点力。反之，由于单元刚阵奇异，不存在逆矩阵，因此，当已知单元节点力时不能求出单元上的节点位移。从物体变形的实际情况来说，单元刚阵奇异是必须的，因为单元除了产生变形外，还会产生刚体位移，即仅仅依靠节点力是无法唯一地确定刚体位移的。事实上，当单元的节点力为零时，单元仍可做刚体运动。3.1.4单元刚度矩阵及其特点例如，假定单元产生了x方向的单位刚体位即并假设此时对应的单元节点力为零，则由

得可以得到，在单元刚度矩阵中1，3，5列中对应行的系数相加为零，由行列式的性质可知，。3.1.4单元刚度矩阵及其特点

4.单元的刚度不随单元或坐标轴的平行移动而改变单元的刚度取决于单元的形状、大小、方向和弹性系数，而与单元的位置无关，即不随单元或坐标轴的平行移动而改变。因此，只要单元的形状、大小、方向和弹性系数相同，无论单元出现在任何位置均有相同的单元刚度矩阵，根据对称性，可以减少计算量。，厚度为t。例3.2求图3.4所示等腰直角三角形单元的刚度矩阵，设解根据1（0，a）、2（0，0）、3（a，0）的坐标和式（3.8）可得3.1.4单元刚度矩阵及其特点单元面积，将以上各数据代入几何矩阵式（3.16）为由于，则平面应力与平面应变问题的弹性矩阵相同，为3.1.4单元刚度矩阵及其特点由式（3.18）可得应力转换矩阵，如下则单元刚度矩阵为3.1.4单元刚度矩阵及其特点从以上求得的结果可以看出，单元刚度矩阵是对称的。1、3、5列相加等于零，2、4、6列相加也等于零，单元刚度矩阵是奇异的。另外在几何矩阵中所涉及到的元素均为“坐标差量”，而与坐标的具体值无关，从而也说明了单元刚度矩阵不随单元或坐标轴的平行移动而改变。1）获取单元节点信息；2）计算单元面积及参数、等；3）计算几何矩阵中各元素；4）根据问题类型（平面应力、平面应变）计算弹性矩阵；5）计算应力转换矩阵；6）计算单元刚度矩阵。从单元刚度矩阵的计算公式可以看出，在求解单元刚度矩阵时，可以不必求得应力转换矩阵，但在后期求解单元应力时要用到该矩阵，所以在此处一并求出，也可以根据实际情况具体考虑。3.1.4单元刚度矩阵及其特点实际上作用在物体上的外力可以直接作用在节点上，也可以不作用在节点上，然而在有限元法中要求作用在物体上的各种外力必须用作用在节点上的力表示。这一点体现了有限法中“离散化”这一概念，即将连续体离散成只有在节点处相连的单元体，单元与单元之间的联系只能通过节点。如前文所述，单元内任意点的位移、应力、应变等变量最终都用单元节点位移来表示。同样，作用在物体上的各种外力也必须用作用在节点上的力表示，这一过程称为外力的静力等效移置，所得到的节点力称为等效节点力。

3.1.5等效节点载荷静力等效原则：

对于刚体来说，所谓静力等效原则就是单元上原有的外力系和将外力系向各节点移置所得的等效节点力，二者向同一点简化应具有相同的主失和主矩；对于弹性体来说，所谓静力等效原则就是指单元上的外力系和将该力系向各节点移置后的等效节点力，二者在虚位移上的虚功相等，也即外力作用在单元上所引起的变形能和移置后的等效节点力在单元上引起的变形能相等，在一定的位移模式下这种移置是唯一的。3.1.5等效节点载荷1.弹性体静力等效原则——虚功原理虚功等效的思路是：就一个单元来说，把作用在单元上的外力移置到节点上后，应当与原来的实际外力所作虚功等效。其计算方法是：对任意允许的微小虚位移，原外力所作虚功等于移置后的等效节点力所作虚功。即式中：——节点虚位移；——单元虚位移场；——集中力作用处产生的虚位移。3.1.5等效节点载荷等式左边表示等效节点力在节点虚位移上所作虚功；右边第一项表示作用在单元任意点c处的集中力所作虚功、第二项表示表面力所作的虚功、第三项表示体积力所作的虚功。将代入上式，并考虑的任意性，得式中：——由集中力等效的单元等效节点力；——由表面力等效的单元等效节点力；——由体积力等效的单元等效节点力；——形函数在集中力作用处的值。3.1.5等效节点载荷单元等效节点力、和的表达式分别为

以上三式不但适用三角形单元，同样适用其它单元，具有普遍意义。下面根据以上三式，分别进行讨论。3.1.5等效节点载荷可以看出，当外力作用在节点上时，如作用在i节点上，此时，，则计算所得到的等效节点力与原外力在数值上是相同的。3.1.5等效节点载荷（1）作用在单元任意点的集中力3.1.5等效节点载荷（2）作用在单元上的表面力现举例说明，如图所示的单元e，在ij边上作用有表面力，假设ij边的长度为l，单元厚度t为常数，其上任一点P距节点i的距离为s。有，，，则例3.3如图所示，设31边的长度为l，单元厚度为t，其上作用沿x方向的均布荷载，求其等效节点荷载。解：在31边上，，

，则等效节点荷载为从计算结果来看，等效后的节点力是将原来作用在31边上的均布面力平分到1、3节点上。3.1.5等效节点载荷（3）作用在单元上的体积力设单元所受体积力为，单元厚度t为常数，如所示，则3.1.5等效节点载荷当单元体是均质、等厚、比重为时，则、，有即当单元有均匀自重时，在三节点三角形单元中每个节点沿y方向的等效节点力等于单元总重量的1/3。3.1.5等效节点载荷2.刚体静力等效原则——力系简化（1）作用在单元任意点的集中力若集中力作用在节点上，则可直接将集中力认为是等效节点力即可。若集中力作用于单元的某一边界上，如图所示，则先将分解为、，然后按线段的比例把分别移置到i，j两点上。即3.1.5等效节点载荷（2）作用在单元上的表面力如图所示，已知在ij边受有均布面力，单元厚度t为常数，则移置到i、j节点上的等效节点力为如图所示，当某一边上有三角形分布的面力时，则等效节点力为

3.1.5等效节点载荷（3）作用在单元上的体积力如图所示，如果三角形单元ijm的重心c上受有体积力，单元厚度t为常数，则按刚体静力等效原理可把体积力直接移置到i，j，m三个节点上，即根据圣维南原理可知，按刚体静力等效原则移置后的等效节点力所引起的应力误差只是局部的不会影响到整体，即物体内部距实际外力作用处相当远的各点处应力，与其外力的具体分布情况关系很小。在线性位移模式下，弹性体和刚体按各自的静力等效原则移置外力的结果是一致的。3.1.5等效节点载荷

整体节点载荷列阵：由各节点所受等效节点力按节点号码以从小到大的顺序排列组成的列阵。等效节点力是由集中力、表面力和体积力共同移置构成的，其中集中力包括直接作用在弹性体上的外力和边界约束力，如支座反力。前文对单元体进行了分折，得到了单元刚度方程，但要解决问题，还必须进一步建立整个计算模型的整体刚度方程。完成这一步的关键，在于怎样将单元的刚度矩阵和节点荷载列阵，分别“组装”成整体刚度矩阵和整体节点荷载列阵。这里通过研究任意节点的平衡来建立整体刚度矩阵，该方法不但比较直观、易懂，而且对怎样编写计算机程序是很有帮助的。为了研究整体刚度矩阵的组装过程，先引入两个概念。

整体节点位移列阵：由各节点位移按节点号码以从小到大的顺序排列组成的列阵。3.1.6整体刚度矩阵及其特点式中：，

3.1.6整体刚度矩阵及其特点不失一般性，仅考虑计算模型中有4个单元，如图所示。四个单元的整体节点位移列阵为对每个单元都可以写出相应的单元刚度方程，即单元节点平衡方程。例如，对①号单元，有式中：——①号单元中第i（i=1,2,3）节点所受力。为了便于组装整体刚度矩阵，将上式以整体节点位移表示，即

3.1.6整体刚度矩阵及其特点——①号单元的扩大刚度矩阵或称为单元贡献矩阵。同理，对于②单元，有式中：——②号单元中第i（i=1,3,4）节点所受力；——②号单元的扩大刚度矩阵。3.1.6整体刚度矩阵及其特点对于④单元，有式中：——④号单元中第i（i=3,4,5）节点所受力；——④号单元的扩大刚度矩阵。3.1.6整体刚度矩阵及其特点对于任意一个节点，可能承受两种力的作用，一种是其它单元给予该节点的反作用力；另一种是作用在节点上的等效节点力。对整体而言，前者属于内力，后者属于外力，每个节点在两种力的作用下处于平衡。将各单元刚度方程左边相加，即将各节点所受力相加，由于对于整体而言，单元给予节点的反作用力属于内力，在相加过程中相互抵消，所以各节点所受力相加的结果只有外力，即等效节点力，从而得到整体节点荷载列阵，如下3.1.6整体刚度矩阵及其特点将各单元刚度方程右边相加，从而得到整体刚度矩阵，如下3.1.6整体刚度矩阵及其特点通过以上分析得，整体节点载荷与整体节点位移之间的关系式，即结构整体有限元方程，如下式中：——整体刚度矩阵。3.1.6整体刚度矩阵及其特点整体刚度矩阵组装的基本步骤：1）将单元刚度矩阵中的每个子块放到在整体刚度矩阵中的对应位置上，得到单元的扩大刚度矩阵。注意对于单元刚度矩阵是按照局部编码排列的，即对应单元刚度矩阵中的i、j、m；对于整体刚度矩阵是按照整体编码排列的，即按节点号码以从小到大的顺序排列。在组装过程中，必须知道单元节点的局部编码与该节点在整体结构中的整体编码之间的关系，才能得到单元刚度矩阵中的每个子块在整体刚度矩阵中的位置。将单元刚度矩阵中的每个子块按总体编码顺序重新排列后，可以得到单元的扩大矩阵。例如在图中，单元②的局部编码为i、j、m，对应整体编码为1、3、4，然后将单元②刚度矩阵中的每个子块按总体编码顺序重新排列后，可以得到单元的扩大矩阵。注意有些书籍中将局部编码表示为1、2、3或1，2，3等；2）将全部单元的扩大矩阵相加得到整体刚度矩阵。3.1.6整体刚度矩阵及其特点通过以上组装过程可以得到组装整体刚度矩阵的一般规则：1）结构中的等效节点力是相关单元结点力的叠加，整体刚度矩阵的子矩阵是相关单元的单元刚度矩阵子矩阵的集成；2）当整体刚度矩阵中的子矩阵中r=s时，该节点（节点r或s）被哪几个单元所共有，则就是这几个单元的刚度矩阵中的子矩阵的相加。如应该是单元①－④中对应子矩阵的集成，即3.1.6整体刚度矩阵及其特点3.1.6整体刚度矩阵及其特点3）当中时，若rs边是组合体的内边，则就是共用该边的两相邻单元刚度矩阵中的子矩阵的相加。如13边为单元①和②的共用边，则4）当中r和s不同属于任何单元时，则=0。如节点r＝1和s＝5不同属于任何单元，此时=0。上述组装基本步骤和规则具有普遍意义，对于不同类型、不同形式的单元，只是相应节点的子矩阵的阶数（节点自由度×节点自由度）可能不同，至于组装整体刚度矩阵的规律仍是相同的。正是应为有了这种组装规律，使得有限元法能够很方便地应用电子计算机进行计算。例3.4如图所示有限元模型，弹性模量为，厚度为，为简化计算取，求整体刚度矩阵。3.1.6整体刚度矩阵及其特点单元编号①②③④整体编码1、2、32、4、55、3、23、5、6局部编码i、j、mi、j、mi、j、mi、j、m以整体编码表示的单元刚度矩阵子块解：该模型中共有6个节点，4个单元，各单元的信息如表3.1所示。表3.1各单元信息3.1.6整体刚度矩阵及其特点同例3.2类似的分析，得根据单元刚度矩阵的性质，可知，由于几何矩阵，所以有，整体刚度矩阵中的各子块是对所有单元相应的子块求和得到的（实际只是对相关单元求和），其中各子块矩阵均为2行×2列，整体刚度矩阵用子块矩阵可以表示为3.1.6整体刚度矩阵及其特点上式中任意一子块矩阵均为2行×2列，在计算过程中，无需将每个单元刚度矩阵进行扩大，只需判断整体刚度矩阵子块的下标，然后利用组装整体刚度矩阵的一般规则进行计算，如，由图形可知，节点2由单元①、②和③所共有，则3.1.6整体刚度矩阵及其特点，由图形可知，25边为单元②和③的共用边，则，由图形可知，节点1、5不同属于任何单元，则采用同样的方法进行计算，得到整体刚度矩阵为3.1.6整体刚度矩阵及其特点3.1.6整体刚度矩阵及其特点4.是奇异矩阵，在排除刚体位移后，它是正定阵3.是稀疏矩阵，非零元素呈带状分布用有限元方法分析复杂工程问题时，节点的数目比较多，整体刚度矩阵的阶数通常也是很高的。那么在进行计算时，如果存储整体刚度矩阵的全部元素，将会浪费较大的资源、降低计算效率。如果根据整体刚度矩阵的特点进行编写程序，可以大大节省资源、并提高计算效率。因此有必要了解和掌握整体刚度矩阵的特点，整体刚度矩阵具有以下几个显著的特点：1.是对称矩阵2.中主对角元素总是正的3.1.6整体刚度矩阵及其特点3.1.6整体刚度矩阵及其特点1.是对称矩阵由单元刚度矩阵的对称性和整体刚度矩阵的组装过程，可知整体刚度矩阵必为对称矩阵。利用对称性，在计算机编写程序时，只存储整体刚度矩阵上三角或下三角部分即可。2.中主对角元素总是正的例如，刚度矩阵中的元素表示节点2在x方向产生单位位移，而其它位移均为零时，在节点2的x方向上必须施加的力；表示节点2在y方向产生单位位移，而其它位移均为零时，在节点2的y方向上必须施加的力。很显然在此情况下力的方向应该与位移方向一致，故应为正号。3.1.6整体刚度矩阵及其特点3.是稀疏矩阵，非零元素呈带状分布如果遵守一定的节点编号规则，就可使矩阵的非零元素都集中在主对角线附近呈带状。整体刚度矩阵中的子矩阵只有当下标s等于r或者s与r同属于一个单元时才不为零，这就说明，在第r双行中非零子矩阵的块数，应该等于节点r周围直接相邻的节点数目加1。可见，中元素一般都不是填满的，是稀疏矩阵，且非零元素呈带状分布。以下图所示的单元网格为例，其整体刚度矩阵中的非零子块（每个子块为2行2列）的分布情况如下图所示，图中阴影部分表示该子块不为零，其它子块部位均为零。3.1.6整体刚度矩阵及其特点3.1.6整体刚度矩阵及其特点显然，带状刚度矩阵的带宽取决于单元网格中相邻节点号码的最大差值D。把半个斜带形区域中各行所具有的非零元素的最大个数叫做刚度矩阵的半带宽（包括主对角元），用B表示，如下

B=2(D+1)通常的有限元程序，一般都利用刚度矩阵的对称和稀疏带状的特点，在计算求解中，只存储上半带的元素，即所谓的半带存储。因此，在划分完有限元网格进行节点编号时，要采用合理的编码方式，使同一单元中相邻两节点的号码差尽可能小，以便节省存储空间、提高计算效率。3.1.6整体刚度矩阵及其特点对于同样的有限元单元网格，按照图（a）的结点编码，最大的半带宽为B=2×(10－4+1)=14；按照图（b）的结点编码，最大的半带宽为B=2×(10－2+1)＝18；按照图（c）的结点编码，最大的半带宽为B=2×(10－6+1)＝10。(a)(b)(c)3.1.6整体刚度矩阵及其特点4.是奇异矩阵，在排除刚体位移后，它是正定阵无约束的弹性体（或结构物）的整体刚度矩阵是奇异的，不存在逆矩阵，即关于位移的解不唯一。这是因为弹性体在外力的作用下处于平衡，外力的分量应该满足三个静力平衡方程。这反映在整体刚度矩阵中就意味着存在三个线性相关的行或列，所以是个奇异阵，不存在逆矩阵。例如：设弹性体在外力的作用下处于平衡，这时相应的解为，然后在给予弹性体以刚体位移而相应的节点位移，这时，仍是问题的解，因为刚体位移不会破坏平衡。注：当排除刚体位移后，整体刚度矩阵是正定矩阵。前文已经提到在排除刚体位移后，整体刚度矩阵是正定的，方程才可求得唯一解。排除刚体位移可以通过引入边界约束条件来实现，这里介绍两种比较简单的引入已知位移的方法。1．代入法2．乘大数法3.1.7边界条件3.1.7边界条件1．代入法该方法保持方程组仍为2n×2n系统，仅对整体刚度矩阵和整体载荷列阵进行修正。下面以一个只有四个方程的简单例子加以说明，方程如下假定系统中节点位移、，则当引入这些节点的已知位移之后，方程就变成3.1.7边界条件若，则然后，用这组维数不变的方程来求解所有的节点位移。显然，其解答仍为原方程的解答。在手算时，可直接将零位移约束所对应的整体刚度矩阵中的行和列直接划去，使得整体刚阵的维数变小，更便于手算。3.1.7边界条件2．乘大数法将中与指定的节点位移相对应的主对角元素乘上一个大数，同时将中的对应元素换成指定的节点位移值、该大数与节点位移相对应的主对角元素三者的乘积。若把此方法用于上面的例子，则方程就变成该方程组的第一个方程为解得，这种方法就是使中相应行的修正项远大于非修正项。3.1.7边界条件1．代入法2．乘大数法在以上的两种方法中，代入法接近人工解法，虽然该方法比较直观，但该方法对刚度矩阵改变较多，程序效率不高。乘大数法对刚度矩阵改变较少，工作量较小，但相乘的“大数”若取得过大，求解时会发生“溢出”、若取得太小则会引起较大的误差。根据前面的讨论，现以三角形常应变单元为例来说明应用有限元法求解弹性力学平面问题的具体步骤。1）力学模型的确定。根据工程实际情况确定问题的力学模型，并按一定比例绘制结构图、注明尺寸、载荷和约束情况等；2）结构进行离散化。将弹性体划分为许多单元，并对节点进行编号，确定全部节点的坐标值；对单元进行编号，并列出各单元三个节点的节点号；3）计算等效节点载荷，形成整体载荷列阵；4）计算单元刚度矩阵；5）组装整体刚度矩阵；6）处理约束，引入边界条件；7）求解线性方程组，得到节点位移和节点力；8）计算单元应力；9）整理计算结果（后处理部分）。3.1.8算例例3.5如图所示为一厚度t=1cm的均质正方形薄板，边长为2m，上下受均匀拉力=106N/m2，材料弹性模量为E，泊松比，不记自重，试用有限元法求其应力分量。

3.1.8算例解：1）力学模型的确定。由于此结构长、宽远大于厚度，而载荷作用于板平面内，且沿板厚均匀分布，故可按平面应力问题处理，考虑到结构和载荷的对称性，可取结构的1/4来研究。3.1.8算例表3.2节点坐标值节点坐标1234x0110y00112）结构进行离散化。为了计算简便并能说明整个计算过程，将该1/4结构离散为两个三角形单元。节点编号、单元划分如图所示，各节点的坐标值如表3.2所示，各单元编码如表3.3所示。3.1.8算例表3.3单元编码单元号①②整体编码1、2、33、4、1局部编码i、j、mi、j、m3）计算等效节点载荷，形成整体载荷列阵。按照弹性体静力等效原则，应用公式（3.41），并参考例3.3可得单元②由均布载荷等效的节点力为3.1.8算例若按照刚体静力等效原则，参考公式（3.45），同样可得到以上结果。考虑到支座反力，整体载荷列阵为4）计算单元刚度矩阵。对于单元①有，局部编码i、j、m对应整体编码1、2、3，利用公式（3.34）并参考例3.2，可得单元①的刚度矩阵为3.1.8算例对于单元②有，局部编码i、j、m对应整体编码3、4、1，由单元刚度矩阵的性质可知，单元②的刚度矩阵与单元①的刚度矩阵相同。

5）组装整体刚度矩阵。在组装整体刚度矩阵时，要注意单元局部编码与整体编码之间的关系，按照组装整体刚度矩阵的一般规则和基本步骤，可得整体刚度矩阵为3.1.8算例

6）处理约束，引入边界条件。由以上分析可知，整体有限元方程为3.1.8算例根据边界约束条件：，采用代入法引入边界条件，划去整体刚度矩阵中1、2、4、7的行和列，得3.1.8算例

7）求解线性方程组，得到节点位移。求解上面方程组可得出相应得节点位移为3.1.8算例所以，整体节点为列阵为代入整体有限元方程，便可得整体载荷列阵为若要求得每个单元所受的节点力，则可以从整体节点位移列阵中提取单元所对应的节点位移，然后在利用，便可求得单元节点所受的力。如单元①所对应的节点位移为3.1.8算例求得单元①所对应的节点力为单元②所对应的节点位移为求得单元②所对应的节点力为注意，这里求得的单元节点位移、和单元所对应的节点力、均对应局部编码。很显然，每个单元均满足力的平衡。8）计算单元应力。先求出各单元的应力矩阵、，然后再求得各单元的应力分量，如下3.1.8算例9）整理计算结果（后处理部分）。对于三节点三角形单元，单元内各点的应力值相等，算出的应力一般作为单元形心处的应力。由于单元应力为常数，整个结构的应力场呈阶梯状，在单元之间不连续。而工程上往往更加关心边界和节点上的受力情况，因此，必须对所得到的应力再次进行处理，得到更加合理的应力场，并得到所需点上的应力值。这里介绍两种简单的方法，一种方法称为节点平均法，即把环绕某一节点的各单元的应力加以平均作为该节点的应力值。例如图中节点3的应力为3.1.8算例为了使通过这样平均得来的应力比较接近实际情况，要求环绕节点的各单元尺寸不应相差太大。这种做法，对内节点比较好，对边界点则可能很差。因此，边界节点处的应力不宜直接由单元应力平均来获得，而要根据内节点的应力构造插值函数推算出来。例如图中边界点1的应力，可以先用节点平均法求得节点2、3、4处的应力，在构造相应的插值函数推算边界点1的应力，如常用的抛物线插值公式如下3.1.8算例另一种方法称为单元平均法，即把两相邻单元的应力加以平均，用以表示公共边界中点的应力。为了由这样平均所得到的应力具有较好的精度，两相邻单元的面积不应相差太大。如图中单元②和③边界的中点处的应力为3.1.8算例在不同的有限元软件中均具有各自的后处理方法，但无论后期怎样处理，应力场来源于应力的计算，应力的精度主要依赖于单元的尺寸、单元的类型（位移模式）。对于有限元这种数值计算方法，一般总是希望随着网格的逐步细分所得到的解能够收敛于问题的精确解。根据前面的分析，可知在有限元分析中，一旦确定了单元的形状之后，位移模式的选择将是非常关键的。由于载荷的移置、应力矩阵和刚度矩阵的建立等等，都依赖于单元的位移模式，所以，如果所选择的位移模式与真实的位移分布有很大的差别，那么就很难获得良好的数值解。为了保证解答的收敛性，要求位移模式必须满足以下三个条件，即1）位移模式必须包含单元的刚体位移。2）位移模式必须包含单元的常应变。3）位移模式在单元内要连续、且在相邻单元之间的位移必须协调。3.1.9收敛准则3.1.9收敛准则

1）位移模式必须包含单元的刚体位移。也就是说，当节点位移是由某个刚体位移所引起时，弹性体内将不会产生应变。所以，位移模式不但要具有描述单元本身形变的能力，而且还要具有描述由于其它单元形变而通过节点位移引起单元刚体位移的能力。2）位移模式必须包含单元的常应变。每个单元的应变一般都是包含着两个部分：一部分是与该单元中各点的坐标位置有关的应变（即所谓各点的变应变）；另一部分是与位置坐标无关的应变（即所谓的常应变）。从物理意义上看，当单元尺寸无限缩小时，每个单元中的应变应该趋于常量。因此，在位移模式中必须包含有这些常应变，否则就不可能使数值解收敛于正确解。3.1.9收敛准则3）位移模式在单元内要连续、且在相邻单元之间的位移必须协调。当选择多项式来构成位移模式时，单元内的连续性要求总是得到满足的，单元间的位移协调性，就是要求单元之间既不会出现开裂也不会出现重叠的现象。通常，当单元交界面上的位移取决于该交界面上节点的位移时，就可以保证位移的协调性。在有限单元法中，把能够满足条件1和2的单元，称为完备单元，完备单元是收敛的必要条件。满足条件3的单元，叫做协调单元或保续单元。前面讨论过的三节点三角形单元，均能同时满足上述三个条件，因此属于完备的协调单元，完备的协调单元是收敛的充分条件。在某些梁、板及壳体分析中，要使单元满足条件3比较困难，所以实践中有时也出现一些只满足条件1和2的单元，其收敛性往往也能够令人满意。特别是放松条件3的单元，即完备而不协调的单元，这时要保证其收敛必须通过分片试验的测试，关于分片试验的内容请参见相关书籍。目前，完备而不协调的单元，已获得了很多成功应用。3.1.9收敛准则在有限元分析中，将实际连续体分成许多单元体的组合后，根据线性或非线性位移的假定，人为地选择一个位移场，通过这些措施所得到的模型比实际连续体的刚性要高。因而，近似模型的刚度是实际连续体刚度的上界。若选择不协调单元，那么这种模型可能由于单元分离、叠加或单元之间形成铰而降低刚性，这种影响就有可能使得不协调元比应用协调元所得的结果要好。不过，应用不协调单元事先不能肯定所得的刚度是真实刚度的上界。换句话说，不协调元不一定象协调单元那样刚硬，可能比较柔软，因此有可能会比协调单元收敛得快。3.1.9收敛准则经验证明，根据巴斯卡（Pascal）三角形来选择二维多项式的各项。在二维多项式中，如果包含有对称轴一边的某一项，那么就必须同时包含有另一边的对称项。选择多项式位移模式时，还应该要考虑的一个因素是，多项式中的项数必须等于或稍大于单元边界上的外节点的自由度数。通常是取项数与单元的外节点的自由度数相等，取过多的项数是不恰当的。从节点平衡得到结构整体有限元方程式的方法物理概念清晰，但这种方法在数学上不够严谨。本节从势能原理导出有限元方程，该方法是一种更加严谨、适用性更加广泛的方法。由第2章可知，系统的总势能为如果有集中力作用等式右边还应加一项集中力所作的功，如下3.1.10从能量原理推导刚度矩阵当把连续体离散为有限个单元，并在各单元内假设了位移模式，在单元内应力与应变用节点位移（或单元中）来描述，、，在连续体上的整体积分等于在每个单元上分别积分再求和，如下由势能原理可知，对于处于平衡状态的弹性体，位移真实解应总使总势能取极小值，反之依然，则由，得3.1.10从能量原理推导刚度矩阵上式进一步可写为如果有集中力作用，有上式进一步可写为上式同式（3.55），即为整体结构有限元方程。

3.1.10从能量原理推导刚度矩阵3.2矩形单元虽然三角形单元具有很好的“适应性”，几乎任何复杂边界的弹性体总可以划分为三角形，并且三角形单元计算公式简单，但精度较低。三角形单元间虽然能够保证位移连续，但应力的精度较差，不能很好的反映弹性体内应力的准确分布规律。为了提高计算精度，准确反映弹性体内的应力状态，可以采用一些较精密的单元类型。本节将介绍常用的矩形单元，它采用了比常应变三角形单元次数更高的位移模式，因而可以更好地反映弹性体中的位移状态和应力状态。另外，对一些边界比较规则且呈直线的平面结构的分析，采用矩形单元较合适。这时单元总数可以减少，相应的原始数据准备工作和单元特征计算工作均可节省。3.2.1位移函数如图所示的矩形单元，不失一般性，令矩形单元的长、宽分别为2a、2b。矩形单元有4个节点，共8个自由度，即共有8个节点位移，采用类似三角形单元的分析方法，同样可以完成对矩形单元的力学特性分析。这里引入一个局部坐标系、，这样可以推出比较简洁的结果。如图所示，取矩形单元的形心o为局部坐标系的原点，和轴分别与整体坐标轴x和y平行，两坐标系存在有以下的坐标变换关系式中：、——矩形形心处坐标。矩形形心处坐标以及矩形长、宽可由下式计算3.2.1位移函数在局部坐标系中，节点i的坐标是，其值分别为±1。如节点1在局部坐标系下的坐标为（－1，－1）。由于矩形有4个节点，共8个自由度，可以选择有8个待定参数的位移模式，如下采用类似三角形的分析方法，将各节点的局部坐标代入上式，便可以求得相应的待定参数，然后将带回式位移模式，可得到用节点位移表示的位移模式，如下3.2.1位移函数式中：——矩形单元的形函数，i＝1，2，3，4；——形函数矩阵；——单元节点位移列阵，，i＝1，2，3，4。3.2.1位移函数（i＝1，2，3，4）形函数的表达式为3.2.1位移函数引入符号，，i＝1，2，3，4，则上式可以统一写为可以看出，矩形单元的形函数具有和三角形单元形函数同样的性质，即：形函数在各单元节点上的值，具有“本点是1、它点为零”的性质；在单元内任意点上，四个形函数之和等于1；单元任意一条边上的形函数，仅与该边的两端节点坐标有关。3.2.1位移函数有了单元的位移模式，就可以利用平面问题的几何方程求出单元内任意点的应变，将位移代入几何方程，得式中的应变转换矩阵的子块（i＝1，2，3，4）为3.2.2应变与应力矩阵求得应变之后，再将应变代入物理方程，便可推导出以节点位移表示的应力，如下式中，应力矩阵为其子块（i＝1，2，3，4）为3.2.2应变与应力矩阵由前面的讨论可以发现，四边形单元的位移模式，式（3.77）比常应变三角形单元所采用的线性位移模式增添了项（即相当于xy项），把这种位移模式称为双线性模式。在这种模式下，单元内的应变分量将不再是常量，这一点可以从的表达式中看出。另外四边形单元的位移模式中的与三角形单元相同，它反映了刚体位移和常应变，而且在单元的边界上(=±1或

=±1)，位移是按线性变化的，显然在两个相邻单元的公共边界上，其位移是连续的。3.2.2应变与应力矩阵由单元的应力矩阵表达式还可以看出，矩形单元中的应力分量也都不是常量。正应力、和剪应力均沿、两个方向线性变化，即沿x、y两个方向线性变化。正因为如此，若在弹性体中采用相同数目的节点时，矩形单元的精度要比常应变三角形单元的精度高。但是，矩形单元也有一些明显的缺点，矩形单元不能适应斜交的边界和曲线边界，不便于对不同部位采用不同大小的单元，以便提高有限元分析计算的效率和精度。3.2.2应变与应力矩阵矩形单元刚度矩阵的推导过程与三节点三角形单元类似，同样可以表示为式（3.30）的形式，即，由前文可知的推导过程与形函数的具体表达形式、节点个数均无关，该表达式具有普遍意义。若单元厚度t

为常量，则可以进一步表示为将单元刚度矩阵写成子块的形式，如下3.2.3单元刚度矩阵上式中每一个子块矩阵均为2行×2列，单元刚度矩阵中的子块矩阵的表达式为

（r、s＝1，2，3，4）将应变转换矩阵子块和弹性矩阵，代入上式，得3.2.3单元刚度矩阵式中：（r、s＝1，2，3，4）3.2.3单元刚度矩阵例3.6如图3.22所示，该模型中有两个四边形单元，弹性模量为＝210GPa，厚度为＝0.025m，泊松比＝0.3，＝1kN，求单元所受应力。3.2.4算例解：单元①所对应的节点为2、4、3、1，单元②所对应的节点为4、6、5、3。在有限元分析过程中，首先求解单元①和②的刚度矩阵；然后组装整体刚度矩阵，组装的过程同三角形单元，此处不再给出；最后引入边界条件（，＝0）并结合受力情况，求得整体节点位移列阵＝10-5×{0，0，0，0，0.1162，-0.1674，-0.1149，-0.1628，0.1514，-0.4707，-0.1568，-0.4978}T。为了求解单元应力，需要先求得每个单元所对应的节点位移列阵，结合单元的节点编号，可从整体位移列阵中提取单元①和单元②的位移列阵，如下＝10-5×{0，0，-0.1149，-0.1628，0.1162，-0.1674，0，0}T＝10-5×{-0.1149，-0.1628，-0.1568，-0.4978，0.1514，-0.4707，0.1162，-0.1674}T3.2.4算例注意：整体节点位移列阵是按照节点编号由小到大排列的，而单元位移列阵是按照单元节点编号排列的，如单元①所对应的节点为2、4、3、1，则单元①的位移列阵中的前两个数则表示节点2的x和y方向位移。将单元①和单元②的位移列阵代入应力矩阵，可求得单元①和单元②的应力，如下3.2.4算例通过以上两式便可以求得单元①和单元②内任意点的应力，以单元①为例，若取，，则表示单元①的13边中点处的应力；若取，，则表示24边中点处的应力。若计算单元形心处的应力，则取，为通过分析结果可知，单元内任意点的应力是坐标的函数，若在弹性体中采用相同数目的节点时，矩形单元的精度显然要高于常应变三角形单元的精度。3.2.4算例3.3平面等参元前一节讨论的矩形单元具有精度较高、形状规则、便于计算自动化等优点，对于结构复杂的曲边外形，只能通过减小单元尺寸，增加单元数量进行逐渐逼近。这样，自由度的数目随之增加，计算时间长，工作量大。另外，这些单元的位移模式是线性模式，是实际位移模式的最低级逼近形式，问题的求解精度受到限制。为了克服以上缺点，人们试图找出这样一种单元：一方面，单元能很好地适应曲线边界和曲面边界，准确地模拟结构形状；另一方面，这种单元要具有较高次的位移模式，能更好地反映结构的复杂应力分布情况，即使单元网格划分比较稀疏，也可以得到比较好的计算精度。等参数单元（等参元）就具备了以上两条优点，等参元为任意四边形单元，其网格划分不受边界形状的限制，单元大小可以不相等，是一种精度高而且应用广泛的单元。3.3平面等参元然而直接对任意四边形进行单元分析是困难的，这是由于它的几何形状不规则，没有统一的形状，对各个单元逐个按不同的公式计算，因其工作量大而难以进行。为了解决这一问题，人们采用坐标变换的方法，如图所示，把坐标系ξη（局部坐标系）内形状规则单元（母单元）变换为另一坐标系xy（整体坐标系）中的任意四边形单元（子单元）。然后导出关于局部坐标系形状规整的单元（母单元）的高阶位移模式的形函数，然后利用形函数进行坐标变换，得到关于整体坐标系的复杂形状的单元（子单元）。a）母单元b）子单元3.3平面等参元从图形变换的角度看，坐标系ξη和坐标系xy可以分别看成是母单元和子单元这两个不同单元的坐标系，它们都是直角坐标系。而从另一角度看，坐标系ξη和坐标系xy又可以看成是同一单元（子单元）的两种不同的坐标系。坐标系xy是子单元的直角坐标系，而坐标系ξη可看成是子单元的曲线坐标系。可以看出坐标系xy始终扮演同一角色，即子单元的直角坐标；而坐标系ξη则扮演两种角色，它既是母单元的直角坐标，又是子单元的曲线坐标。在有限元分析中，两者的作用是不同的。直角坐标系xy在整个结构的所有子单元中共同采用，所以称为整体坐标系。而曲线坐标系ξη则只适用于单个独立的子单元，所以称为局部坐标。整体坐标在整体分析中采用，局部坐标则在单元分析中采用。平面问题中常用的等参元有四节点四边形单元、八节点曲边四边形单元和6~8可变节点曲边四边形单元等，本节以八节点曲边四边形等参元为例介绍平面等参元的分析过程。八节点曲边四边形等参元的母单元是二维二次单元。八个节点分别为正方形的四个角点和四个边中点，母单元（边长为2的正方形）采用直角坐标系ξη。由于母单元为16个自由度，因此单元的位移模式为3.3.1等参元刚度矩阵式中：式中的16个常数用8个节点的位移（，）表示后，则上式为3.3.1等参元刚度矩阵采用坐标变换可使母单元的八个节点坐标与等参元的八个节点坐标建立一一对应的关系，整体坐标和局部坐标的变换式为式中：——母单元的形函数。可知，母单元上任意直边，例如2－3边，有ξ＝1，则有关的形函数N2、N3和N6均是η的二次函数，其余的形函数均为零。这样变换就成为二次非线性变换，它可将母单元的直边2－3映射成子单元的曲边2－3。根据等参元的思想，利用上述的形函数，可得等参元的位移模式为3.3.1等参元刚度矩阵有了单元的位移模式，就可以利用平面问题的几何方程求出单元内任意点的应变，将上式代入式（2.12），得等参元的应变矩阵为式中：——i节点位移，（i=1，2，…，8）；3.3.1等参元刚度矩阵——单元几何矩阵，其中为（i=1，2，…，8）在上式中，是ξ、η的函数，因此必须用坐标变换式来转换导数关系，根据复合函数的求导法则，有（i=1，2，…，8）3.3.1等参元刚度矩阵式中的J称为雅可比（Jacobi）矩阵，为则3.3.1等参元刚度矩阵雅可比（Jacobi）矩阵的逆矩阵为将单元应变代入平面问题的物理方程式，就得到平面八节点等参元的应力列阵，为（i=1，2，…，8）3.3.1等参元刚度矩阵等参元刚度矩阵仍可用式（3.30）计算，即在上式中，应变矩阵为局部坐标系ξη下的显式，因此在计算上式时应将微面积dxdy变换成dξ和dη所包围的面积dA，由于所以，有3.3.1等参元刚度矩阵则等参元刚度矩阵为应该指出，上式是对ξ和η的重积分，尽管其积分区域十分简单，但其被积函数却比较复杂，一般该积分没有显式，需要采用数值积分法求解。至于等参元的等效节点力的计算一般也要采用数值积分的方法求解，这里不再讨论，可参考相关书籍。由以上分析可知，在划分单元时，只需确定单元节点的整体坐标值，而不必画出等参元的具体形状，因为在计算中实际使用的只有单元八个节点在整体坐标系下的位移值。等参元变换的条件为，因此在有限元网格划分时，要特别注意这一点。3.3.1等参元刚度矩阵如图所示，在图b中3、4点退化为一个点，在该点，同理在图c中2、3点退化为一个点，在该点，在图d在点1、2、3处>0，在点4处<0，又因为在单元内连续变换，所以单元内肯定存在＝0，这是由于单元过分歪曲造成的。a）正常b）不正常c）不正常