陳莘莘， 肖樹(shù)聰， 周書(shū)濤

(1.華東交通大學(xué) 土木工程國(guó)家實(shí)驗(yàn)教學(xué)示范中心，南昌 330013；2. 北京強(qiáng)度環(huán)境研究所，北京 100076)

作為近年來(lái)興起的一種新的數(shù)值方法，無(wú)網(wǎng)格法[1-3]可以克服有限元等傳統(tǒng)數(shù)值方法對(duì)網(wǎng)格的依賴性，在很大程度上緩解了網(wǎng)格扭曲導(dǎo)致的數(shù)值困難。目前，一系列的無(wú)網(wǎng)格法相繼被提出，如：無(wú)單元Galerkin法[4-5]、無(wú)網(wǎng)格局部Petrov-Galerkin法[6- 7]、邊界無(wú)單元法[8-9]、自然單元法[10-11]和雜交邊界點(diǎn)法[12-13]等。與大多數(shù)無(wú)網(wǎng)格方法不同，無(wú)網(wǎng)格局部Petrov-Galerkin(meshless local Petrov-Galerkin, MLPG)法是基于微分方程的局部弱形式，進(jìn)行積分計(jì)算時(shí)只需要布置局部的積分單元，被譽(yù)為是一種真正的無(wú)網(wǎng)格法。然而，基于移動(dòng)最小二乘近似的MLPG方法不能直接施加本質(zhì)邊界條件。為了克服這一困難，Cai等[14-15]提出的無(wú)網(wǎng)格自然鄰接點(diǎn)Petrov-Galerkin法不僅允許權(quán)函數(shù)和試函數(shù)取自不同的函數(shù)空間，而且以節(jié)點(diǎn)的真實(shí)變量作為未知函數(shù)，從而可以直接施加本質(zhì)邊界條件。在這種無(wú)網(wǎng)格法中，任意節(jié)點(diǎn)的子域可由以該節(jié)點(diǎn)為公共頂點(diǎn)的Delaunay三角形構(gòu)成的多邊形區(qū)域，且積分可轉(zhuǎn)化為構(gòu)成該子域的Delaunay三角形區(qū)域上的積分之和。此外，線性有限元形函數(shù)作為權(quán)函數(shù)可以減少被積函數(shù)的階次，因而具有計(jì)算高效的特點(diǎn)。鑒于這些優(yōu)良特性，近年來(lái)無(wú)網(wǎng)格自然鄰接點(diǎn)Petrov-Galerkin法的在很多領(lǐng)域均得到了廣泛的應(yīng)用[16-19]。

軸對(duì)稱問(wèn)題在工程中非常常見(jiàn)，其幾何形狀、約束條件及作用的載荷都對(duì)稱于某一固定軸。 若利用其軸對(duì)稱特點(diǎn)，可將其轉(zhuǎn)化為二維問(wèn)題求解，進(jìn)而減少未知量。近年來(lái)，陳莘莘等[19-20]都致力于無(wú)網(wǎng)格法求解復(fù)雜軸對(duì)稱問(wèn)題的研究，并相繼取得了一些進(jìn)展。其中，陳莘莘等采用無(wú)網(wǎng)格自然鄰接點(diǎn)Petrov-Galerkin法求解了軸對(duì)稱彈性動(dòng)力學(xué)問(wèn)題。然而，在結(jié)構(gòu)的動(dòng)力響應(yīng)過(guò)程中，通?？偸羌扔袕椥宰冃?，又有塑性變形，而且這兩種變形以及它們之間的分界面都隨時(shí)間變化[21]。因此，本文在采用自然鄰接點(diǎn)插值構(gòu)造徑向位移和軸向位移近似函數(shù)的基礎(chǔ)上，采用局部Petrov-Galerkin法在多邊形子域上推導(dǎo)軸對(duì)稱結(jié)構(gòu)動(dòng)力彈塑性分析的離散方程，并采用預(yù)校正形式的Newmark法進(jìn)行迭代求解。最后，數(shù)值算例驗(yàn)證了本文建立的軸對(duì)稱結(jié)構(gòu)動(dòng)力彈塑性分析方法的有效性。

1 控制方程

對(duì)于軸對(duì)稱結(jié)構(gòu)，通常采用圓柱坐標(biāo)系(r,θ,z)，以對(duì)稱軸作為z軸，所有應(yīng)力、應(yīng)變和位移都與θ無(wú)關(guān)。任一點(diǎn)的位移只有兩個(gè)方向的分量，即沿r方向的徑向位移ur和沿z方向的軸向位移uz。在軸對(duì)稱面Ω上，Γu和Γt分別為問(wèn)題域的本質(zhì)邊界和自然邊界。如果不考慮阻尼，軸對(duì)稱結(jié)構(gòu)動(dòng)力彈塑性問(wèn)題的平衡方程可寫(xiě)為

(1)

(2)

(3a)

(3b)

ui(x,t)|t=0=ui(x,0)

(4a)

(4b)

材料進(jìn)入塑性狀態(tài)后，繼續(xù)加載時(shí)的應(yīng)力應(yīng)變關(guān)系為

dσ=Depdε

(5)

式中：dσ和dε分別為應(yīng)力增量和應(yīng)變?cè)隽?；Dep為彈塑性矩陣，且有

(6)

(7)

Dep=De-Dp

(8)

式中，De和Dp分別為彈性矩陣和塑性矩陣，且有

(9)

式中：σs為屈服應(yīng)力；G和Ep分別為剪切模量和塑性模量；sr,sz和sθ為應(yīng)力偏量。

2 動(dòng)力彈塑性問(wèn)題的離散化實(shí)現(xiàn)

2.1 自然鄰接點(diǎn)插值

對(duì)軸對(duì)稱面Ω上的任一節(jié)點(diǎn)x1，其一階Voronoi結(jié)構(gòu)TI可定義為

TI={x∈R2∶d(x,xI)

(10)

式中，d(x,x1)為點(diǎn)x和xI之間的距離。顯然，TI為以節(jié)點(diǎn)xI為最近離散點(diǎn)的空間點(diǎn)位置的集合。

為了對(duì)某點(diǎn)進(jìn)行插值，還需要引入插值點(diǎn)的二階Voronoi結(jié)構(gòu)TIJ,其數(shù)學(xué)表達(dá)式為

TIJ={x∈R2∶d(x,xI)d(x,xK),?K≠I(mǎi),J}

(11)

實(shí)際上，TIJ是以節(jié)點(diǎn)xI為最近點(diǎn)，節(jié)點(diǎn)xJ為第二近點(diǎn)的空間點(diǎn)位置的集合。圖1所示為平面中7個(gè)節(jié)點(diǎn)的Voronoi結(jié)構(gòu)以及插值點(diǎn)x的二階Voronoi結(jié)構(gòu)。

圖1 點(diǎn)x的一階和二階Voronoi結(jié)構(gòu)Fig.1 First-order and second-order Voronoi cell about x

設(shè)AI(x)和A(x)分別為插值點(diǎn)x的二階Voronoi結(jié)構(gòu)TxI和一階Voronoi結(jié)構(gòu)Tx的面積，則自然鄰接點(diǎn)形函數(shù)及其導(dǎo)數(shù)可以表達(dá)為

φI(x)=AI(x)/A(x)

(12)

φI,j(x)=(AI,j(x)-φI(x)AI,j(x))/A(x)

(13)

定義了各節(jié)點(diǎn)的插值函數(shù)后，點(diǎn)x的位移函數(shù)可寫(xiě)為

(14)

式中：uI(I=1,2,…,n)為點(diǎn)x周?chē)鷑個(gè)自然鄰接點(diǎn)的節(jié)點(diǎn)位移；φI(x)為對(duì)應(yīng)節(jié)點(diǎn)的形函數(shù)。

2.2 平衡離散方程

(15)

(16)

式中，vI為加權(quán)函數(shù)。

利用格林公式對(duì)式(15)中積分的前兩項(xiàng)進(jìn)行分部積分，可得

(17)

(18)

類似地,由式(16)可得

(19)

為了便于數(shù)值計(jì)算，聯(lián)立式(18)和式(19)可得

(20)

其中，

u=[ur,uz]T,σ=[σr,σz,σθ,τrz]T

(21a)

(21b)

(21c)

由于只對(duì)空間域進(jìn)行離散，軸對(duì)稱面Ω內(nèi)的位移u(x,t)可由式(14)寫(xiě)為

(22)

圖2 局部多邊形子域Fig.2 The local polygonal sub-domains

將式(22)代入式(20)，可得軸對(duì)稱結(jié)構(gòu)動(dòng)力彈塑性分析的離散控制方程為

(23)

其中，

(24a)

(24b)

(24c)

其中，

(25)

3 時(shí)間積分方案

求解運(yùn)動(dòng)方程式(23),本文采用預(yù)校正形式的Newmark方法。采用Newton-Raphson迭代，時(shí)間t+Δt的運(yùn)動(dòng)方程可由式(23)改寫(xiě)為

(26a)

t+Δtu(k+1)=t+Δtu(k)+Δu(k)

(26b)

式中,KT為切線剛度矩陣。

(27)

(28a)

(28b)

式中，α和β為按積分精度和穩(wěn)定性要求決定的參數(shù)。聯(lián)立式(26b)和式(28b)，得到

(29)

將式(29)代入式(26a)，則可得到如下的靜力等效問(wèn)題

K*Δu(k)=t+Δtf(k)

(30)

其中，

K*=KT+M/(α(Δt)2)

(31)

(32)

靜力等效問(wèn)題式(30)的具體迭代過(guò)程可參見(jiàn)文獻(xiàn)[22]。

4 數(shù)值算例

為了驗(yàn)證所提數(shù)值方法的有效性，本文對(duì)軸對(duì)稱結(jié)構(gòu)動(dòng)力彈塑性分析的一些典型算例進(jìn)行計(jì)算和對(duì)比分析。如果不作特別說(shuō)明，材料均取為服從Von Mises屈服條件的理想彈塑性材料。

4.1 受突加內(nèi)壓的厚壁圓筒

無(wú)限長(zhǎng)受突加內(nèi)壓P=185 Pa的厚壁圓筒，內(nèi)半徑a=100 m，外半徑b=200 m，彈性模量E=2.1×105Pa，泊松比v=0.3，屈服應(yīng)力σs=355 Pa，質(zhì)量密度ρ=7.85×10-6kg/m3。截取一段h=100 m作為計(jì)算模型，其節(jié)點(diǎn)布置方案如圖3所示，且時(shí)間步長(zhǎng)取為Δt=1.0×10-5s。圖4給出了厚壁圓筒內(nèi)表面A點(diǎn)的徑向位移隨時(shí)間的變化曲線，圖5給出了厚壁圓筒r=150 m處的彈塑性徑向應(yīng)力σr隨時(shí)間的變化曲線。從圖4和圖5可以看出本文數(shù)值解與有限元軟件Abaqus的計(jì)算結(jié)果吻合很好。

圖3 厚壁圓筒節(jié)點(diǎn)布置Fig.3 Nodal distribution for the thick-walled cylinder

圖4 厚壁圓筒內(nèi)表面的徑向位移Fig.4 Radial displacement history at the inner surface of the thick-wallled cylinder

圖5 r=150 m處的彈塑性徑向應(yīng)力σr Fig. 5 Elastoplastic radial stress σr when r=150 m

4.2 受突加內(nèi)壓的厚壁球殼

考慮一受突加均布內(nèi)壓P=75 MPa的厚壁球殼，其外壁半徑R1=0.2 m，內(nèi)壁半徑R2=0.16 m。彈性模量E=2.1×1011Pa，泊松比v=0.3，屈服應(yīng)力σs=200 MPa，質(zhì)量密度ρ=7.85×103kg/m3。采用如圖6所示的節(jié)點(diǎn)布置方案，時(shí)間步長(zhǎng)取為Δt=3.0×10-7s。圖7給出了厚壁球殼內(nèi)表面A點(diǎn)的徑向位移隨時(shí)間的變化曲線。顯然，本文數(shù)值解與有限元軟件Abaqus的計(jì)算結(jié)果吻合很好，這進(jìn)一步驗(yàn)證了本文所提方法的有效性。

圖6 厚壁球殼節(jié)點(diǎn)布置Fig. 6 Nodal distribution for the thick-walled spherical shell

圖7 厚壁球殼內(nèi)表面的徑向位移Fig. 7 Radial displacement history at the inner surface of the thick-walled spherical shell

5 結(jié) 論

本文將無(wú)網(wǎng)格自然鄰接點(diǎn)Petrov-Galerkin法進(jìn)一歩推廣求解軸對(duì)稱結(jié)構(gòu)動(dòng)力彈塑性問(wèn)題。在采用自然鄰接點(diǎn)插值構(gòu)造徑向位移和軸向位移近似函數(shù)的基礎(chǔ)上，采用線性有限元形函數(shù)作為權(quán)函數(shù)詳細(xì)推導(dǎo)了軸對(duì)稱結(jié)構(gòu)動(dòng)力彈塑性分析的無(wú)網(wǎng)格法。這種方法不僅有效地避免了復(fù)雜的網(wǎng)格劃分和網(wǎng)格畸變的影響，而且能夠方便地通過(guò)加密節(jié)點(diǎn)提高計(jì)算精度，有利于發(fā)展相關(guān)的自適應(yīng)算法。另外，這種方法沒(méi)有任何人為參數(shù)的選擇問(wèn)題，可以直接準(zhǔn)確地施加本質(zhì)邊界條件，子域構(gòu)造方式簡(jiǎn)單，計(jì)算效率高。本文的數(shù)值算例結(jié)果表明，采用無(wú)網(wǎng)格自然鄰接點(diǎn)Petrov-Galerkin法求解軸對(duì)稱結(jié)構(gòu)動(dòng)力彈塑性問(wèn)題是可行的，而且具有較高的計(jì)算精度。