軸對稱結(jié)構(gòu)動力彈塑性分析的無網(wǎng)格自然鄰接點(diǎn)Petrov-Galerkin法

陳莘莘， 肖樹聰， 周書濤

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

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

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

1 控制方程

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

(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)變增量；Dep為彈塑性矩陣，且有

(6)

(7)

Dep=De-Dp

(8)

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

(9)

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

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

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

對軸對稱面Ω上的任一節(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)位置的集合。

為了對某點(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,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ù)可寫為

(14)

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

2.2 平衡離散方程

(15)

(16)

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

利用格林公式對式(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)

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

(22)

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

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

(23)

其中，

(24a)

(24b)

(24c)

其中，

(25)

3 時(shí)間積分方案

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

(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)，則可得到如下的靜力等效問題

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

(30)

其中，

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

(31)

(32)

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

4 數(shù)值算例

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

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

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