有限元第五章

1、5-5 数值积分数值积分1. 问题的提出问题的提出(i) 单元刚度矩阵 ddJtBEBtdxdyBEBkTTedet1111 (5-4-5)(ii) 体积力的等效结点力(iii) 边界力的等效结点力(iv) 温升载荷的等效结点力 ddJtffNtdffNryxTyxTeVdet1111 (5-4-6) dyyfdxxftdsqpNrT1111,(5-4-7) dJdtBTdet01111T (5-4-8)式（5-4-5）（5-4-8）分别归结为计算以下两种形式的积分 111111),()(dxdxyxfdxxf、2. 数值积分的基本概念数值积分的基本概念 任何积分工作取决于三个要素：任何积分工

2、作取决于三个要素： (1) 给定的积分区间给定的积分区间; (2) 给定的被积函数给定的被积函数; (3) 具体的积分方法。具体的积分方法。下面以一维情况为例介绍数值积分的基本概念。下面以一维情况为例介绍数值积分的基本概念。(i) 梯形法梯形法abxixi+1f(xi)f(xi+1)hf(x)xf(x) 函数 f(x) 在区间 (a,b) 的积分可以表达为 iniibaxfWdxxfI1iWix ixf：权系数；：积分样点；：积分样点的函数值。梯形法的求积公式为 niiibahxfxfdxxfI1121nabhabWnii1(ii) 当被积函数为当被积函数为n-1次多项式次多项式Pn-1(x)

3、时，则由时，则由n个样点及其样点值个样点及其样点值（xi, Pn-1(xi),i=1,n）可以确定这个多项式，从而可以得到精确解。）可以确定这个多项式，从而可以得到精确解。（注意：用以确定多项式的样点不必刻意选取）对多项式形式的被积函数（注意：用以确定多项式的样点不必刻意选取）对多项式形式的被积函数进行积分可以采用高斯求积法。进行积分可以采用高斯求积法。 3. Legendre 多项式多项式L0=11-1xL1= x1-1xL2-1L3-1xx1xx2-1/211110)(11xxxL根33,33)31(23)(212xxxxL根515, 0,515)53(25)(32123xxxxxxL根3

4、512015,3512015)3512015)(3512015(835)(3241234、根xxxxxL1)(0 xLLegendre多项式的定义域为-1,1 零阶Legendre多项式三阶Legendre多项式四阶Legendre多项式二阶Legendre多项式一阶Legendre多项式一般n阶Legendre多项式的定义为nnnnnxdxdnx) 1(!21)(2n 阶Legendre多项式是二阶变系数齐次微分方程在区间-1,1上的有界解。01212 ynnyxyxLn(x) 在区间（-1, 1）上有n个相异实根（零点）若再补充定义 1)(0 x、)()()()(0 xxxxn则得一个定义

5、在-1, 1上的多项式系列)() 1(4) 1() 1(2) 1() 1()() 1(2) 1() 1(22222212222112122xPxxnnxnxxdxdxpxxnxxdxdnnnnnnn)() 1() 1(22xPxxdxdkknnkk01) 1(2xxdxdnkk多项式的性质，涉及下面的关系：对于任何kn ）当m=n时则有122)1 () !(21) 1() 1() !(21)(1122211222112mdxxmdxxmdxxLmmmmmm0)()(111dxxPxLmm(2) 若在 (1) 的证明中将 Ln (x) 换成任何次数不超过m-1次的多项式 P m(x) 则有这表明

6、：Lm(x) 与任何一个次数不超过 m-1的多项式正交。(3) 若q(x) 是（，）上平方可积的函数，则可将q(x)展开成)()(0 xLCxqmmm11)()(212dxxqxLmCmm对于n次多项式Pn(x)有 )()(0 xLCxPmnmmn4. 一维情况一维情况设需要计算积分11)(dxxfI)0()0(2)0(111fWfdxfIf(x)xx00 x2x1 (x1 )-(a )(b )图5-24f(x)()()()(112121xxxxxfxfxf dxxxxxxfxfxfI111121211)()()()()()(2)(2221121211122xfWxfWxfxxxxfxxx取x

7、1=0为积分点以常量f(0) 代替f(x) 进行积分，作为的近似值 也可以任取两个点x1、x2 为积分点用一个线性插值函数 当当f(x)是一次函数时显然可得到的精确值（取一个积分点已经能做到是一次函数时显然可得到的精确值（取一个积分点已经能做到这一点），当这一点），当f(x)是二次、三次函数时又将如何？是二次、三次函数时又将如何？ 一般若取任取一般若取任取n个积分点个积分点x1、x2、 xn，作，作n-1次插值多项式，积分的次插值多项式，积分的近似值可表示为近似值可表示为)(1iniixfWI(5-5-1) 当积分点取为当积分点取为n 阶阶Legendre多项式的零点多项式的零点时，如果时，如

8、果被积函数被积函数f(x) 为次为次数不超过数不超过2n1次的多项式次的多项式，（，（5-5-1）将给出积分的）将给出积分的精确值精确值。这就是高斯。这就是高斯求积分法，上述积分点又称为求积分法，上述积分点又称为高斯点高斯点。高斯点的个数又称为。高斯点的个数又称为积分阶数积分阶数，有，有限元分析中一般限元分析中一般n=24。可以证明如下。可以证明如下 一般来说，对于一个2n-1阶的多项式，需用2n个样点及其样点值才能精确重构该多项式，或者说，需用2n个积分样点才能给出精确积分。 xLxRxQxLxPnnnnn1112 xRxLxQxPnnnn1112 1111111111112dxxRdxxR

9、xdxLxQdxxPnnnnn ininxRxP112ni1 xLn若取 的n个零点为积分样点，则 结论：用结论：用n个个Legendre 多项式的零点作为多项式的零点作为积分样点，式积分样点，式(5-5-1)可以给出精确的积分值，可以给出精确的积分值，这种方式的积分称为高斯积分。可见，高斯这种方式的积分称为高斯积分。可见，高斯积分方法即减少了积分样点数也优化了积分积分方法即减少了积分样点数也优化了积分样点。样点。 积分阶数n高斯点坐标权系数2n1x1=0W1=21X 1，2= 0.5773502692W1,2=13x 1，3 = 0.7745966692 x2=0W1,3=5/9W1=8/9

10、5X 1,4 = 0.8611363116X 2,3 = 0.3399810436W1,4= 0.3478548451W2,3= 0.65214515497常见的高斯点坐标和权系数5. 二维情况二维情况 111111),(),(ninjjijiyxfWWdxdyyxfI 积分阶数积分阶数n是对每一个自变量而言的积分阶数，而积分点总数在二维情是对每一个自变量而言的积分阶数，而积分点总数在二维情况下为况下为n2，在三维情况下为，在三维情况下为n3。11223344一种经常采用的（并非唯一可能的）选择方式是：沿x、y方向取同样个数的积分点 单元刚度公式（5-4-5）的数值积分形式 ddJtBEBtd

11、xdyBEBkTTedet1111 ),(det11jininjTjiJtBEBWWk(5-5-3)(5-4-5)(3)对于每个积分点都必须计算detJ 的值。 (1) 其中(i，j )为积分点的自然坐标。注意：(2) 矩阵在积分点上的函数值是确定的，因为它与形函数相关，而形函数是一个确定的函数。(4)八结点单元而言k共有162=256个元素，利用对称性仍需对其中的136个元素进行数值积分。. 有限元解的误差有限元解的误差有限元解的误差是四个方面因素综合影响的结果。有限元解的误差是四个方面因素综合影响的结果。 可以加大计算机的字长（例如用双精度量）使其减少。(i) 插值误差插值误差 这是在单元

12、内用多项式代替真实解（在整个求解域内则表现为用有限自由度代替了无限自由度）所引起的。这个因素在协调位移单元中将导致“刚度偏大”。（ii）边界形状以及边界条件的误差边界形状以及边界条件的误差 即使采用了曲边单元，单元边界仍有它本身的特点（例如是某个自然坐标的二次函数），不可能作到与实际曲线边界完全吻合。边界形状的误差使得实际边界条件不能得到精确满足。这种误差一般只在边界附近影响较大（奇点除外）。 (iii) 数值积分误差数值积分误差 这种“误差”如果利用得好，可以与（i）引起的误差“抵消”，但处理不当也将影响解的精度。（iv）截断误差截断误差5-6 积分阶数的选择积分阶数的选择1. 积分精度积分

13、精度(1) B接近于常量（与常量之差为高阶小量）；(2) detJ则反映单元的畸变程度在网格加密的过程中单元畸变未必会减轻。则反映单元的畸变程度在网格加密的过程中单元畸变未必会减轻。在单元 内要尽量得到detJ积分的精确值；(3) 对于四结点单元，detJ 是， 的双线性函数，四结点单元的积分阶数至少为；(4) 对于八结点单元detJ中出现，的最高次数是3，3。，八结点单元的积分阶 数至少为。当网格分得很细时当网格分得很细时:在一般情况下：在一般情况下： 为保证积分精度，对于四结点单元一般取积分阶数为。对八结点单元积分阶数则为。对于畸变严重的单元积分阶数可以取为，但这种单元应尽量对于畸变严重的

14、单元积分阶数可以取为，但这种单元应尽量少用。少用。 .人为的选择较低的积分阶数人为的选择较低的积分阶数 协调位移单元的刚度偏大，而采用较低的积分阶数，却可能使刚度下降，从而改善了解的质量。例证如下： 在一维情况下，设单元e为一长度为的区间。设位移场（试探函数）为三次多项式，几何矩阵 B 的元素将是 x 的二次多项式,BTEB的元素将是 x 的四次函数。取积分阶数为即可得到数值积分的精确值。从另一个角度看，每个B的元素bi是 x 的二次函数，可以展成)()()(22110 xLCxLCCxbiiii)()()()(2211022110 xLCxLCCxLCxLCCEjjjiii)5131(222

15、1100jijijiijCCCCCCEk)31(2)33)(33()33)(33(110010101010jijijjiijjiiijCCCCECCCCCCCCEkijijkkBTEB的第i行第j列的被积表达式将是利用Legendre多项式的正交性，精确的积分值将是(取积分阶数为3)若取L2(x)的两个零点为积分点，即取积分阶数为，得到的数值积分为. 零变形能零变形能(位移位移)模式模式一个精确的单元刚度矩阵k是一个半正定阵，单元变形能 ukuVTe21 仅对刚体位移模式才为零。在积分点个数过少，刚度阵为近似值的情况下，零特征值的数目可能会多于独立刚体位移的个数。以至会出现这样一种或几种位移模

16、式：它们不是刚体位移模式，但对应的变形能为零，这样的位移模式称为零变形能模式零变形能模式。 （a）（b）（c）图5-2e(0,0)(0,0)(a)(b)图5-27（a）（b）图5-28零变形能模式要出现下结构中，必须具备三个条件：零变形能模式要出现下结构中，必须具备三个条件： (i) 对单元刚度矩阵进行数值积分时积分点个数较少（或者说积分阶数较低）。积分点个数 每点应变分量数单元结点自由度总数单元独立刚体位移模式个数 (ii) 对于可能出现的零变形能位移模式，单元之间是相容的。 (iii) 加在区域上的强制边界条件不能约束零变形位移模式。 e图5-29 对于八结点矩形单元，取 个积分点即可得到

17、单元刚度矩阵的精确积分。当取 个积分点时332231634 存在着一个零变形能模式（图5-29）在这种零变形能模式中，单元外凸边受拉，内凹边受压，相邻单元之间不相容（即相互制约）。因而不会在结构中出现。 x,u z,wy,v图5-30当用这样的单元分析一个悬臂梁（图5-30）而在单元分析中又采用个积分点时，有限元解将受到其中一个零变形能模式的严重干扰。 2224. 积分点的其他几个问题积分点的其他几个问题（）应力计算（）应力计算 第四章中给出的有限元解的误差估计是平均意义下的（能量模意义），有些点上的精度比这个估计差，有些点上则比这个估计好。对等参数单元而言，高斯积分点上应力的精度比其他点好，

18、当单元接近矩形时将获得收敛性。 为了提高输出应力的精度。可以先在高斯点上计算应力，再用插值的方法为了提高输出应力的精度。可以先在高斯点上计算应力，再用插值的方法（内插和外推）得出其他应力输出点上的应力。（内插和外推）得出其他应力输出点上的应力。 （2） 不可压缩条件与减缩积分不可压缩条件与减缩积分 “不可压缩“条件 0yx（5-6-2） 条件（5-6-2）必须在区域的每一点得到满足，因而是一个很强的约束条件。实际问题是无限维的，附加这一条件不会发生什么问题。有限元空间则是一个有限维的子空间，要求单元内每点满足条件（5-6-2）将使单元自由度减少，或者说使有限元空间的容量降低。以至会了生这样的情

19、况：在施加了位移强制边界条件后，结构将丧失全部自由度，所有结点不能参生位移，出现所谓自锁现象自锁现象。 在对类似橡胶和一些高分子材料进行有限元分析时可能出现这类现象。在对类似橡胶和一些高分子材料进行有限元分析时可能出现这类现象。 出现自锁现象的原因不在于实际问题本身不合理，而在于有限元空间的容量不够，克服这一困难的方法之一就是使用所谓减缩积分减缩积分：在对某些矩阵进行数值积分时有意识地减少积分点个数。这样做将削弱刚度矩阵对位移场的约束能力，在一般的情况下将导致出现零变形能位移模式。但在我们现在所讨论的问题中却在一般的情况下将导致出现零变形能位移模式。但在我们现在所讨论的问题中却可以使有限元空间

20、容量不足的矛盾得到解决。可以使有限元空间容量不足的矛盾得到解决。（3） 试探函数的缺项问题试探函数的缺项问题 例如例如，八结点等参数单元的试探函数为 、 的二次多项式：282726524321u229282726524321u、 的双二次多项式（对于、每一个都是完全二次多项式）少一项（）。 在试探函数有在试探函数有“缺项缺项“时，考虑到插值本身的误差，完全可以用较时，考虑到插值本身的误差，完全可以用较少的积分点而保持同等的精度。少的积分点而保持同等的精度。 5-7 其它形式的等参数单元其它形式的等参数单元 对于等参数单元而言，起决定作用的是它的母体单元（又称为参考单元）和对于等参数单元而言，起

21、决定作用的是它的母体单元（又称为参考单元）和形函数。形函数。1. 一维等参数单元一维等参数单元-3131图5-31 母体单元：长的直线段结点个数为24。二结点情况下形函数为)21()1 (21)(、iNii 其他情况下的形函数，可用“修改“的方法得到。这种单元除用来描述直杆外，还可用来描述“曲的“受拉物体，例如悬索和悬链以及它们的几何非线性问题。（图5-31）. 三角形二维等参数单元三角形二维等参数单元图5-32母体单元：母体单元：直角边长度为的等腰直角三角形。自然坐标：自然坐标： 、有时还引入第三个坐标有时还引入第三个坐标 1， 但但、之间不独立。之间不独立。结点个数：结点个数： 6（或更多

22、）。结点个数为3时，形函数 ）（31)(iNii 实际单元即为常应变三角形。 当结点个数为6时，形函数可用“修改“的方法得到，实际单元为六结点直边或曲边三角形。. 奇异单元奇异单元 等参数单元有可能很方便地构造奇异单元，以便在一些特殊的结点处（例如，前面讲过的角点处）描述应力集中问题。4(a)(b)(4, 2)(4, 0)(1, 2)(1, 0)(0,0)(0, 2)156432图5-33xy 六结点单元（图5-33）。形函数 :)1)(1 (21),()1)(1 (21),()1)(1 (41),()1)(1 (41),()1)(1 (41),()1)(1 (41),(26254321NNN

23、NNN形函数是的线性函数，是的二次函数。当实际单元的坐标如图所示时，其坐标变换：1),(2),(2),(2)1 ()1 ()1 (2),(),(),(4),(4643226532NNNYNNNNx11yxxyyxyxxu65432126542321逆变换试探函数(b)(4, 2)(4, 0)(1, 2)(1, 0)(0,0)(0, 2)156432图5-33xy)(21xoxu在 x=0 的线上 出现 的奇异性，即可以描述应力集中现象。图5-344(a)(b)(a, b)(a, 0)(a/2, b)(a/4, 0)(0,0)(0, b)15743286(a, b/2)(0, b/4)yx1,4

24、,862y(a+c)/2,d/2)x753(a,o)(c,d)(a/4,o)(c/4,d/4)(c) 当八结点单元的结点如图5-34(b) 配置时，将在（x, y）=(0, 0) 点处，有的奇异性， 有 的奇异性。xu)(21xoyu)(21yoru)(21ro 而当实际单元为图5-35(b) 所示的三角形单元时，在（x, y）=(0, 0) 处将使有 的奇异性, 其中r为点M到 O的距离。 5-8 三维等参元的讨论三维等参元的讨论图5-36 1. 母体单元母体单元是一边长为的正立方体是一边长为的正立方体 自然坐标：、。 （图5-36） 结点个数：827。其中18为角点。920为各棱中点， 21 为体心，2227为各面面心。结点个数至少为8。形函数：当结点个数为8时）（81)1)(1)(1 (81),(iNiiii 实际单元的棱为直棱，面可以是双曲面。增加结点个数时，要补充新的形函数，同时对原有的形函数进行“修改“。当结点个数加到20时，每条棱都允许是曲棱，每个面都可以是曲面。2. 三维等参元的数值积分三维等参元的数值积分6186138故，故，27 结点的等参元不能采用缩减积分。结点的等参元不能采用缩减