彈性動(dòng)力學(xué)的邊界無(wú)單元法

中國科學(xué)G輯物理學(xué)力學(xué)天文學(xué)2005,35(4):435-448435彈性動(dòng)力學(xué)的邊界無(wú)單元法程玉民彭妙娟上海大學(xué)上海市應用數學(xué)和力學(xué)研究所,上海200072)上海大學(xué)土木工程系,上海200072)摘要討論了 Hilbert空間上的改進(jìn)的移動(dòng)最小二乘法,并將彈性動(dòng)力學(xué)的邊界積分方程方法與改進(jìn)的移動(dòng)最小二乘法結合,提出了彈性動(dòng)力學(xué)的邊界無(wú)單元法.改進(jìn)的移動(dòng)最小二乘法避免了求解病態(tài)方程組,既提高了效率,又提高了精度.邊界無(wú)單元法是邊界積分方程的無(wú)網(wǎng)格方法的直接列式法,容易引入邊界條件,且具有更高的精度.最后給出了彈性動(dòng)力學(xué)的邊界無(wú)單元法的數值實(shí)現方法和具體的算例關(guān)鍵詞移動(dòng)最小二乘法邊界積分方程無(wú)網(wǎng)格方法彈性動(dòng)力學(xué)改進(jìn)的移動(dòng)最小二乘法邊界無(wú)單元法無(wú)網(wǎng)格方法近年來(lái)得到了較大發(fā)展,是目前計算力學(xué)研究的熱點(diǎn)之無(wú)網(wǎng)格方法采用基于點(diǎn)的近似,在處理諸如大變形、裂紋擴展等問(wèn)題時(shí)不需進(jìn)行網(wǎng)格重構,具有較為明顯的優(yōu)越性.目前發(fā)展的無(wú)網(wǎng)格方法有擴散單元法( Diffuse element method,即DEM)、無(wú)單元 Galerkin法( lement-Free GalerkinMethod,即EFG)、Hp- -clouds無(wú)網(wǎng)格方法、有限點(diǎn)法( The Finite Point method,即FPM)、無(wú)網(wǎng)格局部 Petrov-Galerkin方法( Meshless local petrov- Galerkin Method,即MLPG)、多尺度重構核粒子方法(Mult- Scale Reproducing Kernel ParticleMethod,即MRKP)、小波粒子方法( Wavelet Particle Method,即WPM)、徑向基函數法( Radial basis Functions,即RBF)、無(wú)網(wǎng)格有限元法( Meshless finite elementMethod,即MFEM)、移動(dòng)粒子有限元法( Moving particle Finite Element Method,即 MPFEM)以及邊界積分方程的無(wú)網(wǎng)格方法等邊界積分方程的無(wú)網(wǎng)格方法是將邊界積分方程方法與無(wú)網(wǎng)格方法中的移動(dòng)最小二乘法結合而形成的.目前這一方面的研究主要有 Mukherjee等人提出的邊界點(diǎn)法( Boundary Node method,即BNM)3、Atui等人提出的局部邊界積分方004-08-02收稿,2005-04-22收修改稿中國煤化工*上海市重點(diǎn)學(xué)科建設項目資助(批準號:Y0103)*sE-mail:yecheng@sh163net:yecheng@staff.shu.edu.cnCNMHGSCIENCE IN CHINA Ser. G Physics, Mechanics Astronomy436中國科學(xué)G輯物理學(xué)力學(xué)天文學(xué)第35卷程方法 Local Boundary Integral Equation Method,即LBI)、姚振漢等人提出的雜交邊界點(diǎn)法( Hybrid Boundary Node Method)9、以及本文作者提出的邊界無(wú)單元法( Boundary Element-Free Method,即BEFM)等與邊界點(diǎn)法和局部邊界積分方程方法相比,邊界無(wú)單元法在形成的邊界積分方程中采用節點(diǎn)變量的真實(shí)解為未知量,而邊界點(diǎn)法和局部邊界積分方程方法均采用節點(diǎn)變量的近似解為未知量.邊界無(wú)單元法是邊界積分方程無(wú)網(wǎng)格方法的直接解法.邊界無(wú)單元法不僅直觀(guān)而且具有較高的精度,還可方便地引入邊界條件本文在 Hilbert空間上討論了改進(jìn)的移動(dòng)最小二乘法,并將彈性動(dòng)力學(xué)的邊界積分方程方法與改進(jìn)的移動(dòng)最小二乘法結合,提出了彈性動(dòng)力學(xué)的邊界無(wú)單元法,最后給出了彈性動(dòng)力學(xué)的邊界無(wú)單元法的數值實(shí)現方法和具體的算例1移動(dòng)最小二乘法在移動(dòng)最小二乘法( Moving least- square approximation)中,取試函數n(x)=∑n(x)a(x)=p(x)a(x)(x∈9)(1)為函數u(x)的逼近函數.這里m是基函數的個(gè)數,P(x)是基函數,a1(x)是相應的系數.基函數的通常形式為線(xiàn)性基:p1=(1,x)(對一維區域)p=(,x,x2)(對二維區域二次基p2=(1,x1,x2)(對一維區域p=(1,x1,x2,x,xx2,x2)(對二維區域)對應于(1)式的整體逼近,在點(diǎn)x的鄰域內的局部逼近定義為n'(x,)=∑(x)a1(x)=p(x)a(x),(6)其中系數a1(x)根據加權最小二乘法來(lái)確定,它使得對函數n(x)的局部逼近誤差最小定義J=2w(r-x,u(r,x)-u(xn)w(x-x1)∑P(x7)a中國煤化工其中wx-x1)是具有緊支集特性的權函數,xCNMH點(diǎn)SCIENCE IN CHINA Ser. G Physics, Mechanics Astronomy第4期程玉民等:彈性動(dòng)力學(xué)的邊界無(wú)單元法437(7)式可用矩陣形式表示為J=(Pa-u) w(r(pa-u)(8)其中(r,) p2(n)l pm(n)P=|()P2()LPm(與)(10)p,(An) P2(r,n)L Pm(n)L為了得到a(x),對J取極值,即得A(x)a(x)-B(r)u=0(12)A(x)a(r)=b(x )u(13)其中矩陣A(x)和B(x)分別為A(r)=P w(r)P(14)B(x)=Pw(x)(15)由(13)式,可得a(r)=A (xB(x)u(16)這樣,逼近函數u(x)的表達式為x)=φ(x)m=∑(17)其中Φ(x)為形函數Φ(x)=閩(x,2(x),,2(x)=p(x)A(x)B(x)(18)2改進(jìn)的移動(dòng)最小二乘法在移動(dòng)最小二乘法中,當m較大時(shí),方程(13)有時(shí)是病態(tài)的,甚至是奇異的這樣,方程(13就難以求解或獲得正確的解.若選取正交函數作為基函數,則所得到的方程即不病態(tài)也不奇異,而且不需求中國煤化工方程組的解THCNMHG中國科學(xué)G輯物理學(xué)力學(xué)天文學(xué)第35卷21 Hilbert空間span(p)為了討論正交基函數,我們先來(lái)證明pam(p)是一個(gè) Hilbert空對vf(x),g(x)∈spmn(p),定義I=1那么(f,g)是一個(gè)內積,spmn(p)是一個(gè)內積空間.下面我們給出證明(a)由(19)式可得(f,g)=(g,f)(b)和β為常數,那么(f+βg,h)=(f,h)+β(g,h),(21)(c)對wf(x)∈spmn(p),可得1)f(x1)2≥0(22)若(f,f)=0,那么w(x-x1)=0或f(x1)=0,I=1,2,…,N.在移動(dòng)最小二乘法中(x-x1)≠0、所以我們有f(x1)=0,I=1,2,…,N所以(f,g)是一個(gè)內積由于spm(p)是一個(gè)線(xiàn)性空間,則spmn(p)是一個(gè)內積空間對內積空間 spani(p)來(lái)說(shuō),如果它是完備的,那么它是一個(gè) Hilbert空間.現在我們來(lái)證明spn(p)是完備的設{(x)為span(p)中的任一收斂級數,且lim fn (x)=f(r)(23)將fn(x)表示為fn(r)1,2,L其中an為常數.那么mf(x)=lm∑ankP(x)=∑(imak)(x)由于極限limf1(x)存在,則極限 lim ank存在,即有lim可得中國煤化工CNMHGSCIENCE IN CHINA Ser. G Physics, Mechanics Astronomy第4期程玉民等:彈性動(dòng)力學(xué)的邊界無(wú)單元法439lim fn(xf(x)=∑akPk(x)(28)因此f(x)∈spon(p).這就證明了spn(p)是完備的,所以span(p)是一個(gè) Hilbert空間2.2帶權的正交函數族對于點(diǎn)集{x;}和權函數{w},若一組函數q1(x),q2(x,L,qn(x)滿(mǎn)足如下條件0 k(k,9)=∑呢92(x)(x)(k,j=1,2,L,m)則稱(chēng)φ1(x),q2(x),L,φn(x)是關(guān)于點(diǎn)集{x;}帶權{w}的正交函數族.當q(x),q2(x,L,φn(x)是多項式時(shí),就稱(chēng)q1(x),q2(x),L,pn(x)是關(guān)于點(diǎn)集{x;}帶權{w}的正交多項式2.3改進(jìn)的移動(dòng)最小二乘法方程(13)可寫(xiě)成(P1,P1)(P1,P2)L(1,Pmn)a(x)|(D1,4)(P2,P1)(P2P2)L(P2Pm)ax2(x)_(p2,a)(30)(Pm, P1)(Pm, P2)L (Pm, Pm)Lam(x)L(pm, u,)若{P(x),i=1,2,L,m,為 Hilbert空間span(p)上的關(guān)于點(diǎn)集{x}的帶權的正交基函數族,即Pi, Pj(≠j)則方程(30)可寫(xiě)成0(P2,P20p2,ur(32)這樣我們可以直接得到系數a(x),即.1TYH魏中國科學(xué)G輯物理學(xué)力學(xué)天文學(xué)第35卷寫(xiě)成矩陣形式a(x)=A(rb(x)u(34)其中(P1,P1)(P2,P2)AO(35)將(33)代入(1)式可得)=(x)=∑(x(36)其中Φ(x)為形函數Φ(x)=應1(x,2(x),,(x)=P(x)A(x)B(x)(37)這樣,系數a(x)可以簡(jiǎn)單、直接地得到,不需要求矩陣A(x)的逆,避免了求解病態(tài)或奇異的方程組,既提高了效率,又提高了精度.2.4離散點(diǎn)上帶權的正交基函數的構造利用 Schmidt正交化方法,帶權的正交基函數可構造如下P11,PPki=2.3k= (pk,Pk)也可表示為P1p239)P=(r-a1)P-1-bP1其中rP;-1,P2-1)(P1-1,Pb=(41)中國煤化工對一維問(wèn)題,r=x;對二維問(wèn)題,r=f(YHCNMHG+x2或SCIENCE IN CHINA Ser. G Physics, Mechanics Astronomy第4期程玉民等:彈性動(dòng)力學(xué)的邊界無(wú)單元法r=x+x2等;對三維問(wèn)題,r=f(x,x2,x,例如可取r=√x2+x2+雞或r=x+寺另外,也可利用 Schmidt正交化方法將類(lèi)似于(2)-(5)式的基函數正交化.例如對基函數q=(q1)=(1,x,,,x2,x2…)(42)正交化得到的正交基函數為P;=q;Pk(i=1,2,3,L)若選取(38)或(39)式作為基函數,則在在試函數階次相同的情況下,試函數中的待定常數的數目比原來(lái)要少.對線(xiàn)性基,原來(lái)的待定常數是3個(gè),現在是2個(gè);對二次基,原來(lái)的待定常數是6個(gè),現在是3個(gè).這樣,對任一場(chǎng)點(diǎn)來(lái)說(shuō),其緊支域(影響域中所含的最小節點(diǎn)數就大大減少了,進(jìn)而在整個(gè)求解域中所需選取的節點(diǎn)數目也可以大大減少.所以在精度相同的情況下,改進(jìn)的移動(dòng)最小二乘法形成的無(wú)網(wǎng)格方法比移動(dòng)最小二乘法形成的無(wú)網(wǎng)格方法的節點(diǎn)數目要少得多3彈性動(dòng)力學(xué)的邊界無(wú)單元法對彈性動(dòng)力學(xué)的基本方程進(jìn)行 Fourier變換,得到 Fourier變換域中的基本方程.在 Fourier變換域中應用邊界無(wú)單元法進(jìn)行求解,然后應用數值 Fourier本征反變換即得時(shí)間域中的解彈性動(dòng)力學(xué)的基本方程:運動(dòng)微分方程(C2-C2)1(x,D)+Clu1(x,1)+f(x,)=(x,1)(x∈92),(44)其中+2pC(45)分別是彈性體內膨脹波和畸變波的傳播速度,λ和μ是Lame常數,ρ是體密度物理方程o(x,)=p(C2-2C2)xmm(x,16+C2(n1(x,1)+m1(x,1)(x∈9),(46)其中δ為 Kronecker應力和位移滿(mǎn)足如下的邊界條件t, (r, t)=o(x,)中國煤化工(47)l(x,1)=q1(x,1)CNMHGhina co442中國科學(xué)G輯物理學(xué)力學(xué)天文學(xué)第35卷和初始條件(x∈),其中r和r分別表示已知面力和位移的邊界,r∪I"=,⌒r"=Φ;n(x)為r上點(diǎn)x處的外法線(xiàn)的方向余弦;和分別為已知的位移和面力分量;u0和io分別為彈性體的初始位移和初始速度場(chǎng)將 Fourier變換f(x,o)=FLf(x,D1=6-f(r, t) - dr作用于時(shí)間域中的彈性動(dòng)力學(xué)基本方程,并令l(x,O)=F[L1(x,1)](r, o)=Flo,(r, t)t;(x,0)=F[t1(x,t)]f (x, o)=FIf, (r, t)],Pi(FLP; (x, t)l,即得 Fourier變換域中的彈性動(dòng)力學(xué)基本方程:運動(dòng)微分方程(C2-C2x0(x,0)+C1(x,0)+f(x,0)=02(x,0)(x∈9),(51)邊界條件1(x,0)=G(x,0)m(x)=(x,0)(x∈r(52)u(x,0)=q;(x,O)物理方程0(x,0)=p[(C1-2C2)mm(x,O)6+C2(1(x,)+:(x,0)(x∈9)(53)為方便起見(jiàn),我們假設初始條件和體力均為零.這樣,由加權殘數法可得彈性動(dòng)力學(xué)問(wèn)題在 Fourier變換域中的邊界積分方程lCn(q)x(q)=U(p,9)(p)dI-」T(P,q)(p)dr其中q為邊界點(diǎn)源點(diǎn),P為場(chǎng)點(diǎn),Cn(q)為自由項,U和T為 Fourier變換域中彈性動(dòng)力學(xué)方程的基本解.將邊界r離散為N個(gè)邊界子域,這些邊于,亞數無(wú)關(guān),僅供積分時(shí)使用,則邊界積分方程(54)變?yōu)門(mén)HE中國煤化工CNMHGSCIENCE IN CHINA Ser. G Physics, Mechanics Astronomy第4期程玉民等:彈性動(dòng)力學(xué)的邊界無(wú)單元法C(o(=∑J1(p(p-∑∫Tk(p, q)u, (p)dr在各邊界子域中配置一定數量的點(diǎn),并分別選取相應的緊支域,這些緊支域可以相交,但其并集需覆蓋整個(gè)邊界.由改進(jìn)的移動(dòng)最小二乘法得插值公式∑Φ(x∑(x)其中u;l=u(r)L r=t;(x將插值公式(56和(57代入(55)式,則有C()(q)=∑∫(,9∑畫(huà)(p)qmR(P,41)∑(P)(q(60)其中q為節點(diǎn),n1個(gè)q的緊支域覆蓋p,所選用的權函數是Gaus函數(d≤1)0其中c是常數,d是緊支域大小的度量.通常對二維問(wèn)題的邊界區域,緊支域為一維區域;對三維問(wèn)題的邊界區域,緊支域可選規則的圓域或矩形域對離散化的邊界積分方程(60)進(jìn)行數值積分后,寫(xiě)成矩陣形式即為([C]+[H])U]=[G][T其中U]=[1,12,l13,21,21T中國煤化x(660將邊界條件代入后,對線(xiàn)性代數方程組(HYH中邊界節CNMHG44中國科學(xué)G輯物理學(xué)力學(xué)天文學(xué)第35卷點(diǎn)的位移和面力變換域中內點(diǎn)x的位移可由內點(diǎn)的離散的邊界積分方程∑∫rU(pq∑畫(huà)(p)∑nr(p9,(p()d(67)得到,然后由幾何方程和物理方程即可得到變換域中內點(diǎn)的應力.在求得變換域中的解后,利用數值 Fourier本征反變換即可得到時(shí)間域中的解.這就是彈性動(dòng)力學(xué)的邊界無(wú)單元法4數值 Fourier本征反變換ll與 Fourier變換(49)式對應的反變換為f(x,)=FI(,0)=6_f(x,o)e'o'do(68)由 Fourier本征變換2可知,(68)式可以表示為f(x,1)=∑7a1甲(),其中(7(x,0)9n「f(x.)n(o)doO)qn()為本征基函數9()=/~1(71)由(68)式可得數值 Fourier本征反變換的表達式為(2r2f(x,1)=∑"anh1(e2=∑anhn()e2(72)其中dh, (t)∑4(f(x,0)中國煤化工CNMHGSCIENCE IN CHINA Ser. G Physics, Mechanics Astronomy第4期程玉民等:彈性動(dòng)力學(xué)的邊界無(wú)單元法5彈性動(dòng)力學(xué)邊界無(wú)單元法的數值實(shí)現對上述導出的彈性動(dòng)力學(xué)邊界無(wú)單元法,其數值實(shí)現過(guò)程如下(i)在問(wèn)題所在域的邊界配置n個(gè)節點(diǎn)(ⅱi)確定各節點(diǎn)的緊支域rn的大小d,使得UIn=In(i)選取基函數和權函數;(iv)計算形函數(ⅴ)將邊界離散為N個(gè)積分子域In,得到離散的邊界積分方程(ⅵi)通過(guò)數值積分,并代入邊界條件,得到線(xiàn)性代數方程組;(ⅶi)求解此線(xiàn)性代數方程組,即得各邊界節點(diǎn)的位移和面力(ⅷⅲ)由改進(jìn)的移動(dòng)最小二乘法的插值公式得到其他邊界點(diǎn)的位移和面力(ⅸx)由內點(diǎn)的離散的邊界積分方程計算域內任意一點(diǎn)的位移(x)由幾何方程和物理方程計算域內任意一點(diǎn)的應力(xi)由數值 Fourier本征反變換求時(shí)間域中的解6數值算例以下采用本文提出的彈性動(dòng)力學(xué)的邊界無(wú)單元法對受突加荷載作用的中心圓孔板和中心裂紋板進(jìn)行計算61中心圓孔板長(cháng)為36cm、寬為20cm的矩形板,中央有一直徑是10cm的圓形孔洞,在平面應力情況下兩側當t=0時(shí)作用突加荷載p=75MNm2,如圖1所示.材料的彈性模量E=2.1×105MNm2, Poisson比=0.3,質(zhì)量密度p=0785kNm2x2中國煤化工圖1受突加荷載作用哨YHCNMHGhina. com446中國科學(xué)G輯物理學(xué)力學(xué)天文學(xué)第35卷我們用邊界無(wú)單元法對本例進(jìn)行了計算.由于對稱(chēng)性,我們僅考慮此矩形裂紋板的四分之一情形.邊界無(wú)單元法的邊界節點(diǎn)設置如圖2所示計算所得的點(diǎn)A的位移值與有限元法計算結果的比較如圖3所示.可以看出,本文討論的邊界無(wú)單元法的計算結果和圖2中心圓孔板的邊界節點(diǎn)配置有限元法的計算結果吻合較好.0.00400020.008邊界無(wú)單元法0.012有限元法0.0010.0020.0030.004圖3點(diǎn)A處的位移62中心裂紋板受突加荷載作用的矩形裂紋板如圖4所示.裂紋長(cháng)度為2a,a=0.24cm.材料的彈性模量E=2.0×103MNm2, Poisson比v=0.3,質(zhì)量密度p=0.005MN.s2/n我們用邊界無(wú)單元法對本例進(jìn)行了計算.由于對稱(chēng)性,我們僅考慮此矩形裂紋板的四分之一情形.邊界無(wú)單元法的邊界節點(diǎn)設置如圖5所示計算所得的正則動(dòng)態(tài)應力強度因子(K()/K1,K1為靜態(tài)應力強度因子)與邊界元法的計算結果的比較如圖6所示由圖6可以看出,邊界無(wú)單元法的計算結果與邊界元法的計算結果吻合得很好中國煤化工CNMHGSCIENCE IN CHINA Ser. G Physics, Mechanics Astronomy第4期程玉民等:彈性動(dòng)力學(xué)的邊界無(wú)單元法圖4中心裂紋板圖5中心裂紋板的邊界節點(diǎn)配置3.02.50.50.0邊界元單元法圖6正則應力強度因子7結論和討論本文討論了 Hilbert空間上的改進(jìn)的移動(dòng)最小二乘法.改進(jìn)的移動(dòng)最小二乘法不會(huì )形成奇異或病態(tài)的方程組,可提高計算精度和效率改進(jìn)的移動(dòng)最小二乘法中所含的待定常數較少,這樣只需較少的節點(diǎn)的變量的值就可通過(guò)逼近函數得到任意場(chǎng)點(diǎn)的變量的值,那么由改進(jìn)的移動(dòng)最小二乘法形成的無(wú)網(wǎng)格方法求解問(wèn)題時(shí)可大大減少中國煤化工本文建立了求解彈性動(dòng)力學(xué)問(wèn)題的邊界eHCNMHG公式,提hina co中國科學(xué)G輯物理學(xué)力學(xué)天文學(xué)第35卷出了具體的數值實(shí)現方法.數值算例表明,本文方法是有效的本文方法可推廣到其他邊界元法可以解決的問(wèn)題.參考文獻I Belytschko T, Krongauz Y, Organ D, et al. Meshless methods: An overview and recent developments.Computer Methods in Applied Mechanics and Engineering, 1996. 139: 3-472 Li S F, Liu W K. Meshfree and particle methods and their applications. Applied Mechanics Review2002,55:1-343 Mukherjee Y X, Mukherjee S The boundary node method for potential problems. International Journal foNumerical Methods in Engineering, 1997, 40: 797-8154 Kothnur V S, Mukherjee S, Mukherjee Y X. Two dimensional linear elasticity by the boundary nodemethod, International Journal of Solids and Structures, 1999: 36: 1129~1145 Chati M K, Mukherjee S, Mukherjee Y x. The boundary node method for three-dimensional linearelasticity, International Journal for Numerical Methods in Engineering, 1999. 46: 1163-11846 Chati M K, Mukherjee S. The boundary node method for three-dilal problems in potential theorInternational Journal for Numerical Methods in Engineering, 2000, 47: 1523-1547 Zhu T L, Zhang J D, Atluri S N, A meshless numerical method based on the local boundary integralequation (LBIE)to solve linear and non-linear boundary value problems. Engineering Analysis withBoundary Elements, 1999, 23: 375-3898 Atluri S N, Sladek J, Sladek V, et al. The local boundary integral equation(LBIE)and it's meshlessiplementation for linear elasticity. Computational Mechanics, 2000, 25: 180-199 Zhang J M, Yao Z H, Li H. A hybrid boundary node method. International Journal for NumericalMethods in Engineering, 2002, 53: 751-76310程玉民,陳美娟.彈性力學(xué)的一種邊界無(wú)單元法.力學(xué)學(xué)報,2003,35(2):181-1861嵇醒,臧躍龍,程玉民.邊界元法進(jìn)展及通用程序.上海:同濟大學(xué)出版社,199712陳紹汀,韓海潮.傳立葉本征變換FI.應用力學(xué)學(xué)報,1987,4(1):33~38中國煤化工CNMHGSCIENCE IN CHINA Ser. G Physics, Mechanics Astronomy