二維粘彈性不可壓分析在擴展有限元中的實現①

顧志旭，鄭 堅，彭 威，殷軍輝

(軍械工程學院火炮工程系，石家莊 050003)

0 引言

裂紋作為影響固體火箭發動機藥柱結構完整性的重要因素，長期以來受到廣泛關注。傳統有限元法在處理斷裂問題時存在著固有缺陷，如裂紋面和裂尖必須位于單元邊和單元節點上，裂紋只能沿單元邊界擴展，擴展時需重新劃分網格等［1－2］。此外，為獲得較為精確的結果，往往需要對裂紋局部的網格進行加密處理，大大增加了計算成本。針對上述弊端，1999年，Belytschko和M?es基于單位分解思想對傳統有限元進行了重要的改進，提出了一種新的計算方法－擴展有限元法(XFEM)［3－4］。該方法是在常規有限元位移模式中，加入了能夠反映裂紋面不連續性的階躍函數和裂尖局部特性的漸進位移場函數，使得非連續體的幾何構型完全獨立于網格，無需對裂紋局部的網格加密處理，同時也避免了傳統有限元重復劃分網格的不便。應用擴展有限元法，Belytschko等［5］成功模擬了裂紋分叉和擴展，Sukumar等［6］模擬了含孔洞和夾雜的非均勻材料，Gracie等［7］首次實現了在細觀尺度上二維和三維固體材料中的位錯的有限元模擬。方修君等［8］提出了基于擴展有限元法的粘聚裂紋模型，并用于重力壩地震開裂過程。Zhuang Z等［9］基于連續體的殼單元建立了新的殼體擴展有限元格式。此外，在剪切帶演化［10］、多相流［11］、納米界面力學［12］、多尺度模擬［13］等方面也有廣泛的應用。經十多年的發展完善，擴展有限元法逐漸成為了一種處理非連續場、局部變形和斷裂等復雜力學問題的功能強大、極具應用前景的新方法［14］。

固體推進劑的泊松比接近0.5，是一種近似不可壓縮的粘彈性材料。采用基于Herrmann泛函建立的不可壓和近似不可壓粘彈性增量有限元法［15－16］分析時，可避免計算中出現體積“自鎖”的現象，且計算準確度較好。然而，該方法是在傳統有限元位移模式建立的，在分析斷裂問題時，仍然會遇到網格劃分不便的問題。

本文基于單位分解思想，將由Herrmann泛函建立的不可壓和近似不可壓粘彈性增量有限元法推廣至擴展有限元模式下，充分結合了前者計算精度高和后者網格劃分(前處理)簡潔的優勢。最后給出一個算例，結果表明了本文方法的有效性和便捷性。

1 擴展有限元與粘彈性不可壓增量有限元

1.1 擴展有限元的位移模式

擴展有限元法分析斷裂問題時，先不考慮裂紋面的位置，直接劃分網格，然后應用單位分解法的思想，對裂紋周圍的節點自由度進行加強。對于平面裂紋，裂紋面的近似位移為

式中 Sc為所有常規單元節點的集合;Sh為完全被裂紋貫穿的單元節點的集合(圖1中小方塊表示的節點);St為含裂尖單元節點的集合(圖1中小圓圈表示的節點);{u0i，v0i}T、{aui，avi}T和{buα，bvα}T分別為常規單元節點、貫穿單元節點和裂尖單元節點的位移;Nl(x)(l=i，j，k)為常規有限元中的形函數;H(f(x))為階躍函數;φα(x)為裂尖漸進位移函數。

階躍函數H(f(x))用來反映裂紋面位移的不連續性，定義為

式中 f(x)為符號距離函數構造的水平集函數。

裂尖漸進位移函數φα(x)用來反映裂尖區域的局部特性，常取裂紋位移解析解的各項，即

式中 r、θ為裂尖局部極坐標。

圖1 網格中的任意裂紋Fig.1 An arbitrary crack placed on a mesh

在擴展有限元中，節點自由度對應著局部近似空間的基函數系數［17］。St、Sh、Sc集中的節點在每個常規自由度方向分別增加4、1、0個自由度，則節點的廣義自由度分別為 10、4、2。

1.2 粘彈性不可壓有限元法

Herrmann通過引入平均應力函數，將經典彈性本構方程修改為方程(4)和方程(5)，解決了泊松比等于0.5時方程出現奇異的問題。

利用上述本構關系建立有限元列式時，因平均應力函數R的引入，單元內增加了一個單自由度的節點［15］。對于平面問題，單元的廣義節點位移列陣U和單元位移列陣u分別為

此時，形函數矩陣為

對于粘彈性材料，將松弛模量表示為Prony級數形式:

利用Herrmann變分原理，粘彈性對應原理及虛功原理，推導得出粘彈性不可壓增量有限元法的支配方程［15－16］:

其中

平均應力函數R的引入，使單元內增加了一個單自由度的節點，可視為對傳統有限單元進行的增強處理，但這不同于上文介紹的擴展有限元法。相比于傳統有限元法，擴展有限元法只是增加了部分節點的自由度，并未增加節點的數目。上述不可壓粘彈性有限元法則相反，不改變原有節點的自由度，而在全體單元內增加一個節點。前者是局部增強，實現便捷地劃分網格，后者是全局增強，實現不可壓材料的準確分析。

2 擴展有限元中的粘彈性不可壓法

2.1 有限元列式與相關矩陣

在擴展有限元模式下，基于本構方程(4)、(5)建立粘彈性不可壓法的支配方程的方法同文獻［15－16］。難點在于確定剛度矩陣K、形函數矩陣N、應變矩陣B及載荷列陣F等矩陣的具體形式。為便于說明，對圖1中的單元進行編號。圖中不含自由度增強節點的單元，如25號單元，相關矩陣見文獻［15］。本文只討論含增強型節點的單元，如9、12、15、16、21 號單元。

首先，分析裂尖所在的12號單元。依據式(1)，單元內節點k(k=i，j，m，n為單元節點局部編號，下同)的位移向量為

考慮到引入的平均應力函數R，單元廣義節點位移列陣應為

單元內任一點的廣義位移同式(7)。

形函數矩陣整體形式同式(8)，但因節點自由數的增加，結合式(1)為

其中

Nl=(1+ξlξ)(1+ηlη)/4 為常規矩形單元的形函數。

應變矩陣依據式(12)、式(13)及幾何方程得出:

則，單元內任一點的應變為

可見，因節點及節點自由度數的增加，形函數矩陣N和應變矩陣B的規格相比于傳統有限元中的規格有所增加。為便于下文推導，依據B矩陣的特點，對彈性矩陣D進行分塊處理，即

將式(16)代入式(10)中剛度矩陣公式，得

其中

單元廣義節點載荷列陣為

其中，Fk(k=i，j，m，n)為廣義等效節點載荷，由式(19)確定:

式中 f=［fu，fv，0］T、p=［pu，pv，0］T分別稱為廣義分布體力和廣義分布面力，0分量對應于廣義位移中的R分量［15］。

由式(17)、式(12)及式(18)求得裂尖12號單元的剛度矩陣K、廣義節點位移列陣U及廣義節點載荷列陣F后，利用式(10)可對此單元進行分析。

其他含有增強型節點的單元處理方法相似，關鍵在于確定廣義節點列陣U、形函數矩陣N及應變矩陣B，下文以裂紋穿透的15單元為例補充說明。

15號單元為裂面增強型單元，根據式(1)，節點k的位移向量為

單元廣義節點位移列陣同(12)式，形函數矩陣N與應變矩陣B的整體形式同式(13)、式(14)，但其中的分量不同，如下式所示:

在求解應變矩陣時需要對形函數求偏導，這里直接給出結果［18］:

其中，φα，x、φα，y通過坐標變換，由裂尖局部坐標系下的偏導獲得，即

式中 φ為裂尖切線與全局坐標系x軸之間的夾角。

φα，x'、φα，y'由以下分量確定:

2.2 應力強度因子的求解

在擴展有限元中直接求解應力強度因子K時，需要借助于相互作用積分［4］:

在裂尖區域，相互作用積分I(1，2)與真實變形場和附加變形場的應力強度因子K有如下關系:

式中 平面應力時E*=E(彈性模量)，平面應變時E*=E/(1－ν2)，ν為泊松比。

3 算例

為了驗證本文方法的合理性，以一受拉的帶單邊裂紋的矩形平板為例進行試算，見圖2。

圖2 含單邊裂紋的受拉板Fig.2 An edge-cracked plate under tension

相關參數為:σ =10 Pa，L/W=16/7，a/W=1/2，W=0.07 m。板的粘彈性由Burgers模型描述，見圖3。其中，E1=0.495 MPa，E2=4.854 MPa，η1=1 × 1015MPa·s，η2=1.55 ×103MPa·s。

圖3 Burgers模型Fig.3 Burgers model

前文建立有限元列式時，用到了松弛模量的Prony級數形式。對此，由四參量的Burgers模型求得級數形式的松弛模量［19］:

此問題的解析解為

式中 KI(t)為蠕變斷裂應力強度因子［20］為線彈性應力強度因子，由文獻［21］給出;fig(t)為能量釋放率的時間因子［20］。

不同泊松比(不同程度的近似不可壓性)下，裂紋的蠕變應力強度因子的計算結果見圖4。其中，圖4(a)給出了ν=0.495時的解析解、本文解以及常規有限元解。可知，蠕變開始時刻t=0時，三者分別為9.386 2、9.255 2、9.257 6 Pa·m1/2較吻合。隨著蠕變的進行，本文解和解析解趨近于15.916 8 Pa·m1/2，而常規有限元解趨于15.697 7 Pa·m1/2，比本文解誤差大。圖4(b)用解析法和本文法計算了 ν分別為0.400、0.450、0.500 時的蠕變應力強度因子。結果表明，在蠕變初期(0＜t＜120)，不同泊松比下應力強度因子相近。隨著時間的推移，泊松比的影響逐漸凸顯，最終蠕變應力強度因子趨于不同值。圖4(c)計算了ν分別為0.490和0.499時的蠕變應力強度因子，兩者變化規律與圖4(a)和(b)相同。從結果來看，泊松比相差1.84%，最終應力強度因子相差2.75%。可見，在此種工況下，粗略計算時泊松比精確到小數點后第二位即可。此外，圖4(a)～(c)中均表現出本文解的誤差，在蠕變起始階段和穩定階段均小于中間階段。這主要是由Burgers模型轉換得到的松弛模量中級數項數較少，粘彈性表征不足所致。

綜上所述，對于彈性材料泊松比對應力強度因子的影響較小(t=0 s時刻，蠕變未發生時，不同泊松比下所得應力強度因子近似為同一值，曲線起始于同一點)，而對于粘彈性材料，隨著蠕變的發生，不同泊松比的影響逐漸放大。因此，對粘彈性斷裂問題，如果僅關注短時響應，粗略計算時，可忽略泊松比變化的影響。如果關注長時響應，就必須考慮泊松比大小的影響。

圖4 蠕變載荷下的應力強度因子Fig.4 Variation of SIF under creep loading

4 結論

(1)嘗試將基于Herrmann泛函建立的不可壓和近似不可壓粘彈性增量有限元法，推廣至擴展有限元模式。這樣在有限元分析粘彈性材料的斷裂問題時，將不受泊松比和網格劃分的限制，既保證了計算的準確度，同時又降低了有限元前處理的難度。算例結果表明了該方法的有效性和便捷性。

(2)以四節點矩形單元為例建立了有限元列式，但該方法并不受單元類型的限制，適用于任意節點的二維單元。在選擇裂尖增強單元時，本文只選取了1層，而在實際應用時，為提高計算精度，往往選擇多層單元，文中方法仍可推廣。

(3)本文的工作只表明該方法可行，所建有限元列式僅能用于二維分析。如要對推進劑藥柱結構分析，鑒于其結構和載荷的復雜性，必須要建立三維的有限元列式，進行三維分析，這將是下一步的工作。

［1］茹忠亮，朱傳銳，張友良，等.斷裂問題的擴展有限元法研究［J］.巖土力學，2011，32(7):2171－2176.

［2］郭歷倫，陳忠富，羅景潤，等.擴展有限元方法及應用綜述［J］.力學季刊，2011，32(4):612－625.

［3］Belytschko T，Black T.Elastic crack growth in finite elements with minimal re－meshing［J］.Int.J.Numer.Meth.Engng.，1999，45(5):601－620.

［4］M?es N，Dolbow J，Belytschko T.A finite element method for crack growth without remeshing［J］.Int.J.Numer.Meth.Engng.，1999，46:131－150.

［5］Belytschko T，et al.Dynamic crack propagation based on loss of hyperbolicity and a new discontinuous enrichment［J］.Int.J.Numer.Meth.Engng.，2003，58:1873－1905.

［6］Sukumar N，et al.Modeling holes and inclusions by level sets in the extended finite－element method［J］.Computer Methods in Applied Mechanics and Engineering，2001，190(46－47):6183－6200.

［7］Gracie R，Oswald J，Belytschko T.On a new extended finite element method for dislocations:Core enrichment and nonlinear formulation［J］.Journal of the Mechanics and Physics of Solids，2008，56(1):200－214.

［8］方修君，金峰，王進挺.基于擴展有限元法的Koyna重力壩地震開裂過程模擬［J］.清華大學學報，2008，48(012):2065－2069.

［9］Zhuang Z，Cheng B B.A novel enriched CB shell element method for simulating arbitrary crack growth in pipes［J］.Science China Physics，Mechanics ＆ Astronomy，2011，54(8):1520－1531.

［10］Samaniego E，Belytschko T.Continuum－discontinuum modelling of shear bands［J］.Int.J.Numer.Meth.Engng.，2005，62(13):1857－1872.

［11］Chess J，Belytschko T.An extended finite element method for two－phase fluids:Flow simulation and modeling［J］.Journal of Applied Mechanics，2003，70(1):10－17.

［12］Farsad M，Vernerey F J，Park H S.An extended finite element/level set method to study surface effects on mechanical behavior and properties of nanomaterials［J］.Int.J.Numer.Meth.Engng.，2010，84(12):1466－1489.

［13］Belytschko T，Gracie R，Ventura G.A review of extended/generalized finite element methods for material modeling［R］.Modeling and Simulation in Materials Science and Engineering，2009，17(043001).

［14］莊茁，柳占成，成斌斌，等.擴展有限單元法［M］.北京:清華大學出版社，2012.

［15］張海聯，周建平.固體推進劑藥柱的近似不可壓縮粘彈性增量有限元法［J］.固體火箭技術，2001，24(2):36－40.

［16］田四朋，雷勇軍，李道奎，等.固體火箭發動機藥柱不可壓和近似不可壓三維分析［J］.固體火箭技術，2006，29(6):395－399.

［17］金峰，方修君.擴展有限元法及與其他數值方法的聯系［J］.工程力學，2008，25(1):1－16.

［18］Soheil Mohammadi.Extended Finite Element Method［M］.Blackwell Publishing Ltd.，2008.

［19］潘曉明，余俊，楊釗，等.一種將線性粘彈性微分型本構方程應用到ABAQUS的方法［J］.華僑大學學報，2010，31(5).

［20］張淳源.粘彈性斷裂力學［M］.武漢:華中理工大學出版社，1992.

［21］趙建生.斷裂力學及斷裂物理［M］.武漢:華中科技大學出版社，2003:126.