蘇毅, 王生楠, 劉儉輝

(西北工業(yè)大學(xué) 航空學(xué)院, 陜西 西安 710072)

為了更好地了解非均質(zhì)材料的力學(xué)性能,常用有限元法對(duì)這種材料進(jìn)行仿真分析。但是,常規(guī)有限元法網(wǎng)格依賴性使之在處理連續(xù)-非連續(xù)問題時(shí)表現(xiàn)出明顯的局限性。擴(kuò)展有限元(XFEM)在處理這類不連續(xù)問題表現(xiàn)了極大的優(yōu)勢(shì)。因?yàn)閄FEM[1-3]所使用的網(wǎng)格與結(jié)構(gòu)內(nèi)部幾何或物理界面無關(guān),從而克服了裂紋尖端等高應(yīng)力和變形集中區(qū)內(nèi)網(wǎng)格劃分的困難,使得在模擬裂紋生長(zhǎng)過程中無需對(duì)網(wǎng)格進(jìn)行重新劃分。網(wǎng)格劃分也無需與材料細(xì)觀結(jié)構(gòu)的內(nèi)部邊界相協(xié)調(diào)。XFEM在處理不連續(xù)問題上的優(yōu)點(diǎn)使其從問世以來,在國(guó)際上得到了很快的發(fā)展和廣泛的應(yīng)用。Zhang H.H.等人[4]用XFEM模擬了粘彈性材料的夾雜問題。2006年,文獻(xiàn)[5]用擴(kuò)展有限元對(duì)正交各向異性材料中的裂紋采用了新的加強(qiáng)法。用算例證明了他的準(zhǔn)確性。王志勇等[6]運(yùn)用不含裂尖增強(qiáng)函數(shù)的擴(kuò)展有限元法,分別研究了靜態(tài)及動(dòng)態(tài)載荷作用下顆粒增強(qiáng)復(fù)合材料的斷裂行為。假定基體和顆粒都為線彈性材料并且綁定完好,研究了不同顆粒位置、數(shù)量對(duì)基體裂紋尖端斷裂參數(shù)的影響。數(shù)值模擬結(jié)果清晰地顯示了含剛性顆粒和柔性顆粒時(shí)不同的失效機(jī)制。

本文的主要內(nèi)容:①建立了含夾雜的帶裂紋問題的位移場(chǎng)給出了均勻材料的相互作用積分和非均質(zhì)材料的相互作用積分;②介紹了本文的相互作用積分的離散問題；③用3個(gè)算例模擬了帶顆粒夾雜加強(qiáng)的裂紋板,軟硬顆粒對(duì)應(yīng)力強(qiáng)度因子或能量釋放率的影響。

1 XFEM的基本原理

1.1 單位分解法

1996年Melenk和Babuska[7]及Duarte和Oden[8]先后提出了單位分解法(PUM),其基本思想是任意函數(shù)φ(x)都可以用域內(nèi)一組局部函數(shù)Ni(x)φ(x)表示,即

(1)

式中：Ni(x)為有限元形函數(shù),它形成一個(gè)單位分解

(2)

基于此,可以對(duì)有限元形狀函數(shù)根據(jù)需要進(jìn)行改進(jìn)。XFEM是基于PUM對(duì)單元的形函數(shù)加以改進(jìn),從而考慮所研究問題的不連續(xù)、奇異性和邊界層等特性。

1.2 位移逼近方程的建立

符合單位分解概念的向量函數(shù)u的逼近具有下列形式:

(3)

式中:Ni為有限元形函數(shù),φα為改進(jìn)函數(shù),M表示改

進(jìn)函數(shù)個(gè)數(shù),根據(jù)(3)式,有限元空間(φ1≡1,其他φα≡0)將是改進(jìn)空間的子空間。

XFEM的重點(diǎn)就集中在用最少的加強(qiáng)函數(shù)來描述不連續(xù)問題如裂紋位移逼近函數(shù)。用不連續(xù)函數(shù)和二維漸進(jìn)裂紋尖端位移場(chǎng)來建立各向同性線彈性的裂紋模型。用水平集函數(shù)來建立材料界面的局部改進(jìn)函數(shù),這使得在裂紋擴(kuò)展時(shí)不需要重新劃分網(wǎng)格。

圖1 裂紋和顆粒需要改進(jìn)的自由度的節(jié)點(diǎn)集

擴(kuò)展有限元的位移逼近為:

當(dāng)節(jié)點(diǎn)為常規(guī)單元的節(jié)點(diǎn)時(shí),只有第1項(xiàng);當(dāng)節(jié)點(diǎn)為被裂紋貫穿單元的節(jié)點(diǎn)時(shí),有第1項(xiàng)和第2項(xiàng);當(dāng)節(jié)點(diǎn)為裂紋尖端所在單元的節(jié)點(diǎn)時(shí),有第1項(xiàng)和第3項(xiàng);節(jié)點(diǎn)同時(shí)處于裂紋貫穿單元和裂紋尖端所在單元,應(yīng)優(yōu)先屬于裂紋尖端所在單元,加強(qiáng)方式選擇裂紋尖端所在單元節(jié)點(diǎn)加強(qiáng)方式。當(dāng)節(jié)點(diǎn)屬于夾雜單元節(jié)點(diǎn)時(shí),有第1項(xiàng)和第4項(xiàng)。也就是說當(dāng)板中只有裂紋的話只需要前3項(xiàng),當(dāng)板中只有夾雜的話需要第1項(xiàng)和第4項(xiàng),當(dāng)板中同時(shí)有裂紋和雜質(zhì)時(shí),4項(xiàng)都需要。

1.3 XFEM離散方程的建立

和常規(guī)有限元方法一樣,將(4)式代入彈性體的虛功方程,可導(dǎo)出XFEM的支配方程為

Kδ=F

(5)

式中：δ是節(jié)點(diǎn)位移列陣,K和F分別為結(jié)構(gòu)總剛度矩陣和節(jié)點(diǎn)載荷列陣,K按常規(guī)步驟由單元?jiǎng)偠染仃嚱M集而成。單元層次上的k和f分別為

(6)

單元荷載列陣f:

(7)

子矩陣和單元荷載力矢量分量表示如下:

(8)

(9)

(10)

(11)

(12)

常規(guī)應(yīng)變矩陣:

(13)

裂紋貫穿單元附加應(yīng)變矩陣:

(14)

裂紋尖端單元附加應(yīng)變矩陣:

(15)

夾雜單元附加的應(yīng)變矩陣:

(16)

1.4 裂紋尖端應(yīng)力強(qiáng)度因子計(jì)算

傳統(tǒng)的J-積分形式為:

(17)

式中:Γ是積分路徑,nj是該路徑上弧元素外法線的方向余弦。加入輔助場(chǎng)后得到J-積分形式為:

(18)

選取輔助場(chǎng)形式為:

(19)

(20)

將(18)式展開,交叉部分為相互作用積分可以表示為:

(21)

式中：q為權(quán)函數(shù),mi為單位法向量。當(dāng)裂尖的積分區(qū)域包含夾雜顆粒區(qū)域,也同時(shí)包含了2種材料的界面。因此積分區(qū)域是非均質(zhì)材料的積分。此時(shí)如果材料在空間上具有光滑變化的屬性時(shí),同時(shí)將線積分轉(zhuǎn)化為面積積分。相互作用積分可表示為:

(22)

Ih為均勻材料的相互積分,Inonh為材料非均勻性衍生出的項(xiàng)。

εij=Sijkl(x)σkl(i,j,k,l=1,2)

(23)

從上面的定義發(fā)現(xiàn),輔助應(yīng)力場(chǎng)(20)式由于不涉及材料屬性,是平衡的。但除了裂尖點(diǎn)之外,輔助位移場(chǎng)和輔助應(yīng)變場(chǎng)之間的關(guān)系部再成立。定義裂尖處的應(yīng)變場(chǎng)。

(i,j,k,l=1,2)

(24)

(25)

Ih為均勻材料的相互積分,Inonh也可以表示為:

(26)

考慮不計(jì)體力的平衡狀態(tài)及應(yīng)力張量的對(duì)稱性,有

此時(shí)上式可簡(jiǎn)化為:

εjkδ1j),iqdA

(27)

因?yàn)?/p>

(28)

經(jīng)過變換包含材料界面的相互作用積分表達(dá)式為:

(29)

(30)

至此,可以看出當(dāng)裂紋尖端距離顆粒很近時(shí),(29)式的積分仍是可行的。但如果裂尖距離顆粒非常近或已經(jīng)落在材料界面上時(shí),應(yīng)力不再滿足反平方根奇異性,則需選取相對(duì)應(yīng)的輔助場(chǎng)計(jì)算。當(dāng)積分區(qū)域包含材料界面,如圖2所示。

圖2 裂尖局部坐標(biāo)系及相互作用積分區(qū)域

對(duì)于各向同性材料相互作用積分與應(yīng)力強(qiáng)度因子的關(guān)系可以表示為:

(31)

式中：

(32)

1.5 裂紋擴(kuò)展方向

判斷裂紋擴(kuò)展方向常用的準(zhǔn)則有:最大周向拉應(yīng)力準(zhǔn)則、最大周向拉應(yīng)變準(zhǔn)則、最小應(yīng)變能密度強(qiáng)度準(zhǔn)則、等應(yīng)變能密度線上最大周向拉應(yīng)力準(zhǔn)則、能量釋放率準(zhǔn)則等[10],其中以最大周向拉應(yīng)力準(zhǔn)則計(jì)算最為簡(jiǎn)便。在最大周向拉應(yīng)力準(zhǔn)則里,裂紋尖端擴(kuò)展方向θc由下式確定:(非均質(zhì)有顆粒)

(33)

為了模擬裂紋的擴(kuò)展過程,可采用裂紋增量形式。

2 相互作用積分的數(shù)值離散

在XFEM中為了很好地運(yùn)用相互作用積分,(29)式可以被離散為:

(34)

式中：eA是積分區(qū)域A的的單元的數(shù)量;pe在單元中積分點(diǎn)的數(shù)量;|J|p為雅克比行列式;wp為積分點(diǎn)p的權(quán)重。實(shí)際位移的導(dǎo)數(shù)為:

(35)

對(duì)于非均質(zhì)材料,建立剛度矩陣時(shí),材料屬性按照積分點(diǎn)的材料屬性選取。由于在形函數(shù)的變分矩陣中含有不連續(xù)的函數(shù),對(duì)裂紋分割單元,夾雜單元如果直接進(jìn)行高斯積分,計(jì)算的結(jié)果是不準(zhǔn)確的。文獻(xiàn)[11-13]給出了積分策略。在本文中,對(duì)于常規(guī)單元采用2階Gauss積分,對(duì)于混合單元采用6階Gauss積分,對(duì)于裂紋貫穿和夾雜的單元,進(jìn)一步細(xì)分成子三角形單元,每個(gè)子三角形單元采用3階Gauss積分,對(duì)于裂尖單元,進(jìn)一步細(xì)分成子三角形單元,每個(gè)子三角形單元采用7階Gauss積分。

3 數(shù)值算例

算例1 無限大板中心裂紋與單個(gè)圓形顆粒間的相互作用

圖3 帶圓形夾雜顆粒的中心裂紋板

表1 Wang和Chau[14](E2/E=∞,υ=0.35)

表2 本文結(jié)果(E2/E1=104,υ=0.35)

通過比較發(fā)現(xiàn),本文結(jié)果和Wang和Chau[14]的結(jié)果相對(duì)誤差很小,證明了合理性。當(dāng)b/a增大時(shí),KI(-a)/K0和KI(a)/K0隨之增大,說明當(dāng)裂紋與顆粒的位置越遠(yuǎn)時(shí),顆粒加強(qiáng)的作用會(huì)減小。

算例2 單邊裂紋擴(kuò)展與單個(gè)圓形顆粒間的相互作用

如圖4所示,一個(gè)二維平板尺寸為H×L,受到單向拉伸載荷σ作用,其內(nèi)部含一個(gè)長(zhǎng)度為a的邊裂紋,圓形顆粒半徑為R。顆粒中心與裂紋尖端垂直距離為b,H=15;L=15;b/R=1;a/L=1/15;R=a;EP/E=100,10,5或是EP/E=0.01,0.1,0.2下標(biāo)P代表板的彈性模量υ1=υ2=0.35;σ=1;KI0為沒有顆粒夾雜時(shí)裂紋尖端的應(yīng)力強(qiáng)度因子,裂紋擴(kuò)展60次,每次步長(zhǎng)為0.2。

圖4 帶圓形夾雜顆粒的單邊裂紋板

圖5表示在軟硬顆粒下,裂紋擴(kuò)展軌跡。在軟顆粒的時(shí)候,裂紋擴(kuò)展偏向顆粒,顆粒彈性模量越小最后裂紋越靠近顆粒。這證明了裂紋的擴(kuò)展偏向于強(qiáng)度小的地方,在無顆粒的情況下,裂紋的路徑基本沒有什么改變,在勻質(zhì)板的情況下,裂紋路徑?jīng)]什么改變。同樣硬顆粒的時(shí)候,裂紋擴(kuò)展遠(yuǎn)離顆粒。顆粒的彈性模量越大最終越偏離原來軌跡。

圖6表示無量綱量KI/KI0隨著x/R變化趨勢(shì), 對(duì)于硬顆粒而言,應(yīng)力強(qiáng)度因子一直比無顆粒的時(shí)候小,說明了硬顆粒的確對(duì)板有加強(qiáng)作用,顆粒的彈性模量越大對(duì)板的加強(qiáng)效果越明顯。而對(duì)軟顆粒而言,應(yīng)力強(qiáng)度因子比無顆粒加強(qiáng)時(shí)大,軟顆粒降低了板的強(qiáng)度,從圖中可以看出在x/R在-2到0的區(qū)間,顆粒對(duì)板的加強(qiáng)或減弱作用明顯。

圖5 裂紋擴(kuò)展軌跡圖6 裂紋擴(kuò)展時(shí),KI/KI0的變化

算例3 孔邊單裂紋擴(kuò)展與單個(gè)圓形顆粒間的相互作用

如圖7所示,一個(gè)二維平板尺寸為L(zhǎng)×H,受到單向拉伸載荷σ作用,其中心有一個(gè)孔洞半徑為R孔邊有一長(zhǎng)度為a的裂紋,圓形顆粒半徑為R。其離孔洞中心的水平距離為b=4.5,H=15;L=15;a/W=1/30;R=a;Ep/E=100,10或是Ep/E=0.01,0.1;下標(biāo)p代表顆粒的彈性模量;υ1=υ2=0.35;σ=1;G0為沒有顆粒夾雜時(shí)裂紋尖端的能量釋放率,裂紋擴(kuò)展16次,每次步長(zhǎng)為0.2.c=0.2,0.5。

圖8表示無量綱量G/G0隨著x變化趨勢(shì),對(duì)于硬顆粒而言,能量釋放率一直比無顆粒的時(shí)候小,說明了硬顆粒的確對(duì)板有加強(qiáng)作用,顆粒的彈性模量越大對(duì)板的加強(qiáng)效果越明顯。G/G0隨著x增大越來越小,也就是隨著裂紋尖端越靠近硬顆粒,硬顆粒對(duì)裂紋尖端的能量釋放率的影響越明顯。距離c對(duì)G/G0影響也很明顯,c越小影響越明顯。而對(duì)軟顆粒而言,能量釋放率比無顆粒加強(qiáng)時(shí)大,軟顆粒降低了板的強(qiáng)度,隨著裂紋尖端越靠近軟顆粒,硬顆粒對(duì)裂紋尖端能量釋放率的影響越明顯。 因?yàn)榱鸭y板是不對(duì)稱的,裂紋的擴(kuò)展路徑是偏離水平方向的。

圖7 帶圓形夾雜顆粒的孔邊單裂紋板 圖8 裂紋擴(kuò)展時(shí),G/G0的變化 圖9裂紋擴(kuò)展的軌跡(放大)

4 結(jié) 論

1) 應(yīng)力強(qiáng)度因子是基于2個(gè)相容許的場(chǎng)實(shí)際場(chǎng)和輔助場(chǎng)。

2) 用擴(kuò)展有限元計(jì)算裂紋擴(kuò)展時(shí)很方便,不用在裂紋尖端重新劃分網(wǎng)格。

3) 顆粒對(duì)材料的強(qiáng)度的增加或減弱有很明顯的作用。

參考文獻(xiàn)：

[1] Daux C, Mo?s N, Dolbo W J, Sukumar N, Belytschko T. Arbitrary Branched and Intersecting Cracks with the Entended Finite Element Method[J]. International Journal for Numerical Methods in Engineering, 2000, 48: 1741-1760

[2] Mo?s N, Dolbow J, Belytschko T. A Finite Element Method for Crack Growth without Remeshing[J]. International Journal for Numerical Methods in Engineering, 1999, 46:131-150

[3] Sukumar N, Srolovitz D J, Baker T J, Prevost J H. Brittle Fracture in Polycrystalline Microstructures with the Extended Finite Element Method[J]. International Journal for Numerical Methods in Engineering, 1999, 45: 601-620

[4] Zhang H H, Li L X. Modeling Inclusion Problems in Viscoelastic Materials with the Extended Finite Element Method[J]. Finite Elements in Analysis and Design, 2009, 45: 721-729

[5] Asadpoure A, Mohammadi S. Developing New Enrichment Functions for Crack Simulation in Orthotropic Media by the Extended Finite Element Method[J]. International Journal for Numerical Methods in Engineering, 2007, 69:2150-2172

[6] 王志勇,馬力,吳林志,于紅軍. 基于擴(kuò)展有限元法的顆粒增強(qiáng)復(fù)合材料靜態(tài)及動(dòng)態(tài)斷裂行為研究[J]. 固體力學(xué)學(xué)報(bào),2011,32(6): 566-573

Wang Zhiyong, Ma Li, Wu Linzhi, Yu Hongjun. Investigation of Static and Dynamic Fracture Behavior of Particle-Reinforced Composite Materials by the Extended Finite Element Method[J]. Chinese Journal of Solid Mechanics, 2011, 32(6):566-573 (in Chinese)

[7] Melenk J M, Bubska I. The Partition of the Unity Finite Element Method: Basic Theory and Applications[J]. Computer Methods in Applied Mechanics and Engineering, 1996, 139: 289-314

[8] Duarte G A, Oden J T. An H-P Adaptive Method Using Clouds[J]. Computer Methods in Applied Mechanics and Engineering, 1996, 139: 237-262

[9] Yu H J, Wu L Z, Guo L C, Du S Y, He Q L. Investigation of Mixed-Mode Stress Intensity Factors for Nonhomogeneous Materials Using an Interaction Integral Method[J]. International Journal of Solids and Structures, 2009, 46:3710-3724

[10] 張行. 斷裂力學(xué)[M]. 北京:中國(guó)宇航出版社, 1990

Zhang Hang. Fracture Mechanics[M]. Beijing: China Aerospace Press, 1990 (in Chinese)

[11] Laborde P, Pommier J, Renard Y, Salaun M. High-Order Extended Finite Element Method for Cracked Domains[J]. International Journal for Numerical Methods in Engineering, 2005, 64: 354-381

[12] Mousavi S E, Xiao H, Sukumar N. Generalized Gaussian Quadrature Rules on Arbitrary Polygons[J]. International Journal for Numerical Methods in Engineering, 2010, 82: 99-113

[13] Natarajan S, Bordas S, Mahapatra D R. Numerical Integration over Arbitrary Polygonal Domains Based on Schwarz-Christoffel Conformal Mapping[J]. International Journal for Numerical Methods in Engineering, 2009, 80:103-134

[14] 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