第6章曲边等参单元和数值积分

1、计算力学第六章单元和数值本章目录单元基本概念6.16.2 参数曲线坐标6.3 雅可比矩阵6.4 数值单元矩阵计算6.5 数值阶次的选择*6.6 应力结果的处理6.1单元基本概念 1引言前面介绍了C0插值获得一些常用单元，这类单元在实际应用中有一个无法回避的缺点：在不增加计算工作量的 情况下，采用增加节点数方法提高单元计算精度，必然 减少单元数量，描述相对复杂的几何形状的能力也随之 下降。这是一个实际问题，也即仅采用简单的三角形和四边形单元是不够的，因此，本章讨论由简单单元转换 为其他复杂形状单元的问题。1单元基本概念如下图所示，二维单元可被为扭曲的单元形状。(-1,1)(1,1)yx(-1,-

2、1)(1,-1)畸变单元(实际单元)基本单元变换函数x 总体坐标局部坐标 x y = f h z z 可以看出，局部坐标系（x ,h )或者（L1,L2,L3)，在笛卡尔（x,y)空间中被为新的曲线形状。6.2参数曲线坐标1坐标变换形函数N1很明显，基本单元和畸变单元从形状上就可以看出不是一样的单元，要在两者之间建立，首先应进行坐标变换，建立坐标变换最方便的办法就是使用我们最熟悉的，用以描述 位移函数（或其他场函数）变化的标准C0形状函数。局部坐标与整体坐标的变换局部坐标系(基准坐标系)整体坐标系(物理坐标系)1x = N x+ N x+ N x+ N x11223344(6-1)y = N

3、y+ N y+ N y+ N y11223344mmmx =y =z =N xN yN ziiiiiii =1i =1i =1式中N是基于局部坐标给出的标准形状函数,于是可以马上得到需要的变换关系.可以看出,坐标(x1,y1).的点正好就是单元边界上的对应点(从标准形状函数的定义可知,在对应点上函数值为1,而在其他点为零)，可在这些点上预先建立节点。1将坐标变换式用于任意四边形单元，可得：12在四个结点处给出结点的整体坐标。在四条边上的整体坐标是线性变化的。只要给出任意四边形单元四个结点的整体坐标，用（6-1）式就可以建立局部坐标系中的正方形单元和整体坐标系中的任意四边形单元之间的坐标变换关系

4、。关键在于，首先要解决进行任意形状单元和规则单元之间的几何变换。变换的关键是什么？1mmmx =y =z =N xN yN ziiiiiii =1i =1i =1局部坐标下的每一点都对应整个笛卡尔坐标的一个点，而且一般说来变换是唯一的。但若出现对应的非唯一性，则将引起单元的剧烈畸变一个变换对应有一个雅可比（Jacobian）行列式，可由下式定义：D(x, y, z)D(x ,h,z )J = 即双线性单元不能有内角大于或等于或接近180度情况。雅可比行列式，具体定义可见6.3节。 0J2sin(dx , dh) 01 n 亚参数单元：m 0J194 关于等参数单元说明总结3） 等参单元的优点是

5、当单元边界呈二次以上的曲线时， 容易用很少的单元去逼近曲线边界。将不规则单元转换为规则母单元后，容易构造位移函数和形函数。4） 上述等参单元的理论公式可适应三次以上的曲线型等 参元，只是阶次提高，单元自由度相应增加，计算更复杂，更，实际中，很少超过3次曲线型。5） 上述推导要求：保持坐标变换中几何模式阶次与描述 单元位移函数中形函数的阶次相同。如取坐标变换的几何模式阶次较单元的位移函数的阶次高，则称此单元为超参单元，反之，为亚参单元。这两类单元的收敛性也可得到满足。6） 当然，也可取描述单元几何形状的几何模式不是形函 数的，如p-element6.3雅可比矩阵计算单元矩阵1.雅可比矩阵有限元分

6、析中，建立求解方程，需要进行各个单元体积内和面积内的，其求解的形式：V GdV =G(x, y, z)dxdydz体：面：SgdS = g(x, y, z)dSG、g的形式取决于不同的单元，通常包含Ni对x, y, z的导数NiNiNi= ?= ?= ?xyz1.雅可比矩阵求导的计算(x, y, z)(x ,h,z )Ni = Nix + Niy + Ni zNih= Ni x + Niy + Nizhxxxz xx hhxyyzNi = Nix + Niyz+ Nizzzz zxy写成矩阵形式：Ni xyxz NiNi x xx x x Ni = xz NiyhyzN= Ji h h y y

7、h xNiz NiNi zz z zz J - -雅可比矩阵1.雅可比矩阵Ni Ni x xNiN-1= Ji y h NiNi z zymmmi=1引入: x =z =N xN yN ziiiiiii=1i=1N N N N N N mmmi i i xyz 1 2x m xixixi xx xyzi=1i=1i=1111 xN N N N N N yzmmmJ = i i i z = 1 m 222xy 2h hihihi hh i=1i=1i=1 xN N N N N N yzmmm mm i i i z 1mxy 2z m zzizizi z i=1i=1i=12.单元矩阵的计算(局部

8、坐标下的变换)1 体为对过程中的变量和区域进行变换，需要相应的包括雅可比行列式在内的标准处理方法。例如，体积元的计算。设 x, y, z坐标系的轴为i, j, kdx = x dx i + y dx j + z dx kxxxdh = x dhi + y dh j + z dh khhhdz = x dz i + y dzj + zdz kzzz中的微元体积 dxdydz = dx (dh dz ) =dx dhdz则体J111： Gdxdydz = G(x ,h,z ) Jdx dhdz六面体单元体-1-1-1V区间虽然简单，但却无法给出G的显示表达。除最简上式的单的单元外，大部分单元的都无

9、法给出。2.单元矩阵的计算(局部坐标下的变换) 如何计算三维问题中的面？发生在z =常数的表面上设面ijkixxjyxkzxx dxy dxz dxdx dhdx dh则 dS =xxxz =constx dhy dhz dhxhyhzhhhh= Adx dh其中 A =+ b2 + c2a2a = y z - y zb = x z - x zc = x y- xyx hh xh xx hx hh xA并非面积2.单元矩阵的计算(局部坐标下的变换)11(z =常数的表面上)： gdS = g(x ,h) Adx dh面-1-1S二维单元下的面积:z = zx y - x y= 0A =xhx

10、hh x2 面2.单元矩阵的计算(局部坐标下的变换) 如何计算二维问题的线？例如：在x =常数的曲线上dl = dh = x dhi + y dh jhh( x )2+ ( y )2 dhdl =dlhh3.面积坐标、体积坐标下的矩阵计算z= L- L34雅(L3cons体积 ： GdV =V4.例题 例1：，边界为直线的三角形二维6节 3x = L1h = L2-1)-1)1 - x -h = L35= L (2L= L (2L -1)N N 1112226= L (2L= 4L LyN N 3334122N = 4L LN = 4L L5236134N N N N 1=-= 4L - 1=

11、 0 1 1 1 2xxx1LL13N N N N N N = 1 - 4L= 4L=-= 0= 4L-1 3x 4x 1 1 1 2h32h2LL2330点单元的边中点正处在中间，试证明此单元的Jacobi矩阵是常数矩阵。4.例题由于边中点正好处x4 =12之后314. 例题例2：的平面8节点单元，试计1 (40,50)5(22.5,45)Q2 (5,40)8 (35,35)提示：由于所有边中点正好处在中点位置，故可以用4节点定义几何形状！6 (7.5,25)y4 (30,20)7(20,15)3(10,10)x算N1 , N2 在自然坐标Q(1 , 1)点的值。xx2 24.例题N = 1

12、 (1+ x )(1+h)(x +h -1)N= 1 (1- x )(1+h)(-x +h -1)1244N1N2= 1 (1+h)(2x +h)= - 1 (1+h)(-2x +h) x x414NN11 h2 h=(1+ x )(2h + x )=(1- x )(2h - x )44形状变换函数：= 1 (1+ x )(1+h)= 1 (1- x )(1+h)N N 1244= 1 (1- x )(1-h)= 1 (1+ x )(1-h)N N 34444.例题N N 1 2 x = xJ N N Q点1 2 h h- 32- 12- 124 3 2121 -70 0-=J5-112 11

13、24.例题NN1 -1 x x NQ点N1 h yNN2-1 x x NQ点N2 y h 4.例题1 (40,50) 例3：，1-4边作 4(22.5,45)2 (5,40)18 (35,35)荷载作用在1 - 4边上，故节点等效力只与1,4,8号节点有关6 (7.5,25)N = 1 (1 + x )(1 + h)(x + h - 1)14y4(30,20)N= 1 (1 + x )(1 -h)(x -h -1)7(20,15)443(10,10)N= 1 (1 + x )(1 -h 2 )x82用有x方向的均布荷载，线密度为1，求节点等效力。4.例题在x = 1边上计算NiNx = Nx

14、+ N+ N x+ Nxx1 1223344在y = N1 y1 + N2 y2 + N3 y3 + N4 y4 x 2 y 2qN T?s= 5 10dhx dsds =+R=e h h s4.例题1=N q ds =R1x1x-s11N q ds =R44x-s11N q ds =R8x-s1注意：1. 本例题但是，2. 如果边5.边中点位置对Jacobi矩阵的影响正方形单元中间节点在1/2处，Jacobi行列式为常数1/4距离变化10.250.2510.250.50.500-0.5-0.5-1-1-1-0.500.51-1-0.500.511 0.250.2510.50.500-0.5-

15、0.5-1-1-1-0.500.51-1-0.500.51= 1/ 4 + x (1-h) / 61/6处J= 1/ 4 + 3x (1-h) / 20J1/5处0.250.250.250.250.250.250.250.250.250.250.250.251/4处 J = 1/ 4 + x (1-h) / 81/3处 J = 1/ 4 + x (1-h) /12110.250.50.500-0.5-0.5-1-1-1-0.500.51-1-0.500.51110.50.500-0.5-0.5-1-1-1-0.500.51-1-0.500.51= 1/ 4 + 7x (1-h) / 301/3

16、0处J= 1/ 4 + 9x (1-h) / 40J1/20处0.250.250.250.250.250.250.250.250.25 0.250.250.251/10处 J = 1/ 4 + x (1-h) / 51/7处 J = 1/ 4 + 5x (1-h) / 285.边中点位置对Jacobi矩阵的影响长方形单元单元长宽比10:1距离变化在1/2处Jacobi行列式为常数5/2110.50.500-0.5-0.5-1-1-1-0.500.51-1-0.500.5101-0.20.5-0.40-0.6-0.5-0.8-1-1-1-0.8-0.6-0.4-0.20-1-0.500.511/

17、5处 J = 5/ 2 + 3x (1-h) / 21/6处 J = 5/ 2 + 5x (1-h) / 31/4处J= 5/ 2 + 5x (1-h) / 41/3处 J = 5/ 2 + 5x (1-h) / 600-0.2-0.2-0.4-0.4-0.6-0.6-0.8-0.8-1-1-1-0.8-0.6-0.4-0.20-1-0.8-0.6-0.4-0.2000-0.2-0.2-0.4-0.4-0.6-0.6-0.8-0.8-1-1-1-0.8-0.6-0.4-0.20-1-0.8-0.6-0.4-0.202.42.41/20处 J = 5/ 2 + 9x (1-h) / 41/30处

18、 J = 5/ 2 + 7x (1-h) / 31/10处 J = 5/ 2 + 2x (1-h)1/7处 J = 5/ 2 + 25x (1-h) /145.边中点位置对Jacobi矩阵的影响长方形单元单元长宽比10:1距离变化结果相当于把前一种的和对调在1/2处Jacobi行列式为常数5/2110.50.500-0.5-0.5-1-1-1-0.500.51-1-0.500.51110.50.500-0.5-0.5-1-1-1-0.500.51-1-0.500.511/5处 J = 5/ 2 + 3h(1- x ) / 21/6处 J = 5/ 2 + 5h(1- x ) / 31/4处 J

19、 = 5/ 2 + 5h(1- x ) / 41/3处 J = 5/ 2 + 5h(1- x ) / 6110.50.500-0.5-0.5-1-1-1-0.500.51-1-0.500.51110.50.500-0.5-0.5-1-1-1-0.500.51-1-0.500.512.41/20处 J = 5/ 2 + 9h(1- x ) / 41/30处 J = 5/ 2 + 7h(1- x ) / 31/10处 J = 5/ 2 + 2h(1- x )1/7处 J = 5/ 2 + 25h(1- x ) /146.4数值1.Newton-Cotes 基本思想bf (x )dx;I =求a构造

20、一个多项式y（x），使在xi (i = 1, 2, n)上有y（xi）=f (xi )，然后bb利用近似函数y（xi）的y (x )dx来近似原被积函数x )f (x )dx ,f (的aaxi称为点或者取样点。点xi的数目和位置决定了y (x )近似f (x )的程度，因而决定了数值的精度。491.Newton-Cotes1f (x )dxI =-1，x1、x2xn在-1,1区间等间隔分布点的Newton - Cotes对n个ni=1y (x ) =(n-x ) f (x )1)l(被积函数的近似多项式：(n -1次多项式）iiy(xi ) =1f (xi )x )dx =显然nI =yH

21、f (xi )(i-1i=11x )dxH =(n-1)其中l(ii-11.Newton-Cotes当：n = 2 :n = 3 :I = f (-1) + f (1)I = 1 f (-1) + 4 f (0) + f (1)3I = 1 f (-1) + 3 f (- 1) + 3 f ( ) + f (1)1n = 4 :43误差：O(Dn )32.高斯点位置的Newton - Cotes点的精度可以达到2n -1次，误差为：O(D2n )优化了n个布点：xi(i = 1, 2, n)确定点的方法：1x kP(x )dx = 0(k = 0,1, 2,n -1)-1n其中 P(x ) =

22、 (x - xi )i=1P(x )的特点：P(xi ) = 0P(x )与x 0,x1,x n-1在-1,1区间内正交2.高斯被积函数用2n-1次多项式y (x )近似代替n-1ni=1k =0y (x ) =x ) f (x ) +b x P(x )(n-1)kl(iikn -1次2n -1次显然满足：y (xi ) = f (xi )n1y (x )dx =H f (xi )I =i-1i=11(n-x )dxH =1)l(其中ii-12.求两点高斯高斯的例：点位置与权值P(x ) = (x -x1)(x -x2 )11x )dx = 0x P(x )dx = 0P(由-111x = -

23、x= -0.57735 = -123x - x2x -x111dx = 1dx = 1H =H=x - xx - x12-1112212.高斯特点：点xi不等间隔(1)(2) y(x )是2n -1次多项式，精度达到2n-1阶nn11f (x ,h)dx dh = Hi H j f (xi ,h j )i=1 j=1 二维： I =-1-1nnnf (x,h,z )dxdhdz = Hi H j Hk f (xi ,h j ,z k )i=1 j=1 k =1111 三维： I =-1-1-1二维高斯公式Hi H j f (x j ,hi )nn11f (x ,h)dxdh = -11i=1i

24、=1式中点和权函数仍按上表采用。ll单元刚度矩阵k = H BT DB | J | hHmnm=1n=1qVx llm =1n =1体积力等效结点力：RV i=Hm Hn Ni qh | J|Vy q2 x 2 y lm =1Sxh=+ 面力等效结点力：R HN SimihhqSy 等参元数值中一般取2-3就可取得足够精度3.高斯的另一种推导例1：3点高斯，应该精确到5阶，而在-1,1上奇数阶的都等于0-1权：-AW10W2A1：1x2x4只需精确即可高斯的结果精确的结果1+W2=1 dx = 22W1 A = 0.774596669241483W1 = 0.55555555556-1x2dx

25、 = 231+ 0=22 A W1-1W= 0.8888888888922 A4W251=x4dx =1-13.高斯的另一种推导例2：平面单元的8点高斯(1,1)到5阶）（可精确(B,0)对应权为W2(-1,-1)(A,-A) 对应权为W1需精确以下多项式项：1xy x3y2x2xyy2 y5x3x2yxy2y3x4xy3x2y2x3yy4x5xy4x2y3x4y在-1,1上带奇次项的结果都等于0：1x2x4x2y2只需精确4.高斯的另一种推导的结果高斯的结果精确4W + 4W=12W9=4 A W + 2B2W=2112497+=4W4W A =234 A4W=18斯的最9点高斯的最4点高斯

26、的最8点高斯最容6.5数值阶次的选择1.保证收敛精度的阶次1 完全（精确）,必然要引起额外的误差,因此的计算量还是很大的,因此采用数值来代替精确要尽量减少这一误差值.采用数值需要做到:1)在保证收敛的情况下,如何尽量减少数值;2)在逼近精确解时,用什么条件来保证其收敛速度.工作量插值函数N的完全多项式阶次：p微分算子L的导数阶次：m这种并不实用收敛性要求：能再现m阶导数的任一常数值对C0型单元，m = 1，只要求V dV被正确即可,因此只要单点即可保证收敛性;在曲线坐标系中,则需要对VJ dx dhdV dV进行正确.1.保证收敛精度的阶次精度可知，要使P阶被积函数达成完全（精确）根据高斯点的

27、数目n为:2n -1 P,即n P +1，需要2插值函数包含1、 项，刚度矩阵被积项包含1、 2、2、 ，n P +1 = 2 +1 = 1.5, 所以n = 222点数目为2 2;因为为平面问题，所以例如：二维四节点平面线性单元，若J=const，其刚度矩阵完全则需要多少高斯点?1.保证收敛精度的阶次2 减缩（降阶）中,能够保证不降低收敛速度的条件下求解各种有限元问题的在数值最小阶次,比精确低阶的可称为减缩1。dxBT DB J例如,在一维问题中：-1P=2(p-m)次多项式= const,只要求达到2( p - m)阶精度即可点 n = p - m +1，精度可达2(p-m)+1有xy项存

28、在若 J高斯然而在二维或者三维单元中，微分算子L并不能使被积函数的阶次降低m次，BT DB的阶次比2( p - m)高。三维情况更复杂一些。不管一维、二维、三维，都采用n=p-m+1来确定积而减缩分阶次1.保证收敛精度的阶次例：二维四节点平面线性单元，p-m+1=1，即采用11阶的单点高斯采用减缩（降阶）往往可以得到更好的解！是由非完全项决定的，而有限元的离原因：(1) 精确的散精度是由完全多项式决定的。精确常常是由插值函数中非完全项的最高方次要求，而决定有限元精度的是完全多项式的方次。这些非完全的最高方次项往往不能很好地提高精度，反而可能带来不好的影响。取较低阶的高斯，使精度正好保证完全多项

29、式方次的要求，而不包括更高次的非完全多项式的要求，其实质是相当用一种新的插值函数替代原来的插值函 数，从而一定情况下改善了单元的精度。1.保证收敛精度的阶次(2) 有限元解偏刚，而减缩两者由互补关系基于最小位能原来基础上建立的位移有限元，其解答具有下限性质。即有限元的计算模型具有较实际结构偏大的整（降阶）使模型刚度降低，体刚度。选取缩减方案将使有限元计算模型的刚度有所降低，因此可能有助于提高计算精度。另外，这种缩减方案对于泛函中包含罚函数的情况也常常是必须的，用以保证和罚函数相应的矩阵的奇异性（见相应），否则将可能导致完全歪曲了的结果。2.由数值引起的矩阵奇异性容易理解：若未知量的数目大于所有点处提供的关系的数目，则矩阵K必定奇异。（奇异的充分条件）关系的数目等于D的维数（一个高斯当于一组方程）从物理意义上容易理解，当然也可以用代数的秩的关系 进行严格的证明：点就相Kd = P已引入边条件2.由数值引起的矩阵奇异性ngK e = HBDB JTiiiii=1点的个数，D为d d的矩阵，秩一般为d，B是d nf 的矩阵ng为一