模擬動斷裂問題的多尺度擴(kuò)展有限元方法

張凝,尹碩輝,趙子衡

(湘潭大學(xué) 機(jī)械工程學(xué)院,湖南湘潭 411105)

部件受外部動載荷影響,其內(nèi)部微觀層面的缺陷會出現(xiàn)應(yīng)力集中現(xiàn)象,以致產(chǎn)生裂紋萌芽、拓展,進(jìn)而產(chǎn)生宏觀層面的裂縫,隨著裂縫的增多最終造成整個機(jī)構(gòu)失穩(wěn)[1],因此動斷裂的研究對機(jī)構(gòu)的安全性和穩(wěn)定性是非常重要。

動斷裂裂紋存在不穩(wěn)定擴(kuò)展的特性,通過控制和引導(dǎo)裂紋擴(kuò)展成為了工程中最有效、快速的結(jié)構(gòu)分離模式[2],在結(jié)構(gòu)安全性設(shè)計、大型機(jī)構(gòu)產(chǎn)品定期維修、某些材料分割加工等具有很大的實用價值。

研究動斷裂裂紋擴(kuò)展問題,國內(nèi)外主要采用試驗研究和數(shù)值理論分析兩種手段。試驗方面,李清等[3]利用落錘沖擊試驗以及動焦散線試驗方法研究了沖擊下預(yù)制裂紋梁的動態(tài)斷裂行為;岳中文等[4-5]利用對巖體和梁進(jìn)行動態(tài)焦散線三點彎曲沖擊試驗來研究梁和巖體的動斷裂行為;楊仁樹等[6]利用動態(tài)焦散線三點彎曲沖擊試驗研究裂紋角度對相互貫通裂紋動斷裂特性的影響;鐘紅等[7]利用軸拉試驗研究混凝土和花崗巖的動態(tài)斷裂性能;Sharma[8]利用氣槍沖擊試驗結(jié)合人工神經(jīng)網(wǎng)絡(luò)分析長徑比對顆粒聚合物復(fù)合材料動態(tài)斷裂韌性的影響。

在數(shù)值理論研究方面,因動裂紋存在不穩(wěn)定擴(kuò)展,而擴(kuò)展有限元方法因不需要預(yù)制裂紋路徑優(yōu)勢,廣受學(xué)術(shù)界青睞。Song等[9]分別對擴(kuò)展有限元方法、單元刪除法和界面裂紋法對脆性材料的有限元分析,結(jié)果表明XFEM更有優(yōu)勢;Menouillard等[10]利用無網(wǎng)格技術(shù)對擴(kuò)展有限元法進(jìn)行改進(jìn)來模擬動斷裂問題;Asareh等[11]結(jié)合顯式時間積分和內(nèi)聚規(guī)律改進(jìn)擴(kuò)展有限元法,來預(yù)測脆性材料動斷裂問題;Zeng等[12]利用嵌入式有限元法和擴(kuò)展有限元法模擬動態(tài)裂紋擴(kuò)展并與實驗作比較;Yin等[13]利用多尺度擴(kuò)展等幾何方法對彈性固體進(jìn)行靜態(tài)和動態(tài)斷裂分析;Teng等[14]在擴(kuò)展有限元基礎(chǔ)上利用自適應(yīng)網(wǎng)格方法研究微觀缺陷對裂紋擴(kuò)展影響。

工程中動斷裂問題數(shù)值模擬需要很大的計算量,網(wǎng)格尺度是計算量大的主要成因之一,本文基于此問題提出多尺度擴(kuò)展有限元法,通過MATLAB進(jìn)行數(shù)值計算,并通過與經(jīng)典算例數(shù)值計算結(jié)果進(jìn)行對比,驗證該方法的有效性和可靠性。

1 基本方程

二維裂紋問題示意圖如圖1所示,邊界Γ=Γu∪Γt∪Γc,其中,Γu、Γt、Γc分別為模型的位移邊界、應(yīng)力邊界、裂紋面邊界。模型所有區(qū)域的集合用Ω表示。

圖1 二維裂紋問題示意圖

動態(tài)斷裂模型滿足以下平衡方程和邊界條件:

(1)

(2)

模型的本構(gòu)方程為

(3)

2 多尺度擴(kuò)展元法

2.1 位移模式

擴(kuò)展有限元(XFEM)是傳統(tǒng)有限元法的擴(kuò)展,相對于傳統(tǒng)有限元法,在求解不連續(xù)斷裂問題上具有更高的精度和效率[15]。XFEM的動態(tài)裂紋的位移模式可表示為

(4)

式中:Ni、Nj、Nk為標(biāo)準(zhǔn)有限元結(jié)點形函數(shù);ui為結(jié)點位移向量的連續(xù)部分;Ω為所有結(jié)點的集合;aj、bαk為富集結(jié)點變量;Ω1、Ω2分別為Heaviside富集函數(shù)和裂尖分支富集函數(shù)的結(jié)點的集合;H(x)、χα(x)分別為Heaviside函數(shù)和裂尖分支函數(shù)。Heaviside函數(shù)和裂尖分支函數(shù)結(jié)點分布如圖2所示。

圖2 裂紋結(jié)點分布圖

Heaviside函數(shù)表達(dá)式示為

(5)

式中:x為參考點;x*為裂紋上離參考點x最近的一點;en為在裂紋x*處的外法線單位向量。

裂尖分支函數(shù)能反映裂紋尖端附近的奇異性和不連續(xù)性,對于二維各向同性的彈性體問題,裂紋分支函數(shù)的表達(dá)式為

(6)

式中(r,θ)為裂紋尖端局部的極坐標(biāo)。其中,裂紋分支函數(shù)的第一項可反映裂紋面的不連續(xù)性,在裂紋尖端附近的裂紋面可表現(xiàn)為斷裂;其余3項可改善裂紋尖端奇異性。

2.2 控制方程

在忽略阻尼的影響下,動態(tài)斷裂的控制方程為

(7)

式中:δ=[u,a,b]T為未知列向量系數(shù);K為整體剛度矩陣;M為質(zhì)量矩陣以及;R為載荷矩陣。

(8)

式中:i、j為節(jié)點編號;u、a、b分別為標(biāo)準(zhǔn)結(jié)點、Heaviside函數(shù)結(jié)點、裂尖分支函數(shù)結(jié)點對剛度矩陣分量的影響上標(biāo)。

相應(yīng)單元質(zhì)量矩陣元素為:

(9)

(10)

(11)

(12)

(13)

(14)

(15)

(16)

(17)

(18)

(19)

(20)

單位載荷矩陣元素為:

(21)

(22)

(23)

式中:Fb為體力;Fs為面力;F為集中力。

2.3 時間積分方案

控制方程求解方式有隱式積分和顯示積分兩大類,本文采用隱式積分中的Newmark積分對控制方程時間積分。離散方程在第n個計算步可以表示為[13]:

(24)

式中:α、β為Newmark法參數(shù);Δt為計算時間的步長增量。

2.4 可變結(jié)點

計算量與網(wǎng)格規(guī)模有關(guān),通過適當(dāng)?shù)脑黾訂卧?guī)模,可以提高計算精度。但同尺度加密網(wǎng)格單元會造成計算量的浪費(fèi),因此只針對重要部位采用局部網(wǎng)格加密是更好的方法,如圖3所示。

圖3 多尺度網(wǎng)格示意圖

采用可變結(jié)點單元[16]連接不同尺度網(wǎng)格,如圖3交叉部分網(wǎng)格單元所示。可變結(jié)點單元四周邊上的結(jié)點數(shù)在允許的范圍內(nèi)可以是任意值,如圖4所示,雖然結(jié)點數(shù)增加,但任一兩個相鄰點的插值特性仍是線性的。

圖4 (4+k+l+m+n)結(jié)點單元示意圖

可變結(jié)點的位移模式可表示為

(25)

(26)

插值點表示為

uh(ξ)=aTp(ξ)=UTq-1p(ξ)

(27)

式中:q=[p(ξ1),…,p(ξ4+k+l+m+n)];U=[u1,…,u4+k+l+m+n]T。通過聯(lián)立式(25)與式(27)可以得到

[N1,…,N4+k+l+m+n]T=q-1p(ξ)

(28)

2.5 動應(yīng)力強(qiáng)度因子

本文計算應(yīng)力強(qiáng)度因子是利用互作用積分[17]形式推導(dǎo)。因考慮二維裂紋問題,因此裂紋可以用J積分描述,如圖5所示,用一條周線Γ0圍繞裂尖。

圖5 J積分輪廓圖

利用散度定理,并將平衡方程和幾何方程代入,可以得到J積分表達(dá)式[18]為

(29)

J(1,2)=J(1)+J(2)+I(1,2)

(30)

互作用積分I(1,2)為

(31)

(32)

因此,第三平衡狀態(tài)可以改寫成

(33)

通過式(30)與式(33)可以得到

(34)

令平衡狀態(tài)2為Ⅰ型和Ⅱ型的漸近場,可以得到

(35)

3 數(shù)值算例

3.1 算例1:半無限板Ⅰ型裂紋動態(tài)斷裂

研究半無限板上Ⅰ型裂紋在拉力作用下的動裂紋擴(kuò)展。如圖6所示,在半無限板選取L×2H=10×4 m有限區(qū)域為研究對象,裂紋長度a=5 m。板的彈性模量E=210 GPa,泊松比ν=0.3,密度ρ=8 000 kg/m3。有限區(qū)域頂端受到一分布載荷σ0=1 500 MPa,為避免應(yīng)力波影響,計算時間為t=1×10-3s≤3tc=H/cd=1.009×10-3s[18],其中,cd為應(yīng)力波的傳播速度,tc為應(yīng)力波第一次到達(dá)裂尖的時間。

圖6 Ⅰ型裂紋計算模型示意圖

為了方便比較,將動應(yīng)力強(qiáng)度因子歸一化,在不考慮裂紋擴(kuò)展的情況下半無限板Ⅰ型裂紋動態(tài)應(yīng)力強(qiáng)度因子計算公式為[19]

(36)

定義歸一化動態(tài)應(yīng)力強(qiáng)度因子的相對誤差為

(37)

3.1.1 數(shù)值算例驗證與比較

為驗證多尺度擴(kuò)展有限元方法的優(yōu)劣性,分別對算例的解析解、標(biāo)準(zhǔn)擴(kuò)展有限元法以及多尺度擴(kuò)展有限元法的數(shù)值解作比較,并計算相對誤差,如圖7和圖8所示。

圖7 不同方法歸一化動強(qiáng)度因子結(jié)果比較

圖8 不同方法歸一化動強(qiáng)度因子結(jié)果比較

數(shù)值計算中的網(wǎng)格模型分為65×25、多尺度網(wǎng)格、79×41、195×75這4種,如圖9所示。

圖9 半無限板網(wǎng)格示意圖

多尺度網(wǎng)格模型是由原始網(wǎng)格模型65×25中的裂紋附近網(wǎng)格進(jìn)行3×3細(xì)化,共3 239個網(wǎng)格單元,計算自由度為6 812。與之作比較的網(wǎng)格單元規(guī)模為65×25(自由度3 592)、79×41(自由度6 908)、195×75(自由度30 212),對應(yīng)優(yōu)化比較、同等自由度比較、同等細(xì)網(wǎng)格尺寸比較。

從圖7和圖8可以看出多尺度擴(kuò)展有限元法數(shù)值解與同等細(xì)網(wǎng)格尺寸下的標(biāo)準(zhǔn)擴(kuò)展有限元法數(shù)值解能很好的吻合解析解,但同等自由度下的標(biāo)準(zhǔn)擴(kuò)展有限元法精度沒有多尺度擴(kuò)展有限元法精度高。結(jié)果表明,多尺度擴(kuò)展有限元法比標(biāo)準(zhǔn)擴(kuò)展有限元法有更高精度。

3.1.2 計算效率

如表1所示,對各個模型進(jìn)行計算所消耗的時間進(jìn)行記錄比較,可以看出,多尺度擴(kuò)展有限元法比195×75網(wǎng)格模型的標(biāo)準(zhǔn)擴(kuò)展有限元法計算的時間短很多,即在得到相似精度數(shù)值解情況下,多尺度擴(kuò)展有限元法相較于標(biāo)準(zhǔn)擴(kuò)展元法計算效率提高188%左右。

表1 算例運(yùn)算時間記錄表

3.2 算例2:邊界斜裂紋板動斷裂

為驗證混合裂紋模式下多尺度擴(kuò)展有限元法的適用性和可靠性,如圖10所示,對帶有斜裂紋的矩形板進(jìn)行計算。

圖10 邊界斜裂紋板示意圖

板的幾何尺寸L×H=44 mm×32 mm,裂紋距離左邊界的距離D=16 mm,裂紋長度a=22.63 mm,與下邊界的角度成45°銳角。板的彈性模量E=29.4 GPa,泊松比ν=0.286,密度ρ=2 450 kg/m3。右邊界受一均布拉伸沖擊載荷σ0H(t),其中H(t)為Heaviside階躍函數(shù),板的上、左、下邊界受法向位移約束。設(shè)定計算時間為3×10-5s。

網(wǎng)格劃分如圖11所示,用于標(biāo)準(zhǔn)擴(kuò)展有限元法數(shù)值計算的網(wǎng)格模型為45×21;用于多尺度擴(kuò)展有限元法數(shù)值計算網(wǎng)格模型如圖11b)所示,原始網(wǎng)格模型為15×7,裂紋附近采用3×3網(wǎng)格細(xì)化。

圖11 邊界斜裂紋板網(wǎng)格示意圖

采用同等細(xì)網(wǎng)格比較,通過數(shù)值計算驗證多尺度擴(kuò)展有限元的實用性和可靠性,計算結(jié)果如圖12所示。

圖12 不同方法歸一化動強(qiáng)度因子結(jié)果比較

從圖12可看出,標(biāo)準(zhǔn)擴(kuò)展有限元法和多尺度擴(kuò)展有限元法兩種不同方法計算出來的數(shù)值解吻合較好,表明多尺度擴(kuò)展有限元也適用于分析混合裂紋模式的動裂紋擴(kuò)展。

4 結(jié)論

通過給出多尺度擴(kuò)展有限元數(shù)值計算結(jié)果與解析解、不同網(wǎng)格模型的標(biāo)準(zhǔn)擴(kuò)展有限元數(shù)值計算結(jié)果進(jìn)行對比。結(jié)果表明:

1) 多尺度擴(kuò)展有限元數(shù)值計算結(jié)果與解析解及標(biāo)準(zhǔn)擴(kuò)展有限元數(shù)值計算結(jié)果一致。

2) 多尺度有限元法比標(biāo)準(zhǔn)擴(kuò)展有限元法具有更高的精度和更低的計算成本。

3) 多尺度擴(kuò)展有限元法適用于混合裂紋下動斷裂裂紋擴(kuò)展。