軸對稱結構動力彈塑性分析的無網格自然鄰接點Petrov-Galerkin法

陳莘莘， 肖樹聰， 周書濤

(1.華東交通大學 土木工程國家實驗教學示范中心，南昌 330013；2. 北京強度環境研究所，北京 100076)

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

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

1 控制方程

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

(1)

(2)

(3a)

(3b)

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

(4a)

(4b)

材料進入塑性狀態后，繼續加載時的應力應變關系為

dσ=Depdε

(5)

式中：dσ和dε分別為應力增量和應變增量；Dep為彈塑性矩陣，且有

(6)

(7)

Dep=De-Dp

(8)

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

(9)

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

2 動力彈塑性問題的離散化實現

2.1 自然鄰接點插值

對軸對稱面Ω上的任一節點x1，其一階Voronoi結構TI可定義為

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

(10)

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

為了對某點進行插值，還需要引入插值點的二階Voronoi結構TIJ,其數學表達式為

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

(11)

實際上，TIJ是以節點xI為最近點，節點xJ為第二近點的空間點位置的集合。圖1所示為平面中7個節點的Voronoi結構以及插值點x的二階Voronoi結構。

圖1 點x的一階和二階Voronoi結構Fig.1 First-order and second-order Voronoi cell about x

設AI(x)和A(x)分別為插值點x的二階Voronoi結構TxI和一階Voronoi結構Tx的面積，則自然鄰接點形函數及其導數可以表達為

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

(12)

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

(13)

定義了各節點的插值函數后，點x的位移函數可寫為

(14)

式中：uI(I=1,2,…,n)為點x周圍n個自然鄰接點的節點位移；φI(x)為對應節點的形函數。

2.2 平衡離散方程

(15)

(16)

式中，vI為加權函數。

利用格林公式對式(15)中積分的前兩項進行分部積分，可得

(17)

(18)

類似地,由式(16)可得

(19)

為了便于數值計算，聯立式(18)和式(19)可得

(20)

其中，

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

(21a)

(21b)

(21c)

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

(22)

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

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

(23)

其中，

(24a)

(24b)

(24c)

其中，

(25)

3 時間積分方案

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

(26a)

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

(26b)

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

(27)

(28a)

(28b)

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

(29)

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

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

(30)

其中，

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

(31)

(32)

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

4 數值算例

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

4.1 受突加內壓的厚壁圓筒

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

圖3 厚壁圓筒節點布置Fig.3 Nodal distribution for the thick-walled cylinder

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

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

4.2 受突加內壓的厚壁球殼

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

圖6 厚壁球殼節點布置Fig. 6 Nodal distribution for the thick-walled spherical shell

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

5 結 論

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