4第4章 平面问题有限元法

第四章 平面问题有限元法第一节 弹性力学有关知识在有限元法中经常要用到弹性力学的基本方程，关于它们的详细 推导可从弹性力学的有关教材中查到。现将它们连同相应的矩阵和张 量表达式形式综合引述于后。一、弹性力学基本方程的矩阵形式弹性体在载荷作用下，体内任意一点的应力状态可由 6 个应力分 量Q Q Q ,T ,T ,T来表示。其中Q Q Q为正应力；T ,T ,T为剪应x y z xy yz zx x y z xy yz zx力。应力分量的正负号规定如下：如果某一面的外法线方向与坐标轴 的正方向一致，这个面上的应力分量就以沿坐标轴正方向为正，与坐 标轴反向为负；相反，如果某一个面的外法线方向与坐标轴的负方向 一致，这个面上的应力分量就以沿坐标轴负方向为正，与坐标轴同向 为负。应力分量及其正方向见图 4.1。图 4-1应力分量的矩阵表示称为应力列阵或应力向量，即r 、cxcyczTxyTyzTzxc c T T T 1yz xyyz zx4-1)弹性体在载荷作用下，还将产生位移和变形，即弹性体位置的移动和形状的改变。弹性体内任一点的位移可由沿直角坐标轴方向的 3 个位移分量u, v, w来表示，它的矩阵形式是=L v wT4-2)称作位移阵或位移向量。弹性体内任意一点的应变,可以由6个应变分量L ,l ,l ,Y ,Y ,Yx y z xy yz zx 为剪应变。应变的正负号为正应变；来表示。其中L ,L ,L为正应变；Y ,Y ,Yx y zxy yz zx与应力的正负号相对应，即应变以伸长时为正，缩短为负；剪应变是以两个沿坐标轴正方向的线段组成的直角变小为正，反之为负。应变的矩阵形式是CALxLyL zYxyYyzYzxLLLxyY y Y xyyzzx4-3)称作应变列阵或应变向量。1. 平衡方程弹性体V域内任一点沿坐标轴x, y, z方向的平衡方程为氏溯氏+ f = 0Xx + -Qx产+ QyzxQzQcQcQcxy +y +QxQyQzQcQcQc+QxQyQz+ f = 0 y+ f = 0 z4-4)其中f , f , f为单位体积的体积力在x, y, z方向的分量。且有 xyzT =T ,T =T ， P。xy yx yz zy zx xz平衡方程的矩阵形式为Ao+ f =0 (在V内）4-5)其中， A 是微分算子，ddxQQydQydQxQQzddzQQzQQyQQx4-6)f是体积力向量，f =yz2. 几何方程应变-位移关系在微小位移和微小变形的情况下，略去位移导数的高次幂，则应变向量和位移向量间的几何关系为QuQv8 , 8 -xQxyQyQw8 -zQzQuY -+xyQy=Yyx4-7)Qv 丫 yz 飞 +=YzyQu 丫 zx 飞 +=Yxz几何方程的矩阵形式为8x8y8zVxyVyzVzx8= LuV=TxyG xy1V=TyzG yz1V=TzxG zx（4-10）用矩阵形式表示为11 |LX1 |LX0O10x1 |LX1 |LXay100-VazE（1卩）s .1 |LX1 |LXT（1 + 卩）（1 - 2 卩）0001 2 |LXxyT2（1-卩）yzT0000zx0000000001 2p2（1-卩）000001 2p2（1-卩）VxyVyzVzxy8V z 尸 D 8（4-11）其中E (1-v)(1 +v)(l - 2v)1v1 -vv1 -v0v11 -v01000001 - 2v002 (1 -v)1 - 2v02(1-v)1 - 2v2(1-v)00(4-12)D称为弹性矩阵。它完全取决于弹性体材料的弹性模量E和泊松比v。表征弹性体的弹性，也可以采用拉梅(Lam e)常数g和九，它们和 E， v 的关系如下EvX = (1 +v)(1 - 2v)G 也称为剪切弹性模量。注意到2X+g=(i+v)(i - 2v)物理方程中的弹性矩阵D亦可表示为九+ 2G 九九+ 2GFl九九+ 2G000G0000G00000G(4-13)物理方程的另一种形式是e =c其中：C是柔度矩阵。C二D-1，它和弹性矩阵是互逆关系。弹性体V的全部边界为S。一部分边界上已知外力T,T,T称为力的边 xyz界条件，这部分边界用S表示；另一部分边界上已知位移，称为几何 a边界条件或位移边界条件，这部分边界用 S 表示。这两部分边界构成 u弹性体的全部边界，即S + S 二 Sau4、力的边界条件 弹性体在边界上单位面积的内力为T,T,T，在边界S上已知弹性体单x y za位面积上作用的面积力为T ,T ,T，根据平衡应有xyzT 二 T T = T T 二 Tx x y y z z设边界外法线的方向余弦为n ,n ,n，则边界上弹性体的内力可由下式xyz确定T 二 n a + n t + n txx xy yxz zxT 二 n t + n a + n tyx xyy yz zyx xzy yz z z上 式的矩阵形式为T二T （在S上）a其中0n0nyz0nn0xzn0nnzyx二 na0ny05、几何边界条件在S上弹性体的位移已知为u,v,w， u即有v =v w= w用矩阵形式表示：u = u（在S上）u以上是三维弹性力学问题中的一组基本方程和边界条件，同样，对于平面问题，轴对称问题等也可以得到类似的方程和边界条件。把弹性力学方程记作一般形式，它们是：平衡方程Ab + f 二 0 （在 V 内）几何方程Lu(在V内)物理方程b = De（在 V 内）边界条件nb = T（在 S 内）bu = u（在S上）u并有S + S二S，S为弹性体全部边界。b u对于不同类型问题，几何方程和物理方程的有关矩阵符号的意义汇集 于表。板与壳的基本方程将分别在本书有关章节中给出。6、弹性体的应变能和余能 单位体积的应变能(应变能密度)为U (e )二& TDE(4-14)2应变能是个正定函数，只有当弹性体内所有的点都没有应变时，应变 能才为零。单位体积的余能(余能密度)为V (b) = - b tCo(4-15)2余能也是个正定函数。在线性弹性力学中弹性体的应变能等于余能 二 弹性力学的平面应力问题与平面应变问题在实际问题中，任何一个弹性体严格地说都是空间物体，它所受 的外力一般都是空间力系。但是，当所考察的弹性体的形状和受力情 况具有一定特点时，只要经过适当的简化和力学的抽象处理，就可以 归结为弹性力学平面问题。平面问题分为平面应力问题和平面应变问题。（一）、平面应力问题 等厚度薄板，板边承受平行于板面并且不沿厚度变化的面力，同时体 力也平行于板面并且不沿厚度变化（如图4-2）。图42特点：1） 长、宽尺寸远大于厚度2） 沿板边受有平行板面的面力，且沿厚度均布，体力平行于板面且不沿厚度变化，在平板的前后表面上无外力作用。y图 4-3（二）、平面应变问题 很长的柱体，在柱面上承受平行于横截面并且不沿长度变化的面力，同时体力也平行于横截面并且不沿长度变化，如图4-4。E = 0 T = 0 T = 0z zx zy如：水坝、受内压的圆柱管道和长水平巷道等。xy图 4-4平衡微分方程为dodr_X +yx + X = 0dxdy(4-16)do drL +XL + Y = 0 dydx几何方程为du = _xdxdvyEyxYxy2(1 +卩)TE xyLl且有： = -(C +C )z E x y平面应变问题的物理方程(cx1l2(c1l(4-20)Yxy2(1 +1) TE xy第二节 平面问题有限元法一、有限元法的基本思想假想的把一连续体分割成数目有限的小体（单元），彼此间只在数目有限的指定点（结点）出相互连结，组成一个单元的集合体以代替原来的连续体（如图 4-5 所示），再在结点上引进等效力以代替实 际作用于单元上的外力。选择一个简单的函数来近似地表示位移分量 的4.2 有限元法的基本概念 点力之间的关系。离散连续体人为分割组合体人为假想分网meshing1* 1节点（n单元(Element、mesh)Finite Element MethodFEM图 4-5 有限元离散思想有限元法的实质是：把有无限个自由度的连续体，理想化为只有 有限个自由度的单元集合体，使问题简化为适合于数值解法的结构型 问题。二、 三角形常应变单元（一）、 离散化在运用有限单元法分析弹性力学平面问题时，第一步就是要对弹性体进行离散化，把一个连续的弹性体变换为一个离散的结构物。对 于平面问题，三角形单元是最简单、也是最常用的单元，在平面应力 问题中，单元为三角形板，而在平面应变问题中，则是三棱柱。假设采用三角形单元，把弹性体划分为有限个互不重叠的三角 形。这些三角形在其顶点（即节点）处互相连接，组成一个单元集合 体，以替代原来的弹性体。同时，将所有作用在单元上的载荷（包括 集中载荷、表面载荷和体积载荷），都按虚功等效的原则移置到节点 上，成为等效节点载荷。由此便得到了平面问题的有限元计算模型， 如图 4-6 所示。图 4-6 弹性体和有限元计算模型y图 4-7 平面三角形单元二)、 位 移首先，我们来分析一下三角形单元的力学特性，即建立以单元节 点位移表示单元内各点位移的关系式。设单元e的节点编号为i、j、 m,如图4-7所示。由弹性力学平面问题可知，每个节点在其单元平 面内的位移可以有两个分量，所以整个三角形单元将有六个节点位移 分量,即六个自由度。用列阵可表示为：匕e 二 L T 8 T 8二 U v u v u v (4 21)i j m i i j j m m其中的子矩阵=L vli i i式中u、v是节点i在x轴和y轴方向的位移。ii在有限单元法中,虽然是用离散化模型来代替原来的连续体,但 每一个单元体仍是一个弹性体,所以在其内部依然是符合弹性力学基 本假设的,弹性力学的基本方程在每个单元内部同样适用 从弹性力学平面问题的解析解法中可知,如果弹性体内的位移分量函 数已知,则应变分量和应力分量也就确定了。但是,如果只知道弹性 体中某几个点的位移分量的值,那么就不能直接求得应变分量和应力 分量。因此,在进行有限元分析时,必须先假定一个位移模式。由于 在弹性体内,各点的位移变化情况非常复杂,很难在整个弹性体内选 取一个恰当的位移函数来表示位移的复杂变化,但是如果将整个区域 分割成许多小单元,那么在每个单元的局部范围内就可以采用比较简 单的函数来近似地表示单元的真实位移,将各单元的位移式连接起 来,便可近似地表示整个区域的真实位移函数。这种化繁为简、联合 局部逼近整体的思想,正是有限单元法的绝妙之处。基于上述思想，我们可以选择一个单元位移模式，单元内各点的 位移可按此位移模式由单元节点位移通过插值而获得。线性函数是一 种最简单的单元位移模式，故设u 二a +a x+a y123v = a + a x+a y（4-22）456式中a、a、a是待定常数。因三角形单元共有六个自由度，且1 2 6位移函数 u、 v 在三个节点处的数值应该等于这些点处的位移分量的 数值。假设节点i、j、m的坐标分别为（x , y ）、（x , y ）、iijj（x , y ），代入 （4-22） 式，得：mmu a +a x +a yi 1 2 i 3 iu a +a x +a y j 1 2 j 3 ju a +a x +a ym 1 2 m 3 mv a +a x +a yj 45 i 6 i (4-23)v a +a x +a yj 4 5 j 6 jv a +a x +a ym 4 5 m 6 m由 （4-23） 式左边的三个方程可以求得1uixiyi1a uxy,a i 2Ajjj22Auxym m m1uy1 x u1ii1iiuy,a 1 x u1jj32Aj juy1 x umm/血、m(4-24)其中1 xi2A 1 xjyiyjym1 x从解析几何可知，式中的A就是三角形i、j、m的面积。为保证求 得的面积为正值，节点i、j、m的编排次序必须是逆时针方向，如图 4-3 所示。将 (4-24) 式代入 (4-22 式的第一式，经整理后得到12A12a+ bx + cy)U +C + bx + cy 1 +(a + bx + c yL i ii ij jj jmmmm+ (ajm+ b x + c y)v + a + b x + c yi i i i j j j+ b x + c y)vm m m12A(ai+bx+cy i(4-25)i,j,m 轮换 ) ( 4-26)xa = ji xmyjym1b =-i 11c =i 1yjym二 y - yjm-xjm4-27)u = Nu + N u + Nui i j j m mv = N v + N v + N vi i j j m m4-28)也可写成矩阵形式NIjN I=|N fcem4-29)式中I是二阶单位矩阵；Ni 、Nj 、Nm是坐标的函数，它们反 映了单元的位移状态，所以一般称之为形状函数，简称形函数。矩阵N 叫做形函数矩阵。三节点三角形单元的形函数是坐标的线性函 数。单元中任一条直线发生位移后仍为一条直线，即只要两单元在公 共节点处保持位移相等。则公共边线变形后仍为密合。(三) 应 变有了单元的位移模式，就可以利用平面问题的几何方程x = yYJ 2求得应变分量。将 (4-25)式代入上式，即得：b0b0s=1i0cj0c2Aijcibicjbjbm0cm0CMembm4-30)可简写成S=b 毗(4-31)其中可写成分块形式B 矩阵叫做单元应变矩阵b = t B B而子矩阵b Libi02Aci0cibi(4-32)由于A和bi、匕、bm 、ci 、 cj 、 cm 等都是常量，所以矩阵ij cmB中的诸元素都是常量，因而单元中各点的应变分量也都是常量,通常称这种单元为常应变单元。（四）应 力求得应变之后，再将（4-31）式代入物理方程便可推导出以节点位移表示的应力。即fc= Id 聪令s =d JB D脸则（4-33）其中S叫做应力幷阵，若写成分块形式，有sms Lb JB B B=k Si j m i j对于平面应力问题，弹性矩阵D为1对Id =E1 |LX 2卩1 称1 - u0 0 HL2所以，S的子矩阵可记为b|HctsLbBL( e)血C4-34)11 2、-卩2厶 ii_也c b2 i 2 i由（4-33）、（4-34）式不难看出，S中的诸元素都是常量，所以每个单元中的应力分量也是常量。可见，对于常应变单元，由于所选取的位移模式是线性的，因而 其相邻单元将具有不同的应力和应变，即在单元的公共边界上应力和 应变的值将会有突变，但位移却是连续的。（五）有限元理论模型1. 单元有限元理论模型为了推导单元的节点力和节点位移之间的关系，可应用虚位移原 理对图 4-2 中的单元 e 进行分析。单元 e 是在等效节点力的作用下处于平衡的，而这种节点力可采用列阵表示为RTjRTmUVjj4-35)假设在单元 e 中发生有虚位移，则相应的三个节点 i、 j、 m 的虚位移为uSvSuSvSuiijjmSv 1m4-36)且假设单元内各点的虚位移为f *,并具有与真实位移相同的位移 模式。所谓虚功原理，即如果内力与外力在物体发生任意位移或虚位移过程，所作的内功和外功相等，则物体上所有的点处于平衡状态。故有4-37)参照(4-31)式，单元内的虚应变* *为4-38)于是，作用在单元体上的外力在虚位移上所做的功可写为( *e)T r(4-39)而单元内的应力在虚应变上所做的功为4-40)这里我们假定单元的厚度t为常量。把(4-33 )式及(4-38)式代 入上式，并将提到积分号的前面，则有(8*e)t口 bt dBetdxdy(4-4D根据虚位移原理，由(4-39)和(4-41)式可得到单元的虚功方程， 即(8*e)tRe = (8*e)t口 b 1 dBbhdxdy注意到虚位移是任意的，所以等式两边与相乘的项应该相等，即得Reb 1 dBdxdy (4-42)记k 】=b d B dxdy(4-43)则单元的有限元方程有Re = k e 8 e(4-44)上式就是表征单元的节点力和节点位移之间关系的刚度方程， ke 就是单元刚度矩阵。如果单元的材料是均质的，那么矩阵D中的元 素就是常量，并且对于三角形常应变单元，B矩阵中的元素也是常 量。当单元的厚度也是常量时，因，所以（4-43）式可以简化为ke =BT DBtA（4-45）2、整体有限元方程K5=R（4-46）有限元计算方法是将位移作为基本未知量，通过边界条件和总体有限元理论模型计算位移，再计算应变和应力。其解题思路如图4-8所示。边界条件V有限元总体理 论模型方程应变应力图 4-8 有限元求解思路（六）CAE有限元分析流程有限元分析流程分三个阶段：前处理计算求解后处理6.1 有限元建模概述图 4-9 有限元分析流程图一有限元分析的三个阶段