《有限元分析第五章.docx》由会员分享，可在线阅读，更多相关《有限元分析第五章.docx（15页珍藏版）》请在三一办公上搜索。
1、5-5数值积分1、问题的提出在上一节中对等参元进行单元分析时要进行下列积分：(i) 单元刚1E止dxdye=j j B k E 】B t| det J|d&d门(5-4-5)-1 -1(ii) 体积力的等效结点力r = JJljVy J f tdoV e1U=jj N:fX jt |det J|d&d 门-1 -1I yJ(iii)边界力的等效结点力rr = j InL jpjtds=j f Qxjj f (y )iy-111)(5-4-7)(iv)温升载荷的等效结点力玄=jjB上& tdetJdRd(5-4-8)T0-1 -1式(5-4-5)(5-4-8)分别归结为计算以下两种形式的积分j
2、f (x)dx、-1j j f (x, y )dxdx- 1 - 1对于上述积分仅在单元的形状十分规则的情况下才能得到解析的结果(精确值)，一般 情况只能用数值积分方法(主要是高斯求积法)求近似值。虽然数值积分是“被迫“采用的， 但后来发现：有选择地控制积分点的个数和位置，可以方便地实现我们的某些特殊意图。 这样一来，数值积分就成为有限元分析的一个重要组成部分，以至本来可以精确积分的三角 形单元也常常采用数值积分。2、数值积分的基本概念任何积分工作取决于三个要素：给定的积分区间，给定的被积函数，具体的积分方法。 下面以一维情况为例介绍数值积分的基本概念(i)梯形法函数f (x)在区间(a,b)
3、的积分可以表达为I = j f (xdx Ew f (x )i ii=1W :权系数;i七:积分样点；f G)：积分样点的函数值。i梯形法的求积公式为i =，G d 火竺&2 h ai=1其中，h = -a，而UW = b - a n 一 1ii=1axih xi+1b图 5-22(ii)当被积函数为n-1次多项式Pn1(x)时，则由n个样点及其样点值(必提习川二项)可以 精确重构这个多项式，从而可以得到精确解。(注意：用以确定多项式的样点不必刻意选取)对多项式形式的被积函数进行积分可以采用高斯求积法，高斯求积法中要利用Legendre多项式的性质，下面先介绍这种多项式及相关的性质。3、Leg
4、endre 多项式Legendre多项式的定义域为-1,1L (x) = x一阶 Legendre多项式二阶 Legendre多项式31l2(x) = 2(x2 - 3)三阶 Legendre 多项式v15x = 0,15四阶 Legendre多项式1、4)(x2 -15 -、1203535，15 - 顽x = + -2、335一般n阶Legendre多项式的定义为Ln =A读者可以验证：n阶Legendre多项式是二阶变系数齐次微分方程 (-x 2* 2 W + n(n + 1)y = 0dn (x2 1)n dxn在区间-1,1上的有界解。Ln(x)在区间(-1,1)上有n个相异实根(零点
5、) 若再补充定义L (x) = 1则得一个定义在-1, 1上的多项式系列L0(x)、Lx) L(x)、.L (x)、.图5-23给出了 Legendre多项式函数列中的前四个元素。为了进一步讨论这个多项式 的性质，先注意下面一个事实：d (x 2 - 1) n = (x 2 1) n 一1 n - 2 x = (x 2 - 1) n-1 p (x)d 2(x2 -1)n = (x2 -1)n-1 n 2 + (x2 -1)n-2 n(n 1) 4x2 dx 2=(x2 1)n-2 P (x)对于任何kn)j L (x) L (x)dx =j- (x2 -1) n m (x2 -1) mdxm
6、n2m+n m!n! dxndxm_1dm-1dn=(x 2 - 1) m ( x 2 - 1) n2 m+n m!n! dxm-1dxnd m-1-11 d m - 1-(x 2 - 1)dx m-1-11d m - n - 12m+n mn(T)瓦二1(x2-1)mdn+1 ,m (x 2 - 1) ndxdxn+1当m=n时则有J L 2 (x)dx = (-1) m f (x 2 - 1) mdx-1-1J (1 - x2)mdx22 m (m!)2-12(2) 则有2m +1若在(1)的证明中将Ln(x)换成任何次数不超过m-1次的多项式Pm_(x)1 L (x)P 1 (x)dx
7、= 0-1这表明：Lm(x)与任何一个次数不超过m-1的多项式正交。(3)若q(x)是(一1,1)上平方可积的函数，则可将q(x)展开成q (x) = C L (x)其中系数m=0C = 2f L (x)q(x)dx-1特别，对于n次多项式Pn(x)有P (x) = c -L (x)4、一维情况设需要计算积分m mm=0I = J f (x)dx-1我们可以取x1=0为积分点(图5-2 4 (a),以常量f(0)代替f(x)进行积分，作为I的近似值I J f (0)dx = 2 f (0) = Wf (0)-1当川;J是一次函数时，可得到I的精确值。也可以任取两个点、为积分点(图5-2 4(b
8、)，用一个线性插值函数f (x ) - / 3 ) Z 、f (x ) +2(X - X )21代替f(x)进行积分，作为I的近似值I J f(X )+ f (X2) f (X1)dxX X。2 x 。f (X )f (X )X X 1 X X 2=Wf (X；) + W2 f (x：) 112X2当f(x)是一次函数时显然可得到I的精确值(取一个积分点已经能做到这一点)，当f(x)是二 次、三次函数时又将如何？一般若取任取n个积分点X、X2、 Xn，作n-1次插值多项式，积分I的近似值可表示 为I n W f (X )(5-5-1)i i i=1其中f(x)为积分点上的函数值，W|为权系数。
9、当积分点取为n阶Legendre多项式的零点时， 如果被积函数f(x)为次数不超过2n1次的多项式，(5-5-1)将给出积分的精确值。这就 是高斯求积分法，上述积分点又称为高斯点。高斯点的个数又称为积分阶数，有限元分析 中一般n=24。可以证明如下一般来说，对于一个2n-1阶的多项式P2n1(X)，需用2n个样点及其样点值才能精确重 构该多项式，或者说，需用2n个积分样点才能给出精确积分。若用、(X)除以多项式P2n 1(X)则pG)八 r)r G)L C)= q + 七 &)即 P2 1(X)= Q 1 (x)L (x)+ R 1(X)所以有j P (X X = j Q Ql(X )dx +
10、 j R (x 认=j R (x 认1 nn1n1111Rn1为同阶多项式，2 n1一一一因此仅需n个样点及其样点值即可精确重构该多项式，进而给出精确的积分值，若取Ln(X)的n个零点为积分样点，则P2 1(X )= R 1(X ) (i = 1 n)结论：用n个Legendre多项式的零点作为积分样点，式(5-5-1河以给出精确的积分值，这种方式的积分称为高斯积分。可见，高斯积分方法即减少了积分样点数也优化了积分样点。现将常见的高斯点坐标和权系数列成下表：积分阶数n高斯点坐标权系数2n 11x1=0W1=212x1 2= + 0.5773502692W12=133x1 3= + 0.7745
11、966692 x2=0W13=5/9W1=8/954x1 4= + 0.3611363116 x2，3= + 0.3399810436W14=0.3478548451W；3=0.652145154975、二维情况二维情况（三维情况与此类似）积分点的选择有更大的灵活性。一种经常采用的（并非 唯一可能的）选择方式是：沿心J方向取同样个数的积分点，积分的近似表达式为(5-5-2)I = j j f (尤,y)dxdy 工切 WW f (x , y )-1 -1i= j=1积分阶数n是对每一个自变量而言的积分阶数，而积分点总数在二维情况下为n2,在三维情况下为n3。图 5-2 5单元刚度公式（5-4-
12、5）的数值积分形式将是(5-5-3)k=空 n ww. &压B-1 - det J L m )i=1 j=1其中（与n ）为积分点的自然坐标。对于每个积分点都必须将（5-4-1）（5-4-5）执行一遍！ 注意:B矩阵在积分点上的函数值是确定的，因为它与形函数相关，而形函数是一个确定 的函数。此外，上述公式均是矩阵运算对八结点单元而言k共有162=256个元素，利用对 称性仍需对其中的136个元素进行数值积分。在单元结点个数增加，积分阶数也随之提高的 情况下，计算量是相当可观的。对三维情况尤其如此。6、有限元解的误差有限元解的误差是四个方面因素综合影响的结果。（i）插值误差 这是在单元内用多项式
13、代替真实解（在整个求解域内则表现为用有限自 由度代替了无限自由度）所引起的。这个因素在协调位移单元中将导致“刚度偏大”。（ii）边界形状以及边界条件的误差即使采用了曲边单元，单元边界仍有它本身的特点（例如是某个自然坐标的二次函数）， 不可能作到与实际曲线边界完全吻合。边界形状的误差使得实际边界条件不能得到精确满 足。这种误差一般只在边界附近影响较大（奇点除外）。(iii) 数值积分误差 这种“误差”如果利用得好，可以与(i)弓I起的误差“抵消”，但 处理不当也将影响解的精度。(iv) 截断误差可以加大计算机的字长(例如用双精度量)使其减少。5-6积分阶数的选择积分阶数的选择要考虑几方面的因素。
14、1、积分精度单元刚度矩阵公式(5-4-5)的被积函数中起主要影响的是两个因素：几何矩阵B和 Jacobi行列式detJ。暂时认为E、t在单元内为常量。当网格分得很细时，每个单元将接近常应变状态，B接近于常量(与常量之差为高阶 小量)。detJ则反映单元的畸变程度在网格加密的过程中单元畸变未必会减轻。为了使数值 积分得到的变形能逼近真实解的变形能，必须在单元内得到detJ积分的精确值。对于四结 点单元，detJ是& n的双线性函数，对于八结点单元detJ中出现& n的最高次数是督，网 故从收敛到真实解这一要求出发，四结点单元的积分阶数至少为1，八结点单元的积分阶数 至少为2。对于矩形单元detJ
15、为常数，单元内B的元素是小 的多项式。对于四结点单元BtEB 的元素是4, n的二次多项式，积分阶数取为2即可得到精确的数值积分、对于八结点单元 BteB的元素是4, n的四次多项式，积分阶数取为3即可得到精确的数值积分。在厚度 t为4, n线性函数的情况下，不必增加积分点个数仍可保证数值积分的精确性。在一般情况下，单元不会很小(单元内不接近常应变状态)，单元也不一定为矩形(detJ 在单元内不是常量)。为保证积分精度，对于四结点单元一般取积分阶数为2。对八结点单 元积分阶数则为23。对于畸变严重的单元积分阶数可以取为4,但这种单元应尽量少用。 随着积分阶数的提高，积分点个数增加很快，计算量也
16、随之迅速增加。2、选择较低的积分阶数，得到的有限元解有时会更接近精确解。这一现象可以解释 为：协调位移单元的刚度偏大，而采用较低的积分阶数，却可能使刚度下降，从而改善了解 的质量。在一维情况下，这种现象可以给于严格的论证。设单元e为一长度为2的区间(不是这 种情况时可以通过简单的坐标变换化为这种情况)。设位移场(试探函数)为三次多项式， 几何矩阵B的元素将是x的二次多项式，BTEB的元素将是x的四次函数。取积分阶数 为3即可得到数值积分的精确值。从另一个角度看，每个B的元素b.是x的二次函数，可 以展成1b (x) = Ci + CiL (x) + CiL (x)i01 12 2BTEB的第i
17、行第j列的被积表达式将是ECi + CiL (x) + CiL (x)Cj + CjL (x) + CjL (x)01 12 20112 2利用Legendre多项式的正交性，精确的积分值将是k = 2 . E . (CiCj + 1 CiCj + 1 CiCj)ij0 031152 2若取L2(x)的两个零点为积分点，即取积分阶数为2,得到的数值积分为k = E . (Ci - Ci .)(Cj - C j . A + (Ci + Ci )(Cj + C j .ij01 3013013013=2 . E. (C C 0 + 3 qq)显然在二维情况下上述证明过程不能加以套用。这是因为：一个X
18、、y的多项式一般不能表 示成两个多项式P（x）和P（y）的乘积。但经验表明：对于矩形单元（detJ为常量），通过降 低积分阶数来改善解的质量常常是有效的。3、零变形能模式一个精确的单元刚度矩阵k是一个半正定阵，单元变形能仅对刚体位移模式才为零。在积分点个数过少，刚度阵为近似值的情况下，零特征值的数目 可能会多于独立刚体位移的个数。以至会出现这样一种或几种位移模式：它们不是刚体位移 模式，但对应的变形能为零，这样的位移模式称为零变形能模式。下面先以四结点单元为例 加以说明。为了直观上清晰，设单元e为矩形，边与x、y轴平行（图5-26（a）。如果取2 x 2个积 分点，可得到刚度阵的精确积分值。现
19、考虑只取一个积分点的情况。积分点坐标为 n）=（0, 0）。若B 提=(0,0) 一8.8 x8 yi xy)则数值积分给出的变形能(0,0)其中A为单元面积。图 5-2 6单元的总体自由度= uv u v u vv 4*共有8个未知量，而（5-6-1）只有三个方程。故齐次方程组（5-6-1 ）有五个独立的非零解， 它们代表了五种使单元变形能（数值积分值）为零的位移模式。其中三种为刚体位移模式（两 个平移，一个旋转），其余两种则为两种零变形能模式，分别示于图5-2 6（b）和（c）。如果 对采用一个积分点得到的单元刚度矩阵进行特征值检查（只进行三角分解，不必求特征值），可以发现零特征值的个数是
20、五而不是三。这两种零变形能模式在相邻的单元之间是相容的， 或者说单元之间不会制约这种位移模式的出现（图5-27）(a)（b）（a）(b)图 5-27加在区域边界上的强制边界条件有时可以约束零变形能模式（图5-28（a）。有时则不能（图图 5-28综上所述，零变形能模式要出现下结构中，必须具备三个条件：（i）对单元刚度矩阵进行数值积分时积分点个数较少（或者说积分阶数较低）。具体条 件是积分点个数x每点应变分量数单元结点自由度总数一单元独立刚体位移模式个数（ii）对于可能出现的零变形能位移模式，单元之间是相容的。（iii）加在区域上的强制边界条件不能约束零变 形位移模式。如果具备了这些条件，零变形
21、能模式 就可能以相当大的数值（齐次方程的解乘任何倍数 仍然是方程的解）迭加在正常的有限元解上，以至 将后者完全淹没。对于八结点矩形单元，聊x3个积分点即可得 到单元刚度矩阵的精确积分。当取2x 2个积分点时图 5-294 x 3 = 1216 - 3 = 1313 -12 = 1 存在着一个零变形能模式（图5-29）在这种零变形能模式中，单元外凸边受拉，内凹边受 压，相邻单元之间不相容（即相互制约）。因而不会在结构中出现。对于这样的单元，若采 用2x2个积分点可以给出相当好的结果。由于对于任何一个具体的单元而言，独立的刚体位移模式总是有限的。所以，随着一个 单元自由度总数的增多，出现零变形能模
22、式的危险也越大。对于长方形的三维20结点等参 元，当采用2 x 2 x 2个积分点时，每个积分点有六个应变分量，刚体位移个数为6。8 x 6 = 4820 x 3 - 6 = 5454 - 48 = 6有六个零变形能模式。当用这样的单元分析一个悬臂梁（图5-30 ）而在单元分析中又采用2 x 2 x 2个积分点时，有限元解将受到其中一个零变形能模式的严重干扰。这种零变形能模 式只有沿X、z两个方向的位移和w，它们是j的线性函数，在）=常数的截面上的分布同 图5-29，很易识别。z,w图 5-304、积分点的其他几个问题意义（1）应力计算 第四章中给出的有限元解的误差估计是平均意义下的（能量模意
23、 义），有些点上的精度比这个估计差，有些点上则比这个估计好。对等参数单元而言，高斯 积分点上应力的精度比其他点好，当单元接近矩形时将获得收敛性。为了提高输出应力的精 度。可以先在高斯点上计算应力，再用插值的方法（内插和外推）得出其他应力输出点上的 应力。这项技术在70年代中期以后已为一些应用程序所采用。（2）不可压缩条件与减缩积分 在有些问题中要考虑“不可压缩“条件 + = 0（5-6-2）条件（5-6-2）必须在区域6的每一点得到满足，因而是一个很强的约束条件。实际问题是 无限维的，附加这一条件不会发生什么问题。有限元空间则是一个有限维的子空间，要求单 元内每点满足条件（5-6-2）将使单元
24、自由度减少，或者说使有限元空间的容量降低。以至 会了生这样的情况：在施加了位移强制边界条件后，结构将丧失全部自由度，所有结点不能 参生位移，出现所谓自锁现象。出现自锁现象的原因不在于实际问题本身不合理，而在于有 限元空间的容量不够，克服这一困难的方法之一就是使用所谓减缩积分：在对某些矩阵进行 数值积分时有意识地减少积分点个数。这样做将削弱刚度矩阵对位移场的约束能力，在一般 的情况下将导致出现零变形能位移模式。但在我们现在所讨论的问题中却可以使有限元空间 容量不足的矛盾得到解决。（3） 试探函数的缺项问题例如，八结点等参数单元的试探函数为4、n的二次多项式：u = a +a g+a+a g2 +
25、 a gn+aq2 + a g2+a &门2 12345678与4、n的双二次多项式（对于4、n每一个都是完全二次多项式）u = a +a g +a n +a g2 +a gn + a n2 +a g2 +a gn2 +a g2门21 23456789相比少一项（4。在三维情况下，这种“缺项“现象要更严重一些。有一种20结点三维 等参元，它的试探函数是4、n、Z的二次多项式，共20项，如果要实现对4、n、Z每一个 都是完全二次多项式则要有27项。在4、n、Z的每一个自变量的方向上取相同的积分阶数 （二维情况下积分点个数为n2,三维情况下积分点个数为n3）可以保证被积函数对4、n、Z 每一个都是
26、完全多项式时有较好的积分精度。但在试探函数有“缺项“时，考虑到插值本身 的误差，完全可以用较少的积分点而保持同等的精度。对于前述的20结点的三维等参元有 一种14样点的积分公式，积分的精度与2 X 2 X 2个积分点的结果相同，但是计算量几乎 减少了一半。还可以用更加灵活的方式选择积分点。例如有时我们不希望剪切变形能进入总变形能 （纯弯曲精确解的剪切变形能为零）只要取剪切应变为零的点为积分点即可达到上述目的。5-7其它形式的等参数单元1e本章着重讲座了以四边形为基础的等参 数单元，以它为代表讨论了等参数单元的理论 和一些技术问题。而实际应用的等参数单元并 非只此一种，本节将其他类型的等参数单元
27、加-101e1 &以介绍。对于一种等参数单兀而言，起决定作-1016用的是匕的母体单兀（又称为参考单元）和形 函数。1、一维等参数单元1e母体单元：长2的直线段结点个数为24。二结点情况下形函数为-11 031316.1N （g）= 2（1+gig） （i = 1、2）图 5-31图 5-32其他情况下的形函数，可用修改“的方法得到。这种单元 除用来描述直杆外，还可用来描述“曲的“受拉物体，例 如悬索和悬链以及它们的几何非线性问题。（图5-31）2、三角形二维等参数单元母体单元：直角边长度为1的等腰直角三角形。（图5-32）自然坐标：、入2有时还引入第三个坐标 入3 = 1X 1 X 2，但人
28、、入2、入3之间不独立。结点个数：36（或更多）。结点个数为3时，形函数N （人）=人（i = 13）实际单元即为常应变三角形。当结点个数为6时，形函数可用“修改“的方法得到，实际单元为六结点直边或曲边三 角形。3、奇异单元等参数单元有可能很方便地构造奇异单元，以便在一些特殊的结点处(例如，前面讲过 的角点处)描述应力集中问题，下面以二维情况为例。考虑一个六结点单元(图5-33)为了书写方 便，实际单元选择了比较规则的形状。形函数可 以用“修改“的方法得到：n (&,n) = 1(1&)(1门)&1 4n (&，门)=1(1+&)(1n)&2 4n (&,n) = 1(1+&)(1+n)&3
29、4n (&,n) =1(1&)(1+n)&4 4n (&,n) = 1(1-&2)(1n)5 2n (&,n) = 1(1&2)(1+n)6 2它们是n的线性函数，是4的二次函数。当实际单元的坐标如图所示时，坐标变换为尤=4n (&,n)+4n (&,n) + n (&,n) + n (&,n)2 356=2(1 + &)& + (1-&2)=(1 + & )2r = 2N (&,n) + 2N (&,n) + 2n (&,n)3 46=1+n(a)图 5-33逆变换试探函数u =以+以g +以& 2 +以n +以gn+以& 2n12_345 _6=以+以t x+以x + 以y +a/xy+以
30、xy123456在x=0的线上2u 1 . k出现 (x-2)的奇异性，即可以描述应力集中现象。 ox(a)(b)图 5-34(a)图 5-35当八结点单元的结点如图5-34(b)配置时，将在(x, y) =(0, 0)点处， 普 有O(x-2)的1奇异性，有O(J-2的奇异性。而当实际单元为图5-35(b)所示的三角形单元时，在(x, 世 湖 .y)=(0, 0)处将使 有O(r-2)的奇异性，其中r为点M到O的距离。dr5-8三维等参元的讨论1、母体单元e是一边长为2的正立方体。(图5-36)自然坐标：& n、Z。结点个数：827。其中18为角点。920为各棱中点，21 为体心，2227为
31、各面面心。结点个数至少为8。形函数：当结点个数为8时一 -1N (& ,n, S = (1+孑)(1+叩)(1+。)(i = 1 8)实际单元的棱为直棱，面可以是双曲面。增加结点个数时，要补充新的形函数，同时对原有的形函数进图 5-36行“修改“。当结点个数加到20时，每条棱都允许是曲棱，每个面都可以是曲面。2、三维等参元的数值积分(i)对于8结点的等参元，2x 2x 2给出精确积分，1 x 1 x 1为缩减积分。但8 x 3 -1 x 6 = 18 6故，8结点的等参元不能采用缩减积分。(ii)对于27结点的等参元而言，应变为八n、Z的二次函数。精确积分的点数为27, 即3 x 3 x 3，
32、2 x 2 x 2为缩减积分，但27 x 3 - 8 x 6 = 33 6故，27结点的等参元不能采用缩减积分。(iii)对于20结点的三维等参元而言， 精确积分的点数为27，即3x 3x 3， 缩减积分为2 x 2 x 2，但20 x 3 - 8 x 6 = 12 6故，20结点的等参元不能采用缩减积分。有一种14积分点公式(图5-37)，精 度与3x 3x3积分点相同，但计算量几乎 减少了一半。20结点的等参元的积分点坐标及权系数如下：AA :山E 土俾土,四 -V 333333 J图 5-37(A1A8位于立方体的形心与各角点之间)例如积分点入1和坐标.19 ,19 ;19V 33. 3
33、333r 0,0, u I 130(BB6位于面心与 各角点之间)121气：吗W8 =商320B - B6 :W9 W14 = 361弁W = 121 x 8 + 320 x 6 = 8(母体单元的体积)i 361361i=120 x 3 -14 x 6 = -24 6因此，这种积分不会出现零变形能模式。5-9结束语等参数单元是目前应用最广的一类单元，它的边可直可曲，精度可高可低。由于采用数 值积分，处理材料非线性问题不会遇到新的困难。鉴于这些优点，在一些通用分析程序中， 等参元成为处理二阶问题的主要单元。但等参数单元也有它的不足之处：精度和计算量之间在存着矛盾。在二维情况下，结点 个数取到20单元才有较满意的适应能力，而这时计算量往往相当可观。解决这一矛盾的办 法是探索新型单元，下一章将介绍的非协调元就是其中的一类。
