陳小翠，杜成斌，江守燕

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

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

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

1 Reduced XFEM的位移模式簡介

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

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

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

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

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

采用文獻(xiàn)［9］中裂尖改進(jìn)函數(shù)，裂尖改進(jìn)節(jié)點(diǎn)的自由度個(gè)數(shù)從8個(gè)減少為4個(gè)，并且該裂尖改進(jìn)函數(shù)仍保留了裂紋尖端區(qū)應(yīng)力奇異性項(xiàng))和位移不連續(xù)性(sin(*)函數(shù)項(xiàng))。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*——考察點(diǎn)P在裂紋面上的投影點(diǎn)坐標(biāo);xk——第k個(gè)裂縫尖端的坐標(biāo);n——裂紋面的單位外法向;t——第k個(gè)裂紋尖端處的單位切向矢量;sign(x)——符號函數(shù)，x＞0時(shí)sign(x)=1，x=0時(shí)sign(x)=0，x＜0時(shí)sign(x)= －1。

2.2 改進(jìn)單元類型的判斷

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

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

3 應(yīng)力強(qiáng)度因子的計(jì)算

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

其中

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

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

4 數(shù)值算例

4.1 改進(jìn)的XF EM計(jì)算KⅠ

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

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

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

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

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

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

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

本算例分析有限尺寸拉伸板兩同側(cè)裂紋之間的影響，計(jì)算的拉伸板寬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)力強(qiáng)度因子的影響，保持C1長度a1=0.5 m不變，分別改變C2的長度a2及2條裂紋之間的距離d，用改進(jìn)的XFEM來計(jì)算C1相應(yīng)的應(yīng)力強(qiáng)度因子值。

計(jì)算結(jié)果如圖4和圖5所示，C1單裂紋情況下KⅠ由文獻(xiàn)［15］公式計(jì)算為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左右時(shí)，C2長度增量對C1的KⅠ影響最大;2條裂紋間的距離d對C1的KⅠ也有一定影響，如圖5所示，C1的應(yīng)力強(qiáng)度因子隨裂紋間距離增加呈“勺”形分布，且這種趨勢在a2/a1≥0.86時(shí)消失，呈遞增趨勢;從圖5左右兩端差異可知，右端無限向單裂紋狀態(tài)下應(yīng)力強(qiáng)度因子接近，這主要是由于當(dāng)C2無限遠(yuǎn)離C1時(shí)，C1更接近單裂紋狀態(tài);左端裂紋越長，C1的應(yīng)力強(qiá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)力強(qiáng)度因子的影響，保持C1長度a1=0.4 m不變，分別改變C2的長度a2及2條裂紋之間距離d，再用改進(jìn)的XFEM來計(jì)算C1相應(yīng)的應(yīng)力強(qiáng)度因子值。

計(jì)算結(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é) 語

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

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

［1］李錄賢，王鐵軍.擴(kuò)展有限元法(XFEM)及其應(yīng)用［J］.力學(xué)進(jìn)展，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ì)材料的擴(kuò)展有限元數(shù)值模擬［J］.河海大學(xué)學(xué)報(bào):自然科學(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］江守燕，杜成斌.一種擴(kuò)展有限元斷裂分析的裂尖單元新型改進(jìn)函數(shù)［J］.力學(xué)學(xué)報(bào)，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.