汪必升 李毅波,2 廖雅詩 李 劍

1.中南大學(xué)機(jī)電工程學(xué)院，長沙，4100832.中南大學(xué)高性能復(fù)雜制造國家重點實驗室，長沙，4100833.中南大學(xué)輕合金研究院，長沙，410083

0 引言

裂紋擴(kuò)展中的應(yīng)力強(qiáng)度因子不僅體現(xiàn)了裂紋尖端的應(yīng)力奇異性，還反映了應(yīng)力與幾何缺陷之間的關(guān)系，是表征裂紋尖端應(yīng)力場與位移場的關(guān)鍵參數(shù)[1]。研究動態(tài)裂紋擴(kuò)展中應(yīng)力強(qiáng)度因子的計算方法，分析影響計算精度的關(guān)鍵參數(shù)并確定其參考范圍，對提高應(yīng)力強(qiáng)度因子計算精度與效率、實現(xiàn)關(guān)鍵結(jié)構(gòu)裂紋擴(kuò)展剩余壽命預(yù)測具有重要的科學(xué)意義與應(yīng)用價值。

計算應(yīng)力強(qiáng)度因子的常用方法有解析法、實驗法和數(shù)值法。解析法和實驗法受條件所限，一般用在簡單幾何結(jié)構(gòu)和靜載條件；有限元法[2]、無網(wǎng)格法[3]及邊界元法[4]等數(shù)值法因不受幾何結(jié)構(gòu)或載荷情況的影響，可模擬出一般帶裂紋件，受到了國內(nèi)外的關(guān)注。

常規(guī)的有限元方法求解非連續(xù)位移場時依賴于單元邊界，進(jìn)行裂紋擴(kuò)展模擬時往往需要預(yù)定裂紋擴(kuò)展路徑或進(jìn)行網(wǎng)格重劃分。BELYTSCHKO等[5]在有限元法基礎(chǔ)上提出的擴(kuò)展有限元法(extended finite element method，XFEM)，不僅保留了傳統(tǒng)有限元法的優(yōu)點，還很好地解決了因材料和幾何結(jié)構(gòu)而帶來的不連續(xù)問題，故XFEM得到了迅速發(fā)展，已成為裂紋擴(kuò)展分析最主要的方法之一。由于存在裂紋尖端坐標(biāo)、裂紋與單元切割點難以精確確定以及裂紋在單元內(nèi)會發(fā)生轉(zhuǎn)折等問題，采用有效方法精確計算XFEM模型中的應(yīng)力強(qiáng)度因子一直是國內(nèi)外研究的熱點。NEHAR等[6]通過XFEM來模擬脆性和雙材料界面的裂紋擴(kuò)展，利用位移跳躍法計算應(yīng)力強(qiáng)度因子，該計算值與理論值能很好地吻合 。周博等[7]采用基于XFEM的應(yīng)力外推法對平板裂紋的應(yīng)力強(qiáng)度因子進(jìn)行計算，他們發(fā)現(xiàn)，純Ⅰ型裂紋的精度較高且易于實現(xiàn)。GUO等[8]通過相互作用積分法計算非均質(zhì)材料受熱載荷下的應(yīng)力強(qiáng)度因子，其結(jié)果與理論解有較好的一致性。ZHENG等[9]利用相互作用積分法與XFEM結(jié)合計算了帶孔裂紋平板準(zhǔn)靜態(tài)下裂紋尖端的應(yīng)力強(qiáng)度因子，并引入再分析技術(shù)，顯著地減少了計算過程中的迭代次數(shù)，但缺乏實驗驗證。艾書民等[10]利用Franc3D軟件求解裂紋擴(kuò)展過程中的應(yīng)力強(qiáng)度因子，通過疲勞壽命預(yù)測得到了驗證，但其基本原理還是基于有限元法，需對網(wǎng)格進(jìn)行加密及重復(fù)網(wǎng)格劃分。

目前，國內(nèi)外針對XFEM模型應(yīng)力強(qiáng)度因子的計算大多集中在某型斷裂模式下的準(zhǔn)靜態(tài)裂紋，有關(guān)在不同斷裂模式下動態(tài)裂紋擴(kuò)展應(yīng)力強(qiáng)度因子計算方面的報道較少，有學(xué)者雖探討了計算參數(shù)對應(yīng)力強(qiáng)度因子計算結(jié)果的影響，但并未提出明確的計算參數(shù)范圍。為此，本文基于ABAQUS軟件的XFEM模塊分析，采用相互作用積分法，通過Fortran語言開發(fā)urdfil用戶子程序，分別對中心裂紋平板Ⅰ型斷裂模式、三點彎曲Ⅱ型斷裂模式下動態(tài)裂紋擴(kuò)展過程的應(yīng)力強(qiáng)度因子進(jìn)行了計算，研究了網(wǎng)格密度和積分半徑對應(yīng)力強(qiáng)度因子計算精度的影響規(guī)律，確定了較優(yōu)的計算參數(shù)范圍；在對單邊帶孔疲勞擴(kuò)展試樣的應(yīng)力強(qiáng)度因子計算的基礎(chǔ)上，基于Paris公式預(yù)測了單邊帶孔試樣的剩余壽命，并進(jìn)行了試驗驗證。

1 擴(kuò)展有限元法的基本理論

擴(kuò)展有限元法(XFEM)在傳統(tǒng)有限元法的框架上引入富集函數(shù)，以反映裂紋擴(kuò)展帶來的非連續(xù)位移場，只有少數(shù)裂紋附近的單元被富集，大部分單元都保持傳統(tǒng)有限元單元位移模式。為此，需在有限元位移函數(shù)的表達(dá)式中增加反映局部特性的附加函數(shù)來間接模擬局部不連續(xù)性[11]，其描述裂紋含裂紋面的近似位移插值矢量函數(shù)可表示為

(1)

對于各項同性材料，裂紋尖端漸進(jìn)位移場附加函數(shù)φα的向量表達(dá)形式如下[12]:

(2)

式中，s、θ分別為積分點在極坐標(biāo)系下裂紋尖端到節(jié)點的距離和角度。

圖1所示為裂紋尖端采用單層加強(qiáng)時(即對裂紋尖端附近的一層節(jié)點進(jìn)行加強(qiáng))的加強(qiáng)節(jié)點分布，其中，“○”表示裂紋尖端節(jié)點，“□”表示貫穿節(jié)點，其余為常規(guī)節(jié)點(圖中未標(biāo)出)。

圖1 加強(qiáng)節(jié)點分布Fig.1 Enriched nodes distribution

對于含有兩種或兩種以上節(jié)點的混合單元，在式(1)位移模式下得到的節(jié)點處計算值不等于真實節(jié)點位移，需通過平移加強(qiáng)函數(shù)來糾正，相應(yīng)的位移插值矢量函數(shù)可表示為

(3)

2 基于相互作用積分法的應(yīng)力強(qiáng)度因子計算

2.1 相互作用積分法基本理論

相互作用積分法通過建立輔助應(yīng)力場和輔助位移場來獲取和分離真實場中的Ⅰ型和Ⅱ型應(yīng)力強(qiáng)度因子[13]。輔助應(yīng)力場、輔助位移場是滿足平衡方程、物理方程和幾何方程關(guān)系的任一應(yīng)力場和位移場；真實場則是物體待求裂紋尖端的應(yīng)力場和位移場。

以二維裂紋擴(kuò)展問題為例，如圖 2所示，裂紋用一條線段表示，在其尖端建立一個局部正交坐標(biāo)系，其中，Г、C0分別為繞裂紋尖端的內(nèi)圓和外圓周線；nj、mj分別為路徑Г和C0外法線的方向余弦；x1、x2分別表示方向1和方向2。

圖2 局部正交坐標(biāo)系Fig.2 Local orthogonal coordinate system

定義J積分如下[14]：

(4)

式中，w為應(yīng)變能密度；δij為Kronecker函數(shù)；σij為應(yīng)力分量；μi,1為方向i對方向1的位移偏導(dǎo)數(shù)。

為了分離混合模式下的應(yīng)力強(qiáng)度因子KⅠ和KⅡ，利用疊加原理，選兩個獨立的平衡態(tài)，狀態(tài)1為真實場狀態(tài)，狀態(tài)2為輔助場狀態(tài)。將真實狀態(tài)和輔助狀態(tài)疊加后可得

(5)

將式(5)轉(zhuǎn)換為

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

(6)

(7)

式中，I(1,2)為真實場和輔助場的交互積分；w(1,2)為真實場和輔助場的交互應(yīng)變能密度。

各向同性材料的J積分與應(yīng)力強(qiáng)度因子之間的關(guān)系可表示為

(8)

式中，E′為等效模量；E為材料的彈性模量，ν為材料的泊松比。

將兩個平衡態(tài)線性疊加可得第三個平衡態(tài)，并將第三個平衡態(tài)代入式(8)可得

(9)

取狀態(tài)2純Ⅰ型、純Ⅱ型的漸進(jìn)應(yīng)力場和位移場分別為f1、f2，由式(6)和式(9)分別可得

(10)

(11)

2.2 應(yīng)力強(qiáng)度因子計算流程

采用ABAQUS軟件建立各試樣裂紋擴(kuò)展的有限元模型，通過用戶子程序urdfil讀取運(yùn)算后的單元編號以及對應(yīng)的節(jié)點信息，包括當(dāng)前增量步、節(jié)點編號、坐標(biāo)和位移等，裂紋尖端坐標(biāo)可根據(jù)單元節(jié)點的水平集插值得到；通過指定J積分區(qū)域，判斷單元是否在積分區(qū)域內(nèi)，并對落于積分區(qū)域內(nèi)的單元節(jié)點賦予權(quán)函數(shù)；基于相互作用積分理論，構(gòu)造輔助場并計算輔助場的應(yīng)力應(yīng)變，將真實場與輔助場的應(yīng)力應(yīng)變代入式(5)～式(11)分別計算得到J(1)、J(2)和KⅠ、KⅡ，并更新增量步直至計算完成。交互計算的流程見圖3，用戶子程序以及相互作用積分過程通過Fortran語言編寫相應(yīng)代碼來完成。

圖3 算法流程Fig.3 Algorithm flow

積分區(qū)域的指定是計算流程的一個關(guān)鍵步驟，圖4所示陰影部分為擬指定的積分區(qū)域，其中R為繞裂紋尖端圓的半徑(即積分半徑)，定義相對積分半徑為

Rk=R/r

(12)

式中，r為單元特征長度；S為單元面積。

圖4 積分區(qū)域Fig.4 Integration area

當(dāng)單元落于指定積分區(qū)域內(nèi)時，為明確積分點的權(quán)函數(shù)q(x)，定義：當(dāng)積分區(qū)域內(nèi)節(jié)點到裂紋尖端的距離大于或等于半徑R時，節(jié)點權(quán)函數(shù)q(x)=1(x為方形節(jié)點)；當(dāng)積分區(qū)域內(nèi)節(jié)點到裂紋尖端的距離小于R時，其他節(jié)點權(quán)函數(shù)q(x)=0(x為圓形節(jié)點)。

3 應(yīng)力強(qiáng)度因子計算精度與影響因素分析

3.1 不同裂紋擴(kuò)展模式下應(yīng)力強(qiáng)度因子的計算精度

受軸向拉伸載荷作用的中心裂紋平板試樣 和受集中載荷作用的三點彎曲試樣常用于Ⅰ型和Ⅱ型斷裂模式下的裂紋擴(kuò)展規(guī)律的研究。分別制作滿足解析解獲得條件的中心裂紋平板標(biāo)準(zhǔn)試樣和三點彎曲標(biāo)準(zhǔn)試樣，并建立裂紋擴(kuò)展規(guī)律分析的XFEM模型，如圖5所示。

圖5a所示為中心裂紋平板試樣，高度H=400 mm，寬度W=200 mm，上下兩端承受均布載荷σ=30 MPa的作用，裂紋長度定義為2a。圖5c所示為三點彎曲試樣，寬度W=20 mm，支撐點距作用點距離L=40 mm，受集中力FP=1 N的作用，裂紋長度定義為a。兩種試樣材料均為2024鋁合金，材料屬性如下[14]：彈性模量E=68 GPa，泊松比ν=0.33，主應(yīng)力σmax=280 MPa。建立XFEM模型時，均采用四邊形平面應(yīng)力單元、結(jié)構(gòu)化網(wǎng)格劃分，兩種試樣的XFEM模型分別見圖5b和圖5d。

(a)中心裂紋平板幾何尺寸 (b)中心裂紋平板擴(kuò)展有限元模型

(c)三點彎曲試樣幾何尺寸

(d)三點彎曲試樣擴(kuò)展有限元模型圖5 不同斷裂模式下應(yīng)力強(qiáng)度因子的計算模型Fig.5 The calculation model of stress intensity factors under different fracture modes

由文獻(xiàn)[15]可知，受軸向拉伸載荷作用的中心裂紋平板的應(yīng)力強(qiáng)度因子的理論計算表達(dá)式如下：

(13)

當(dāng)L=2W時，三點彎曲試件的應(yīng)力強(qiáng)度因子計算表達(dá)式如下：

(14)

式中，F(xiàn)V為支撐點的支持力；B為板厚。

兩種試樣基于相互作用積分法的應(yīng)力強(qiáng)度因子計算值與理論值對比結(jié)果以及計算誤差分別見表1和表2。由表1可知，在中心裂紋平板試樣的裂紋動態(tài)擴(kuò)展過程中，程序計算值與理論值的偏差不超過3%；由表2可知，在三點彎曲試樣的裂紋動態(tài)擴(kuò)展過程中，程序計算值與理論值的偏差不超過2%。這表明本文所提方法和程序合理，可準(zhǔn)確地進(jìn)行不同裂紋擴(kuò)展模式下應(yīng)力強(qiáng)度因子的計算。

表1 中心裂紋平板計算誤差

表2 三點彎曲試樣計算誤差

3.2 網(wǎng)格密度與積分半徑對計算精度的影響

積分半徑和網(wǎng)格密度是影響計算精度的主要因素，為確定合適的積分半徑與網(wǎng)格密度，研究了不同參數(shù)對應(yīng)力強(qiáng)度因子計算精度的影響規(guī)律。

中心裂紋平板試樣和三點彎曲試樣在不同裂寬比(裂長比)、不同積分半徑下，應(yīng)力強(qiáng)度因子計算值與理論值的相對誤差分別見圖6a和圖6b。由圖6可以看出，對于兩種試樣，當(dāng)相對積分半徑Rk< 3時，計算誤差波動較大；當(dāng)相對積分半徑Rk≥3時，中心裂紋平板試樣和三點彎曲試樣的應(yīng)力強(qiáng)度因子計算誤差分別控制在3%和2%以內(nèi)，滿足計算要求。

(a)中心平板裂紋試樣

(b)三點彎曲試樣圖6 積分半徑對應(yīng)力強(qiáng)度因子的影響Fig.6 Effect of integral radius on stress intensity factors

對于圖6a所示的裂寬比為0.1的中心裂紋平板試樣，當(dāng)相對積分半徑Rk=4.5時，對應(yīng)的積分半徑R=20 mm， 積分邊界正好經(jīng)過裂紋另一端點，裂紋尖端的應(yīng)力奇異性導(dǎo)致應(yīng)力強(qiáng)度因子過大，因此出現(xiàn)了應(yīng)力強(qiáng)度因子突然增大的現(xiàn)象，但其計算精度仍然控制在3%以內(nèi)。

為研究網(wǎng)格密度對計算精度的要求，定義網(wǎng)格密度因子為

m=r/l

(15)

式中，l為平板的特征長度。

中心裂紋平板試樣和三點彎曲試樣在不同裂寬比(裂長比)、不同網(wǎng)格密度下，應(yīng)力強(qiáng)度因子計算值與理論值的相對誤差分別見圖7a和圖7b。由圖7可以看出，當(dāng)m≥0.012時，隨著網(wǎng)格越密(即m越小)，計算誤差越小，求解精度越高；當(dāng)m<0.012時，隨著網(wǎng)格越密，計算誤差越大。這是因為當(dāng)網(wǎng)格過密時，積分區(qū)域是繞裂紋尖端很小的一個區(qū)域，計算值受裂尖奇異性的影響，誤差偏大。當(dāng)網(wǎng)格密度因子m在0.012～0.016范圍內(nèi)時，可得到較為準(zhǔn)確的計算結(jié)果。

(a)中心裂紋平板試樣

(b)三點彎曲試樣圖7 網(wǎng)格密度對應(yīng)力強(qiáng)度因子的影響Fig.7 Effect of mesh density on stress intensity factors

4 典型試樣動態(tài)裂紋擴(kuò)展應(yīng)力強(qiáng)度因子計算與剩余壽命預(yù)測

非標(biāo)準(zhǔn)試樣動態(tài)裂紋擴(kuò)展過程應(yīng)力強(qiáng)度因子沒有解析解，無法直接驗證應(yīng)力強(qiáng)度因子計算結(jié)果的準(zhǔn)確性。將應(yīng)力強(qiáng)度因子代入Paris公式實現(xiàn)各類試樣剩余壽命的預(yù)測，是應(yīng)力強(qiáng)度因子的一個典型應(yīng)用，同時也可間接地驗證動態(tài)裂紋擴(kuò)展應(yīng)力強(qiáng)度因子計算的準(zhǔn)確性。

Paris公式計算疲勞裂紋擴(kuò)展的壽命模型可表示為[16]

da/dN=C(ΔK)n

(16)

ΔK=Kmax-Kmin

式中，N為循環(huán)次數(shù)；C、n為材料常數(shù)，本文所用參數(shù)由文獻(xiàn)[17]給出，C=1.484×10-10mm/周次、n=2.761；ΔK為應(yīng)力強(qiáng)度因子的變化，MPa·mm1/2；Kmax、Kmin分別為循壞載荷下應(yīng)力強(qiáng)度因子的最大值和最小值。

本研究的對象為非標(biāo)準(zhǔn)帶孔單邊裂紋擴(kuò)展試樣，其形狀和尺寸見圖8a。其中，試樣高度H=200 mm、寬度W=100 mm，圓孔直徑為10 mm，試樣預(yù)制初始裂紋長度a=17 mm。試樣受軸向應(yīng)力均值為36 MPa、應(yīng)力比為0.1的循環(huán)載荷σ作用，試驗在MTS 810材料試驗機(jī)上進(jìn)行，裂紋擴(kuò)展長度采用顯微鏡直讀法來獲取，測量誤差不超過0.01 mm，見圖8b。在ABAQUS軟件中利用本文所提方法及程序建立非標(biāo)試樣裂紋擴(kuò)展的XFEM模型(圖8c)并計算應(yīng)力強(qiáng)度因子，相關(guān)材料參數(shù)同前文2024鋁合金的材料屬性參數(shù)。建模時的網(wǎng)格密度因子為0.012，相對積分半徑為3。

(a)試樣幾何尺寸 (b) 裂紋擴(kuò)展試驗

(c)XFEM模型圖8 單邊帶孔平板模型Fig.8 Single side hole plate model

圖9所示為裂紋擴(kuò)展過程的仿真計算結(jié)果與試驗實測結(jié)果，可以看出，仿真計算與試驗測定的裂紋擴(kuò)展路徑基本一致。圖10所示為計算得到的動態(tài)裂紋擴(kuò)展應(yīng)力強(qiáng)度因子變化規(guī)律，可以看出，應(yīng)力強(qiáng)度因子KⅠ隨著裂紋的擴(kuò)展不斷增大；應(yīng)力強(qiáng)度因子KⅡ要小于KⅠ兩個數(shù)量級以上，這表明單邊裂紋帶孔平板疲勞擴(kuò)展以Ⅰ型為主，與實際情況相符(圖9)。

(a)疲勞仿真結(jié)果 (b)疲勞試驗結(jié)果圖9 疲勞仿真與試驗的擴(kuò)展路徑Fig.9 Extended path of fatigue simulation and test

圖10 應(yīng)力強(qiáng)度因子Fig.10 Stress intensity factors

(a)疲勞裂紋擴(kuò)展a-N曲線

(b)疲勞裂紋擴(kuò)展速率da/dN曲線圖11 單邊帶孔平板試樣的疲勞壽命預(yù)測Fig.11 Fatigue life prediction of single side hole plate specimens

單邊帶孔平板試樣的疲勞壽命預(yù)測結(jié)果見圖11。圖11a為依據(jù)應(yīng)力強(qiáng)度因子變化規(guī)律所計算得到的試樣a-N曲線，可以看出，理論預(yù)測與試驗檢測結(jié)果具有較好的一致性，試樣疲勞裂紋擴(kuò)展的預(yù)測壽命與實測壽命分別為42 085周次和39 857周次，兩者相差5.3%，在許用誤差范圍內(nèi)。圖11b為疲勞裂紋擴(kuò)展速率da/dN曲線 ，可以看出，由應(yīng)力強(qiáng)度因子計算得到的擴(kuò)展速率與試驗得到的疲勞擴(kuò)展速率基本一致，進(jìn)一步驗證了應(yīng)力強(qiáng)度因子的計算精度能夠滿足工程要求。

5 結(jié)論

(1)在ABAQUS軟件中建立了典型試樣動態(tài)裂紋擴(kuò)展的XFEM模型；基于相互作用積分方法,通過Fortran語言編制用戶子程序，實現(xiàn)了動態(tài)裂紋擴(kuò)展XFEM模型應(yīng)力強(qiáng)度因子的計算。

(2)中心裂紋平板標(biāo)準(zhǔn)試樣和三點彎曲標(biāo)準(zhǔn)試樣的動態(tài)裂紋擴(kuò)展中應(yīng)力強(qiáng)度因子計算值與理論值的對比分析結(jié)果表明：所提方法與所編制的程序可用于Ⅰ型、Ⅱ型裂紋擴(kuò)展應(yīng)力強(qiáng)度因子的準(zhǔn)確計算，計算誤差分別小于3%和2%。

(3)探討了網(wǎng)格密度和積分半徑對計算精度的影響規(guī)律，結(jié)果表明：當(dāng)網(wǎng)格密度因子為0.012～0.016、相對積分半徑Rk=3時，計算結(jié)果收斂于穩(wěn)定值，計算誤差控制在3%內(nèi)。

(4)利用所提方法與程序?qū)Ψ菢?biāo)帶孔單邊裂紋擴(kuò)展試樣的動態(tài)應(yīng)力強(qiáng)度因子進(jìn)行了計算，基于Paris公式預(yù)測了非標(biāo)準(zhǔn)試樣的疲勞裂紋擴(kuò)展壽命。疲勞裂紋擴(kuò)展實測結(jié)果驗證了模型與所提方法的正確性,仿真值與試驗值的偏差為5.3%，在許用誤差范圍內(nèi)。