1、在线教务辅导网：在线教务辅导网：http:/教材其余课件及动画素材请查阅在线教务辅导网教材其余课件及动画素材请查阅在线教务辅导网QQ:349134187 或者直接输入下面地址：或者直接输入下面地址：http:/本节主要研究三节点三角形单元，首先建立以单元节点位移表示单元内各点位移的关系式。设单元e的节点编号为i、j、m，如图所示。由弹性力学平面问题可知，每个节点在其单元平面内的位移有两个分量，所以整个三角形单元将有六个节点位移分量，即六个自由度。用列阵可表示为TTeTTTijmiijjmmuvuvuv 式中的子向量为Tiiiuv （i、j、m 轮换）式中：iuiv节点i在x轴方向的位移分量；节
2、点i在y轴方向的位移分量。在有限单元法中，虽然是用离散化模型来代替原来的连续体，但每一个单元体仍是一个弹性体，所以在其内部依然是符合弹性力学基本假设的，弹性力学的基本方程在每个单元内部同样适用。从弹性力学平面问题的解析解法中可知，如果弹性体内的位移分量函数已知，则应变分量和应力分量也就确定了。但是，如果只知道弹性体中某几个点的位移分量的值，那么就不能直接求得应变分量和应力分量。因此，在进行有限元分析时，必须先假定一个位移模式。由于在弹性体内，各点的位移变化情况非常复杂，很难在整个弹性体内选取一个恰当的位移函数来表示位移的复杂变化，但是如果将整个区域分割成许多小单元，那么在每个单元的局部范围内就
3、可以采用比较简单的函数来近似地表示单元的真实位移，将各单元的位移模式连接起来，便可近似地表示整个区域的真实位移函数。这种化繁为简、联合局部逼近整体的思想，正是有限单元法的绝妙之处。基于上述思想，可以选择一个单元位移模式，单元内各点的位移可按此位移模式由单元节点位移通过插值而获得。假设三节点三角形单元内任意点的位移是x、y的线性函数，如下123456uxyvxy 是待定常数，因三角形单元共有六个自由度，且位移函数 在三个节点处的数值应该等于这些点处的位移分量的数值。假设节点i、j、m的坐标分别为（16uv,iix y,jjx y,m mx y式中、）、（）、（），代入位移函数，得12312312
4、3iiijjjmmmuxyuxyuxy从上式中可以解得12311111 ,1 ,122211iiiiiiijjjjjjjmmmmmmmuxyuyxuuxyuyxuuxyuyxu式中：三角形i、j、m的面积。为保证求得的面积为正值，节点i、j、m的编排次序必须是逆时针方向，如图所示。的计算公式，如下1211iijjmmxyxyxy 将式所求得得系数代入位移函数，得12iiiijjjjmmmmuabxc y uab xc y uab xc y u式中：1111jjijmmjmmjijmmjijmmxyax yx yxyybyyyxcxxx （i、j、m 轮换）同理可以求出单元内任意点的y方向位移，
5、如下12iiiijjjjmmmmvabxc y vab xc y vab xc y v令12iiiiNabxc y （i、j、m 轮换）则单元内任意点的位移可以表示成如下形式000000iijijmeijmjmmuvuNNNuNNNvvuv u式中：eu单元内位移矢量位移函数可以进一步写成如下形式ieeeeeeijmjm uNNNN式中：000000ijmeeeeijmijmNNNNNNNNNN1001eiiiNNNI=（i、j、m 轮换）是坐标的函数，它们反映了单元的位移状态，所以一般称之为形状函数，简称形函数。显然，形函数(,)ex yN(,)eix yN决定了单元内的“位移模式”，反映了
6、i节点位移对单元内任意点位移的贡献。有了单元的位移模式，就可以利用平面问题的几何方程求出单元内任意点的应变，由物理方程得00010002iiiijmjeeeeeeijmijmjjiijjmmmmmuvbbbucccvcbcbcbuv BBBBe eB式中：单元应变列阵；单元应变转换矩阵，或称单元几何矩阵。单元应变转换矩阵eeeeijmBBBB，为36矩阵，其子矩阵为0102ieiiiibccbB（i、j、m 轮换）可以看出，表示i节点位移对单元应变的贡献率。当单元确定后，也就确定了，此时单元内的应变仅依赖于节点的位移。对于三节点三角形单元，由于面积和 等都是常量，所以单元应变转换矩阵 中的诸元
7、素都是常量，因而单元中各点的应变分量也都是常量，即单元内各点的应变相同，通常称这种单元为常应变单元。eiBeiBibiceB、求得应变之后，再将应变代入物理方程 ，即公式（2.19），便可推导出以节点位移表示的应力，如下 Deee DB式中：e 单元应力矩阵。令eeeeeijmSDBSSS则eeeS式中：eS单元应力转换矩阵。22 11122iieeiiiiiibcEbccbSDB的子矩阵为eS单元应力转换矩阵可以看出，表示i节点位移对单元应力的贡献率。当单元确定后，也就确定了，此时单元内的应力仅依赖于节点的位移。对于三节点三角形单元，由于面积和 等都是常量，所以单元应力转换矩阵 中的诸元素都
8、是常量，因而单元中各点的应力分量也都是常量，即单元内各点的应力相同，通常称这种单元为常应力单元。eiSeiSibiceS、可见，对于三节点三角形单元，由于所选取的位移模式是线性的，因而导致单元内各点的应变、应力相同，但相邻单元将具有不同的应力和应变，即在单元的公共边界上应力和应变的值将会有突变，但位移却是连续的。123456uxyvxy00010002iiiijmjeeeeeeijmijmjjiijjmmmmmuvbbbucccvcbcbcbuv BBBBeee DB 由前文可知，形函数 决定了单元内的“位移模式”，反映了i节点位移对单元内任意点位移的贡献，形函数具有如下性质：(,)eix y
9、N1.形函数在各单元节点上的值，具有“本点是1、它点为零”的性质即1(,)0ikkkiN xyki （i、j、m 轮换）证明如下：在节点j、m上1,02ijjiijijNxyabxc y1,02immiimimNxyabxc y类似有,0 ,1 ,0 ,0 ,0 ,1jiijjjjmmmiimjjmmmNxyNxyNxyNxyNxyNxy在节点i上1,12iiiiiiimNxyabxc y2.在单元内任意点上，三个形函数之和等于1即1ijmNNN证明：将单元内任意点的坐标代入上式，得 ,1212ijmiiijjjmmmijmijmijmNxyNxyNxyab xc yab xc yab xc
10、yaaabbbxcccy ，因此容易得到此结论。由此可见，三个形函数中只有二个是独立的。2ijmaaa 0ijmbbb0ijmccc可以证明 ，例如：假设单元各节点的位移均相同并等于 ，则根据此性质可得单元内任意点的位移为0u0v，0000(,)(,)iijjmmijmi ijjm mijmu x yN uN uN uNNNuuv x yN vN vN vNNNvv 显然，该性质反映了单元刚体位移，即当单元做刚体运动时，单元内任意点的位移均等于刚体位移。3.三角形单元任意一条边上的形函数，仅与该边的两端节点坐标有关。例如在ij边上，有0,1,yxNxxxxyxNxxxxyxNmijijijii
11、事实上，因i j 边的直线方程方程为iimmiiijijyxxcbyxxxxyyy(,)mNx y(,)jNx y代入及，得1,2102mmmmmiimmmimibNx yab xcxxycab xc y1,21212mjjjjiimmjjjijijiimjmmjimbNx yab xcxxycb cab xc ybxxxxcb cb cxxc则,ijjixxNx yxx另外，由1ijmNNN可以求得,11iijmjixxNx yNNxx 利用形函数的这一性质可以证明，相邻单元的位移分别进行线性插值之后，在其公共边上将是连续的。例如，对图3.2所示的单元具有ij公共边，可知，在ij边上有0,0
12、,yxNyxNnm这样，不论按哪个单元来计算，公共边ij上的位移均由下式表示iijji ijjuNuN uvN vN v由此可见，在公共边上的位移u、v将完全由公共边上的两个节点i、j的位移所确定，因而相邻单元的位移是连续的。例例3.1 求如图3.3所示等边三角形的形函数。23002ijmmjjmiimmijjiax yx yaax yx yax yx y33022ijmjmimijbyyabyyabyy 22ijmjmimijaacxxcxxcxxa 三角形的面积234a 将以上数据代入式（3.10），得22123322223iiiiaaNabxc yaxya21232223jjjjaaNa
13、b xc yxya21223mmmmNab xc yaya 将各节点坐标代入以上三式，有3(0,0)1(,0)0(,)022iiiaNN aNa3(0,0)0(,0)1(,)022jjjaNNaNa3(0,0)0(,0)0(,)122mmmaNNaNa即满足“本点是1、它点为零”的性质。即满足“在单元内任意点上，三个形函数之和等于1”的性质。将三个形函数相加，得1ijmNNN ，即在ij边上任意位置的位移，将完全由节点i、j的位移所确定。取其它边同样可以得出其结论，即满足“三角形单元任意一条边上的形函数，仅与该边的两端节点坐标有关”的性质。0y 0mN 以ij边为例，在ij边上，则 单元在节点
14、处受力，单元会发生变形，也就是说单元在节点处受到的力与单元节点位移之间有必然的联系。单元间正是通过与节点间的相互作用力连接起来成为整体，而一个单元仍是一个弹性体，如果将每一个单元的受力位移关系找到，则整体的受力位移关系也容易清楚。为了推导单元的节点力和节点位移之间的关系，可应用虚位移原理对单元进行分析。单元在节点处受到的力称为单元节点力，单元在节点力的作用下处于平衡，节点力可以表示为TeixiyjxjymxmyFFFFFFF 式中：ixF节点i沿x方向的节点力分量；iyF节点i沿y方向的节点力分量。单元节点的虚位移为相应的单元内的虚位移场为*Teiijjmmuvuvuv*eeeuN 单元内的虚
15、应变为*eee B 设单元仅在节点上受力，则单元节点力在节点虚位移上的虚功为*1TeeW F 单元内的应力在虚应变上所做的功为*2Teee TeTeeWdVdVB DB 式中：dV在体积上积分。eeTeedVFB DB 令单元刚度矩阵为由虚功原理可知，式（3.27）和（3.28）相等，再由虚位移 的任意性，可得*e eeTedVKB DB上式为单元刚度矩阵 的定义式，其推导过程与形函数的具体表达形式、节点个数均无关，这说明该表达式具有普遍意义，则 eKeeeFK 上式是单元节点力与单元节点位移之间的关系式，即单元节点平衡方程，或称单元刚度方程，具有普遍意义。对于三节点三角形单元，由于 、中的元
16、素都是常量，这里假设单元厚度均匀，则BDeeTetKBDBt式中：单元面积；单元厚度。上式进一步可以表示为TeeeiiiijimeTeeejijmjijjjmTeeemmimjmmBBBt BKKKKBDKKKBKKK将几何矩阵式（3.16）和弹性矩阵式（2.19）代入，可得每个子块的矩阵表达式，如下,21122114 122rx sxrx syeeTersrsry sxry syrsrsrsrsrsrsrsrsKKtKKb bc cb cc bE tc bb cc cb b KBDB,r si j m单元刚度矩阵具有如下特点：1.中的每一个元素都是一个刚度系数eK假设单元的节点位移为 ，由
17、，得到节点力为 100000Te eeeFK,ixix ixiyiy ixjxjx ixejyjy ixmxmx ixmymy ixFKFKFKFKFKFKF 由此可得：表示i节点在水平方向产生单位位移时，在节点i的水平方向上需要施加的节点力；表示i节点在水平方向产生单位位移时，在节点i的垂直方向上需要施加的节点力。选择不同的单元节点位移，可以得到单元刚度矩阵中每个元素的物理含义，如下：表示s节点在垂直方向产生单位位移时，在节点r的垂直方向上需要施加的节点。表示s节点在水平方向产生单位位移时，在节点r的垂直方向上需要施加的节点力；表示s节点在垂直方向产生单位位移时，在节点r的水平方向上需要施加
18、的节点力；表示s节点在水平方向产生单位位移时，在节点r的水平方向上需要施加的节点力；ixixK,ixiyK,sxrxK,sxryK,syrxK,syryK,eKeK因此，中的每一个元素都是一个刚度系数，表示单位节点位移分量所引起的节点力分量，中的每一个子块表示单元某个节点位移矢量对单元某个节点力矢量的贡献率。2.是对称矩阵eK 由式（3.32）可知，是对称矩阵。根据对称性，可以减少计算量和存储量。eK3.是奇异矩阵，即eK0eK 根据这一性质，当已知单元节点位移时，可以从式（3.31）中求出单元节点力。反之，由于单元刚阵奇异，不存在逆矩阵，因此，当已知单元节点力时不能求出单元上的节点位移。从物
19、体变形的实际情况来说，单元刚阵奇异是必须的，因为单元除了产生变形外，还会产生刚体位移，即仅仅依靠节点力是无法唯一地确定刚体位移的。事实上，当单元的节点力为零时，单元仍可做刚体运动。例如，假定单元产生了x方向的单位刚体位即并假设此时对应的单元节点力为零，则由 1 0 1 0 1 0Te eeeFK 得010001000100e K 可以得到，在单元刚度矩阵中1，3，5列中对应行的系数相加为零，由行列式的性质可知，。0eK 4.单元的刚度不随单元或坐标轴的平行移动而改变 单元的刚度取决于单元的形状、大小、方向和弹性系数，而与单元的位置无关，即不随单元或坐标轴的平行移动而改变。因此，只要单元的形状、
20、大小、方向和弹性系数相同，无论单元出现在任何位置均有相同的单元刚度矩阵，根据对称性，可以减少计算量。，厚度为t。0例例3.2 3.2 求图3.4所示等腰直角三角形单元的刚度矩阵，设 解解 根据1（0，a）、2（0，0）、3（a，0）的坐标和式（3.8）可得12313223121331232100byycxxabyyacxxabyyacxx 单元面积 ，将以上各数据代入几何矩阵式（3.16）为 22a 20000001010110000010100000001101eaaaaaaaaaB由于 ，则平面应力与平面应变问题的弹性矩阵0D相同，为21010010010(1)1000.5002EED00
21、20200202002101101eeEaSDB1011010202001031211213014002020101101eeTeEtt KBDB由式（3.18）可得应力转换矩阵，如下则单元刚度矩阵为 从以上求得的结果可以看出，单元刚度矩阵是对称的。1、3、5列相加等于零，2、4、6列相加也等于零，单元刚度矩阵是奇异的。另外在几何矩阵中所涉及到的元素均为“坐标差量”，而与坐标的具体值无关，从而也说明了单元刚度矩阵不随单元或坐标轴的平行移动而改变。1）获取单元节点信息；2）计算单元面积 及参数 、等；3）计算几何矩阵 中各元素；4）根据问题类型（平面应力、平面应变）计算弹性矩阵 ；5）计算应力转
22、换矩阵 ；6）计算单元刚度矩阵 。ibiceBDeeSDBeeTetKBDB 从单元刚度矩阵的计算公式 可以看出，在求解单元刚度矩阵时，可以不必求得应力转换矩阵 ，但在后期求解单元应力时要用到该矩阵，所以在此处一并求出，也可以根据实际情况具体考虑。eeTetKBDBeeSDB 实际上作用在物体上的外力可以直接作用在节点上，也可以不作用在节点上，然而在有限元法中要求作用在物体上的各种外力必须用作用在节点上的力表示。这一点体现了有限法中“离散化”这一概念，即将连续体离散成只有在节点处相连的单元体，单元与单元之间的联系只能通过节点。如前文所述，单元内任意点的位移、应力、应变等变量最终都用单元节点位移
23、来表示。同样，作用在物体上的各种外力也必须用作用在节点上的力表示，这一过程称为外力的静力等效移置，所得到的节点力称为等效节点力。静力等效原则：对于刚体来说，所谓静力等效原则就是单元上原有的外力系和将外力系向各节点移置所得的等效节点力，二者向同一点简化应具有相同的主失和主矩；对于弹性体来说，所谓静力等效原则就是指单元上的外力系和将该力系向各节点移置后的等效节点力，二者在虚位移上的虚功相等，也即外力作用在单元上所引起的变形能和移置后的等效节点力在单元上引起的变形能相等，在一定的位移模式下这种移置是唯一的。1.弹性体静力等效原则虚功原理 虚功等效的思路是：就一个单元来说，把作用在单元上的外力移置到节
24、点上后，应当与原来的实际外力所作虚功等效。其计算方法是：对任意允许的微小虚位移，原外力所作虚功等于移置后的等效节点力 所作虚功。即eFTeee Te Te TccdSdV FuFuquG*e*euecu式中：节点虚位移；单元虚位移场；集中力作用处产生的虚位移。等式左边表示等效节点力在节点虚位移上所作虚功；右边第一项表示作用在单元任意点c处的集中力 所作虚功、第二项表示表面力 所作的虚功、第三项表示体积力 所作的虚功。cFqG将 代入上式，并考虑 的任意性，得*eeeuN e eeTeTeTeeecccqGdSdVFNFNqNGFFFecFeqFeGFecN式中：由集中力等效的单元等效节点力；由
25、表面力等效的单元等效节点力；由体积力等效的单元等效节点力；形函数在集中力作用处的值。ecFeqFeGF单元等效节点力 、和 的表达式分别为eeTcccFNFeeTqdSFNqeeTGdVFNG 以上三式不但适用三角形单元，同样适用其它单元，具有普遍意义。下面根据以上三式，分别进行讨论。000000iicxjeeTccccyjmmcNNFNFNNN FNF 可以看出，当外力作用在节点上时，如作用在i节点上，此时 ，则计算所得到的等效节点力与原外力在数值上是相同的。1iN 0jmNN（1）作用在单元任意点的集中力TccxcyFFF（2）作用在单元上的表面力 现举例说明，如图所示的单元e，在ij边上
26、作用有表面力 ，假设ij边的长度为l，单元厚度t为常数，其上任一点P距节点i的距离为s。q有 ，则1iNs l jNs l0mN 0010lileeTqjmstdslN tdssdSNtdstdslNtdsqqFNqqqq例例3.3 如图所示，设31边的长度为l，单元厚度为t，其上作用沿x方向的均布荷载 ，求其等效节点荷载。q解解：在31边上，20N 0Tqq=，则等效节点荷载为123012300000001000102TleeTqTNNNqdStdsNNNqlt FNq 从计算结果来看，等效后的节点力是将原来作用在31边上的均布面力平分到1、3节点上。000000iixjeGyjmmNNGN
27、tdxdyGNNN F（3）作用在单元上的体积力 设单元所受体积力为 ，单元厚度t为常数，如所示，则TxyGGG=当单元体是均质、等厚、比重为 时，则 、，有0 xG yG 00010000130001iijeGjmmNNNtdxdytNNN F即当单元有均匀自重时，在三节点三角形单元中每个节点沿y方向的等效节点力等于单元总重量的1/3。2.刚体静力等效原则力系简化（1）作用在单元任意点的集中力 若集中力 作用在节点上，则可直接将集中力认为是等效节点力即可。若集中力作用于单元的某一边界上，如图所示，则先将 分解为 、，然后按线段的比例把 分别移置到i，j两点上。即cFcFcxFcyFcxFcy
28、F221100TeccxcycxcyllllFFFFllllF（2）作用在单元上的表面力 如图所示，已知在ij边受有均布面力 ，单元厚度t为常数，则移置到i、j节点上的等效节点力为q0 1 0 1 002Teqqlt F如图所示，当某一边上有三角形分布的面力时，则等效节点力为0120000233Teqq lt F （3）作用在单元上的体积力 如图所示，如果三角形单元ijm的重心c上受有体积力 ，单元厚度t为常数，则按刚体静力等效原理可把体积力直接移置到i，j，m三个节点上，即t000333TeGtttF 根据圣维南原理可知，按刚体静力等效原则移置后的等效节点力所引起的应力误差只是局部的不会影响
29、到整体，即物体内部距实际外力作用处相当远的各点处应力，与其外力的具体分布情况关系很小。在线性位移模式下，弹性体和刚体按各自的静力等效原则移置外力的结果是一致的。整体节点载荷列阵：由各节点所受等效节点力按节点号码以从小到大的顺序排列组成的列阵。等效节点力是由集中力、表面力和体积力共同移置构成的，其中集中力包括直接作用在弹性体上的外力和边界约束力，如支座反力。前文对单元体进行了分折，得到了单元刚度方程 ，但要解决问题，还必须进一步建立整个计算模型的整体刚度方程。完成这一步的关键，在于怎样将单元的刚度矩阵和节点荷载列阵，分别“组装”成整体刚度矩阵和整体节点荷载列阵。这里通过研究任意节点的平衡来建立整
30、体刚度矩阵，该方法不但比较直观、易懂，而且对怎样编写计算机程序是很有帮助的。为了研究整体刚度矩阵的组装过程，先引入两个概念。eeeFK 整体节点位移列阵：由各节点位移按节点号码以从小到大的顺序排列组成的列阵。式中：，12345TTTTTT Tiiiuv(1,2,5)i 不失一般性，仅考虑计算模型中有4个单元，如图所示。四个单元的整体节点位移列阵为 对 每 个 单 元 都 可 以 写 出 相 应 的 单 元 刚 度 方程 ，即单元节点平衡方程。例如，对号单元，有 eeeFK(1)(1)(1)(1)11112131(1)(1)(1)(1)22122232(1)(1)(1)(1)33132333FK
31、KKFKKKFKKK 式中：(1)iF号单元中第i（i=1,2,3）节点所受力。为了便于组装整体刚度矩阵，将上式以整体节点位移 表示，即 (1)(1)(1)(1)11111213(1)(1)(1)(1)22212223(1)(1)(1)(1)(1)3331323345000000000000000000FKKKFKKKKFKKK (1)K号单元的扩大刚度矩阵或称为单元贡献矩阵。同理，对于单元，有(2)(2)(2)(2)111113142(2)(2)(2)(2)(2)33313334(2)(2)(2)(2)444143445000000000000000000FKKKKFKKKFKKK (2)i
32、F(2)K式中：号单元中第i（i=1,3,4）节点所受力；号单元的扩大刚度矩阵。对于单元，有12(4)(4)(4)(4)(4)33334353(4)(4)(4)(4)44344454(4)(4)(4)(4)55354555000000000000000000FKKKKFKKKFKKK (4)iF(4)K式中：号单元中第i（i=3,4,5）节点所受力；号单元的扩大刚度矩阵。对于任意一个节点，可能承受两种力的作用，一种是其它单元给予该节点的反作用力；另一种是作用在节点上的等效节点力。对整体而言，前者属于内力，后者属于外力，每个节点在两种力的作用下处于平衡。将各单元刚度方程左边相加，即将各节点所受力
33、相加，由于对于整体而言，单元给予节点的反作用力属于内力，在相加过程中相互抵消，所以各节点所受力相加的结果只有外力，即等效节点力，从而得到整体节点荷载列阵，如下(1)(2)111(3)(1)222(3)(4)(1)(2)33333(4)(2)444(3)(4)55500000000FFFFFFF=FFFFFFFFFFF 将各单元刚度方程右边相加，从而得到整体刚度矩阵，如下(1)(2)(3)(4)11121314152122232425313233343541424344455152535455(1)(2)(1)(1)(2)(2)111112131314(1)(1)(3)(1)(3)(3)2122
34、22232325(1)3100KKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKKK(2)(1)(3)(1)(2)(3)(4)(2)(4)(3)(4)3132323333333334343535(2)(2)(4)(2)(4)(4)414343444445(3)(3)(4)(4)(3)(4)52535354555500KKKKKKKKKKKKKKKKKKKKKK 通过以上分析得，整体节点载荷与整体节点位移之间的关系式，即结构整体有限元方程，如下FK 式中：K整体刚度矩阵。4()1iiKK整体刚度矩阵组装的基本步骤：1）将单元刚度矩阵中的每个子块放到在整体刚度矩阵
35、中的对应位置上，得到单元的扩大刚度矩阵。注意对于单元刚度矩阵是按照局部编码排列的，即对应单元刚度矩阵中的i、j、m；对于整体刚度矩阵是按照整体编码排列的，即按节点号码以从小到大的顺序排列。在组装过程中，必须知道单元节点的局部编码与该节点在整体结构中的整体编码之间的关系，才能得到单元刚度矩阵中的每个子块在整体刚度矩阵中的位置。将单元刚度矩阵中的每个子块按总体编码顺序重新排列后，可以得到单元的扩大矩阵。例如在图中，单元的局部编码为i、j、m，对应整体编码为1、3、4，然后将单元刚度矩阵中的每个子块按总体编码顺序重新排列后，可以得到单元的扩大矩阵。注意有些书籍中将局部编码表示为1、2、3或1，2，3
36、等；2）将全部单元的扩大矩阵相加得到整体刚度矩阵。通过以上组装过程可以得到组装整体刚度矩阵的一般规则：1）结构中的等效节点力是相关单元结点力的叠加，整体刚度矩阵的子矩阵是相关单元的单元刚度矩阵子矩阵的集成；2）当整体刚度矩阵中的子矩阵 中r=s时，该节点（节点r或s）被哪几个单元所共有，则 就是这几个单元的刚度矩阵中的子矩阵 的相加。如 应该是单元中对应子矩阵的集成，即 rsKrsKersK33K(1)(2)(3)(4)3333333333KKKKK 3）当 中 时，若rs边是组合体的内边，则 就是共用该边的两相邻单元刚度矩阵中的子矩阵 的相加。如13边为单元和的共用边，则 r srsKrsK
37、ersK(1)(2)131313K=KK 4）当 中r和s不同属于任何单元时，则 =0。如节点r1和 s5不同属于任何单元，此时 =0。rsKrsK15K 上述组装基本步骤和规则具有普遍意义，对于不同类型、不同形式的单元，只是相应节点的子矩阵的阶数（节点自由度节点自由度）可能不同，至于组装整体刚度矩阵的规律仍是相同的。正是应为有了这种组装规律，使得有限元法能够很方便地应用电子计算机进行计算。例例3.4 如图所示有限元模型，弹性模量为 ，厚度为 ，为简化计算取 ，求整体刚度矩阵。t0E单元编号整体编码1、2、32、4、55、3、23、5、6局部编码i、j、mi、j、mi、j、mi、j、m以整体编
38、码表示的单元刚度矩阵子块解：解：该模型中共有6个节点，4个单元，各单元的信息如表3.1所示。表3.1 各单元信息(1)(1)(1)111213(1)(1)(1)212223(1)(1)(1)313233KKKKKKKKK(2)(2)(2)222425(2)(2)(2)424445(2)(2)(2)525455KKKKKKKKK(3)(3)(3)555352(3)(3)(3)353332(3)(3)(3)252322KKKKKKKKK(4)(4)(4)333536(4)(4)(4)535556(4)(4)(4)636566KKKKKKKKK同例3.2类似的分析，得(2)(2)(2)222425(
39、2)(2)(2)(2)424445(2)(2)(2)5254551011010202001031211213014002020101101EtKKKKKKKKKK 根据单元刚度矩阵的性质，可知 ，由于几何矩阵 ，所以有 ，整体刚度矩阵 中的各子块是对所有单元相应的子块求和得到的（实际只是对相关单元求和），其中各子块矩阵均为2行2列，整体刚度矩阵用子块矩阵可以表示为(2)(1)(4)KK=K(3)(2)BB(2)(3)KK111213141516212223242526313233343536414243444546515253545556616263646566KKKKKKKKKKKKKKKK
40、KKK=KKKKKKKKKKKKKKKKKK上式中任意一子块矩阵均为2行2列，在计算过程中，无需将每个单元刚度矩阵进行扩大，只需判断整体刚度矩阵子块的下标，然后利用组装整体刚度矩阵的一般规则进行计算，如 ，由图形可知，节点2由单元、和所共有，则 22K(1)(2)(3)22222222KKK+K25K(2)(3)252525KKK，由图形可知，25边为单元和的共用边，则15K150K，由图形可知，节点1、5不同属于任何单元，则采用同样的方法进行计算，得到整体刚度矩阵为1011010000000202000000001061411101001216120210000041610021011012
41、160014000010003121004001200130100000121206121001014111601000000002020000010001101EtK4.是奇异矩阵，在排除刚体位移后，它是正定阵 K3.是稀疏矩阵，非零元素呈带状分布K 用有限元方法分析复杂工程问题时，节点的数目比较多，整体刚度矩阵的阶数通常也是很高的。那么在进行计算时，如果存储整体刚度矩阵的全部元素，将会浪费较大的资源、降低计算效率。如果根据整体刚度矩阵的特点进行编写程序，可以大大节省资源、并提高计算效率。因此有必要了解和掌握整体刚度矩阵的特点，整体刚度矩阵具有以下几个显著的特点：1.是对称矩阵K2.中主对角
42、元素总是正的K1.是对称矩阵K 由单元刚度矩阵的对称性和整体刚度矩阵的组装过程，可知整体刚度矩阵必为对称矩阵。利用对称性，在计算机编写程序时，只存储整体刚度矩阵上三角或下三角部分即可。2.中主对角元素总是正的K 例如，刚度矩阵 中的元素 表示节点2在x方向产生单位位移，而其它位移均为零时，在节点2的x方向上必须施加的力；表示节点2在y方向产生单位位移，而其它位移均为零时，在节点2的y方向上必须施加的力。很显然在此情况下力的方向应该与位移方向一致，故应为正号。K33K44K3.是稀疏矩阵，非零元素呈带状分布K 如果遵守一定的节点编号规则，就可使矩阵的非零元素都集中在主对角线附近呈带状。整体刚度矩
43、阵中的子矩阵 只有当下标s等于r或者s与r同属于一个单元时才不为零，这就说明，在第r双行中非零子矩阵的块数，应该等于节点r周围直接相邻的节点数目加1。可见，中元素一般都不是填满的，是稀疏矩阵，且非零元素呈带状分布。rsKK 以下图所示的单元网格为例，其整体刚度矩阵中的非零子块（每个子块为2行2列）的分布情况如下图所示，图中阴影部分表示该子块不为零，其它子块部位均为零。显然，带状刚度矩阵的带宽取决于单元网格中相邻节点号码的最大差值D。把半个斜带形区域中各行所具有的非零元素的最大个数叫做刚度矩阵的半带宽（包括主对角元），用B表示，如下 B=2(D+1)通常的有限元程序，一般都利用刚度矩阵的对称和稀
44、疏带状的特点，在计算求解中，只存储上半带的元素，即所谓的半带存储。因此，在划分完有限元网格进行节点编号时，要采用合理的编码方式，使同一单元中相邻两节点的号码差尽可能小，以便节省存储空间、提高计算效率。对于同样的有限元单元网格，按照图（a）的结点编码，最大的半带宽为B=2(104)=14；按照图（b）的结点编码，最大的半带宽为B=2(102)18；按照图（c）的结点编码，最大的半带宽为B=2(106)10。(a)(b)(c)4.是奇异矩阵，在排除刚体位移后，它是正定阵 K 无约束的弹性体（或结构物）的整体刚度矩阵是奇异的，不存在逆矩阵，即关于位移的解不唯一。这是因为弹性体在外力的作用下处于平衡，
45、外力的分量应该满足三个静力平衡方程。这反映在整体刚度矩阵中就意味着存在三个线性相关的行或列，所以是个奇异阵，不存在逆矩阵。例如：设弹性体在外力的作用下处于平衡，这时相应的解为 ，然后在给予弹性体以刚体位移而相应的节点位移 ，这时，仍是问题的解，因为刚体位移不会破坏平衡。1 2 12 注：当排除刚体位移后，整体刚度矩阵是正定矩阵。前文已经提到在排除刚体位移后，整体刚度矩阵是正定的，方程才可求得唯一解。排除刚体位移可以通过引入边界约束条件来实现，这里介绍两种比较简单的引入已知位移的方法。1代入法2乘大数法1代入法 该方法保持方程组仍为2n2n系统，仅对整体刚度矩阵 和整体载荷列阵进行修正。下面以一
46、个只有四个方程的简单例子加以说明，方程如下1112131411212223241231323334234142434424KKKKuFKKKKvFKKKKuFKKKKvF 假定系统中节点位移 、，则当引入这些节点的已知位移之后，方程就变成11uc22uc1112131411212223241231323334234142434424KKKKuFKKKKvFKKKKuFKKKKvF11221 123 22224122441 143 242442100000001000cuFK cK cKKvcuFK cK cKKv若 ，则120cc1222412242442410000000010000uKKv
47、FuKKvF 然后，用这组维数不变的方程来求解所有的节点位移。显然，其解答仍为原方程的解答。在手算时，可直接将零位移约束所对应的整体刚度矩阵中的行和列直接划去，使得整体刚阵的维数变小，更便于手算。2乘大数法 将 中与指定的节点位移相对应的主对角元素乘上一个大数，同时将 中的对应元素换成指定的节点位移值、该大数与节点位移相对应的主对角元素三者的乘积。若把此方法用于上面的例子，则方程就变成KF151511112131411112122232421515231323334233241424344410101010uKKKKc KvKKKKFuKKKKc KvKKKKF该方程组的第一个方程为15151
48、1112 11321421111010KuK vK uK vc K解得 ，这种方法就是使 中相应行的修正项远大于非修正项。11uc1代入法2乘大数法 在以上的两种方法中，代入法接近人工解法，虽然该方法比较直观，但该方法对刚度矩阵改变较多，程序效率不高。乘大数法对刚度矩阵改变较少，工作量较小，但相乘的“大数”若取得过大，求解时会发生“溢出”、若取得太小则会引起较大的误差。根据前面的讨论，现以三角形常应变单元为例来说明应用有限元法求解弹性力学平面问题的具体步骤。1）力学模型的确定。根据工程实际情况确定问题的力学模型，并按一定比例绘制结构图、注明尺寸、载荷和约束情况等；2）结构进行离散化。将弹性体划
49、分为许多单元，并对节点进行编号，确定全部节点的坐标值；对单元进行编号，并列出各单元三个节点的节点号；3）计算等效节点载荷，形成整体载荷列阵；4）计算单元刚度矩阵；5）组装整体刚度矩阵；6）处理约束，引入边界条件；7）求解线性方程组，得到节点位移和节点力；8）计算单元应力；9）整理计算结果（后处理部分）。例例3.5 3.5 如图所示为一厚度t=1cm的均质正方形薄板，边长为2m，上下受均匀拉力 =106N/m2，材料弹性模量为E，泊松比 ，不记自重，试用有限元法求其应力分量。q1 3 解：解：1）力学模型的确定。由于此结构长、宽远大于厚度，而载荷作用于板平面内，且沿板厚均匀分布，故可按平面应力问
50、题处理，考虑到结构和载荷的对称性，可取结构的1/4来研究。表3.2 节点坐标值节点坐标1234x0110y00112）结构进行离散化。为了计算简便并能说明整个计算过程，将该1/4结构离散为两个三角形单元。节点编号、单元划分如图所示，各节点的坐标值如表3.2所示，各单元编码如表3.3所示。表3.3 单元编码单元号整体编码1、2、33、4、1局部编码i、j、mi、j、m 3）计算等效节点载荷，形成整体载荷列阵。按照弹性体静力等效原则，应用公式（3.41），并参考例3.3可得单元由均布载荷 等效的节点力为q(2)01010020101002FTqTqqltF6232101 1020.5 10qFql
