3平面问题的有限元法-Read.docx

《3平面问题的有限元法-Read.docx》由会员分享，可在线阅读，更多相关《3平面问题的有限元法-Read.docx（36页珍藏版）》请在三一办公上搜索。

1、3弹性力学平面问题的有限元法本章包括以下的内容：3.1弹性力学平面问题的基本方程3.2单元位移函数3.3单元载荷移置3.4单元刚度矩阵3.5单元刚度矩阵的性质与物理意义3.6整体分析3.7约束条件的处理3.8整体刚度矩阵的特点与存储方法3.9方程组解法3.1弹性力学平面问题的基本方程弹性力学是研究弹性体在约束和外载荷作用下应力和变形分布规律的一门学科。在弹性 力学中针对微小的单元体建立基本方程，把复杂形状弹性体的受力和变形分析问题归结为偏 微分方程组的边值问题。弹性力学的基本方程包括平衡方程、几何方程、物理方程。弹性力学的基本假定如下：1）完全弹性，2）连续，3）均匀，4）各向同性，5）小变形

2、。3.1.1基本变量弹性力学中的基本变量为体力、面力、应力、位移、应变，各自的定义如下。体力体力是分布在物体体积内的力，例如重力和惯性力。面力面力是分布在物体表面上的力，例如接触压力、流体压力。应力物体受到约束和外力作用，其内部将产生内力。物体内某一点的内力就是应力。图3.1如图3.1假想用通过物体内任意一点p的一个截面mn将物理分为1、11两部分。将部 分II撇开，根据力的平衡原则，部分II将在截面mn上作用一定的内力。在mn截面上取包含 p点的微小面积AA，作用于AA面积上的内力为kQ。令AA无限减小而趋于p点时，AQ的极限S就是物体在p点的应力。AQlin SAAt0 AA应力S在其作用

3、截面上的法向分量称为正应力，用。表示；在作用截面上的切向分量 称为剪应力，用T表示。显然，点p在不同截面上的应力是不同的。为分析点p的应力状态，即通过p点的各个 截面上的应力的大小和方向，在p点取出的一个平行六面体，六面体的各楞边平行于坐标轴。图3.2将每个上的应力分解为一个正应力和两个剪应力，分别与三个坐标轴平行。用六面体表 面的应力分量来表示p点的应力状态。应力分量的下标约定如下：第一个下标表示应力的作用面，第二个下标表示应力的作用方向。T.，第一个下标x表示剪应力作用在垂直于X轴的面上，第二个下标y表示剪应力指 向Y轴方向。正应力由于作用表面与作用方向垂直，用一个下标。x表示正应力作用于

4、垂直于X轴的面上，指向X轴方向。应力分量的方向定义如下：如果某截面上的外法线是沿坐标轴的正方向，这个截面上的应力分量以沿坐标轴正方向 为正；如果某截面上的外法线是沿坐标轴的负方向，这个截面上的应力分量以沿坐标轴负方向 为正。剪应力互等：T -T具 -T具 -Txy yx yz zy zx xz物体内任意一点的应力状态可以用六个独立的应力分量。、。、。、T 、T 、Tx y z xy yz zx来表示。位移位移就是位置的移动。物体内任意一点的位移，用位移在x, y, z坐标轴上的投影u、 v、w表示。应变物体的形状改变可以归结为长度和角度的改变。各线段的单位长度的伸缩，称为正应变，用表示。两个垂

5、直线段之间的直角的改变，用弧度表示，称为剪应变，用Y表示。物体内任意一点的变形，可以用、y、y、y六个应变分量表示。x y z xy yz zx3.1.2平面应力和平面应变问题弹性体在满足一定条件时，其变形和应力的分布规律可以用在某一平面内的变形和应力 的分布规律来代替，这类问题称为平面问题。平面问题分为平面应力问题和平面应变问题。 1）平面应力问题设有很薄的等厚薄板，只在板边上受到平行于板面并且不沿厚度变化的面力，体力也平 行于板面且不沿厚度变化。图3.3设板的厚度为t，J匚M，Z 2在板面上：G ）上=0z=士 2由于平板很薄，外力不沿厚度变化，因此在整块板上有，b= 0， T = 0，

6、T . = 0剩下平行于XY平面的三个应力分量。y、Txy未知。2）平面应变问题设有很长的柱形体，支承情况不沿长度变化，在柱面上受到平行于横截面而且不沿长度 变化的面力，体力也如此分布。图3.4以柱体的任一横截面为XY平面，任一纵线为Z轴。假定该柱体为无限长，则任一截 面都可以看作对称面。由对称性，T = 0，T . = 0，W = 0由于没有Z方向的位移，Z方向的应变 z = 0。未知量为平行于XY平面的三个应力分量。，、T.，物体在Z方向处于自平衡状 态。3.1.3平衡方程弹性力学中，在物体中取出一个微小单元体建立平衡方程。平衡方程代表了力的平衡关 系，建立了应力分量和体力分量之间的关系。

7、对于平面问题，在物体内的任意一点有，（3-1）x + + X = 0 dxdy8。8t+ x + Y = 0 dy8x3.1.4几何方程由几何方程可以得到位移和变形之间的关系。对于平面问题，在物体内的任意一点有，dudx（3-2）dvdy=du + dvxydydx刚体位移由位移u=0，v=0可以得到应变分量为零，反过来，应变分量为零则位移分量不为零。 应变分量为零时的位移称为刚体位移。刚体位移代表了物体在平面内的移动和转动。坐=0dxdv 八由 H = 0dydu dv 八+= 0dy dx可以得到刚体位移为以下形式，u = u -Jv = v +必du 八8v 八 由瓦= dy=可得将f

8、f代入疽虱=可得，df1( y) df2( x)1=2= w dy dx积分后得到，f1( y) = u 0 wyf2 (x) = v + wx得到位移分量，u = u wyv = v + wx当u 0。0, v0 = 0, w= 0时，物体内任意一点都沿x方向移动相同的距离，可见u 0代表物体在x方向上的刚体平移。当u0 = 0, v0。0,w= 0时，物体内任意一点都沿y方向移动相同的距离，可见v0代表 物体在y方向上的刚体平移。当u0 = 0, v0 = 0,。0时，可以假定w 0，此时的物体内任意一点P(X，y)的位移分量为，u = wy, v = wxP点位移与y轴的夹角为a对=竺=

9、z=妙 wx xO O I z _ x O /一、勺 I 勺 c u 2 + v 2 =.(Yy )2 + (x)2 =； x2 + y 2 = srr为P点到原点的距离，可见代表物体绕z轴的刚体转动。3.1.5物理方程弹性力学平面问题的物理方程由广义虎克定律得到。1)平面应力问题的物理方程1L)b 呻/1 ()& 一 日。/2(1+) T(3-3)Yxy平面应力问题有，告 +b y)2)平面应变问题的物理方程1 一 |L12E1 一 |L12E(3-4)2(1+1) T平面应变问题有，Yxy+ b)在平面应力问题的物理方程中，将替换为左、口替换为己，可以得到平面应变问题的物理方程；在平面应变

10、问题的物理方程中，将E替换为E （1 + 2 日）（1 +日）2、日替换为赤，可以得到平面应力问题的物理方程。图3.5求解弹性力学平面问题，可以归结为在任意形状的平面区域Q内已知控制方程、在位 移边界su上约束已知、在应力边界S上受力条件已知的边值问题。然后以应力分量为基本 未知量求解，或以位移作为基本未知量求解。如果以位移作为未知量求解，求出位移后，由几何方程可以计算出应变分量，得到物体 的变形情况；再由物理方程计算出应力分量，得到物体的内力分布，就完成了对弹性力学平 面问题的分析。3.2单元位移函数根据有限元法的基本思路，将弹性体离散成有限个单元体的组合，以结点的位移作为未 知量。弹性体内

11、实际的位移分布可以用单元内的位移分布函数来分块近似地表示。在单元内 的位移变化可以假定一个函数来表示，这个函数称为单元位移函数、或单元位移模式。对于弹性力学平面问题，单元位移函数可以用多项式表示，u = a + a x + a y + a x2 + a xy + a y2 +.123456（3-5）v = b + b x + b y + b x2 + b xy + b y2 +.123456多项式中包含的项数越多，就越接近实际的位移分布，越精确。具体取多项，由单元形式来确定。即以结点位移来确定位移函数中的待定系数。图3.6如图3.6所示的3结点三角形单元，结点I、J、M的坐标分别为区,七）、（

12、x_, y）、（x , y ），结点位移分别为u、v、u、v、u、v 。六个节点位移只能确定六个多 m mi i j j m m项式的系数，所以3结点三角形单元的位移函数如下，u = a + a x + a y 123 卜（3-6）v = a + a x + a y J将3个结点上的坐标和位移分量代入公式（3-6）就可以将六个待定系数用结点坐标和位移分量表示出来。将水平位移分量和结点坐标代入（3-6）中的第一式，u = a + a x + a yi 12 i 3 iu = a + a x + a yj 1 2 j 3 ju = a + a x + a ym1 2 m3 m（3-7）写成矩阵形式

13、，r 1u1 x yr 、aiii11 u =1x y1a jjj2u1 x yam mml 3JxiXjXmyiyjym=Et ，则有（3-8）a1 11 = V 3 JT -1T*|T| = 2A，A为三角形单元的面积。T的伴随矩阵为，x yx yyyxxj mmjjmmjx yxyyyxxm iimmiimxyx yyyxxi jjiijjijmaibiciTaiT*=abc=bjjjiabccmmmir r1r 、aaaau11_ijmia卜=bbb 0；当三个结点i，j11X yN = 2A (a + b x + c y)=(a 2 - ax - ay) = 1 - 形态矩阵为x 0

14、 y 0 1 -X-z0N =aaa a0 x 0 z 01 -X-y_ aaa a _三角形面积的计算公式可得，1Xy1a0,11ii_ 1101A =Xya=a 22jj=221Xy100m m如果把三个结点按顺时针方向排列，即i (a，0), j (0, 0), m (0, a)1 x y1 a 0,1ii_ 11A =1 x y1 0 0=a22jj=221 x ymm1 0 a有限元法的求解对象是单元的组合体，因此作用在弹性体上的外力，需要移置到相应的 结点上成为结点载荷。载荷移置要满足静力等效原则。静力等效是指原载荷与结点载荷在任 意虚位移上做的虚功相等。单元的虚位移可以用结点的虚

15、位移*表示为，(3-15)N *X、iYi令结点载荷为耘。XjYjXmYm1）集中力的移置如图3.7所示，(3-16)（Xp，yp）上作用有集中载荷。由虚功相等可得，（8 牛）Re =（8 牛）N T P由于虚位移是任意的，则Re = NT Pr x.iN.10 一Y0N.iXN0J j =Y0NXmNm0Ym0NmPxPy(xp，yp)例题1 ：在均质、等厚的三角形单元ijm的任意一点p2）体力的移置*、，I x令单兀所受的均匀分布体力为 p = j）由虚功相等可得，(8 *e)Re =JJ(8 *)NTptdxdyRe = jj N t ptdxdy3）分布面力的移置设在单元的边上分布有面

16、力P= X,Yt，同样可以得到结点载荷，Re =j NT P tdss(3-17)(3-18)XiYI X =uiviujvjumvm（3-19）记为 = B8 e，B矩阵称为几何矩阵。B矩阵可以表示为分块矩阵的形式B= B BjBm Ibi00cicibi叫2A由物理方程，可以得到单元的应力表达式?= D R= D B e（3-20）（3-21）D称为弹性矩阵，对于平面应力问题，1旦0旦10001 -D=昂定义S = D &为应力矩阵。将应力矩阵分块表示为，S =根据单元的位移函数，uu 1N 0 N 0 N 0 iI七uJ=i j m0 N 0 N 0 Nj vijmmv .m由几何方程可

17、以得到单元的应变表达式，S=仿止=Ei i2 A(1 - 2)bipbi1 pc2 icici或b2 i(3-22)应用虚功原理可以建立单元结点位移与结点力的关系矩阵，单元刚度矩阵。虚功原理：在外力作用下处于平衡状态的弹性体，如果发生了虚位移，则所有外力在虚 位移上做的虚功等于内应力在虚应变上做的虚功。单元的结点力记为F上=U 匕 uj七 Um 匕1单元的虚应变为* =旧 *单元的外力虚功为，I ） F L单元的内力虚功为，JJ *T j tdxdy由虚功原理可得，t * / F e 品* T Jtdxdy * I =而氏 * ) T 4 )附3=糖e =D B 俱*)Fe =*)JJB)DB

18、tdxdyd e fe = JJB)DBtdxdyE e(3-23)e（3-24）定义K1 =JJB）DBtdxdy为单元刚度矩阵。在3结点等厚三角形单元中B和D的分量均为常量，则单元刚度矩阵可以表示为，（3-25）K1 = B ) D B tA单元刚度矩阵表示为分块矩阵:Ks = B )DB Kr =Krx ,sxKKrx ,syK（3-26）（3-27）ry ,sx3.5单元刚度矩阵的性质与物理意义（一）单元刚度矩阵的物理意义假设单元的结点位移如下：= 1 0 0 0 0 0rU iVI U： 宇 jUI俨 m由F e = K e 5 e，得到结点力如下:ix ,ixK .iy，仄K（3-

19、28）K .jy，ixKmx ,ixKI my ,ixK表示i结点在水平方向产生单位位移时，在结点i的水平方向上需要施加的结点力。K y 表示i结点在水平方向产生单位位移时，在结点i的垂直方向上需要施加的结点力。选择不同的单元结点位移，可以得到单元刚度矩阵中每个元素的物理含义：K表示s结点在水平方向产生单位位移时，在结点r的水平方向上需要施加的结点力。Ky以表示s结点在水平方向产生单位位移时，在结点r的垂直方向上需要施加的结点力。K y表示s结点在垂直方向产生单位位移时，在结点r的水平方向上需要施加的结点力。Ky 表示s结点在垂直方向产生单位位移时，在结点r的垂直方向上需要施加的结点力。因此单

20、元刚度矩阵中每个元素都可以理解为刚度系数，即在结点产生单位位移时需要施 加的力。（二）单元刚度矩阵的性质1）对称性利用分块矩阵的性质证明如下：Kr = B tDB K=B tDB Kt= (B tDB )t =BtDtB = BtDB = Ksrsrrsrsrs即Ke = (Ke )T2）奇异性即单元刚度矩阵的行列式为零，|K|e = 0。将定单元产生了 X方向的刚体移动，8 = 1 0 10 1 0T，此时对应的单元结点力为零。01000113 = K e i100100由行列式的性质可0 1T，可以得可以得到，在单元刚度矩阵中1，3, 5列中对应行的系数相加为零，知，Ke=0。同样如果假定

21、单元产生了 y方向上的刚体位移8e = 0 1 0 1 到，在单元刚度矩阵中2, 4, 6列中对应行的系数相加为零。得到了单元刚度矩阵后，要将单元组成一个整体结构，根据结点载荷平衡的原则进行分 析，即整体分析。整体分析包括以下4个步骤：1）建立整体刚度矩阵，2）根据支承条件修改整体刚度矩阵，3）解方程组，求出结点的位移，4）根据结点位移，求出单元的应变和应力。在这里把结点位移作为基本未知量求解。如何得到整体刚度矩阵？基本方法是刚度集成法，即整体刚度矩阵是单元刚度矩阵的集 成。图3.9如图3.9所示，一个划分为6个结点、4个单元的结构。得到了每个单元的单元刚度矩 阵后，要集成为整体刚度矩阵。3.

22、6.1刚度集成法的物理意义由单元刚度矩阵的物理意义可知，单元刚度矩阵的系数是由单元结点产生单位位移时引 起的单元结点力。在如图3.9所示的结构中，使结点3产生单位位移时，在单元（1）中的结点2上引起 结点力。由于结点2、3同时属于单元（1）、（3），在单元（2）中的结点2上同样也引起结 点力，因此，在整体结构中当结点3产生位移时，结点2上的结点力应该是单元（1）、（2） 在结点2上的结点力的叠加。刚体集成法即结构中的结点力是相关单元结点力的叠加，整体刚度矩阵的系数是相关单 元的单元刚度矩阵系数的集成。结点3在整体刚度矩阵的对应系数，应该是单元（1）、（3）、 （4）中对应系数的集成。3.6.2

23、刚度矩阵集成的规则1）将单元刚度矩阵中的每个分块放到在整体刚度矩阵中的对应位置上，得到单元的扩 大刚度矩阵。单元刚度矩阵系数取决于单元结点的局部编号顺序，必须知道单元结点的局部编号与该 结点在整体结构中的总体编号之间的关系，才能得到单元刚度矩阵中的每个分块在整体刚度 矩阵中的位置。将单元刚度矩阵中的每个分块按总体编码顺序重新排列后，可以得到单元的 扩大矩阵。假定单元结点的局部编号与整体的对应关系如下：单元编号单元结点局部编号单元结点整体编号1i31j11m22i52j22m43i53j33m24i34j54m6单元（2）的单元扩大矩阵K）的分块矩阵形式如下，只列出非零的分块:局部编号jm1i1

24、整体编号123456123456K (2)K(2)K (2)K (2)Kmm 3K (2)K J(2)K (2)K (2)2）将全部单元的扩大矩阵相加得到整体刚度矩阵。K = K + K （2）+ K （3）+ K （4）整体刚度矩阵如下所示：整体编号123456123456K (1) jjK(i)jmK_ (1)Kd)Km皿+ K + jKJ+ K (3)K m (2)K (2)+ K(3)K_ (i)Km (DK m (3)K (1)+ K )(3)+ K：(4)K (3)+ K (4)K (4)K (2)孔m)K .(2)K (2)+ K(3)K (3)+ K (4)K (2)K .(2

25、) + K + K (4)K . )(4)K (4)K (4)Km叫图3.9所示的结构的约束和载荷情况，如图3.10所示。结点1、4上有水平方向的位移 约束，整体刚度矩阵K求出后，结构上的结点力可以表示为：F = K出根据力的平衡，结点上的结点力与结点载荷或约束反力平衡。用P表示结点载荷和支杆反力，则可以得到结点的平衡方程：K G= P（3-29）这样构成的结点平衡方程组，在右边向量P中存在未知量，因此在求解平衡方程之前， 要根据结点的位移约束情况修改方程（3-29）。先考虑结点n有水平方向位移约束，与n结点 水平方向对应的平衡方程为：K u + K v +. + K u + K v +. =

26、 P（3-30）2 n-1,1 12 n-1,2 12 n-1,2 n-1 n2 n-1,2 n n2 n-1根据支承情况，方程（3-30）应该换成下面的方程：u = 0（ 3-31）对比公式（3-30）和（3-31），在式（3-29）中应该做如下修正：在K矩阵中，第2n-1行的对角线元素K2n-12n-1改为1，该行中全部非对角线元素改为0；在P中，第2n-1个元素改为0。为了保持K矩阵的对称性，将第2n-1列的全部非对角 元素也改为0。同理，如果结点n在垂直方向有位移约束，则（3-29）中的第2n个方程修改为，v = 0在K矩阵中，第2n行的对角线元素改为1，该行中全部非对角线元素改为0；

27、在P一0 0一r u1 P 110 0vP120 0uP230 0VP240 00 0卜=4000000100000Pn2 n-1000000010000VPn2n0 00 00 0_0 0_I JI J中，第2n个元素改为0。为了保持K矩阵的对称性，将第2n列的全部非对角元素也改为0。(3-32)100000000000-*00*0*00*0*00*0*00*0*00*0Et100000210000对称*0*0*01_对图3.9所示结构的整体刚度在修改后可以得到以下的形式，(3-33)如果结点n处存在一个已知非零的水平方向位移匕，这时的约束条件为，u = u*（3-34）在K矩阵中，第2n-

28、1行的对角线元素K2ni2n1乘上一个大数A，向量P中的对应换成AKtmt：，其余的系数保持不变。方程改为，K u + K v +. + AK u + K v +. = AK u: （3-35）2 n-1,112 n-1,212n-1,2 n-1n2 n-1,2nn2 n-1,2 n-1nA的取值要足够大，例如取1010只有这样，方程（3-35）才能与方程（3-34）等价。如果结点n处存在一个已知非零的垂直方向位移匕，这时的约束条件为，vv*。也可以采用同样的方法修改整体刚度矩阵。3.8整体刚度矩阵的特点与存储方法用有限元方法分析复杂工程问题时，结点的数目比较多，整体刚度矩阵的阶数通常也是 很

29、高的。那么，是否在进行计算时要保存整体刚度矩阵的全部元素？能否根据整体刚度矩阵 的特点提高计算效率？整体刚度矩阵具有以下几个显著的特点：对称性，稀疏性，非零系数带形分布。1）对称性由单元刚度矩阵的对称性和整体刚度矩阵的集成规则，可知整体刚度矩阵必为对称矩阵。利用对称性，只保存整体矩阵上三角部分的系数即可。2）稀疏性单元刚度矩阵的多数元素为零，非零元素的个数只占较小的部分。如图3.11所示的结构， 结点2只和通过单元联接的1、3、4、5结点相关，结点5只和通过单元联接的2、3、4、6、8、9结点相关。由单元刚度矩阵的物理意义和整体刚度矩阵的形成方式可知，相关 结点2、3、4、6、8、9及结点5本

30、身产生位移时，才使结点5产生结点力，其余结点产生位移时不在该结点处引起结点力。在用分块形式表示的整体矩阵中，与相关结点对应的分块矩阵具有非零的元素，其它位置上的分块矩阵的元素为零，如图3.12所示。图3.12整体刚度矩阵的分块矩阵示意3）非零元素带形分布在图3.12中，明显可以看出，整体刚度矩阵的非零元素分布在以对角线为中心的带形 区域内，这种矩阵称为带形矩阵。在包括对角线元素的半个带形区域内，每行具有的元素个数叫做半带宽，用hbd表示。d =（相邻结点的编码的最大差值+1） x 2图3.11所示结构的相邻结点编码的最大差值为4,所以半带宽为10。二维等带宽存储设整体刚度矩阵K为一个nx n的

31、矩阵，最大半带宽为d。利用带形矩阵的特点和对称 性，只需要保存以d为固定带宽的上半带的元素，称为二维等带宽存储。进行存储时，把整 体刚度矩阵K每行中的上半带元素取出，保存在另一个矩阵K*的对应行中，得到一个 n x d 矩阵K*。把元素在K矩阵中的行、列编码记为r、s，在矩阵K*中的行、列编码记为r*、s*，对 应关系如下：r*=r如图3.13(a)所示的最大半带宽为d的整体刚度矩阵K，采用二维等带宽存储后得到如 图3.13(b)所示的矩阵K*。用新的方法存储后，K矩阵中的对角线元素保存在新矩阵中的 第1列中，K矩阵中的r行元素仍然保存在新矩阵的r行中，K矩阵中的s列元素则按照 新的列编码保存在新矩阵的不同列中。采用二维等带宽存储，需要保存的元素数量与K矩阵中的总元素数量之比为-。所存 n储的元素数量取决于最大半带宽d的值，d的值