《有限元分析(第二部分)》由会员分享,可在线阅读,更多相关《有限元分析(第二部分)(14页珍藏版)》请在金锄头文库上搜索。
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)其中为单元面积。单元的总体自由度共有个未知量,而(5-6-1)只有三个方程。故齐次方程组(5-6-1)有五个独立的非零解,它们代表了五种使单元变形能(数值积分值)为零的位移模式。其中三种为刚体位移模式(两个平移,一个旋转),其余两种则为两种零变