有限元基本理论及工程应用：第五章(第二部分)

1、5-5数值积分、 问题的提出在上一节中对等参元进行单元分析时要进行下列积分：(i) 单元刚度矩阵(5-4-5)(ii)体积力的等效结点力(5-4-6)(iii)边界力的等效结点力(5-4-7)(iv)温升载荷的等效结点力(5-4-8)式（5-4-5）（5-4-8）分别归结为计算以下两种形式的积分对于上述积分仅在单元的形状十分规则的情况下才能得到解析的结果（精确值），一般情况只能用数值积分方法（主要是高斯求积法）求近似值。虽然数值积分是“被迫“采用的，但后来发现：有选择地控制积分点的个数和位置，可以方便地实现我们的某些特殊意图。这样一来，数值积分就成为有限元分析的一个重要组成部分，以至本来可以精

2、确积分的三角形单元也常常采用数值积分。、 数值积分的基本概念任何积分工作取决于三个要素：给定的积分区间，给定的被积函数，具体的积分方法。下面以一维情况为例介绍数值积分的基本概念(i) 梯形法 函数在区间(a,b)的积分可以表达为abxixi+1f(xi)f(xi+1)hf(x)xf(x)：权系数；：积分样点；：积分样点的函数值。梯形法的求积公式为图5-22其中，而(ii) 当被积函数为n-1次多项式Pn-1(x)时，则由n个样点及其样点值（xi, Pn-1(xi),i=1,n）可以精确重构这个多项式，从而可以得到精确解。（注意：用以确定多项式的样点不必刻意选取）对多项式形式的被积函数进行积分可

3、以采用高斯求积法，高斯求积法中要利用Legendre多项式的性质，下面先介绍这种多项式及相关的性质。、 Legendre 多项式Legendre多项式的定义域为-1,1一阶Legendre多项式二阶Legendre多项式三阶Legendre多项式四阶Legendre多项式L0=11-1xL1=x1-1xL2-1L3-1xx1xx2-1/21111图5-23一般n阶Legendre多项式的定义为读者可以验证：n阶Legendre多项式是二阶变系数齐次微分方程在区间-1,1上的有界解。Ln(x) 在区间（-1, 1）上有n个相异实根（零点）若再补充定义则得一个定义在-1, 1上的多项式系列图5-2

4、3 给出了Legendre 多项式函数列中的前四个元素。为了进一步讨论这个多项式的性质，先注意下面一个事实：对于任何kn）当m=n时则有（）若在（）的证明中将Ln(x)换成任何次数不超过m-1次的多项式 Pm(x) 则有这表明：Lm(x) 与任何一个次数不超过 m-1的多项式正交。（）若q(x) 是（，）上平方可积的函数，则可将q(x)展开成其中系数特别，对于n次多项式Pn(x)有、 一维情况设需要计算积分我们可以取x1=0为积分点（图5-2(a)），以常量f(0) 代替f(x) 进行积分，作为的近似值f(x)f(x)xx00x2x1 (x1 )(a )(b )图5-24当f(x)是一次函数时

5、，可得到的精确值。也可以任取两个点x1、x2为积分点（图5-2（b），用一个线性插值函数代替f(x) 进行积分，作为的近似值当f(x)是一次函数时显然可得到的精确值（取一个积分点已经能做到这一点），当f(x)是二次、三次函数时又将如何？一般若取任取n个积分点x1、x2、 xn，作n-1次插值多项式，积分的近似值可表示为(5-5-1)其中f(xi)为积分点上的函数值，Wi为权系数。当积分点取为n阶Legendre多项式的零点时，如果被积函数f(x) 为次数不超过2n1次的多项式，（5-5-1）将给出积分的精确值。这就是高斯求积分法，上述积分点又称为高斯点。高斯点的个数又称为积分阶数，有限元分析中

6、一般n=24。可以证明如下一般来说，对于一个2n-1阶的多项式，需用2n个样点及其样点值才能精确重构该多项式，或者说，需用2n个积分样点才能给出精确积分。若用除以多项式则即所以有为n-1阶多项式，因此仅需n个样点及其样点值即可精确重构该多项式，进而给出精确的积分值，若取的n个零点为积分样点，则结论：用n个Legendre 多项式的零点作为积分样点，式(5-5-1)可以给出精确的积分值，这种方式的积分称为高斯积分。可见，高斯积分方法即减少了积分样点数也优化了积分样点。现将常见的高斯点坐标和权系数列成下表：积分阶数n高斯点坐标权系数2n1x1=0W1=21x1，2=0.5773502692W1,2

7、=13x1，3=0.7745966692x2=0W1,3=5/9W1=8/95x1，4=0.3611363116x2，3=0.3399810436W1,4=0.3478548451W2,3=0.65214515497、 二维情况二维情况（三维情况与此类似）积分点的选择有更大的灵活性。一种经常采用的（并非唯一可能的）选择方式是：沿x、y方向取同样个数的积分点，积分的近似表达式为(5-5-2)积分阶数n是对每一个自变量而言的积分阶数，而积分点总数在二维情况下为n2，在三维情况下为n3。图5-2单元刚度公式（5-4-5）的数值积分形式将是(5-5-3)其中(i，j )为积分点的自然坐标。对于每个积分

8、点都必须将（5-4-1）（5-4-5）执行一遍！注意：矩阵在积分点上的函数值是确定的，因为它与形函数相关，而形函数是一个确定的函数。此外，上述公式均是矩阵运算对八结点单元而言k共有162=256个元素，利用对称性仍需对其中的136个元素进行数值积分。在单元结点个数增加，积分阶数也随之提高的情况下，计算量是相当可观的。对三维情况尤其如此。、有限元解的误差有限元解的误差是四个方面因素综合影响的结果。(i)插值误差这是在单元内用多项式代替真实解（在整个求解域内则表现为用有限自由度代替了无限自由度）所引起的。这个因素在协调位移单元中将导致“刚度偏大”。（ii）边界形状以及边界条件的误差即使采用了曲边单

9、元，单元边界仍有它本身的特点（例如是某个自然坐标的二次函数），不可能作到与实际曲线边界完全吻合。边界形状的误差使得实际边界条件不能得到精确满足。这种误差一般只在边界附近影响较大（奇点除外）。 (iii) 数值积分误差这种“误差”如果利用得好，可以与（i）引起的误差“抵消”，但处理不当也将影响解的精度。（iv）截断误差 可以加大计算机的字长（例如用双精度量）使其减少。5-6积分阶数的选择积分阶数的选择要考虑几方面的因素。、 积分精度单元刚度矩阵公式（5-4-5）的被积函数中起主要影响的是两个因素：几何矩阵B和Jacobi行列式detJ。暂时认为E、t在单元内为常量。当网格分得很细时，每个单元将接

10、近常应变状态，B接近于常量（与常量之差为高阶小量）。detJ则反映单元的畸变程度在网格加密的过程中单元畸变未必会减轻。为了使数值积分得到的变形能逼近真实解的变形能，必须在单元内得到detJ积分的精确值。对于四结点单元，detJ是，的双线性函数，对于八结点单元detJ中出现，的最高次数是3，3。故从收敛到真实解这一要求出发，四结点单元的积分阶数至少为，八结点单元的积分阶数至少为。对于矩形单元detJ为常数，单元内B的元素是，的多项式。对于四结点单元BTEB的元素是，的二次多项式，积分阶数取为即可得到精确的数值积分、对于八结点单元BTEB的元素是，的四次多项式，积分阶数取为即可得到精确的数值积分。

11、在厚度t为，线性函数的情况下，不必增加积分点个数仍可保证数值积分的精确性。在一般情况下，单元不会很小（单元内不接近常应变状态），单元也不一定为矩形（detJ在单元内不是常量）。为保证积分精度，对于四结点单元一般取积分阶数为。对八结点单元积分阶数则为。对于畸变严重的单元积分阶数可以取为，但这种单元应尽量少用。随着积分阶数的提高，积分点个数增加很快，计算量也随之迅速增加。、选择较低的积分阶数，得到的有限元解有时会更接近精确解。这一现象可以解释为：协调位移单元的刚度偏大，而采用较低的积分阶数，却可能使刚度下降，从而改善了解的质量。在一维情况下，这种现象可以给于严格的论证。设单元e为一长度为的区间（不

12、是这种情况时可以通过简单的坐标变换化为这种情况）。设位移场（试探函数）为三次多项式，几何矩阵B的元素将是x的二次多项式,BTEB的元素将是x的四次函数。取积分阶数为即可得到数值积分的精确值。从另一个角度看，每个B的元素bi是x的二次函数，可以展成BTEB的第i行第j列的被积表达式将是利用Legendre多项式的正交性，精确的积分值将是若取L2(x)的两个零点为积分点，即取积分阶数为，得到的数值积分为显然在二维情况下上述证明过程不能加以套用。这是因为：一个x、y的多项式一般不能表示成两个多项式P(x) 和P(y) 的乘积。但经验表明：对于矩形单元（detJ为常量），通过降低积分阶数来改善解的质量

13、常常是有效的。、零变形能模式一个精确的单元刚度矩阵k是一个半正定阵，单元变形能仅对刚体位移模式才为零。在积分点个数过少，刚度阵为近似值的情况下，零特征值的数目可能会多于独立刚体位移的个数。以至会出现这样一种或几种位移模式：它们不是刚体位移模式，但对应的变形能为零，这样的位移模式称为零变形能模式。下面先以四结点单元为例加以说明。为了直观上清晰，设单元e为矩形，边与x、y轴平行（图5-2(a)）。如果取22个积分点，可得到刚度阵的精确积分值。现考虑只取一个积分点的情况。积分点坐标为(，)=(0, 0)。若（5-6-1）则数值积分给出的变形能（a）（b）（c）图5-2e(0,0)(0,0)其中为单元

14、面积。单元的总体自由度共有个未知量，而（5-6-1）只有三个方程。故齐次方程组（5-6-1）有五个独立的非零解，它们代表了五种使单元变形能（数值积分值）为零的位移模式。其中三种为刚体位移模式（两个平移，一个旋转），其余两种则为两种零变形能模式，分别示于图5-2（b）和(c)。如果对采用一个积分点得到的单元刚度矩阵进行特征值检查（只进行三角分解，不必求特征值），可以发现零特征值的个数是五而不是三。这两种零变形能模式在相邻的单元之间是相容的，或者说单元之间不会制约这种位移模式的出现（图5-27）(a)(b)图5-27（a）（b）图5-28加在区域边界上的强制边界条件有时可以约束零变形能模式（图5-

15、28(a)）。有时则不能（图5-28(b)）。综上所述，零变形能模式要出现下结构中，必须具备三个条件：(i)对单元刚度矩阵进行数值积分时积分点个数较少（或者说积分阶数较低）。具体条件是积分点个数每点应变分量数单元结点自由度总数单元独立刚体位移模式个数 (ii)对于可能出现的零变形能位移模式，单元之间是相容的。e图5-29(iii)加在区域上的强制边界条件不能约束零变形位移模式。如果具备了这些条件，零变形能模式就可能以相当大的数值（齐次方程的解乘任何倍数仍然是方程的解）迭加在正常的有限元解上，以至将后者完全淹没。对于八结点矩形单元，取个积分点即可得到单元刚度矩阵的精确积分。当取个积分点时存在着一

16、个零变形能模式（图5-29）在这种零变形能模式中，单元外凸边受拉，内凹边受压，相邻单元之间不相容（即相互制约）。因而不会在结构中出现。对于这样的单元，若采用个积分点可以给出相当好的结果。由于对于任何一个具体的单元而言，独立的刚体位移模式总是有限的。所以，随着一个单元自由度总数的增多，出现零变形能模式的危险也越大。对于长方形的三维20结点等参元，当采用个积分点时，每个积分点有六个应变分量，刚体位移个数为。有六个零变形能模式。当用这样的单元分析一个悬臂梁（图5-30）而在单元分析中又采用x,u z,wy,v图5-30个积分点时，有限元解将受到其中一个零变形能模式的严重干扰。这种零变形能模式只有沿x

17、、z两个方向的位移u和w，它们是y的线性函数，在y常数的截面上的分布同图5-29，很易识别。、 积分点的其他几个问题意义（）应力计算第四章中给出的有限元解的误差估计是平均意义下的（能量模意义），有些点上的精度比这个估计差，有些点上则比这个估计好。对等参数单元而言，高斯积分点上应力的精度比其他点好，当单元接近矩形时将获得收敛性。为了提高输出应力的精度。可以先在高斯点上计算应力，再用插值的方法（内插和外推）得出其他应力输出点上的应力。这项技术在年代中期以后已为一些应用程序所采用。（5-6-2）（） 不可压缩条件与减缩积分在有些问题中要考虑“不可压缩“条件条件（5-6-2）必须在区域的每一点得到满足

18、，因而是一个很强的约束条件。实际问题是无限维的，附加这一条件不会发生什么问题。有限元空间则是一个有限维的子空间，要求单元内每点满足条件（5-6-2）将使单元自由度减少，或者说使有限元空间的容量降低。以至会了生这样的情况：在施加了位移强制边界条件后，结构将丧失全部自由度，所有结点不能参生位移，出现所谓自锁现象。出现自锁现象的原因不在于实际问题本身不合理，而在于有限元空间的容量不够，克服这一困难的方法之一就是使用所谓减缩积分：在对某些矩阵进行数值积分时有意识地减少积分点个数。这样做将削弱刚度矩阵对位移场的约束能力，在一般的情况下将导致出现零变形能位移模式。但在我们现在所讨论的问题中却可以使有限元空

19、间容量不足的矛盾得到解决。（） 试探函数的缺项问题例如，八结点等参数单元的试探函数为、的二次多项式：与、的双二次多项式（对于、每一个都是完全二次多项式）相比少一项（）。在三维情况下，这种“缺项“现象要更严重一些。有一种20结点三维等参元，它的试探函数是、的二次多项式，共20项，如果要实现对、每一个都是完全二次多项式则要有27项。在、的每一个自变量的方向上取相同的积分阶数（二维情况下积分点个数为n2，三维情况下积分点个数为n3）可以保证被积函数对、每一个都是完全多项式时有较好的积分精度。但在试探函数有“缺项“时，考虑到插值本身的误差，完全可以用较少的积分点而保持同等的精度。对于前述的20结点的三

20、维等参元有一种14样点的积分公式，积分的精度与个积分点的结果相同，但是计算量几乎减少了一半。还可以用更加灵活的方式选择积分点。例如有时我们不希望剪切变形能进入总变形能（纯弯曲精确解的剪切变形能为零）只要取剪切应变为零的点为积分点即可达到上述目的。图5-315-7其它形式的等参数单元本章着重讲座了以四边形为基础的等参数单元，以它为代表讨论了等参数单元的理论和一些技术问题。而实际应用的等参数单元并非只此一种，本节将其他类型的等参数单元加以介绍。对于一种等参数单元而言，起决定作用的是它的母体单元（又称为参考单元）和形函数。、 一维等参数单元母体单元：长的直线段结点个数为24。二结点情况下形函数为图5

21、-32其他情况下的形函数，可用“修改“的方法得到。这种单元除用来描述直杆外，还可用来描述“曲的“受拉物体，例如悬索和悬链以及它们的几何非线性问题。（图5-31）、三角形二维等参数单元母体单元：直角边长度为的等腰直角三角形。（图5-32）自然坐标：、有时还引入第三个坐标1，但、之间不独立。结点个数：6（或更多）。结点个数为3时，形函数实际单元即为常应变三角形。当结点个数为6时，形函数可用“修改“的方法得到，实际单元为六结点直边或曲边三角形。、奇异单元4(a)(b)(4, 2)(4, 0)(1, 2)(1, 0)(0,0)(0, 2)156432图5-33xy等参数单元有可能很方便地构造奇异单元，

22、以便在一些特殊的结点处（例如，前面讲过的角点处）描述应力集中问题，下面以二维情况为例。考虑一个六结点单元（图5-33）为了书写方便，实际单元选择了比较规则的形状。形函数可以用“修改“的方法得到：它们是的线性函数，是的二次函数。当实际单元的坐标如图所示时，坐标变换为逆变换试探函数在x=0的线上 出现 的奇异性，即可以描述应力集中现象。(b)(a, b)(a, 0)(a/2, b)(a/4, 0)(0,0)(0, b)157432图5-344(a)86(a, b/2)(0, b/4)yx4(a)1,4,862y(a+c)/2,d/2)x753(a,o)(c,d)(a/4,o)(c/4,d/4)图5-35(b)当八结点单元的结点如图5-34(b) 配置时，将在（x, y）=(0, 0) 点处，有的奇异性， 有 的奇异性。而当实际单元为图5-35(b) 所示的三角形单元时，在（x, y）=(0, 0) 处