1、平面问题的有限单元解法南京农业大学工学院机械工程系2022-11-8南京农业大学工学院机械工程系有限元单元法基本思想n有限单元法的思想是将物体（连续的求解域）离散成有限个且按一定方式相互联结在一起的单元组合，来模拟或逼近原来的物体，从而将一个连续的无限自由度问题简化为离散的有限自由度问题求解的一种数值分析法。物体被离散后，通过对其中各个单元进行单元分析，最终得到对整个物体的分析。n有限单元法的分析步骤如下：q物体离散化q单元特性分析q单元组集，整体分析q求解未知节点的位移q由节点的位移求解各单元的位移和应力2022-11-8南京农业大学工学院机械工程系n基本变量 u （位移）（应变）（应力）n
2、基本方程q力的平衡方程q几何方程q物理方程n求解方法q经典解析q半解析q传统数值解法q现代数值解法（计算机硬件、规范化、标准化、规模化）物体变形及受力情况的描述三大方面三大方程即：=E E 弹性模量2022-11-8南京农业大学工学院机械工程系有限元单元模型中几个重要概念n单元q网格划分中每一个小的块体n节点q确定单元形状、单元之间相互联结的点n节点力q单元上节点处的结构内力n载荷q作用在单元节点上的外力（集中力、分布力）n约束p限制某些节点的某些自由度n弹性模量（杨式模量）En泊松比（横向变形系数）n密度单元单元单元单元载荷节点节点力约束2022-11-8南京农业大学工学院机械工程系平面问题
3、有限单元法基本概念n有限单元法(FEM)是20世纪50年代以来随着计算机的广泛应用而发展起来的一种数值解法。简单地说，就是用结构力学方法求解弹性力学问题。n平面问题的有限单元法求解q将连续体变换成为离散化结构。即将连续体划分为有限多个有限大小的单元，这些单元仅在一些结点连接起来，构成一个所谓离散化结构。（对于平面问题，常用的单元是三角形单元）q用结构力学方法进行求解2022-11-8南京农业大学工学院机械工程系有限元单元法分析步骤（一）n结构离散化q 将结构分成有限个小的单元体，单元与单元、单元与边界之间通过节点连接。结构的离散化是有限元法分析地第一步，关系到计算精度和效率，包括以下三个方面：
4、n单元类型的选择。选定单元类型，确定单元形状、单元节点数、节点自由度数等。n单元划分。网格划分越细，节点越多，计算结果越精确，但计算量越大。网格加密到一定程度后计算精度提高就不明显，对应应力变化平缓区域不必要细分网格。n节点编码。注意：有限元分析的结构已不是原有的物体或结构物，而是由同样材料、众多单元以一定方式连接成的离散物体。所以，用有限元分析计算所获得的结果是近似的（满足工程要求即可）。2022-11-8南京农业大学工学院机械工程系有限元单元法分析步骤（二）n单元特性分析q 选择未知量模式n选择节点位移作为基本未知量时，称为位移法；n选节点力作为基本未知量时，称为力法；n取一部分节点位移和
5、一部分节点力作为未知量，称为混合法。q分析单元力学性质n根据单元材料性质、形状、尺寸、节点数目、位置等，找出单元节点力和节点位移关系式，应用几何方程和物理方程建立力和位移的方程式，从而导出单元刚度矩阵。q计算等效节点力n作用在单元边界上的表面力、体积力或集中力都需要等效地移到节点上去，即用等效力来替代所有作用在单元上的力。2022-11-8南京农业大学工学院机械工程系有限元单元法分析步骤（三）n整体分析q集成整体节点载荷矢量 F。结构离散化后，单元之间通过节点传递力，作用在单元边界上的表面力、体积力或集中力都需要等效地移到节点上去，形成等效节点载荷。将所有节点载荷按照整体节点编码顺序组集成整体
6、节点载荷矢量。q组成整体刚度矩阵K，得到总体平衡方程：q引进边界约束条件，解总体平衡方程求出节点位移。通过上述分析可以看出有限单元法的基本思想是“一分一合”，分是为了进行单元分析，合是为了对整体的结构进行综合分析。KF2022-11-8南京农业大学工学院机械工程系弹性力学中的几个基本概念0limVf FVn作用于物体的外力可以分为体积力和表面力。n体力：分布在物体体积内的力，如重力、惯性力。n为了表明物体在某一点P所受体力的大小和方向，在这一点取物体的一小部分，它包含P点，而它的体积为V，作用于其上的体力为F，则体力的平均集度为F/V。当V不断减小，假定体力为连续分布，则F/V将趋于一定的极限
7、f，即:n这个极限矢量f就是该物体在P点所受体力在集度。f的方向就是F的方向，矢量f在坐标轴x,y,z上的投影fx,fy,fz称为该物体在P点的体力分量，以沿坐标轴正方向为正，沿坐标轴负方向为负。2022-11-8南京农业大学工学院机械工程系弹性力学中的几个基本概念n面力：分布在物体表面上的力，如流体压力和接触力。n为了表明物体在某一点P所受面力的大小和方向，在这一点取物体表面的一小部分，它包含P点，而它的面积为S，作用于其上的面力为F，则面力的平均集度为F/S。当S不断减小，假定体力为连续分布，则F/S将趋于一定的极限 ，即:0limSf FSfn这个极限矢量 就是该物体在P点所受面力在集度
8、。的方向就是F的方向，矢量 在坐标轴x,y,z上的投影 称为该物体在P点的面力分量，以沿坐标轴正方向为正，沿坐标轴负方向为负。fff,xyzfff2022-11-8南京农业大学工学院机械工程系弹性力学中应力的方向规定n每一个面上的应力可以分解为一个正应力和两个切应力。n正应力用表示，加上一个下标字母，表示作用面和作用方向。n切应力用表示，并加上两个下标字母，表示作用面和作用方向。前一个字母表示作用面垂直于哪一个坐标轴，后一个字母表示作用方向沿着哪一个坐标轴。2022-11-8南京农业大学工学院机械工程系弹性力学中的基本假定n连续性假定整个物体的体积都被组成这个物体的介质所填满，不留任何空隙。n
9、完全弹性假定物体在引起形变的外力被除去之后能恢复原形，而没有任何剩余形变。n均匀性假定整个物体有同一材料组成的，物体的所有各部分具有相同的弹性。n各向同性假定物体的弹性在所有各个方向都相同。n小变形假定位移和形变是微小的，物体受力之后，整个物体所有各点的位移都远远小于物体原来的尺寸，因而应变和转角都远小于1。2022-11-8南京农业大学工学院机械工程系平面问题的基本理论n任何一个实际的弹性力学问题都是空间问题，但是如果所考察的弹性体具有某种特殊的形状，并且承受的是某些特殊的外力和约束，就可以把空间问题简化为近似的平面问题。n两种典型的平面问题q平面应力问题q平面应变问题2022-11-8南京
10、农业大学工学院机械工程系n由于板很薄，外力不沿厚度变化，应力沿板的厚度又是连续分布的，所以可以认为在整个薄板的所有各点：n只剩下平行于xy面的三个平面应力分量，即：n这种问题成为平面应力问题。平面应力问题n设有很薄的等厚度薄板，只在板边上受有平行于板面并不沿厚度变化的面力或约束。同时，体力也平行于板面不沿厚度变化。n设薄板的厚度为。以薄板的中面为xy面，以垂直于中面的任何一直线为z轴。所以有：222(0,(0,(0zzxzyz=z=z=),xyxyyx0,0,0zzxxzzyyz2022-11-8南京农业大学工学院机械工程系n只剩下平行于xy面的三个形变分量，即：n这种问题成为平面应变问题。n
11、由于z方向的位移处处为0，所以：，由于z方向的伸缩被阻止，一般平面应变问题n设有很长的柱形体，它的横截面不沿长度变化，在柱面上受有平行于横截面而且不沿长度变化的面力或约束。同时，体力也平行于横截面不沿长度变化。n假想该柱体为无限长，以任一横截面为xy面，以任一纵线为z轴，则所有一切应力分量、形变分量和位移分量都不沿z方向变化，而只是xy的函数，所有各点的位移矢量都平行于xy面，这种问题称为平面位移问题。0,0zxxzzyyz0z,xyxyn由对称条件可知：0zxzyn由胡克定律，相应的切应变：0z2022-11-8南京农业大学工学院机械工程系三大基本方程n根据静力学、几何学和物理学三方面条件，
12、建立三套方程。q平面问题中，根据微分体的平衡条件，建立平衡微分方程：(1-1)q根据微分线段上形变与位移之间的几何关系，建立几何方程：(1-2)q根据应力与形变之间的物理关系，建立物理方程：(1-3)(1-3)0yxxxfxy0yxyyfyx,xyxyuuvuxyxy112(1)(),(),xxyyyxxyxyEEE22112(1)(),(),11xxyyyxxyxyEEE2022-11-8南京农业大学工学院机械工程系yyydyyyxyxdyyxxdxxxyyxxCxfyfxyxydxx平衡微分方程n从弹性体中取出一个微分体，根据平衡条件导出应力分量与体力分量之间的关系式，也就是平面问题的平衡
13、微分方程。n从弹性体中取出一个微小的正平行六面体，它在x和y方向的尺寸分别为dx和dy，在z方向的尺寸为一个单位长度。n以x为投影轴，列出投影的平衡方程：0 xF()1xxdx dyxn约简以后，两边除以dxdy，得：0yxxxfxyn同理，以y为投影轴，列出投影的平衡方程，化简得：0yxyyfyx1xdy()1yxyxdy dxy1yxdx1xf dxdy02022-11-8南京农业大学工学院机械工程系n假定已知任一点P处坐标面上的应力分量x，y，x y=y x。求经过该点的，平行于z轴而倾斜于x轴和 y轴的任何倾斜面上应力。n从在P点附近取一个平面AB，它平行于上述斜面，并经过P点划出一个
14、微小的三棱柱PAB。当AB无限小而趋于P点时，平面AB上的应力就成为斜面上的应力。平面问题中一点的应力状态cos(,),cos(,)n xln ymyyxxyxxfyfxpypnnn设斜面AB 的长度为ds，则PB面及A面的长度分别为 lds及mds，而PAB的面积为 ldsmds/2，棱柱的厚度设为1。n由x轴平衡条件，得：02xxxyxldsmdsp dsldsmdsfn其中，fx为体力分量。将上式除以ds，并令ds趋于0（斜面AB趋于P点），即得：xxxyplmn由y轴平衡条件，得：yx yxypmln用n表示斜面AB的外法线方向，其方向余弦为：2022-11-8南京农业大学工学院机械工
15、程系uvvdyyvvdxxvuudxxuudyy几何方程n经过弹性体内的任意一点P，沿x轴和y轴的正方向取两个微小长度的线段PAdx和PBdy。假定弹性体受力后，P,A,B三点分别移动到P,A,B.n线段PA的线应变是：xuudxuuxdxx注:由于位移微小，y方向的位移v引起的PA的伸缩，是高一阶微量，略去不计。n线段PB的线应变是：yvyn线段PA与 PB之间的直角的改变，即切应变n线段PA的转角是：tanvvdxvvxdxxn线段PB的转角是：uyxyvuxy2022-11-8南京农业大学工学院机械工程系物理方程n在理想的弹性体中，形变分量和应力分量之间的关系，在材料力学根据胡克定律导出
16、如下：1()xxyzE,yzxyzxyzzxxyGGG1()yyzxE 2(1)EGn在平面应力问题中，z0，式变为：1()xxyE2(1)xyxyxyGE1()yyxEn在平面应变问题中，只要将上式中的E换为 ，换为 就得到平面应变问题的物理方程。21E12()1xxyE2(1)xyxyE2()1yyxE2022-11-8南京农业大学工学院机械工程系边界条件q若在su部分边界上给定了约束位移分量 和 ，则对于此边界上的每一点，位移函数u和v应满足条件：()u s()v s()()suu s()()svv sq其中（u）s 和（v）s 是位移的边界值，和 在边界上是坐标的已知函数。()u s(
17、)v sn边界条件表示在边界上位移与约束，或应力与面力之间的关系式。它可以分为位移边界条件、应力边界条件和混合边界条件。q位移边界条件：q应力边界条件：q若在su部分边界上给定了面力 和 ，则由平衡条件得出平面应力问题的应力(或面力）边界条件为：()xfs()yfs()()xyxsxlmfs()()yxysymlfs其中，l,m是边界面外法线的方向余弦。cos(,cos(,)n xln ym）2022-11-8南京农业大学工学院机械工程系圣维南原理n在求解弹性力学问题时，应力分量、形变分量和位移分量必须满足区域内的三套基本方程，还必须满足边界上的边界条件。但是，要使边界条件得到完全满足，往往遇
18、到很大的困难。n圣维南原理可为简化局部边界上的应力边界条件提供很大方便。n圣维南原理表明，如果把物体的一小部分边界上的面力，变换为分布不同但静力等效的面力（主矢相同，对同一点的主矩也相同），那么，近处的应力分布将有显著的改变，但是远处所受的影响可以不计。2022-11-8南京农业大学工学院机械工程系圣维南原理的应用n例，设有柱形构件，在两端截面的形心受到大小相等而方向相反的拉力F（a）。如果把一端或两端的拉力变换为静力等效的力，则只有虚线划出的部分的应力分布有显著的改变，而其余部分所受影响是可以不计的。n由于（d）图中，面力连续分布，边界条件简单，应力容易求得。其它三种情况，应力难以求得。把d
19、情况下的应力解答应用到其它三个情况，虽不能满足两端的应力边界条件，但仍然可以表明离杆端较远处的应力状态，而没有显著的误差。n图e，构件右端有位移边界条件，d情况的解答，不能满足位移边界条件，但e图右端的面力，一定是合成为经过截面形心的力F。所以把图d情况的解答应用于图e时，仍然只是在靠近两端处有显著的误差，而在离两端较远之处，误差可以不计。()0,()0ssuuvv2022-11-8南京农业大学工学院机械工程系yfxfxyxxxyxfyfNFMSF圣维南原理的应用(续)n例，厚度=1的梁中，左右两端x=l，的边界面是正、负x面，其上作用有一般分布的面力 。按照严格的应力边界条件，应力分量在边界
20、上满足：n上式要求在边界上y值不同的各点，应力分量与对应的面力分量必须处处相等，这种严格的条件是较难满足的。()(),()()xxlxxyxlyfyfy (),()xyfyfyn当1h时，x=l 是梁的边界的一小部分，可以应用圣维南原理，利用静力等效条件来代替，即，使应力的主矢量和主矩分别等于对应的面力的主矢量和主矩。2022-11-8南京农业大学工学院机械工程系圣维南原理的应用(续)n应力的主矢量和主矩的绝对值分别等于面力的主矢量和主矩的绝对值；n面力的主矢量和主矩的方向就是应力的主矢量和主矩的方向。/2/2/2/2()(),hhxxlxhhdyfy dy/2/2/2/2()(),hhxxl
21、xhhdy yfy dy y/2/2/2/2()()hhxyxlyhhdyfy dy y。/2/2(),hxxlNhdyF/2/2()hxxlhdy yM，/2/2()hxyxlShdyF2022-11-8南京农业大学工学院机械工程系有限单元法中基本量的矩阵表示n有限单元法(FEM)中，为了简洁清晰地表示各个基本量以及它们之间的关系，也为了便于编制程序利用计算机进行计算，广泛采用矩阵表示和矩阵运算。q平面问题中，物体受体力，可用体力列阵表示：(1)q物体受面力，可用面力列阵表示：(2)q3个应力分量的应力列阵表示：(3)q3个形变分量的应变列阵表示：(4)q2个位移分量的位移列阵表示：(5)(
22、)xTxyyfffff()Txyxy()Txyxy()Tduv()xTxyyfffff2022-11-8南京农业大学工学院机械工程系弹性力学中基本方程的矩阵表示q几何方程的矩阵表示为：(6)q物理方程矩阵表示为：(7)q 利用应力列阵和应变列阵（3）、（4）得：(8)其中矩阵 (9)()()TTxyxyuuvuxyxy21010100(1)/2xxyyxyxyED21010100(1)/2ED只于弹性常数E及有关，称为平面问题的弹性矩阵。2022-11-8南京农业大学工学院机械工程系虚位移原理n用u*和v*表示虚位移，用 表示与该虚位移相应的虚应变。n根据虚功方程：在虚位移过程中，外力在虚位移
23、上所做的虚功等于应力在虚应变上所做的虚功。n对于厚度为t的薄板，虚功方程可用矩阵表示为：*,xyxy*()()()(10)TTTxASAdfdxdy tdfds tdxdy t n其中，分别为体力列阵，面力列阵和应力列阵。*(11)Tdu v为虚位移列阵,f f*(12)Txyxyv为虚应变列阵n有限单元法中，作用于弹性体的各种外力常以作用于某些点的等效集中力来代替。在厚度为t的薄板上，设作用于i点的集中力沿x及y方向的分量为Fix,Fiy,作用于j点的力为Fjx,Fjy等。这些集中力以及它们相应的虚位移用列阵表示为：2022-11-8南京农业大学工学院机械工程系*,jyjFv*,jxjFu*
24、,ixiF u*,iyiF v虚位移原理(续)*()()(15)TTxAFdxdy tn代入虚功方程，得：n 上式为集中力作用下的虚功方程。n集中力列阵 （13）(.)TixiyjxjyFFFFFn虚位移列阵 （14）*(.)Tiijjuvuvn外力在虚位移上所做的功为：*()TixiiyijxjjyjF uF vF uF vF2022-11-8南京农业大学工学院机械工程系（1）取三角形单元的结点位移为基本位置量：(a)其中，称为单元的结点位移列阵；（2）应用插值公式，由单元结点位移求出单元的位移函数：(b)其中，N称为形函数矩阵；（3）应用几何方程，由单元的结点位移求出单元的应变：(c)其中
25、，B是表示 与 之间关系的矩阵；三角形单元离散化结构分析步骤()()eTijmTiijjmmuvuvuv(,)(,)eu x ydNv x yeeBe2022-11-8南京农业大学工学院机械工程系 (f)其中，Fe 是单元的结点力，k称为单元劲度列阵；对三角形板单元，节点力为：(e)（5）应用虚功方程，由单元的结点应力求出单元的结点力。假设把单元和节点切开，对右图中的i节点：节点对单元的作用力为节点力，作用于单元上。三角形单元离散化结构分析步骤(续）（4）应用物理方程，由单元的结点位移求出单元的应力：(d)其中，S称为应力转换矩阵；eeFkTTeijmixiyjxjymxmyFFFFFFFFF
26、FeS()TiixiyFFF Fe是作用于单元的外力，此外，单元内部还作用有应力。根据虚功方程，可以将单元的节点力Fe用应力来表示，从而得到节点力的公式：2022-11-8南京农业大学工学院机械工程系（7）列出各结点的平衡方程，组成整个结构的平衡方程组。由于节点i受有环绕节点的单元移置而来的节点载荷 和节点力 因而i节点的平衡方程为：（i=1,2,n）(h)三角形单元离散化结构分析步骤(续）（6）应用虚功方程，将单元中的外力载荷向结点移置，化为结点载荷（即求出单元的节点载荷）：(g)()eTLLiLjLmFFFFiLieeFF()TLixLiyLjxLjyLmxLmyFFFFFF()TLiLi
27、xLiyFFF()TiixiyFFF 将（f）代入（h），整理得：（j)其中，K称为整体刚度矩阵，FL是整体结点载荷列阵，是整体结点位移列阵。LKF 在上述求解步骤中，（2）至（6）是针对每个单元进行的，称为单元分析；（7）是针对整个结构进行的称为整体分析。n对三角形三个结点i,j,m结点，位移函数应当等于该节点的位移值，即：三角形单元的位移模式n对每个单元，只要求得单元中的位移函数，就可以应用几何方程求得应变，再应用物理方程求得应力。有限单元法中常取结点位移为基本未知量，由单元的结点位移求出单元中的位移函数是首先必须解决的问题。n可以假定一个位移模式，来表示单元中的位移函数（即在单元中做出位
28、移插值函数）。三角形单元中，可以假定位移分量只是坐标的线性函数，即假定：123456,uxy vxy111,(,)(17)111jjmmiiijjmmxyxyxyNi j mxyxyxy123456,iiiiiixyuxyv123456,jjjjjjxyuxyv123456,mmmmmmxyuxyvn6个方程解出1-6，代入u，v式整理得：其中：iijjmmuN uN uN uiijjmmvN vN vN v(16)2022-11-8南京农业大学工学院机械工程系三角形单元的位移模式nNi也可以该写成为：()/2,(,)(18)iiiiNab xc yAi j m111(20)21iijjmmx
29、yAxyxy其中系数ai,bi,ci是：11,(,)(19)11jjjjiiimmmmxyyxabci j mxyyx 其中A就等于三角形ijm的面积：n 按照解析几何学，在图示的坐标系中，为了得出的面积A不致成为负值，节点i，j，m的次序必须是逆时针转向的。n Ni，Nj，Nm这三个函数，表明了单元ijm的位移次形态（也就是位移在单元内的变化规律），因而称为形态函数，简称形函数。2022-11-8南京农业大学工学院机械工程系三角形单元的位移模式n位移模式的表示式可用矩阵表示为：iijjmmiijjmmN uN uN uudN vN vN vv 简写为：(23)edN其中()(24)TeTij
30、miijjmmuvuvuv是单元的节点位移列阵。000(25)000ijmijmNNNNNNN是形态函数矩阵或形函数矩阵。n有限单元法中，应力转换矩阵和劲度矩阵的建立以及载荷的移置等，都依赖于位移模式。2022-11-8南京农业大学工学院机械工程系简写为：其中矩阵B可写成分块形式：其子矩阵为：单元的应变列阵和应力列阵n利用几何方程和物理方程，求出单元中的应变和应力，用结点位移表示：将位移函数(16)和(18)代入几何方程(6)，得出用结点位移表示单元应变。(26)eB00010002iiijmjijmjiijjmmmmuuvxbbbuucccvycbcbcbuvuxyv(27)ijmBBBB0
31、10,(,)(28)2iiiiibBci j mcb2022-11-8南京农业大学工学院机械工程系n将D表达式(9)和B表达式(27)代入上式，并写成分块形式，即得到平面应力问中的应力转换矩阵：单元的应力列阵（续）n再将单元的应变式(26)代入物理方程(8)，得出用结点位移表示单元中应力的表达式。(29)eDDB(32)ijmSSSS21122,(,)(33)2(1)iiiiiiibcESbci j mAcb(30)eS(31)SDB其中子矩阵为：简写为：其中，2022-11-8南京农业大学工学院机械工程系n由式(26)引起的虚应变为：n由于结点力在虚位移上的虚功应当等于应力在虚应变上的虚功，
32、即：单元的结点列阵与劲度矩阵n对于任一单元，均假设所受的外力载荷已经被移置到结点上，并且单元已经切开，如右图所示：n单元只受到结点对单元的作用力，即结点力：*()()eTiijjmmuvuvuv(34)TixiyjxjymxmyFFFFFFTeijmFFFF*()eB*()()e TeTAFdxdyt*()e TTeAB DBdxdytn假想在结点i,j,m处发生了虚位移，即：对单元而言，这些结点力是外力，使单元内部产生应力。n从而建立了单元结点力和结点位移之间的关系。对于三角形单元，B中的元素为常量。n并且 ，因此，k可简写为：k称为单元的劲度矩阵。单元的结点列阵与劲度矩阵(续）n由于 中的
33、元素是常量，并且虚位移的值可以是任意的：则*()eeTeAFB DBdxdytAdxdyATkB DBtA(37)iiijimjijjjmmimjmmkkkkkkkkkk21122,(,;,)114(1)22rsrsrsrsrsrsrsrsrsb bc cb cc bEtkri j m si j mAc bb cc cb bn将B和D表达式代入上式，得：(35)TAkB DBdxdyt令eeFk则式可以简写为2022-11-8南京农业大学工学院机械工程系LjyFLjxFLiyFLixFLmyFLmxFPyfPfPxf载荷向节点移置，单元的载荷列阵n设单元ijm在坐标为（x，y）的任意一点M，在
34、单位厚度上受有集中载荷fP，其坐标方向的分量为fPx 和fPy，用矩阵表示为fP=（fPx fPy）T，将此集中力移置到单元的节点处，转换为节点载荷，并且单元节点载荷列阵表示为：*()()eTiijjmmuvuvuv,(39)TLixLiyLjxLjyLmxLmyFFFFFFTeLLiLjLmFFFFn假想单元的各点发生了虚位移：n由位移模式，相应于集中力fP的作用点（x,y）的虚位移为：*(,)(,)()Tedux yv x yNn集中载荷的移置2022-11-8南京农业大学工学院机械工程系载荷向节点移置，单元的载荷列阵（续）*()()eTiijjmmuvuvuvn由于虚位移可以是任意的，所
35、以：(40)eTLPFNf tn把N的表达式（25）代入上式，上式改写为：TeLLixLiyLjxLjyLmxLmyFFFFFFF(41)TiPxiPyjPxjPymPxmPyt N fN fN fN fN fN f其中，Ni,Nj,Nm，为它们在M点的函数值：*()()()()e TeTe Te TTLPPPFdf tNf tNf tn根据静力等效原则，节点载荷在节点虚位移上的虚功，等于原载荷集中力在其作用点的虚位移上的虚功，即：2022-11-8南京农业大学工学院机械工程系载荷向节点移置，单元的载荷列阵（续）n例，设单元ijm的密度为，试求自重的等效节点载荷。分析：因为fx=0，fy=-g
36、，故由式（43）得：0,LixLjxLmxFFF1,(,)3LiyiAFgtN dxdygtAi j m n由设上述单元受有分布的体力f=（fx fy）T，可将微分体积tdxdy上的体力ftdxdy当作集中力，利用（40）式积分，得到：(42)eTLAFtNfdxdy(43)TixiyjxjymxmyAtN fN fN fN fN fN fdxdyn体力的移置q注意单元的自重为-gtA，可见移置到每个节点的载荷均为1/3自重。2022-11-8南京农业大学工学院机械工程系载荷向节点移置，单元的载荷列阵（续）(44)eTLSFtN fds(45)TixiyjxjymxmyStN fN fN fN
37、 fN fN fds,0LixiLjxjLmxmLiyLjyLmyijijijFqtN ds FqtN ds FqtN dsFFFn由设上述单元的某一边上受有分布的面力 ，可将微分面积tds上的面力 当作集中载荷，利用（40）式积分，得到：()TxyfffftdsTeLLixLiyLjxLjyLmxLmyFFFFFFFn面力的移置n例，设在ij边上受有沿x方向的均布面力q，试求等效节点载荷。分析：因为 ，故由式（45）得：,0 xyfq f1,02LixLjxLmxFFqtij F2022-11-8南京农业大学工学院机械工程系注意：式(46)和(48）中的编码i,j,m仅是每个单元的局部编码，
38、对于整个结构，则将结点的平衡方程按整体结点编码1，2，n排列起来，就组成整个结构的结点平衡方程组：整体的结构分析 节点平衡方程组,()iiiiijjimminnn i j mFkkkk (49)LKF(46)iLieeFF,(48)ixLixiyLiyeeeeFFFFn因此，结点i的平衡方程是：n以上几节的分析都是针对单元进行的，即一将单元上的外力载荷都向节点移置而成为节点载荷；另一方面求出节点载荷与单元之间的相互作用力，如左图所示。节点对单元的作用力是节点力，相反，单元对节点的作用力是节点力的负值。于是，作用于结点i上的力，有结点载荷FLi ，和结点力的负值，即：其中，是对环绕结点i的单元求
39、和，写成标量形式：e2022-11-8南京农业大学工学院机械工程系n由整体平衡方程组，解出结点位移，便可由式(23)和(30)求出每个单元的位移函数和应力。整体的结构分析 节点平衡方程组rsrseKk其中，整体结点位移列阵：整体结点载荷列阵：K是整体刚度矩阵，其元素是：12(.)(51)LLLLnFFFF12(.)(50)n(49)LKFn整个结构的结点平衡方程组即整体劲度矩阵的元素,Krs就是按整体节点编码的、同下标rs的单元劲度矩阵元素叠加而得到的。21122114(1)22rsrsrsrsrsrsrsrsrsb bc cb cc bEtkAc bb cc cb b2022-11-8南京农
40、业大学工学院机械工程系平面有限元解法（例）n设有对角受压的正方形薄板（如上图所示），载荷沿厚度均匀分布，为2N/m。试对该结构进行整体分析，建立整体刚度矩阵和整体结点载荷列阵，建立整体结点方程组，通过编程求解出结点的位移，并从而求出各单元的应力。（为简单起见，取板的厚度t=1,弹性常数E=1,泊松比0）2022-11-8南京农业大学工学院机械工程系平面有限元解法划分单元n由于平面薄板沿xz面和yz面均对称，所以只取1/4之一部分作为分析和计算对象。将对象划分成4个单元，共有6个结点，单元和结点上均编上号码，其中结点的整体编码1至6，以及个单元的结点局部编码i,j,m,均示于上图中。单元号局部编
41、码整体编码i3526j1253m24352022-11-8南京农业大学工学院机械工程系平面有限元解法整体劲度矩阵n每个单元，结点的局部编码和整体编码对应关系已经确定，每个单元劲度矩阵中任一子矩阵在整体劲度矩阵中的位置及其力学意义也就明确了。如单元的kii，即k33，它的四个元素表示当结构的结点3沿x或y方向有单位位移时，在结点3的x方向或y方向引起的结点力。n暂时不考虑位移边界条件，把所分析结构的整体结点平衡方程组列出：111213141516112122232425262231323334353633414243444546445152535455565561626364656666LLLL
42、LLKKKKKKFKKKKKKFKKKKKKFKKKKKKFKKKKKKFKKKKKKFn整体劲度矩阵写成66的矩阵，它的每个子块是22的矩阵，实际它是一个1212的矩阵。如K23，它的四个元素表示当结构的结点3沿x或y方向有单位位移时，在结点2的x方向或y方向引起的结点力。2022-11-8南京农业大学工学院机械工程系平面有限元解法整体劲度矩阵续n由于于结点3和结点2在结构中是通过和这两个单元相联系，因而K23应是单元 的k23和单元 的k23之和。同理，可以找到各单元劲度矩阵中所有子矩阵在整体劲度矩阵K中的位置，得到整体劲度矩阵。111112313223113134344222233422
43、344444000000000000jjjmjimjmmjjiimiimjmjiijjiimmiiimmjjmjjmjimjmmmiijjijmmjimiijjmmmiijimiikkkkkkkkkkkkkkkkkkkkkKkkkkkkkkkkkkkkk n式中k的上标1,2,3,4表示是哪一个单元的劲度矩阵中的子矩阵，空白处是22的零矩阵。2022-11-8南京农业大学工学院机械工程系平面有限元解法整体劲度矩阵续n对于单元、，根据公式，可求得A=0.5m2，0.50000.5000.250.2500.250.2500.250.2500.250.250000.500.50.50.250.250
44、0.750.2500.250.250.50.250.75kEn将上式中各子块的具体数值代入整体刚度矩阵K表达式中，得出整体刚度矩阵。1,0,1,0,1,1ijmijmbbbccc n对于单元，根据公式，可求得A=0.5m2，1,0,1,0,1,1ijmijmbbbccc n把0，t1m，代入单元的劲度矩阵，得两种单元的劲度矩阵k都是：(37)2022-11-8南京农业大学工学院机械工程系平面有限元解法整体劲度矩阵续n整体刚度矩阵K0.2500.250.2500.2500000000.500.5000000000.2501.50.2510.250.250.2500.25000.250.50.25
45、1.50.250.500.50.250000010.251.50.25000.50.2500.250.2500.250.50.251.5000.25100000.250000.750.KE250.50.2500000.250.5000.250.7500.25000000.250.50.250.501.50.250.50.25000.2500.2510.250.250.251.500.25000000000.500.5000000.250000.250.2500.25(38)2022-11-8南京农业大学工学院机械工程系平面有限元解法位移边界条件n位移边界条件为：(39)1244560uuuvv
46、vn因此，整体结点的位移列阵就简化为：123356()Tvvuvuun与这6个零位移分量相应的6个平衡方程不必建立，因此，将整体刚度矩阵中，第1、3、7、8、10、12各行以及同序号的各行划去，因而整体劲度矩阵K简化为：0.50.500000.51.50.250.50.25000.251.50.250.5000.50.251.50.25000.250.50.251.50.500000.50.5KE2022-11-8南京农业大学工学院机械工程系平面有限元解法整体结点载荷列阵n确定了每个单元的结点载荷列阵：(38)n根据各单元的结点局部编码与整体编码的关系，确定三个子块FLi,FLj,FLm在FL
47、中的位置。TLixLiyLjxLjyLmxLmyFFFFFFTeLLiLjLmFFFF111232134324234546LjLmLjLiLiLmLjLmLiLjLmLiLLLLLLLFFFFFFFFFFFFFFFFFFFn由于该结构只是在结点1受有向下1N/m的载荷，因而，非零元素子块，只有101LLjFFn在考虑了边界条件后，整体载荷列阵为：(100000)LF 2022-11-8南京农业大学工学院机械工程系平面有限元解法求解整体结点载荷列阵n求解化简后的整体刚度矩阵：1233560.50.5000010.51.50.250.50.250000.251.50.250.50000.50.25
48、1.50.250000.250.50.251.50.5000000.50.50vvuEuuu(39)LKFn求解以后，得结点位移：1233563.2531.2530.08810.3740.1760.176vvuuEuu2022-11-8南京农业大学工学院机械工程系平面有限元解法求解应力转换矩阵n应用单元的应力转换矩阵S，求出各单元中的应力：q根据0，以及已求出的A、b和c的值，再由式(21)和(22)得出应力转换矩阵如下，对于单元、：n对于单元10001000010100.50.500.50.5SE10001000010100.50.500.50.5SE2022-11-8南京农业大学工学院机械
49、工程系平面有限元解法求解各单元中的应力（续）n应用单元的应力转换矩阵S，求出各单元中的应力：331121000100.08800001012.00000.50.500.50.50.4400 xyxyuvEvv Pa单元单元52201000100.17600001011.25300.50.500.50.5000 xyxyuEv Pa2022-11-8南京农业大学工学院机械工程系平面有限元解法求解各单元中的应力n应用单元的应力转换矩阵S，求出各单元、中的应力：2533301000100.0880001010.374000.50.500.50.50.308xyxyvuEuv Pa单元单元633450
50、10001000001010.37400.50.500.50.50.1320 xyxyuuEvu Pa2022-11-8南京农业大学工学院机械工程系平面有限元解法计算机编程解题步骤n 划分单元格，并按照一定的规律将所有的结点和单元格分别编上号码，需要注意单元编码和整体编码的对应关系。n选定一个直角坐标系。按照计算程序的要求，填写各种输入信 息有：每个结点的坐标值，即x1,y1,x2,y2 等；材料的弹性常数值；各种载荷信息，即载荷点的点号及载荷的大小等；约束信息，即哪些结点哪个方向上的位移为零或为某个已知值。将这些信息按照计算机程序规定的格式输入。n计算程序中对输入的各种信息进行加工、运算，一
