張凝,尹碩輝,趙子衡

(湘潭大學 機械工程學院,湖南湘潭 411105)

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

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

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

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

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

1 基本方程

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

圖1 二維裂紋問題示意圖

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

(1)

(2)

模型的本構方程為

(3)

2 多尺度擴展元法

2.1 位移模式

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

(4)

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

圖2 裂紋結點分布圖

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

(5)

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

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

(6)

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

2.2 控制方程

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

(7)

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

(8)

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

相應單元質量矩陣元素為:

(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 可變結點

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

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

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

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

可變結點的位移模式可表示為

(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 動應力強度因子

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

圖5 J積分輪廓圖

利用散度定理,并將平衡方程和幾何方程代入,可以得到J積分表達式[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)斷裂

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

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

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

(36)

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

(37)

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

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

圖7 不同方法歸一化動強度因子結果比較

圖8 不同方法歸一化動強度因子結果比較

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

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

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

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

3.1.2 計算效率

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

表1 算例運算時間記錄表

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

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

圖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ù),板的上、左、下邊界受法向位移約束。設定計算時間為3×10-5s。

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

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

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

圖12 不同方法歸一化動強度因子結果比較

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

4 結論

通過給出多尺度擴展有限元數(shù)值計算結果與解析解、不同網(wǎng)格模型的標準擴展有限元數(shù)值計算結果進行對比。結果表明:

1) 多尺度擴展有限元數(shù)值計算結果與解析解及標準擴展有限元數(shù)值計算結果一致。

2) 多尺度有限元法比標準擴展有限元法具有更高的精度和更低的計算成本。

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