深入浅出的讲清楚有限元法

1、弹簧单元描述：图 1-12个节点：i, j 节点位移：Uj,Uj 节点力：fi , f j 单元自由度：2弹簧的物理特性:IL_ kkujJ(1-2)“有限元法基础及应用”补充讲义（一）顾克秋（2005 年 3 月），、引子一一弹簧单元与弹簧系统目标：掌握离散结构直接刚度法分析的原理和形式。了解有限元位移法列式的形 式和基本概念。1、典型弹簧单元分析已知弹簧力一一位移关系： F二k k-弹簧刚度二Uj _ ui 弹簧伸长量F - 弹簧力，拉伸为正考虑弹簧变形平衡时的条件和弹簧物理特性，得到下列方程:fi _ _ F - - k（uj _ 5） = kui _ ku j(1-1)fj = F 二

2、 k(uj - uj = -kui kuJ写成矩阵形式:写成矩阵符号形式:(1-3)f 二 kd式（1-2）、（1-3）为弹簧单元的刚度方程，反映了单元特性：节点力与节点位移 之间的关系。式中：k弹簧单元的刚度矩阵d单元节点位移列阵f单元节点力列阵（注意：单元节点力是节点对单元的作用力）弹簧单元刚度方程讨论：1）k有何特点？对称、奇异、主对角元素恒正2）k中元素代表什么含义？刚度系数大小等于弹簧刚度；每列元素代表一端固定、另一端产生单位 位移时加在弹簧单元上的节点力。3）上面单元刚度方程可以求解吗？为什么？不可以。刚度方程仅仅表征一个典型单元的弹性特性，单元水平上无法 确定单元节点位移。只有把

3、系统中所有单元特性集成后，在系统水平上才可 能求出所有未知位移和反力。单元水平上，若已知单元的节点位移，可由刚 度方程求出所有单元节点力分量。若节点力已知，单元节点位移不能确定， 单元可作刚体运动（小位移）。这也是单元刚度矩阵奇异性的物理解释。2、弹簧系统整体分析原理ki4Icn以右图的一个弹簧系统为例，研究如何由单元特性集成系统特性并建立对系 统进行求解的控制方程。由前面得到的弹簧单元的刚度方程公式（1-2），分别写出2个弹簧单元的特性方程如下:单元1单元2-ky3图1-3（1-4）(1-5)（注：右端节点力分量的下标1, 2为单元节点的局部编号，上标是单元号）下面按两个方法完成系统特性的装

4、配和控制方程的建立。并在特定条件下求1）由节点平衡方程导出：系统处于平衡时，考虑各节点（1, 2, 3节点）的平衡条件：节点受到的外载荷与节点受到与其连接的所有单元对其作用力（单元节点力的反作用力）之和等于零。因此有下列（节点）平衡方程（组）F111(1-6)把单元特性（1-4），（ 1-5）代入（1-6）得到:(1-7)k1u2F2 二 _k1u1 (kk2)u2 - k2u3F3 二-k2k2U3写成矩阵形式:-右0-何+ k 丿0-他*(1-8)或矩阵符号形式:（1-9）式（1-8 ）, （1-9）就是系统平衡方程，该方程建立了离散系统的外载荷与节 点位移之间的关系，是求解节点位移的控制

5、方程。弹簧系统的结构总刚度矩阵D 系统节点位移列阵F 系统节点载荷列阵讨论：（1）K有那些特点和性质?（2）上述方程能求解吗？2）由单元刚度方程叠加导出将单元1, 2的刚度方程（1-4），（ 1-5）进行增广（扩大到系统规模）：处 一 k 0IL/；q一& & 0 *i* = J 1 fi（1-10）_ 00 ojv00 0 0 _比010 k2 -k2=卜 （1-11）0 -k. k、上述两个矩阵方程叠加，得:k0 广一何A +心-k2 *0（1-12）-k2k2 V J jfl 上式中代入节点力平衡关系（1-6）,就得到与（1-8 ）相同的节点平衡方程。 上述两种方法都必须考虑1）单元特性

6、集成；2）离散结构的节点上外载荷（系 统外力）与节点力（系统内力）的平衡。因此方程（1-8 ）的本质是节点的力平 衡关系，左边是由节点位移表示的（总）节点力，右边是节点所受外载荷。3）给定载荷和约束条件下的求解 设边界条件为：(1-13)比=0F F P23则节点平衡方程（1-8 ）变化为:0 0代十处卜=_ 0Pk J(1-14)该方程组展开后分为2个部分:第1个方程变化为:第2, 3个方程变化为：（1-15）(1-18)(1-16)先后解方程（1-15 ）、（ 1-16）得到:(1-17)卜be卜:=-2 P从而解出了系统的未知位移和未知反力，并可以进一步求弹簧力3、例题k/ki P -f

7、NWMWMNWFWWW- x1234图 1-4图1-4所示一个3个弹簧的系统。=100N/mm,k2 =200N/mm,k3 =100N / mm, P =500N,Ui 二 U4 =0:（a）系统总刚度矩阵(b)节点2, 3的位移(c)节点1、4的反力(d)弹簧2中的力解：（a）:写出各单元刚度矩阵:100-100-100100 200-200-200200_ 100-I0CT-100100(N/inm)(N/nnn)应用叠加法直接得到系统总刚度矩阵:100-1( 00 -100100+200-2(1000200200十1001000000100町112或:10010000100300-20

8、000-200300-10000-100100该总刚度矩阵特点：对称性、奇异性、稀疏、非零元素沿主对角线呈带状分布。(b)参考前面的做法（1-8 ）式）和求出的总刚度矩阵，写出系统节点平衡方程:_ 100-10000-100300-20000-200300-1000-100100f 1%片0（1-19） =斗考虑到位移边界条件：Ui = 5 = 0则平衡方程组（1-19）第2，3方程化为:300-200求解上式得:mm)(c):由(1-19)的方程1, 4得：F =-100划=-200 (N)F4 = -100u3 =-300 (N)(d):弹簧2内力为：F2 二 k2 2 二 k2(U3 -

9、 U2)-200 3 2= 200(N)(拉力)4图1-54、练习题对图示弹簧系统，试用叠加法 求其总刚度矩阵。并根据节点平衡 方程的含义，尝试由各单元刚度矩 阵的元素直接写出总刚度矩阵的非零】、杆单元目标：通过杆单元特性方程的建立，初步掌握有限元法单元分析的过程和原理 了解杆系结构分析的原理。1、等截面杆单元及其刚度矩阵研究2节点等截面杆单元：单元上的力学量和基本关系如下：L 杆长A截面积E弹性模量dug V-dxu = u(x)杆单元位(x)杆单兀应变；-；(x)杆单兀应力应变一位移关系：单元节点位移：应力一应变关系：if血单元节点力：d三A下面研究杆单元的单元特性。1)直接法导出杆单元特

10、性采用材料力学基本知识对单元进行力学分析:杆单元伸长量: II；A,E图 2-1(2-1)(2-2)(2-3)杆应变:杆应力:杆内力:EA.：L(2-4)(2-5)(2-6)杆的轴向刚度:(2-7)由于轴向变形模式下，杆单元的行为与弹簧单元相同,因此可比照弹簧单元的刚度方程（1-2），考虑到（2-7），直接写出杆单元的刚度方程:fiEA 1 -1 Uif = 一I$ 卜L 广11 jUjj写成符号形式：f = kd杆单元刚度矩阵为：_ EA 1-1k = 一 IL 厂11 一(2-9)(2-10)2）公式法导出杆单元特性步骤如下：（1）在单元上定义近似位移场把一个单元上的位移分布假设为简单多项

11、式函数。有限元法中用插值法通过节点位移分量作为待定参数来构造单元位移函数起见引入局部坐标对图2-1的杆单元，方便01由于该杆单兀只有2个未知位移分量，因此单兀上假设的简单位移函数米用 一次多项式。故对单元的节点位移进行线性插值。则容易定义出节点的插值函数如下：忙 =1爲（ 2-11）因此单元上近似位移函数的插值形式为：-:宀一! 一*l-i-（ 2-12）该位移函数也称为单元的位移模式，这里是线性位移模式。式（2-11）中的插值函数又称为形状函数，简称形函数。式（2-12）写成矩阵形式为：上式中N称为单元的形函数矩阵。u = Ni Nj 宀二 Nd（2-13）J lujJ式（2-13）是有限元

12、法中最重要的关系式之一， 通过该式把单元上的近似位 移分布函数用节点位移来表示，为进行单元层次上的分析打下了基础。（2）单元应变和单元应力由杆一维变形的应变位移方程（2-1 ）和单元的位移函数（2-13）求出 单元的应变分布和节点位移的关系：dudx_dx Nd = Bd(2-14)(2-15)式中:B= d Ni( ) Nj( )1- 1-1/L 1/L dxB称为单元的位移应变转换矩阵，简称应变矩阵由一维杆的应力一一应变关系（2-2 ），得单元应力和单元节点位移的关系:(2-16)（3）用弹性体的虚位移原理导出杆单元刚度方程变形体的虚位移:假想在弹性体上发生的，满足位移许可条件（内部连续，

13、边界协调）的微小、 任意位移场。可以理解为某个位移场的微小扰动（变分） 虚位移的特征：1）假想的，与真实位移无关；2）几何上是许可的：连续、协调;3）微小、任意大小。虚位移原理：弹性体受力平衡时，若发生虚位移，则外力虚功等于弹性体内的虚应变能（应 力在虚应变上做的虚功）。下面把虚位移原理应用在所研究的杆单元上。定义杆单元的虚位移：节点虚位移 ，单元虚位移，单元虚应变节点虚位移:单元虚位移:单元虚应变:=旦（孔）=B6ddx那么，节点力虚功:单元虚应变能:、；T；dV = 、dTBtEBddVVV=6d （B E BdV d iV丿对杆单元应用虚功原理，那么上述节点力（外力）虚功等于虚应变能，因

14、此有下列关系:fBT EBdV(2-25)则上式的第2个方程为：字即讣P(2-26)求解该方程后得到系统的位移解:PL1亠1 3EA0计算应力: 单元1E：(2-27)L 13EAE，L匚 u3 一 u2PL0-3EA丿3A单元2提示:1）本例中单元应力的计算采用了材料力学中的方法，与采用有限元单元 应力公式心=E ；二EBd的结果相同，请验证。2）对锥形杆，单元截面积可以用平均值，从而转化为类似本题的问题 求解。3）求应力之前需要先求出节点位移，因此本方法称为有限元位移法。图 2-32-3。先检查杆右端（节点3）与4）如果杆上受连续分布的轴向载荷或节点之间受轴向集中载荷，分析 时可以按照虚功

15、相等的原则先把单元上的载荷等效移置到节点上。已知:尸二 60 x lOg A =2.0x 04N/mm A 250rnm:, = 150mm, A 1.2 mm求：杆两端的支反力。单元、节点的定义如图 墙壁是否接触。计算右端的自由伸长:EA 1-10 _-1 2-1 =V卜;0-11K,J J系统平衡方程为:引入下列载荷和位移边界条件:所以，右端间隙将闭合。-、(2-28)出=(K； = POxOi M_ 1-10 I01fA；lEA2 -1=P L_ 0一 1 1A1 J则有限元平衡方程（2-28）成为:分离出第二个方程:(2-29)(2-30)学2 T炸旧即：(2-31)=L5 mm全部位

16、移解为:广F码卜二斗15卜12LJmm)(2-32)根据上式位移解，从系统平衡方程（2-28 ）的第1, 3个方程分别求出支反 力如下：= -l.Ox IO4 N解毕2、2-D和3-D空间中的杆单元（平面和空间桁架单元）1）2-D空间中 的杆单元2-D空间中建立杆单元的基本思路是根据前面在杆的一维局部坐标系 下建立的单元特性方程通过坐标变换， 转换为2-D总体坐标系下的方程，同时得 到坐标变换后的单元刚度矩阵。而系统整体分析的原理和方法与一维情况相同。图 2-4（1）变换图2-4为一个杆单元及其局部坐标系与 2-D总体坐标的关系。节点的位移分量和节点力分量在 2-D局部坐标系x-y下描述，杆节

17、点i具有2个自由度:位移分量为Ui，vi ；节点力分量为fxi，fyi其中只有x方向的位移分量和节点力分量用来描述单元特性。节点上的位移和节点力向量在2-D局部坐标系与2-D总体坐标系下的变换如下:M- = u. cosd 十叫 sin Q = I/ mvt = -ui sin 6? + 片 cos仪二-in上述变换的矩阵形式:称为方向余弦。1/ m =00/tniV00-mIV.1门k J J(2-36)或:d 二 Td其中T T 0T =丨0 T 一比照（2-37）得到单元节点力的坐标变换式:f = Tf(2-37)(2-38)(2-39)(2-40)（2） 2-D空间刚度矩阵下面可以导出

18、单元刚度方程和单元刚度矩阵的坐标变换式。已经知道杆的一维局部坐标系下的刚度方程为：EA 1L -1把该方程扩充到-j0-1EA000L-101000(2-41)01Wffx：l01Vf 1| yi卜=丿2L 0 2一V J(2-52)解得：:氏：L ；R =电M ： EA f2(2-53)F面按公式（2-48 ）分别求得单元1，2的应力:E 2 1-1L 20“ L 0-.21 1 1一J = ?!(r+p2)EAR 2A(2-54)E和2 r1 LP2i P2=1-1 -1 1 J丿H ( P - p2 )L 2EA1。32A|p2(2-55)试对节点2的载荷分量不同情况，讨论本题解的平衡、对称、反对称等。3）3-D空间中的杆单元3-D空间中建立杆单元的基