基于XFEM的主次裂紋間應(yīng)力強度因子相互作用

陳小翠，杜成斌，江守燕

(河海大學(xué)力學(xué)與材料學(xué)院，江蘇南京 210098)

盡管有限元法已存在數(shù)十年，但是有限元法仍不能有效地模擬工程中的裂紋及裂紋生長問題。對于裂紋問題，有限元法需布置高密度網(wǎng)格，并且在裂紋生長計算時有限元網(wǎng)格需要重剖分［1］。擴展有限元法(extended finite element method，XFEM)就是基于有限元處理不連續(xù)問題的弊端而提出的一種新的數(shù)值計算方法［2］。自1999年美國西北大學(xué)Belytschko教授及其研究組［3］提出擴展有限元法后，因XFEM在計算不連續(xù)問題時的優(yōu)越性，十多年內(nèi)得到迅速發(fā)展:Belytschko等［4］在XFEM中引入新的開裂準(zhǔn)則，來判斷裂紋生長路徑和速度;Song等［5］增強了含裂紋單元的描述;Fries［6］、Cheng等［7］在常規(guī) XFEM 基礎(chǔ)上提出了Corrected-XFEM，以提高數(shù)值求解精度;Liu等［8］通過譜單元與XFEM結(jié)合，改善動態(tài)裂紋擴展的數(shù)值擾動問題。在發(fā)展過程中，XFEM的理論基礎(chǔ)在不斷更新進步，以提高數(shù)值求解的精度。

筆者采用文獻［9］提出的新型裂尖改進函數(shù)，計算分析次裂紋的位置及長度對主裂紋應(yīng)力強度因子的影響。改進的裂尖改進函數(shù)在保留了裂紋尖端場的應(yīng)力奇異性和裂紋上、下表面位移不連續(xù)性的基礎(chǔ)上，減少了裂尖改進單元的附加自由度。

1 Reduced XFEM的位移模式簡介

式中:x——考察點的坐標(biāo);I——求解域中所有結(jié)點的集合;——被裂紋完全貫穿的單元結(jié)點集合;——裂尖改進結(jié)點集合;Ni(x)——結(jié)點i處的常規(guī)有限元插值形函數(shù);(x)——單位分解函數(shù)，其形式可以與Ni(x)相同，也可以不同，文中取二者相同;H(x)——Heaviside改進函數(shù);ui——常規(guī)有限元部分在結(jié)點i處的未知量;ai——與Heaviside改進相關(guān)的結(jié)點未知量;(x)——裂尖單元結(jié)點改進函數(shù)——與裂尖改進相關(guān)的結(jié)點未知量。

二維裂紋體的位移場用常規(guī)XFEM可表示為［10］

常規(guī)擴展有限元方法根據(jù)線彈性斷裂力學(xué)裂尖解析位移場的基本形式，裂尖單元結(jié)點的改進函數(shù)可用式(2)表示［10］:

式中:r、θ——裂尖極坐標(biāo)。

式(3)中π/4前的符號與θ保持一致。

采用文獻［9］中裂尖改進函數(shù)，裂尖改進節(jié)點的自由度個數(shù)從8個減少為4個，并且該裂尖改進函數(shù)仍保留了裂紋尖端區(qū)應(yīng)力奇異性項)和位移不連續(xù)性(sin(*)函數(shù)項)。XFEM位移模式可表示為

2 水平集法

2.1 裂紋的水平集法描述

裂紋用水平集函數(shù) ψ(x，t)和 φk(x，t)(k=1，2)描述［11］，裂紋面的位置通過 ψ(x，t)的零水平集函數(shù)ψ(x，t)=0來描述。用符號距離函數(shù)來構(gòu)造水平集函數(shù)，ψ(x，t)可表示為

波前水平集 φk(x，t)(k=1，2)與 ψ(x，t)正交，可表示為

式中:x*——考察點P在裂紋面上的投影點坐標(biāo);xk——第k個裂縫尖端的坐標(biāo);n——裂紋面的單位外法向;t——第k個裂紋尖端處的單位切向矢量;sign(x)——符號函數(shù)，x＞0時sign(x)=1，x=0時sign(x)=0，x＜0時sign(x)= －1。

2.2 改進單元類型的判斷

應(yīng)用XFEM計算裂紋問題時，共有裂尖改進單元和Heaviside改進單元這2類改進單元。文中分析時先計算出每個單元內(nèi)各節(jié)點到裂尖的值(φk(x，t))i及 ψi(x，t)，再比較得到各個單元的 φmax、φmin和 ψmax、ψmin。裂尖改進單元為單元節(jié)點的水平集函數(shù)滿足如下條件:

Heaviside改進單元為單元節(jié)點的水平集函數(shù)滿足如下條件:

3 應(yīng)力強度因子的計算

應(yīng)力強度因子是衡量裂紋尖端區(qū)應(yīng)力場強度的重要參數(shù)，也是斷裂力學(xué)中裂紋的失效判據(jù)。應(yīng)力強度因子計算方法主要是解析解法和數(shù)值解法。文中采用相互作用積分法計算應(yīng)力強度因子，裂紋尖端相互作用能量積分公式為［12］

其中

裂紋尖端相互作用積分與應(yīng)力強度因子的關(guān)系為

取XFEM計算的數(shù)值解作為真實場，線彈性斷裂力學(xué)裂尖解析場為輔助場(應(yīng)力場、位移場)，由式(11)可得真實場的Ⅰ型應(yīng)力強度因子:

4 數(shù)值算例

4.1 改進的XF EM計算KⅠ

本算例考察有限尺寸拉伸板的裂紋問題，如圖1所示。拉伸板寬b=1 m，高h(yuǎn)=2 m，彈性模量E=1 MPa，泊松比υ=0.3，軸向拉伸應(yīng)力σ=1 kPa。網(wǎng)格劃分為19×39(共741個)均勻網(wǎng)格，如圖1(c)所示。計算時分別取縫長 a為0.1 m、0.2 m、0.3 m、0.4 m、0.5 m。

邊緣裂紋問題應(yīng)力強度因子的解析解為［13］

中心裂紋問題應(yīng)力強度因子的解析解為［14］

圖1 有限尺寸拉伸板Fig.1 A finite stretched plate

采用改進的XFEM來計算圖1(a)和圖1(b)所示裂紋問題的KⅠ，并與解析解的計算結(jié)果進行比較，如圖2所示。由圖2可知，改進的XFEM方法的計算結(jié)果與解析解結(jié)果吻合較好。

圖2 解析解與XFEM數(shù)值解比較Fig.2 Comparison of analytical and XFEM numerical results

4.2 有限尺寸拉伸板同側(cè)裂紋問題

本算例分析有限尺寸拉伸板兩同側(cè)裂紋之間的影響，計算的拉伸板寬b=1 m，高h(yuǎn)=2 m，彈性模量E=1 MPa，泊松比υ=0.3，軸向拉伸應(yīng)力σ=1 kPa。板邊緣設(shè)2條裂紋C1和C2，位置如圖3所示。為研究C2的位置和長度對C1的應(yīng)力強度因子的影響，保持C1長度a1=0.5 m不變，分別改變C2的長度a2及2條裂紋之間的距離d，用改進的XFEM來計算C1相應(yīng)的應(yīng)力強度因子值。

計算結(jié)果如圖4和圖5所示，C1單裂紋情況下KⅠ由文獻［15］公式計算為3.545 kPa·m1/2。由圖4和圖5可知，由于C2存在，C1的KⅠ明顯比單裂紋狀態(tài)結(jié)果(3.545 kPa·m1/2)小;由圖4可知，C1長度a1不變，隨著C2的長度a2從0.1 m遞增到0.6 m，C1的KⅠ呈明顯遞減趨勢;當(dāng)a2/a1=0.9左右時，C2長度增量對C1的KⅠ影響最大;2條裂紋間的距離d對C1的KⅠ也有一定影響，如圖5所示，C1的應(yīng)力強度因子隨裂紋間距離增加呈“勺”形分布，且這種趨勢在a2/a1≥0.86時消失，呈遞增趨勢;從圖5左右兩端差異可知，右端無限向單裂紋狀態(tài)下應(yīng)力強度因子接近，這主要是由于當(dāng)C2無限遠(yuǎn)離C1時，C1更接近單裂紋狀態(tài);左端裂紋越長，C1的應(yīng)力強度因子相對單裂紋的減少量越大。

圖3 拉伸板Fig.3 A stretched plate

圖4 C2長度對C1的KⅠ影響Fig.4 Effect of length of C2 on KⅠ of C1

圖5 裂紋間距離對C1的KⅠ影響Fig.5 Effect of crack distance on KⅠ of C1

4.3 有限尺寸拉伸板異側(cè)裂紋問題

本算例分析有限尺寸拉伸板的異側(cè)2條裂紋之間的影響，仍采用4.2節(jié)的算例，并在板的左右設(shè)2條裂紋，位置如圖6所示。為研究C2的位置和長度對C1的應(yīng)力強度因子的影響，保持C1長度a1=0.4 m不變，分別改變C2的長度a2及2條裂紋之間距離d，再用改進的XFEM來計算C1相應(yīng)的應(yīng)力強度因子值。

計算結(jié)果如圖7所示，C1長度a1不變，隨著裂紋間距離d增大，C2長度變化對C1的KⅠ影響減小，而隨著C2的長度a2從0.1 m遞增到0.4 m，C1的KⅠ值呈明顯遞減趨勢。

圖6 有限尺寸拉伸板(2條裂紋)Fig.6 A finite stretched plate(two cracks)

5 結(jié) 語

用文獻［9］提出的改進XFEM進行次裂紋對主裂紋應(yīng)力強度因子的影響研究，通過斷裂力學(xué)中有限尺寸板的算例，得出該XFEM計算的應(yīng)力強度因子精度較高，與解析解結(jié)果吻合。在此基礎(chǔ)上研究次裂紋位置和長度對主裂紋應(yīng)力強度因子的影響，結(jié)果表明:對于平行裂紋問題，只要結(jié)構(gòu)有另一條裂紋存在，主裂紋的應(yīng)力強度因子會明顯比單裂紋狀態(tài)下的小;同側(cè)裂紋情況，主裂紋的應(yīng)力強度因子隨次裂紋長度增加有明顯的遞減趨勢，當(dāng)d/a1達到1.5時主裂紋的應(yīng)力強度因子變化趨于平緩，主裂紋的應(yīng)力強度因子隨裂紋間距離增大呈“勺”形，且這種趨勢在a2/a1≥0.86時消失，呈遞增趨勢;對于異側(cè)裂紋情況，主裂紋的應(yīng)力強度因子隨次裂紋長度的增加呈整體遞減趨勢，但隨著2條裂紋間距離的增大，次裂紋長度變化對主裂紋應(yīng)力強度因子的影響減小。

圖7 C2長度和位置對C1的影響Fig.7 Effects of length and position of C2 on C1

［1］李錄賢，王鐵軍.擴展有限元法(XFEM)及其應(yīng)用［J］.力學(xué)進展，2005(1):5-20.(LI Luxian，WANGTiejun.The extended finite element method and its applications［J］.Advances in Mechanics，2005(1):5-20.(in Chinese))

［2］應(yīng)宗權(quán)，杜成斌，程麗.含夾雜非均質(zhì)材料的擴展有限元數(shù)值模擬［J］.河海大學(xué)學(xué)報:自然科學(xué)版，2008，36(4):546-549.(YING Zongquan，DU Chengbin，CHENG Li.Application of extended finite element method in heterogeneous materials with inclusions［J］.Journal of Hohai University:Natural Sciences，2008，36(4):546-549.(in Chinese))

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

［4］ BELYTSCHKO T，CHEN Hao，XU Jingxiao，et al.Dynamic crack propagation based on loss of hyperbolicity and a new discontinuous enrichment［J］.Int J Numer Meth Engng，2003，58(12):1873-1905.

［5］SONGJ，BELYTSCHKOT.A method for dynamic crack and shear band propagation with phantom nodes［J］.Int JNumer Meth Engng，2006，67:863-893.

［6］FRIES T.A corrected XFEM approximation without problems in blending elements［J］.International Journal for Numerical Methods in Engineering，2008，75(5):503-532.

［7］CHENGK W，F(xiàn)RIEST.Higher-order XFEM for curved strong and weak discontinuities［J］.International Journal for Numerical Methods in Engineering，2010，82(5):564-590.

［8］LIU Z L，MENOUILLARD T，BELYTSCHKO T.An XFEM/Spectral element method for dynamic crack propagation［J］.Int J of Fracture，2011，169(2):183-198.

［9］江守燕，杜成斌.一種擴展有限元斷裂分析的裂尖單元新型改進函數(shù)［J］.力學(xué)學(xué)報，2013，45(1):134-138.(JIANG Shouyan，DU Chengbin.A novel enriched function of elements containing crack tip for fracture analysis in the framework of extended finite element methods［J］.Chinese Journal of Theoretical and Applied Mechanics，2013，45(1):134-138.(in Chinese))

［10］NICOLAS M，JOHN D，BELYTSCHKO T.A finite element method for crack growth without remeshing［J］.International Journal for Numerical Methods in Engineering，1999，46(1):131-150.

［11］DUFLOT M.A study of the representation of cracks with level sets［J］.2007，70(11):1261-1302.

［12］TOSHIO N，YOUHEI O，SHUICHI T.Stress intensity factor analysis of interface cracks using XFEM［J］.International Journal for Numerical Methods in Engineering，2003，56:1151-1173.

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

［14］ ZAMANI A，GRACIE R，ESLAMI M R.Cohesive and non-cohesive fracture by higher-order enrichment of XFEM［J］.International Journal for Numerical Methods in Engineering，2012，90:452-483.

［15］MOHAMMADI S.Extended finite element method for fracture analysis of structures［M］.Oxford，UK:Blackwell Publishing Ltd，2008.