基于XFEM研究含顆粒夾雜材料的疲勞裂紋行為

束一秀，李亞智，尚海江，姜　薇

(西北工業(yè)大學 航空學院，西安　710072)

?

基于XFEM研究含顆粒夾雜材料的疲勞裂紋行為

束一秀，李亞智，尚海江，姜薇

(西北工業(yè)大學 航空學院，西安710072)

材料中的顆粒夾雜會改變附近的應力分布，對疲勞裂紋的擴展行為有顯著的影響。在XFEM方法的基礎上，借助一種考慮界面影響的交互積分方法以及ABAQUS的二次開發(fā)環(huán)境，模擬疲勞裂紋在非勻質材料中的擴展，研究夾雜顆粒對疲勞裂紋擴展的影響模式。裂紋和夾雜使用擴充函數(shù)結合水平集函數(shù)隱式模擬，使用最大周向拉應力準則判斷裂紋擴展方向。通過模擬疲勞裂紋經過對稱夾雜和非對稱夾雜的情況，分別研究夾雜對裂紋擴展速率以及擴展路徑的影響模式，探討了夾雜的材料性質對結果的影響，考慮了裂紋在特定情況下可能刺入夾雜的情況。將部分數(shù)值結果與更新網(wǎng)格方法以及以往的擴展有限元模擬結果進行了對比，體現(xiàn)了該文方法的優(yōu)點和模擬結果的合理性。

疲勞裂紋擴展；擴展有限元法；交互積分；應力強度因子；材料界面

0　引言

如果材料內存在材料性質與母材不同的夾雜，不管這些夾雜是偶然存在(材料加工過程中的雜質)或是人為添加(顆粒增強復合材料)，都會顯著地影響材料的力學性能，其中夾雜對裂紋擴展的影響是受關注的問題之一。由勻質材料形成的結構中，填充孔和空孔也可看作是夾雜問題的特殊形式。夾雜與裂紋之間的相互作用對于很多工程材料和結構的斷裂以及疲勞壽命等問題具有重要的意義。

解決此類問題有很多種方法。一些學者[1-3]推導了夾雜對裂紋作用的解析解，證明夾雜對其附近的裂紋尖端有明顯的“抑制”或“放大”作用。數(shù)值解法也被用來求解夾雜和裂紋的相互作用，比如有限元方法[4-5]和邊界元方法[6-8]。擴展有限元方法[9-12](XFEM)是一種解決區(qū)域內包含不連續(xù)問題的數(shù)值方案，與傳統(tǒng)有限元方法相比，其優(yōu)點是不依賴于不連續(xù)體的幾何形狀，只要針對不同的間斷類型選擇不同的擴充函數(shù)，即可模擬域內包含裂紋和夾雜的問題，避免了網(wǎng)格更新。于紅軍、王志勇[13-14]使用擴展有限元方法結合一種考慮材料界面的交互積分形式計算了夾雜對裂尖能量釋放率的影響，但他們沒有使用裂尖擴充函數(shù)，而是通過調整裂尖處的網(wǎng)格來提高精度，其優(yōu)點是求解不依賴于裂尖應力場的形式，缺點是需要對裂尖處的網(wǎng)格作特殊處理；Natarajan[15]使用擴展有限元法計算了裂尖前方的圓形夾雜對裂尖能量釋放率的影響，但沒有對受夾雜影響之后的裂紋擴展路徑進行分析。雖然很多學者[7-8，14]模擬了裂紋受夾雜影響下的擴展路徑分析，但沒有考慮裂紋刺入夾雜的情況，而實際上如果夾雜的剛度小于母材，或者夾雜的材料雖然剛度大于母材但差距不大，不能視為剛性夾雜，在多個夾雜綜合作用的情況下，有可能會出現(xiàn)裂紋刺入夾雜的情況。

本文基于擴展有限元方法研究母材中含顆粒夾雜的情況下疲勞裂紋的擴展行為。第一部分介紹了相關數(shù)值過程以及斷裂參數(shù)的計算方法，數(shù)值計算的實現(xiàn)基于ABAQUS軟件的二次開發(fā)；第二部分通過算例分析了規(guī)則分布的夾雜對Ⅰ型疲勞裂紋擴展壽命的影響并與更新網(wǎng)格法作對比；第三部分模擬了不同性質的夾雜對裂紋擴展路徑的影響，探討了裂紋經過夾雜區(qū)的擴展規(guī)律，分析了裂紋刺入夾雜時應力強度因子的變化規(guī)律。算例合理地反映了裂紋經過顆粒夾雜區(qū)域時的擴展行為。

1　數(shù)學模型

1.1擴展有限元離散形式

Belytschko[10]等提出了求解線彈性問題的擴展有限元方法，其核心思想是在標準位移場上增加附加位移場來模擬不連續(xù)問題。如果區(qū)域內包含裂紋和夾雜，則其位移函數(shù)可寫成式(1)的形式：

(1)

式中Ni為形狀函數(shù)；H為階躍函數(shù)，表征被裂紋貫穿單元的節(jié)點的擴充函數(shù)，在裂紋兩側分別取1和-1；φ為表征夾雜的擴充函數(shù)[12]，其形式見式(2)；Ψj(r，θ)為裂尖擴充函數(shù)，其形式見式(3)。

(2)

式中ζi為節(jié)點到夾雜邊界的符號距離函數(shù)。

(3)

單元剛度矩陣和節(jié)點力向量分別為式(4)和式(5)所示：

(4)

(5)

其中

(6)

式中B為幾何矩陣；D為彈性矩陣。

1.2裂紋和夾雜的水平集表征

擴展有限元方法中的裂紋和夾雜不依賴于網(wǎng)格，用水平集方法來表征任意形式的裂紋[16-17]及顆粒夾雜[12]的幾何形狀。

裂紋屬于非閉合曲線，需要使用2個水平集函數(shù)φ和ψ來確定其幾何形狀及位置(圖1)。其中，φ用來表征裂紋面，ψ用來輔助確定裂尖位置，則一條任意位置及形狀的裂紋可描述為φ=0，ψ≤0。

圖1　擴充節(jié)點和水平集函數(shù)的零等值面

顆粒夾雜邊界是封閉曲線，只需要用一個水平集函數(shù)ζ即可精確描述其幾何形狀，ζ=0即為夾雜的輪廓線。以圖1中的橢圓形夾雜為例，節(jié)點的水平集函數(shù)可用式(7)計算：

(7)

2　疲勞裂紋擴展

2.1交互積分

文獻[13]提出了一種考慮材料界面影響的交互積分形式：

(8)

現(xiàn)在考察如果裂紋與材料界面相交，交互積分是否依然滿足式(8)的形式。如圖2所示，積分區(qū)域被分割為A1、A2和A33個材料性質均勻的積分區(qū)域，定義3個區(qū)域的輪廓為Γ01、Γ02和Γ03，則交互積分形式可表示為式(9)的形式：

(9)

(10)

圖2　裂紋穿過材料界面時的積分域

分別對式(9)中的3個環(huán)路積分使用散度定理，將線積分轉換為面積分，則區(qū)域內的交互積分的形式為

(11)

(12)

2.2積分方案

對于裂紋或者材料界面相交的單元，由于單元內部的位移函數(shù)不連續(xù)，無法直接使用高斯積分，一種可行的方法是將單元分割成若干子單元，在每個子單元內位移函數(shù)連續(xù)可微，分別對子區(qū)域積分之后累加即可得到精確的單元積分。對于裂紋與材料界面穿過同一個單元的情況(圖3)，視裂紋與材料界面的相對位置，將單元分為若干個區(qū)域，然后以區(qū)域內的形心、裂尖或交點為基準，將各區(qū)域劃分為三角形子區(qū)域，再進行積分。對于混合單元，通過加密積分點提高計算精度，使用5×5高斯積分，對于劃分之后的單元，試算之后發(fā)現(xiàn)，每個三角形子區(qū)域取3個積分點，已經能獲得較高的精度，因此每個子區(qū)域分配3個高斯積分點。

圖3　裂紋與材料界面經過同一個單元時的積分方案

2.3應力強度因子

交互積分與應力強度因子的關系為

(13)

(14)

式中Etip為裂尖所在位置的彈性模量；νtip為裂尖所在位置的泊松比。

2.4裂紋擴展方向

判定裂紋擴展方向基于最大周向應力準則[19]，即尋找θc使該方向的剪應力為零，得到

(15)

3　 數(shù)值結果和分析

本文數(shù)值方法的實現(xiàn)基于ABAQUS軟件的二次開發(fā)。其中，有限元方程的創(chuàng)建以及求解在ABAQUS環(huán)境下完成，前后處理(模型創(chuàng)建、裂紋及夾雜初始化、水平集更新、斷裂參數(shù)計算以及節(jié)點應力計算等)使用FORTRAN編程實現(xiàn)，并通過ABAQUS子程序UEXTERNALAB進行數(shù)據(jù)交換，使用軟件TECPLOT繪制結果云圖。

3.1網(wǎng)格密度對精度的影響

為了驗證網(wǎng)格密度對本文方法精度的影響，考慮圖4(a)所示的邊緣裂紋拉伸平板，裂紋長度a/L=0.2。使用四邊形網(wǎng)格劃分模型，交互積分區(qū)域半徑取3倍單元特征長度。為了更真實地反應XFEM方法的精度，在水平和豎直方向分別劃分奇數(shù)個單元，保證裂紋面貫穿單元，且裂尖位于單元內部。

3.2夾雜對疲勞裂紋擴展壽命的影響

考慮圖4(b)所示的模型及加載方式，初始裂紋長度a=20 mm，在裂紋行進路線的兩側對稱分布著一對圓形夾雜，研究夾雜材料、尺寸及位置對疲勞裂紋擴展壽命的影響。母板彈性模量為68 320 N/mm2，泊松比為0.3，使用101×151的四邊形網(wǎng)格。

3.2.1夾雜-裂紋距離對裂紋擴展速率的影響

取夾雜半徑R=5 mm，夾雜與母材的彈性模量比EⅠ/E=2，改變夾雜邊緣之間的最短距離S，計算裂紋擴展曲線及速率。

圖6是XFEM和FEM計算模型的網(wǎng)格劃分對比?？煽吹剑褂肵FEM方法，時裂紋和夾雜是完全獨立于網(wǎng)格的。

圖7是無量綱能量釋放率曲線。其中，G0是裂紋在不含夾雜材料中的能量釋放率，裂紋擴展速率的變化正比于能量釋放率的變化。從圖7可見，夾雜性質不變的情況下，對裂紋擴展速率影響的程度與離裂紋路徑的距離成反比。對于硬夾雜，受到前方夾雜的影響，裂紋擴展速率相比于勻質材料中降低了；當裂尖經過夾雜中心位置之后，裂紋擴展速率有小幅的提高，然后逐漸趨于在勻質材料中的擴展速率。注意到結果中沒有給出S=5 mm時的更新網(wǎng)格方法的模擬結果，因為此時裂紋距離夾雜邊緣太近，沒有足夠的空間劃分裂尖奇異網(wǎng)格，S取10、15、20 mm時，XFEM模擬結果與更新網(wǎng)格法模擬結果相差極小。

(a)均質平板　　　(b)含對稱分布夾雜平板

圖5　應力強度因子數(shù)值結果相對于節(jié)點數(shù)的收斂性

(a)XFEM　　　　　　　　(b)更新網(wǎng)格法

圖7　不同S值下的裂紋擴展速率比

當S=5 mm，裂紋擴展到夾雜附近時，裂尖距離夾雜很近，此時積分域內將包含材料界面。圖8是使用標準交互積分以及本文方法計算S=5 mm時，裂紋擴展過程中的無量綱能量釋放率。從圖8可看出，當積分域包含材料界面時，標準的交互積分方法計算的結果出現(xiàn)了很大的波動，而考慮了材料界面影響的交互積分能夠得到合理的結果。

圖8　交互積分方法對數(shù)值結果的影響

3.2.2夾雜材料對裂紋擴展速率的影響

除夾雜與裂紋的距離之外，夾雜對能量釋放率的影響程度還與材料性質相關。圖9為保持夾雜與裂紋距離不變的情況下，改變夾雜材料時裂尖無量綱能量釋放率的變化。由圖9可見，硬夾雜對裂紋擴展速率影響的嚴重程度與EⅠ/E的大小成正比；軟夾雜對裂紋擴展速率的影響的作用規(guī)律與硬夾雜相反，裂紋首先加速擴展，在經過夾雜中心位置之后，又有小幅的減速；然后，逐漸趨于在勻質材料中的擴展速率，影響的嚴重程度與EⅠ/E的大小成反比。

從3.2節(jié)的模擬結果可發(fā)現(xiàn)，位于裂紋擴展路徑兩側的夾雜會影響疲勞裂紋擴展速率，軟夾雜對加速裂紋擴展，而硬夾雜會抑制疲勞裂紋擴展，影響程度與夾雜和母材的彈性模量比、半徑以及到裂紋路徑的距離有關。

圖9　不同EⅠ/E下的裂紋擴展速率比

3.3夾雜對裂紋擴展路徑的影響

3.3.1裂紋經過單個夾雜

考察圖10所示的方形平板(W=L=100 mm)，平板含一條邊緣裂紋以及一個圓形夾雜。初始裂紋長度a0=30 mm，d=5 mm。分別取夾雜和板的彈性模量比EⅠ/E為2和0.5，取裂紋擴展步長為1 mm，計算裂紋擴展路徑以及無量綱化的能量釋放率G/G0，并與文獻[14]的結果對比。

裂紋靠近夾雜時的擴展行為如圖11所示。

圖10　含夾雜邊緣裂紋平板

從圖11(b)可看出，大約x/r=-5～-6的位置開始，裂尖能量釋放率受到了夾雜的影響，當裂紋擴展到x/r=-1的位置時，裂紋才開始偏離原來的擴展方向：硬夾雜迫使裂紋向偏離夾雜的方向擴展，軟夾雜吸引裂紋向夾雜擴展；當裂紋越過夾雜中心點之后，夾雜對裂紋的影響逐漸減弱，裂紋路徑逐漸恢復到Ⅰ型裂紋的形式；裂紋擴展路徑如圖11(a)所示。當裂紋刺入軟夾雜后，裂尖處在一個被較硬基體包圍的軟材料內，因此能量釋放率遠低于裂尖位于夾雜外的情況。從圖11(c)可看到，此時的應力強度因子KI出現(xiàn)了突變，直到裂紋刺穿夾雜后，應力強度因子才恢復到原先的變化趨勢。

(a)裂紋擴展路徑

(b)能量釋放率

(c)含有軟夾雜平板的應力強度因子

改變夾雜與裂紋的相對位置與彈性模量，比較不同情況下的裂紋擴展路徑，對比結果如圖12所示：夾雜離裂紋擴展路徑越近，對裂紋路徑的影響越大；對于硬夾雜，夾雜與母材的彈性模量比EⅠ/E越大，夾雜對裂紋的影響越大；對于軟夾雜，夾雜與母材的彈性模量比EⅠ/E越小，夾雜對裂紋的影響越大。

從本節(jié)算例的結果可看出，在純拉伸狀態(tài)下，裂紋沿著垂直于載荷的方向擴展，當附近存在夾雜時，裂紋會首先向材料剛度小的方向擴展，在越過夾雜之后逐漸恢復到垂直于載荷的方向。

(a)d增大時裂紋擴展路徑

(b)EI增大時裂紋擴展路徑

3.3.2裂紋經過一對夾雜

考慮文獻[14]中裂紋經過一對夾雜的情況，夾雜半徑R=5 mm，EⅠ/E=6.43，v=0.33，v1=0.17，2個夾雜的圓心連線與水平面的夾角θ分別為60°和30°，距離S=3R。母版的尺寸、材料以及加載方式與3.3.1節(jié)相同，裂紋增量為1 mm。模擬結果如圖13所示。對于θ=30°的情況，裂紋在刺入夾雜之前的路徑與文獻[14]的計算結果基本吻合，裂紋跨過第一個夾雜之后，裂尖后方的區(qū)域被加強了，由于位置關系，第二個夾雜對裂紋的抑制作用無法迫使裂紋完全偏離夾雜而只是方向略往上偏，最終刺入夾雜。由于此時無量綱能量釋放率遠大于裂尖在夾雜外的情況，影響了曲線特征的顯示，故沒有顯示該段的數(shù)據(jù)曲線。從圖13(d)的應力強度因子曲線可看出，裂尖的應力強度因子KI在裂紋刺入夾雜之后出現(xiàn)了突變，遠大于裂尖在母材中的情況，直到裂紋刺穿夾雜后才恢復到原先的變化趨勢。圖14是進入夾雜影響區(qū)域后裂尖在幾個典型位置的y方向應力云圖，反映了一對夾雜相互作用下裂尖附近應力場的分布情況。

(a)裂紋擴展路徑

(b)θ=60°能量釋放率

(c)θ=30°能量釋放率

(d)θ=30°應力強度因子

(a)Step 20，當θ=60°　　　　(b) Step 30，當θ=60°

(c)Step 25，當θ=30°　　　　(d) Step 35，當θ=30°

3.3.3裂紋經過一組不規(guī)則分布的夾雜

在實際情況中，夾雜通常是成群存在的。使用圖6所示的模型和加載方式，在板中心的矩形區(qū)域設置圖15所示的一組12個夾雜，圖中H=20 mm，圓形夾雜的半徑R=3 mm，橢圓形夾雜的長軸半徑a=3.75 mm，短軸半徑b=2.5 mm，夾雜與母材的彈性模量比EⅠ/E=10。取初始裂紋長度a0=25 mm，裂紋增量為1 mm，計算裂紋擴展路徑。裂紋經過夾雜區(qū)域時的擴展路徑如圖15所示，裂紋方向只有輕微的變化；從圖16可看出,裂尖處在夾雜區(qū)域時無量綱能量釋放率始終小于1，表明兩側的硬夾雜抑制了疲勞裂紋的擴展。

圖15　裂紋經過一組12個夾雜

圖16　無量綱能量釋放率

4　結論

(1)與更新網(wǎng)格方法的對比結果證明,使用擴展有限元方法結合水平集方法模擬疲勞裂紋在含夾雜體中的擴展壽命具有很好的精度，且在使用較稀疏網(wǎng)格的情況下即可獲得較高的計算精度。

(2)剛度大于母材的夾雜能降低附近的裂紋擴展速率，并能迫使裂紋向偏離自己的方向擴展；雖然夾雜在較遠的位置(大約5～6倍夾雜半徑)就會影響裂尖的能量釋放率，但裂紋直到非?？拷鼕A雜時才會改變擴展方向。

(3)當夾雜剛度小于母材，或者多個夾雜綜合作用時，裂紋可能刺入夾雜，此時應力強度因子以及能量釋放率出現(xiàn)突變，直到裂紋刺穿夾雜，應力強度因子的變化才回復到原先的變化趨勢。

[1]Tamate O.The effect of a circular inclusion on the stresses around a line crack in a sheet under tension[J].International Journal of Fracture Mechanics,1968,4(3):257-266.

[2]Atkinson C.The interaction between a crack and an inclusion[J].International Journal of Engineering Science,1972,10(2):127-136.

[3]Erdogan F,Gupta G D,Ratwani M.Interaction between a circular inclusion and an arbitrarily oriented crack[J].Journal of Applied Mechanics,1974,41(4):1007-1013.

[4]Thomson R D,Hancock J W.Local stress and strain fields near a spherical elastic inclusion in a plastically deforming matrix[J].International Journal of Fracture,1984,24(3):209-228.

[5]Zhang J,Katsube N.A hybrid finite element method for heterogeneous materials with randomly dispersed elastic inclusions[J].Finite Elements in Analysis and Design,1995,19(1):45-55.

[6]Wang Y B,Chau K T.A new boundary element method for mixed boundary value problems involving cracks and holes:Interactions between rigid inclusions and cracks[J].International Journal of Fracture,2001,110:387-406.DOI:10.1023/A:1010853804657.

[7]Knight M G,Wrobel L C,Henshall J L,et al.A study of the interaction between a propagating crack and an uncoated/coated elastic inclusion using the BE technique[J].International Journal of Fracture,2002,114(1):47-61.

[8]Kitey R,Phan A V,Tippur H V,et al.Modeling of crack growth through particulate clusters in brittle matrix by symmetric-Galerkin boundary element method[J].International Journal of Fracture,2006,141(1-2):11-25.

[9]Melenk J M,Babu?ka I.The partition of unity finite element method:basic theory and applications [J].Computer methods in applied mechanics and engineering,1996,139(1):289-314.

[10]Belytschko T,Black T.Elastic crack growth in finite elements with minimal remeshing[J].International Journal for Numerical Methods in Engineering,1999,45(5):601-620.

[11]Dolbow J,Belytschko T.A finite element method for crack growth without remeshing[J].Int.J.Numer.Meth.Engng.,1999,46(1):131-150.

[12]Sukumar N,Chopp D L,Moes 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):6183-6200.

[13]Yu H,Wu L,Guo L,et al.Investigation of mixed-mode stress intensity factors for nonhomogeneous materials using an interaction integral method[J].International Journal of Solids and Structures,2009,46(20):3710-3724.

[14]Wang Z,Ma L,Wu L,et al.Numerical simulation of crack growth in brittle matrix of particle reinforced composites using the XFEM technique[J].Acta Mechanica Solida Sinica,2012,25(1):9-21.

[15]Natarajan S,Kerfriden P,Mahapatra D R,et al.Numerical analysis of the inclusion-crack interaction by the extended finite element method[J].International Journal for Computational Methods in Engineering Science and Mechanics,2014,15(1):26-32.

[16]Duflot M.A study of the representation of cracks with level sets[J].International Journal for Numerical Methods in Engineering,2007,70(11):1261-1302.

[17]Sethian J A.A fast marching level set method for monotonically advancing fronts[J].Proceedings of the National Academy of Sciences,1996,93(4):1591-1595.

[18]Hutchinson J W.Mixed mode cracking in layered materials[J].Advances in Applied Mechanics,1991:63-191.

[19]Kim J H,Paulino G H.Consistent formulations of the interaction integral method for fracture of functionally graded materials[J].Journal of Applied Mechanics,2005,72(3):351-364.

(編輯：薛永利)

Fatigue crack behavior in metallic panels with inclusions using XFEM

SHU Yi-xiu，LI Ya-zhi，SHANG Hai-jiang，JIANG Wei

(College of aeronautics，Northwestern Polytechnical University，Xi’an710072，China)

Fatigue crack propagation behaviour is remarkably affected by inclusions in homogenous material because of the stress redistribution.For the purpose of analyzing the influence pattern of inclusions on crack when a crack propagates in non-homogeneous material, an interaction integral method combined with extended finite element method(XFEM)was presented and implemented in ABAQUS using the secondary development environment.Both crack and inclusion are implicitly simulated using enrichment functions and level set method within XFEM framework.The maximum hoop stress criterion was used for crack path prediction.The situation of a crack which propagates through an area containing symmetric or asymmetric inclusions was simulated to analyze the influence of inclusions and their material properties on crack propagation rate and crack propagation path.The situation of a crack penetrating into an inclusion in certain cases was simulated and analyzed.Some of the simulations results agree well with the available results of XFEM,and are validated by the FEM remeshing approach.

fatigue crack propagation；extended finite element method；interaction integral；stress intensity factor；material interface

2015-05-13；

2015-06-24。

束一秀(1988—)，男，博士，研究方向為結構多部位損傷及疲勞斷裂。E-mail：shuyixiu@mail.nwpu.edu.cn

V250

A

1006-2793(2016)04-0547-08

10.7673/j.issn.1006-2793.2016.04.018