第6章曲邊等參單元和數(shù)值積分

1、計(jì)算力學(xué)第六章單元和數(shù)值本章目錄單元基本概念6.16.2 參數(shù)曲線坐標(biāo)6.3 雅可比矩陣6.4 數(shù)值單元矩陣計(jì)算6.5 數(shù)值階次的選擇*6.6 應(yīng)力結(jié)果的處理6.1單元基本概念 1引言前面介紹了C0插值獲得一些常用單元，這類單元在實(shí)際應(yīng)用中有一個(gè)無法回避的缺點(diǎn)：在不增加計(jì)算工作量的 情況下，采用增加節(jié)點(diǎn)數(shù)方法提高單元計(jì)算精度，必然 減少單元數(shù)量，描述相對(duì)復(fù)雜的幾何形狀的能力也隨之 下降。這是一個(gè)實(shí)際問題，也即僅采用簡單的三角形和四邊形單元是不夠的，因此，本章討論由簡單單元轉(zhuǎn)換 為其他復(fù)雜形狀單元的問題。1單元基本概念如下圖所示，二維單元可被為扭曲的單元形狀。(-1,1)(1,1)yx(-1,-

2、1)(1,-1)畸變單元(實(shí)際單元)基本單元變換函數(shù)x 總體坐標(biāo)局部坐標(biāo) x y = f h z z 可以看出，局部坐標(biāo)系（x ,h )或者（L1,L2,L3)，在笛卡爾（x,y)空間中被為新的曲線形狀。6.2參數(shù)曲線坐標(biāo)1坐標(biāo)變換形函數(shù)N1很明顯，基本單元和畸變單元從形狀上就可以看出不是一樣的單元，要在兩者之間建立，首先應(yīng)進(jìn)行坐標(biāo)變換，建立坐標(biāo)變換最方便的辦法就是使用我們最熟悉的，用以描述 位移函數(shù)（或其他場函數(shù)）變化的標(biāo)準(zhǔn)C0形狀函數(shù)。局部坐標(biāo)與整體坐標(biāo)的變換局部坐標(biāo)系(基準(zhǔn)坐標(biāo)系)整體坐標(biāo)系(物理坐標(biāo)系)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是基于局部坐標(biāo)給出的標(biāo)準(zhǔn)形狀函數(shù),于是可以馬上得到需要的變換關(guān)系.可以看出,坐標(biāo)(x1,y1).的點(diǎn)正好就是單元邊界上的對(duì)應(yīng)點(diǎn)(從標(biāo)準(zhǔn)形狀函數(shù)的定義可知,在對(duì)應(yīng)點(diǎn)上函數(shù)值為1,而在其他點(diǎn)為零)，可在這些點(diǎn)上預(yù)先建立節(jié)點(diǎn)。1將坐標(biāo)變換式用于任意四邊形單元，可得：12在四個(gè)結(jié)點(diǎn)處給出結(jié)點(diǎn)的整體坐標(biāo)。在四條邊上的整體坐標(biāo)是線性變化的。只要給出任意四邊形單元四個(gè)結(jié)點(diǎn)的整體坐標(biāo)，用（6-1）式就可以建立局部坐標(biāo)系中的正方形單元和整體坐標(biāo)系中的任意四邊形單元之間的坐標(biāo)變換關(guān)系

4、。關(guān)鍵在于，首先要解決進(jìn)行任意形狀單元和規(guī)則單元之間的幾何變換。變換的關(guān)鍵是什么？1mmmx =y =z =N xN yN ziiiiiii =1i =1i =1局部坐標(biāo)下的每一點(diǎn)都對(duì)應(yīng)整個(gè)笛卡爾坐標(biāo)的一個(gè)點(diǎn)，而且一般說來變換是唯一的。但若出現(xiàn)對(duì)應(yīng)的非唯一性，則將引起單元的劇烈畸變一個(gè)變換對(duì)應(yīng)有一個(gè)雅可比（Jacobian）行列式，可由下式定義：D(x, y, z)D(x ,h,z )J = 即雙線性單元不能有內(nèi)角大于或等于或接近180度情況。雅可比行列式，具體定義可見6.3節(jié)。 0J2sin(dx , dh) 01 n 亞參數(shù)單元：m 0J194 關(guān)于等參數(shù)單元說明總結(jié)3） 等參單元的優(yōu)點(diǎn)是

5、當(dāng)單元邊界呈二次以上的曲線時(shí)， 容易用很少的單元去逼近曲線邊界。將不規(guī)則單元轉(zhuǎn)換為規(guī)則母單元后，容易構(gòu)造位移函數(shù)和形函數(shù)。4） 上述等參單元的理論公式可適應(yīng)三次以上的曲線型等 參元，只是階次提高，單元自由度相應(yīng)增加，計(jì)算更復(fù)雜，更，實(shí)際中，很少超過3次曲線型。5） 上述推導(dǎo)要求：保持坐標(biāo)變換中幾何模式階次與描述 單元位移函數(shù)中形函數(shù)的階次相同。如取坐標(biāo)變換的幾何模式階次較單元的位移函數(shù)的階次高，則稱此單元為超參單元，反之，為亞參單元。這兩類單元的收斂性也可得到滿足。6） 當(dāng)然，也可取描述單元幾何形狀的幾何模式不是形函 數(shù)的，如p-element6.3雅可比矩陣計(jì)算單元矩陣1.雅可比矩陣有限元分

6、析中，建立求解方程，需要進(jìn)行各個(gè)單元體積內(nèi)和面積內(nèi)的，其求解的形式：V GdV =G(x, y, z)dxdydz體：面：SgdS = g(x, y, z)dSG、g的形式取決于不同的單元，通常包含Ni對(duì)x, y, z的導(dǎo)數(shù)NiNiNi= ?= ?= ?xyz1.雅可比矩陣求導(dǎo)的計(jì)算(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.單元矩陣的計(jì)算(局部

8、坐標(biāo)下的變換)1 體為對(duì)過程中的變量和區(qū)域進(jìn)行變換，需要相應(yīng)的包括雅可比行列式在內(nèi)的標(biāo)準(zhǔn)處理方法。例如，體積元的計(jì)算。設(shè) x, y, z坐標(biāo)系的軸為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區(qū)間雖然簡單，但卻無法給出G的顯示表達(dá)。除最簡上式的單的單元外，大部分單元的都無

9、法給出。2.單元矩陣的計(jì)算(局部坐標(biāo)下的變換) 如何計(jì)算三維問題中的面？發(fā)生在z =常數(shù)的表面上設(shè)面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.單元矩陣的計(jì)算(局部坐標(biāo)下的變換)11(z =常數(shù)的表面上)： gdS = g(x ,h) Adx dh面-1-1S二維單元下的面積:z = zx y - x y= 0A =xhx

10、hh x2 面2.單元矩陣的計(jì)算(局部坐標(biāo)下的變換) 如何計(jì)算二維問題的線？例如：在x =常數(shù)的曲線上dl = dh = x dhi + y dh jhh( x )2+ ( y )2 dhdl =dlhh3.面積坐標(biāo)、體積坐標(biāo)下的矩陣計(jì)算z= L- L34雅(L3cons體積 ： GdV =V4.例題 例1：，邊界為直線的三角形二維6節(jié) 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點(diǎn)單元的邊中點(diǎn)正處在中間，試證明此單元的Jacobi矩陣是常數(shù)矩陣。4.例題由于邊中點(diǎn)正好處x4 =12之后314. 例題例2：的平面8節(jié)點(diǎn)單元，試計(jì)1 (40,50)5(22.5,45)Q2 (5,40)8 (35,35)提示：由于所有邊中點(diǎn)正好處在中點(diǎn)位置，故可以用4節(jié)點(diǎn)定義幾何形狀！6 (7.5,25)y4 (30,20)7(20,15)3(10,10)x算N1 , N2 在自然坐標(biāo)Q(1 , 1)點(diǎn)的值。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形狀變換函數(shù)：= 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點(diǎn)1 2 h h- 32- 12- 124 3 2121 -70 0-=J5-112 11

13、24.例題NN1 -1 x x NQ點(diǎn)N1 h yNN2-1 x x NQ點(diǎn)N2 y h 4.例題1 (40,50) 例3：，1-4邊作 4(22.5,45)2 (5,40)18 (35,35)荷載作用在1 - 4邊上，故節(jié)點(diǎn)等效力只與1,4,8號(hào)節(jié)點(diǎn)有關(guān)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，求節(jié)點(diǎn)等效力。4.例題在x = 1邊上計(jì)算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.邊中點(diǎn)位置對(duì)Jacobi矩陣的影響正方形單元中間節(jié)點(diǎn)在1/2處，Jacobi行列式為常數(shù)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.邊中點(diǎn)位置對(duì)Jacobi矩陣的影響長方形單元單元長寬比10:1距離變化在1/2處Jacobi行列式為常數(shù)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.邊中點(diǎn)位置對(duì)Jacobi矩陣的影響長方形單元單元長寬比10:1距離變化結(jié)果相當(dāng)于把前一種的和對(duì)調(diào)在1/2處Jacobi行列式為常數(shù)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數(shù)值1.Newton-Cotes 基本思想bf (x )dx;I =求a構(gòu)造

20、一個(gè)多項(xiàng)式y(tǒng)（x），使在xi (i = 1, 2, n)上有y（xi）=f (xi )，然后bb利用近似函數(shù)y（xi）的y (x )dx來近似原被積函數(shù)x )f (x )dx ,f (的aaxi稱為點(diǎn)或者取樣點(diǎn)。點(diǎn)xi的數(shù)目和位置決定了y (x )近似f (x )的程度，因而決定了數(shù)值的精度。491.Newton-Cotes1f (x )dxI =-1，x1、x2xn在-1,1區(qū)間等間隔分布點(diǎn)的Newton - Cotes對(duì)n個(gè)ni=1y (x ) =(n-x ) f (x )1)l(被積函數(shù)的近似多項(xiàng)式：(n -1次多項(xiàng)式）iiy(xi ) =1f (xi )x )dx =顯然nI =yH

21、f (xi )(i-1i=11x )dxH =(n-1)其中l(wèi)(ii-11.Newton-Cotes當(dāng)：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.高斯點(diǎn)位置的Newton - Cotes點(diǎn)的精度可以達(dá)到2n -1次，誤差為：O(D2n )優(yōu)化了n個(gè)布點(diǎn)：xi(i = 1, 2, n)確定點(diǎn)的方法：1x kP(x )dx = 0(k = 0,1, 2,n -1)-1n其中 P(x ) =

22、 (x - xi )i=1P(x )的特點(diǎn)：P(xi ) = 0P(x )與x 0,x1,x n-1在-1,1區(qū)間內(nèi)正交2.高斯被積函數(shù)用2n-1次多項(xiàng)式y(tǒng) (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.求兩點(diǎn)高斯高斯的例：點(diǎn)位置與權(quán)值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.高斯特點(diǎn)：點(diǎn)xi不等間隔(1)(2) y(x )是2n -1次多項(xiàng)式，精度達(dá)到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式中點(diǎn)和權(quán)函數(shù)仍按上表采用。ll單元?jiǎng)偠染仃噆 = H BT DB | J | hHmnm=1n=1qVx llm =1n =1體積力等效結(jié)點(diǎn)力：RV i=Hm Hn Ni qh | J|Vy q2 x 2 y lm =1Sxh=+ 面力等效結(jié)點(diǎn)力：R HN SimihhqSy 等參元數(shù)值中一般取2-3就可取得足夠精度3.高斯的另一種推導(dǎo)例1：3點(diǎn)高斯，應(yīng)該精確到5階，而在-1,1上奇數(shù)階的都等于0-1權(quán)：-AW10W2A1：1x2x4只需精確即可高斯的結(jié)果精確的結(jié)果1+W2=1 dx = 22W1 A = 0.774596669241483W1 = 0.55555555556-1x2dx

25、 = 231+ 0=22 A W1-1W= 0.8888888888922 A4W251=x4dx =1-13.高斯的另一種推導(dǎo)例2：平面單元的8點(diǎn)高斯(1,1)到5階）（可精確(B,0)對(duì)應(yīng)權(quán)為W2(-1,-1)(A,-A) 對(duì)應(yīng)權(quán)為W1需精確以下多項(xiàng)式項(xiàng)：1xy x3y2x2xyy2 y5x3x2yxy2y3x4xy3x2y2x3yy4x5xy4x2y3x4y在-1,1上帶奇次項(xiàng)的結(jié)果都等于0：1x2x4x2y2只需精確4.高斯的另一種推導(dǎo)的結(jié)果高斯的結(jié)果精確4W + 4W=12W9=4 A W + 2B2W=2112497+=4W4W A =234 A4W=18斯的最9點(diǎn)高斯的最4點(diǎn)高斯

26、的最8點(diǎn)高斯最容6.5數(shù)值階次的選擇1.保證收斂精度的階次1 完全（精確）,必然要引起額外的誤差,因此的計(jì)算量還是很大的,因此采用數(shù)值來代替精確要盡量減少這一誤差值.采用數(shù)值需要做到:1)在保證收斂的情況下,如何盡量減少數(shù)值;2)在逼近精確解時(shí),用什么條件來保證其收斂速度.工作量插值函數(shù)N的完全多項(xiàng)式階次：p微分算子L的導(dǎo)數(shù)階次：m這種并不實(shí)用收斂性要求：能再現(xiàn)m階導(dǎo)數(shù)的任一常數(shù)值對(duì)C0型單元，m = 1，只要求V dV被正確即可,因此只要單點(diǎn)即可保證收斂性;在曲線坐標(biāo)系中,則需要對(duì)VJ dx dhdV dV進(jìn)行正確.1.保證收斂精度的階次精度可知，要使P階被積函數(shù)達(dá)成完全（精確）根據(jù)高斯點(diǎn)的

27、數(shù)目n為:2n -1 P,即n P +1，需要2插值函數(shù)包含1、 項(xiàng)，剛度矩陣被積項(xiàng)包含1、 2、2、 ，n P +1 = 2 +1 = 1.5, 所以n = 222點(diǎn)數(shù)目為2 2;因?yàn)闉槠矫鎲栴}，所以例如：二維四節(jié)點(diǎn)平面線性單元，若J=const，其剛度矩陣完全則需要多少高斯點(diǎn)?1.保證收斂精度的階次2 減縮（降階）中,能夠保證不降低收斂速度的條件下求解各種有限元問題的在數(shù)值最小階次,比精確低階的可稱為減縮1。dxBT DB J例如,在一維問題中：-1P=2(p-m)次多項(xiàng)式= const,只要求達(dá)到2( p - m)階精度即可點(diǎn) n = p - m +1，精度可達(dá)2(p-m)+1有xy項(xiàng)存

28、在若 J高斯然而在二維或者三維單元中，微分算子L并不能使被積函數(shù)的階次降低m次，BT DB的階次比2( p - m)高。三維情況更復(fù)雜一些。不管一維、二維、三維，都采用n=p-m+1來確定積而減縮分階次1.保證收斂精度的階次例：二維四節(jié)點(diǎn)平面線性單元，p-m+1=1，即采用11階的單點(diǎn)高斯采用減縮（降階）往往可以得到更好的解！是由非完全項(xiàng)決定的，而有限元的離原因：(1) 精確的散精度是由完全多項(xiàng)式?jīng)Q定的。精確常常是由插值函數(shù)中非完全項(xiàng)的最高方次要求，而決定有限元精度的是完全多項(xiàng)式的方次。這些非完全的最高方次項(xiàng)往往不能很好地提高精度，反而可能帶來不好的影響。取較低階的高斯，使精度正好保證完全多項(xiàng)

29、式方次的要求，而不包括更高次的非完全多項(xiàng)式的要求，其實(shí)質(zhì)是相當(dāng)用一種新的插值函數(shù)替代原來的插值函 數(shù)，從而一定情況下改善了單元的精度。1.保證收斂精度的階次(2) 有限元解偏剛，而減縮兩者由互補(bǔ)關(guān)系基于最小位能原來基礎(chǔ)上建立的位移有限元，其解答具有下限性質(zhì)。即有限元的計(jì)算模型具有較實(shí)際結(jié)構(gòu)偏大的整（降階）使模型剛度降低，體剛度。選取縮減方案將使有限元計(jì)算模型的剛度有所降低，因此可能有助于提高計(jì)算精度。另外，這種縮減方案對(duì)于泛函中包含罰函數(shù)的情況也常常是必須的，用以保證和罰函數(shù)相應(yīng)的矩陣的奇異性（見相應(yīng)），否則將可能導(dǎo)致完全歪曲了的結(jié)果。2.由數(shù)值引起的矩陣奇異性容易理解：若未知量的數(shù)目大于所有點(diǎn)處提供的關(guān)系的數(shù)目，則矩陣K必定奇異。（奇異的充分條件）關(guān)系的數(shù)目等于D的維數(shù)（一個(gè)高斯當(dāng)于一組方程）從物理意義上容易理解，當(dāng)然也可以用代數(shù)的秩的關(guān)系 進(jìn)行嚴(yán)格的證明：點(diǎn)就相Kd = P已引入邊條件2.由數(shù)值引起的矩陣奇異性ngK e = HBDB JTiiiii=1點(diǎn)的個(gè)數(shù)，D為d d的矩陣，秩一般為d，B是d nf 的矩陣ng為一