二維粘彈性不可壓分析在擴(kuò)展有限元中的實現(xiàn)①

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

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

0 引言

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

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

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

1 擴(kuò)展有限元與粘彈性不可壓增量有限元

1.1 擴(kuò)展有限元的位移模式

擴(kuò)展有限元法分析斷裂問題時，先不考慮裂紋面的位置，直接劃分網(wǎng)格，然后應(yīng)用單位分解法的思想，對裂紋周圍的節(jié)點自由度進(jìn)行加強(qiáng)。對于平面裂紋，裂紋面的近似位移為

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

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

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

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

式中 r、θ為裂尖局部極坐標(biāo)。

圖1 網(wǎng)格中的任意裂紋Fig.1 An arbitrary crack placed on a mesh

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

1.2 粘彈性不可壓有限元法

Herrmann通過引入平均應(yīng)力函數(shù)，將經(jīng)典彈性本構(gòu)方程修改為方程(4)和方程(5)，解決了泊松比等于0.5時方程出現(xiàn)奇異的問題。

利用上述本構(gòu)關(guān)系建立有限元列式時，因平均應(yīng)力函數(shù)R的引入，單元內(nèi)增加了一個單自由度的節(jié)點［15］。對于平面問題，單元的廣義節(jié)點位移列陣U和單元位移列陣u分別為

此時，形函數(shù)矩陣為

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

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

其中

平均應(yīng)力函數(shù)R的引入，使單元內(nèi)增加了一個單自由度的節(jié)點，可視為對傳統(tǒng)有限單元進(jìn)行的增強(qiáng)處理，但這不同于上文介紹的擴(kuò)展有限元法。相比于傳統(tǒng)有限元法，擴(kuò)展有限元法只是增加了部分節(jié)點的自由度，并未增加節(jié)點的數(shù)目。上述不可壓粘彈性有限元法則相反，不改變原有節(jié)點的自由度，而在全體單元內(nèi)增加一個節(jié)點。前者是局部增強(qiáng)，實現(xiàn)便捷地劃分網(wǎng)格，后者是全局增強(qiáng)，實現(xiàn)不可壓材料的準(zhǔn)確分析。

2 擴(kuò)展有限元中的粘彈性不可壓法

2.1 有限元列式與相關(guān)矩陣

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

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

考慮到引入的平均應(yīng)力函數(shù)R，單元廣義節(jié)點位移列陣應(yīng)為

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

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

其中

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

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

則，單元內(nèi)任一點的應(yīng)變?yōu)?/p>

可見，因節(jié)點及節(jié)點自由度數(shù)的增加，形函數(shù)矩陣N和應(yīng)變矩陣B的規(guī)格相比于傳統(tǒng)有限元中的規(guī)格有所增加。為便于下文推導(dǎo)，依據(jù)B矩陣的特點，對彈性矩陣D進(jìn)行分塊處理，即

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

其中

單元廣義節(jié)點載荷列陣為

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

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

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

其他含有增強(qiáng)型節(jié)點的單元處理方法相似，關(guān)鍵在于確定廣義節(jié)點列陣U、形函數(shù)矩陣N及應(yīng)變矩陣B，下文以裂紋穿透的15單元為例補(bǔ)充說明。

15號單元為裂面增強(qiáng)型單元，根據(jù)式(1)，節(jié)點k的位移向量為

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

在求解應(yīng)變矩陣時需要對形函數(shù)求偏導(dǎo)，這里直接給出結(jié)果［18］:

其中，φα，x、φα，y通過坐標(biāo)變換，由裂尖局部坐標(biāo)系下的偏導(dǎo)獲得，即

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

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

2.2 應(yīng)力強(qiáng)度因子的求解

在擴(kuò)展有限元中直接求解應(yīng)力強(qiáng)度因子K時，需要借助于相互作用積分［4］:

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

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

3 算例

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

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

相關(guān)參數(shù)為:σ =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級數(shù)形式。對此，由四參量的Burgers模型求得級數(shù)形式的松弛模量［19］:

此問題的解析解為

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

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

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

圖4 蠕變載荷下的應(yīng)力強(qiáng)度因子Fig.4 Variation of SIF under creep loading

4 結(jié)論

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

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

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

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

［2］郭歷倫，陳忠富，羅景潤，等.擴(kuò)展有限元方法及應(yīng)用綜述［J］.力學(xué)季刊，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］方修君，金峰，王進(jìn)挺.基于擴(kuò)展有限元法的Koyna重力壩地震開裂過程模擬［J］.清華大學(xué)學(xué)報，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］莊茁，柳占成，成斌斌，等.擴(kuò)展有限單元法［M］.北京:清華大學(xué)出版社，2012.

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

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

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

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

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

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

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