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