第五章 线性三角形单元

1、三角形单元 3:212引 言 杆梁结构：由于有自然的连接关系，可以凭一种直觉将其进行自然的离散。 连续体：它的内部没有自然的连接节点，必须完全通过人工的方法进行离散。三维问题平面问题平面应力平面应变平面问题平面应力平面应变离散3:213三节点平面三角形单元 x, u y, v 1 (x1, y1) (u1, v1) 2 (x2, y2) (u2, v2) 3 (x3, y3) (u3, v3) A fsx fsy 112233 euvuvuvd节点1的位移节点2的位移节点3的位移三节点三角形单元的位移函数可假设为：123456,u x yxyv x yxy“位移函数”也称“位移模式”，是单元内

2、部位移变化的数学表达式，是坐标的函数。有限元分析必须事先给出（设定）位移函数。一般而论，位移函数选取会影响甚至严重影响计算结果的精度。弹性力学中，恰当选取位移函数不是一件容易的事情。有限单元法中当单元划分得足够小时，把位移函数设定为简单的多项式也可得到相当精确的结果。这正是有限单元法具有的重要优势之一。引入位移函数的概念：3:214平面三角形单元显然，三角形三个节点的的位移可由下列方程给出，在各节点上的水平位移方程为： u1=1+2 x1 +3 y1 u2=1+2 x2 +3 y2 u3=1+2 x3 +3 y3解出11112223331,111xyuu x yxyxyuxyu11111222

3、23333111xyuxyuxyu3:215平面三角形单元11112322331111xyNNNxyxyxy假设求得1233223322311331133122112211()()21()()21()()2Nx yx yyy xxxyANx yx yyy xxxyANx yx yyyxxx yA11223311121xyAxyxy其中A是三角形的面积3:216平面三角形单元1 1223 3,v x yN vN vN v式中N1，N2和N3是坐标的函数，反映了单元内近似解的形态，称为单元的形函数，数学上它反应了由节点的场量对单元内任意一点场量的插值，也叫做插值函数。三个函数其实描述的就是单元上近

4、似解的插值关系，它决定了近似解在单元上分布的形状，所以称它为形函数（shape function）。这里值得注意一下的是近似解，前面我们说过，假设位移模式是线性变化的，实际情况并不一定是线性变化的，所以我们通过所做假设得到的结果只能说是近似解，而不能说是精确解。 为什么叫形函数?1 12233,u x yN uN uN u同理( , )( , )ex yx yuNd3:217平面三角形单元12iiiiNab xc yA其中ijki = 1, 2, 3j = 2, 3, 1k = 3, 1, 2ijkmmjjjmmjiyxyxyxyxamjmiiyyyyb11mjjmixxxxc111211ii

5、jjmmxyAxyxy三角形的形函数可统一表示为：3:218形函数的性质在单元任一点上三个形函数之和等于1(单位分解性)1. 三个形函数只有两个是独立的2. 当三角形单元的三个结点的位移相等*ijmuuuu*( , )()iijjmmijmu x yN uN uN uNNNuu( , )( , )( , )1() ()()2ijmijmijmijmN x yN x yNx yaaabbb xcccyA第一列与它的代数余子式乘积之和第一列与第二列的代数余子式乘积之和第一列与第三列的代数余子式乘积之和2A0013:219形函数Ni 在节点i 上的值等于1，在其它节点上的值等于0。0),(0),(1

6、),(mmijjiiiiyxNyxNyxN0),(1),(0),(mmjjjjiijyxNyxNyxN1),(0),(0),(mmmjjmiimyxNyxNyxNNi =1ijmNj =1ijmNm =1ijm形函数的性质3:2110在三角形单元的边界ij上任一点（x，y），有: 形函数的性质0),(),(1),(yxNxxxxyxNxxxxyxNmijijijiixxixjxyNi(xi，yi)j (xj，yj)m (xm，ym)Ni(x、y)1证jijixxxxyxN1),(ijijijiijijixxxxxxxxxxxxxxyxN1),( , )( , )( , )1ijmN x yN

7、x yNx y()jiiijimiimyyyyxxxxbyyxxcij方程3:2111形函数的性质相邻单元的位移在公共边上是连续的ijpm0),(),(1),(yxNxxxxyxNxxxxyxNmijijijii形函数在单元上的面积分和边界上的线积分公式为 ijijiAildlNAdxdyN213式中 为 边的长度。 ijlijxxixjxyNi(xi，yi)j (xj，yj)m (xm，ym)Ni(x、y)1Ni =1ijm3:2112形函数的性质完备性包含常应变项和刚体位移项如果在势能泛函中所出现的位移函数的最高阶导数是m阶，则选取的位移函数至少是m阶完全多项式。协调性相邻单元公共边界保持

8、位移连续如果在势能泛函中所出现的位移函数的最高阶导数是m阶，则位移函数在单元交界面上必须具有直至(m-1)阶的连续导数，即Cm-1连续性。 如果在单元交界面上位移不连续，表现为当结构变形时将在相邻单元间产生缝隙或重叠，这意味着将引起无限大的应变，这时必然会发生交界面上的附加应变能补充到系统的应变能中去，有限元解就不可能收敛于真正解。收敛单元尺寸趋于零时，有限元解趋于真解3:2113形函数的性质当单元的位移函数满足完备性要求时，称单元是完备的（通常较容易满足）。当单元的位移函数满足协调性要求时，称单元是协调的。当势能泛函中位移函数的导数是2阶时，要求位移函数在单元的交界面上具有C1或更高的连续性

9、，这时构造单元的插值函数往往比较困难。在某些情况下，可以放松对协调性的要求，只要单元能够通过分片试验 (Patch test)，有限元分析的解答仍然可以收敛于正确的解。这种单元称为非协调单元。分片试验由B.M.Irons首先提出，已经证明它给出了收敛性的充分条件。 3:2114单元应变和应力矩阵 0( , )( , )0 ( , )( , )( , )( , ) xyeexyxu x yx yx yx yx yv x yyyx uNdBd应变矩阵123123 0 0 0 0( , ) 0 0 0 0 xNNNx yNNNyyx BN3:2115单元应变和应力矩阵123123123112233

10、0 0 01( , )0 0 0 2 bbbx ycccAcbcbcbBBBB312112233112233 0 0 00 0 0 bbbccccbcbcbBBB 由于 与x、y无关，都是常量，因此B矩阵也是常量。单元中任一点的应变分量是B矩阵与单元节点位移的乘积，因而也都是常量。因此，这种单元被称为常应变单元。332211,cbcbcbA3:2116单元应变和应力矩阵21 0( , ) 1 0( , )11-0 0 2xxyyxyxyEx yx yD平面应力：( , )( , )( , )( , )eex yx yx yx yDDBdSd 123123( , ) x y SDBD BBBSS

11、S2 2 (1)1-1- 22iiiiiiiibcEbcAcbSDB应力矩阵平面应变：用平面应变弹性矩阵代入得到类似结果。3:2117单元应变和应力矩阵由于同一单元中的D、B矩阵都是常数矩阵，所以S矩阵也是常数矩阵。也就是说，三角形三节点单元内的应力分量也是常量。当然，相邻单元的E，A和bi、ci(i，j，m)一般不完全相同，因而具有不同的应力，这就造成在相邻单元的公共边上存在着应力突变现象。但是随着网格的细分，这种突变将会迅速减小。3:2118单元分析 6 1ed 3 1 3 6B几何关系位移函数 3 3D本构关系 3 1平衡关系 6 1F 3 6SDB6 6?eK单元刚度矩阵3:2119单

12、元应变能AAxyxyyyxxhdxdyhdxdyUT21)(21AAhdxdyhdxdyUDTT2121eAeAeehdxdyhdxdydDBBdDBdBdTTTT2121单元应变能 U为：TTTT)(DDeBdijmxyh注意到弹性矩阵D的对称性3:2120刚度矩阵引入刚度矩阵K：AehdxdyDBBKTeeeUdKdT21则：注意：hdxdy的实质是任意的微体积dv，于是得 Ke的一般式：VeVdTDBBK3:2121单元外力功单元受到的外力一般包括体积力、表面力和集中力。自重属于体积力范畴。表面力指作用在单元表面的分布载荷，如风力、压力，以及相邻单元互相作用的内力等。 ijmxyqVij

13、mxyqsijmxyfcVyVxVqqqsysxsqqqcycxcfff3:2122单元外力功（1）体积力所做的外力功ijmxyqVAVAVyVxVhdxdyhdxdyvquqWquTeuNdAVeVhdxdyWqNdTTVVeVVWdTTqNd3:2123单元外力功（2）面力所做的外力功lSlSySxShdlhdlvquqWquTeuNdlSeShdlWqNdTTtVeSWdTTqNdijmxyqsqs3:2124单元外力功（3）集中力所做的外力功cecWfdTijmxyfc当结构受到集中力时，通常在划分单元网格时就把集中力的作用点设置为节点。于是单元集中力fc的势能Vc为综合以上诸式，单元

14、外力的总外力功V为 eecVVeceVeVeCSVttttWWWWfdfqNqNdfdqNdqNdTTTTTTTTTddddcVVettfqNqNfddTT3:2125系统势能KddT121NeeUU扩充叠加TT12UW d Kdd fT1NeeWWd fTeeeW d f扩充叠加系统势能eeeeUdKdT213:2126单元刚度矩阵的扩充叠加T12310 0 0 0 0 00 0 0 0 0 00 0 01 0 0 020 0 eeeiiijimenneeejijjjmemimUKKKd KduuuuuKKKKK1231 00 0 0 0 0 0neenjmmuuuuuKmijmij单元编号

15、 ijmeeeeUdKdT213:2127单元等效节点载荷列阵的扩充叠加T12310 0eieennjemWfd fuuuuuffmijTeeeW d f单元编号 ijm3:2128能量原理和系统平衡方程TT12UW d Kdd f系统势能根据弹性力学能量原理：结构处于稳定平衡的必要和充分条件是总势能有极小值。00fKdfKd 上式是从能量原理导出的系统平衡方程。这个方程表达了节点力与节点位移之间的关系。3:2129刚度矩阵单元刚度矩阵：VeVdTDBBK11122122 1 -1( )( ) -1 1 eeeeeeTeeeeelKKE Ax E Ax dxlKKKBBTT eeeiiijim

16、eeeejijjjmAeeemimjmmhdxdytAKKKKB DBB DBKKKKKK1D:2D:系统刚度矩阵：1eNeeKK弹性矩阵D的对称性Ke对称K 对称3:2130刚度矩阵刚度矩阵K的详细内容为：ssNNjNiNNNjNjjjijjiNijiiiiNjiNjiNNKKKKKKKKKKKKKKKKKKKKKKKKKssssssssss21212122222211111211Ki、j是行列号，Ns为系统自由度数。3:2131刚度矩阵(1) 刚度矩阵中每个元素有明确的物理意义。例如，Kij表示当节点位移中第j个元素为1(dj=1)其余元素为零时，引起的单元力中的第i个节点力fi。把平衡方

17、程写开主对角线上元素Kii(i=1, Ns)恒为正值：位移和作用力同向ssssssssssssNjiNjiNNjNiNNNjNjjjijjiNijiiiiNjiNjifffffdddddKKKKKKKKKKKKKKKKKKKKKKKKK2121212121222222111112113:2132刚度矩阵(2) K的每一行或每一列元素之和为零以上式中第i行为例:iNiNjijiiiiifdKdKdKdKdKss2211f2i-1 =0 f2i =0f2j-1=0 f2j=0f2m-1 =0 f2m =0ijmxyijm11当所有节点沿x向或y向都产生单位位移时，单元作平动运动，无应变，也无应力，

18、因而单元结点力为零（不含初应力）。 所以有即，K的每一行元素之和为零。由于对称性，每一列元素之和也为零。021siNijiiiiKKKKK3:2133刚度矩阵(3) 系统刚度矩阵是奇异矩阵（ 即K的行列式为零）(4) 系统刚度矩阵是常量矩阵系统的节点力和节点位移成线性关系是基于弹性理论的结果。刚度矩阵是在系统处于平衡状态的前提下得出的。作用在它上面的外力必定是平衡力系。然而，研究系统平衡时没有引入约束。承受平衡力系作用的无约束系统，其变形是确定的，但位移不是确定的。所以出现性质（2）中的“平动问题”，即可以发生任意的刚体运动。从数学上讲，系统平衡方程的解不是唯一的或不能确定的。由此，系统刚度矩

19、阵一定是奇异的。单元刚度矩阵也一定是奇异的。3:2134位移边界条件的处理系统刚度矩阵是奇异矩阵，其物理原因是结构缺少刚性位移的约束，实际的工程结构都受有足够的支承约束，排除了发生任何刚体位移的可能性，因此，必须引入位移约束。有限元中，位移约束都设置在节点处。这里，只讨论刚性约束情况，即被约束的位移分量为零。设讨论的结构有Nn个节点，每个节点有ndf个自由度。则系统的总自由度为Ns，且dfnsnNN节点总位移列向量d中共包含Ns个分量T1sNidddd3:2135为了引入位移约束，把节点总位移列向量d分成两部分。一部分是不受约束的位移分量，记为df。另一部分是受刚性约束的位移分量，记为dr。不

20、失一般性，设1N号位移分量是不受约束的；N+1N+Nr共Nr个分量是受刚性约束的。即：rNNNNNdddddd2121,rfdd位移边界条件的处理3:2136位移边界条件的处理 111rsNNNrfddd显然不受约束的节点位移的总数N为 N=Ns-Nr 对方程中的刚度矩阵K和节点荷载向量列阵f也作相应分割，则得到rfrfrrrffrffffddKKKK式中，ff是已知力边界，fr是约束反力。rfddd3:2137位移边界条件的处理按矩阵乘法规则得frfrffffdKdK每个受刚性支承约束的位移分量都等于零，即0dr从而得到rfrfrrrffrffffddKKKKrrrrfrffdKdKffff

21、fdKrfrffdK3:2138位移边界条件的处理 Kff 引入约束后的约化的系统刚度矩阵。 这是一个非奇异 矩阵，它的逆矩阵Kff-1是存在的。fffffdK引入约束后的约化的系统平衡方程在分析计算时，从无约束的系统刚度矩阵K中删去与受约束位移号对应的行和列，再将矩阵压缩排列成NN阶方阵，即为约化后的结构刚度矩阵。rfrfrrrffrffffddKKKK3:2139位移边界条件的处理ssssssssssssNjiNjiNNjNiNNNjNjjjijjiNijiiiiNjiNjifffffdddddKKKKKKKKKKKKKKKKKKKKKKKKK21212121212222221111121

22、100010iiud iiuf 显然iNiuddds0101iiud 3:2140位移边界条件的处理111211111212222222121212sssssssssssijNijNiiiiijiNiijjjjijjjNNNNN iN jN NKKKKKdfKKKKKdfKKKKKdfdKKKKKdKKKKKsjNffiiud iiiiuKf显然iiiNiNiiiiuKdKdKdKss11iiud iiiiiiuKdKiiK3:2141节点位移和约束反力f1ffffKd通过求解平衡方程即可解出全部未知的节点位移：把解出的d代入未经修改的平衡方程，即可得到约束反力：关于上述方程的解算方法，一般不

23、采用求逆的方法求解，而是直接采用高斯消元法等求解线性方程组的方法求解求解。施加边界条件后，得到修改后的平衡方程fdK(未约化的)fffffdK(约化的)或fKd1fKd 或rfrffdK3:2142单元应变和应力( , )( , )xyexyx yx yBd( , )( , )( , )xyexyx yx yx yDDBd根据三角形节点的位移，求出单元应力应变为ijk如何求系统应变能和节点应力？3:2143有限元解的收敛性由于在有限元计算中引入了结构离散和位移模式，导致有限元计算结果和真实解的偏差。单元划分越小、位移模式取得越接近真实变形，解答越收敛于真实解。当单元的位移模式采用解析的位移解时

24、，有限元的计算结果和解析结果是相同的。然而，许多情况无真实的位移模式可以借用，只能寻求其近似函数，不可避免带来计算精度问题。实践证明：只要位移模式满足单元的完备性准则和协调性条件，就保证了有限元的解答收敛于真实解。3:2144有限元计算过程框图剖分结构为有限个单元，对节点、单元编号构建单元刚度矩阵和单元等价节点力向量组装系统刚度矩阵并引入约束组装整体等价节点荷载向量从节点平衡方程解未知节点位移计算结构应变、应力3:2145解综合方程Kd=f得结构节点位移d从d中找单元位移de用公式用公式=Bde和和 =D，计算应力计算应力应变应变把单元刚度矩阵组装成系统刚度矩阵K离散结构为若干单元建立单元刚度

25、矩阵Ke形成等价节点荷载f形成单元等价节点力有限元计算流程图3:2146关于三角形单元形函数的一点补充12223322332331111()()() 221xyAxyx yx yyy xxxyxy11eALA 2-3-P: 同样, 3-1-PA2 1-2-PA322eALA33eALA i,1 j, 2 m, 3 x y P A1 面积坐标的定义 在三角形内任意一点P定义 123 1LLL3:2147关于三角形单元形函数的一点补充 i,1 j, 2 m, 3 x y P A1 面积坐标与形函数的关系1111221121212ijjiiimmiiiiiijjjjjjmmmmmmxyxyab xc

26、 yxyLab xc yNLab xc yNLab xcyN1iijjmmiijjmmijmxx Lx Lx Lyy Ly Ly LLLL 面积坐标与直角坐标的关系3:21(0,0)A(0, 2)B(4,0)D(4,1)C1F q 考虑一个平面应力问题如图所示，假设厚度考虑一个平面应力问题如图所示，假设厚度h=1，材料为各项，材料为各项同性，杨氏模量为同性，杨氏模量为E=1，泊松比为，泊松比为=0，相关力和位移边界条，相关力和位移边界条件如图中所示，问题左端为固定约束。试用两个三角形单元件如图中所示，问题左端为固定约束。试用两个三角形单元分析此问题，三角形单元的网格划分如图所示。试求问题各分析

27、此问题，三角形单元的网格划分如图所示。试求问题各节点位移节点位移u、v和应力和应力x，y和和xy。1243 1 2例题3:2112323311212332132111223332231331211200000010000002ebbbyyyyyycccxxxxxxAcbcbcbxxyyxxyyxxyyB对于三角形单元，其对于三角形单元，其B B矩阵的表达式为：矩阵的表达式为：23322313211()()()2eAx yx yyy xxxy1243 1 2123TeehAkB DB2101001001011000.5002EvDFor e=1:1922472216.501215.5204020

28、141024116722492215.501216.5ek110201010400048410241eB1, (2)3, (1)2, (3)1, (2)2, (3)3, (1)例题3:2112323311212332132111223332231331211200000010000002ebbbyyyyyycccxxxxxxAcbcbcbxxyyxxyyxxyyB23322313211()()()2eAx yx yyy xxxy1243 1 2123TeehAkB DB2101001001011000.5002EvDFor e=2:216016404032032001601842414324330116002020404101ek200101010404004404101eB1, (3)3, (1)2, (4)1, (3)2, (4)3, (1)对于三角形单元，其对于三角形单元，其B B矩阵的表达式为：矩阵的表达式为：例题3:21组装刚度矩阵 12112722424217.5215.5410172922400215.5216.50100124202001641641410340