位移外推法和擴(kuò)展有限元法計(jì)算應(yīng)力強(qiáng)度因子的比較

張亞洲, 鐘紅, 王立強(qiáng)

(中國(guó)水利水電科學(xué)研究院, 北京 100048)

混凝土作為非均質(zhì)材料,其表面和內(nèi)部存在大量微裂縫甚至宏觀缺陷[1]?；炷翂沃械牧芽p可能損害大壩的整體性,改變其受力狀態(tài),降低其耐腐蝕性,從而縮短其使用壽命[2-4]。對(duì)混凝土抗裂性能的研究有利于延緩、減輕混凝土的開裂以及對(duì)已開裂混凝土采取止裂措施[5],從而提高混凝土結(jié)構(gòu)的耐久性和安全性[6]。其中混凝土材料中存在的微裂縫可能影響裂尖應(yīng)力強(qiáng)度因子。作為判斷結(jié)構(gòu)失效和裂縫的擴(kuò)展,高效準(zhǔn)確地計(jì)算裂尖應(yīng)力強(qiáng)度因子尤為重要。

自1957年Irwin提出應(yīng)力強(qiáng)度因子(stress intensity factor,SIF)概念后,該值成為斷裂力學(xué)的重要參量,引起廣大學(xué)者的研究[7-8]。常用的計(jì)算方法有根據(jù)定義求解的外推法和間接法,間接法利用應(yīng)力強(qiáng)度因子的相關(guān)參數(shù)求解,例如由Rybicki等[9]提出的虛擬裂縫閉合法,該法通過計(jì)算虛擬裂縫后的節(jié)點(diǎn)力和節(jié)點(diǎn)位移從而避免求解裂縫尖端復(fù)雜的應(yīng)力場(chǎng)和位移場(chǎng)。Rice[10]提出的J積分法,該法與積分路徑無關(guān),因此對(duì)于復(fù)雜的裂縫尖端可以轉(zhuǎn)換成其他路徑進(jìn)行求解,大大降低難度,并且該方法不僅適用線彈性模型[11],在彈塑性斷裂分析中也可廣泛運(yùn)用,Abaqus有限元軟件已集成基于J積分的擴(kuò)展有限元法(extended finite element method,XFEM)計(jì)算應(yīng)力強(qiáng)度因子。Song等[12]基于比例邊界有限元法進(jìn)行應(yīng)力強(qiáng)度因子的計(jì)算,該方法精度較高,是由于其相似中心會(huì)盡可能趨近于裂縫尖端位置,并且該方法不僅適用于均質(zhì)材料,對(duì)于界面材料同樣適用。文獻(xiàn)[13-14]也通過有限元法對(duì)裂尖應(yīng)力強(qiáng)度因子進(jìn)行計(jì)算。此外,學(xué)者們采用聲發(fā)射[15]、光彈性貼片[16]試驗(yàn)與有限元法相結(jié)合的方式計(jì)算裂尖應(yīng)力強(qiáng)度因子。

上述研究對(duì)裂尖應(yīng)力強(qiáng)度因子的求解做出了突出貢獻(xiàn),但是目前尚未有使用兩種不同的方法對(duì)同一模型進(jìn)行求解,并與解析解或文獻(xiàn)解進(jìn)行對(duì)比的研究。基于此,現(xiàn)通過4個(gè)經(jīng)典案例,分別使用位移外推法和擴(kuò)展有限元法計(jì)算裂尖應(yīng)力強(qiáng)度因子,并檢驗(yàn)兩種方法的精度。

1 位移外推法計(jì)算應(yīng)力強(qiáng)度因子

外推法最先由Chan等[17]利用應(yīng)力場(chǎng)和位移場(chǎng)中奇異項(xiàng)進(jìn)行展開,基于最小二乘法繪制K與r曲線擬合所得應(yīng)力強(qiáng)度因子,公式如下。

(1)

(2)

2 擴(kuò)展有限元法計(jì)算應(yīng)力強(qiáng)度因子

相比于傳統(tǒng)有限元法,XFEM以單位分解為理論基礎(chǔ)并在模型中裂紋是可以穿過網(wǎng)格,如圖1所示。XFEM常通過J積分法來計(jì)算應(yīng)力強(qiáng)度因子,J積分在裂縫尖端周圍的曲線曲面積分存在路徑無關(guān)情況,其斷裂韌度指標(biāo)以能量形式呈現(xiàn),而相互積分法作為其擴(kuò)展可以合理、自然地描述總能量的釋放[18],因此可以運(yùn)用在擴(kuò)展有限元法中。通過引入輔助場(chǎng)[式(3)～式(5)中用上標(biāo)(2)表示]與真實(shí)場(chǎng)[式(3)～式(5)中用上標(biāo)(1)表示]進(jìn)行疊加得到J積分分量,稱為相互作用積分,其表示式如式(3)所示。而J為J(1)、J(2)和M的和,因此J的表達(dá)式如式(5)所示。

圖1 擴(kuò)展有限元法的各單元示意圖

(3)

(4)

(5)

3 實(shí)例計(jì)算

3.1 三點(diǎn)彎曲梁模型

含有預(yù)制裂縫的三點(diǎn)彎曲梁模型示意圖如圖2所示,其中長(zhǎng)L=800 mm、跨度S=600 mm、高度H=200 mm、厚度B=60 mm,預(yù)制裂縫長(zhǎng)度a0=20、40、60、80和100 mm,即縫高比α范圍為0.1～0.5,跨高比β=3。模型上表面中間施加F=4 kN壓力荷載,下部左端施加水平向和豎直向邊界條件,右端施加豎直向邊界條件。材料屬性為:彈性模量E=30 GPa,泊松比υ=0.167。

圖2 三點(diǎn)彎曲梁模型

應(yīng)力強(qiáng)度因子的計(jì)算公式[19]如下。

(6)

(7)

(8)

3.1.1 位移外推法

在Abaqus中建立二維模型用于位移外推法計(jì)算應(yīng)力強(qiáng)度因子,裂尖網(wǎng)格加密處理,考慮平面應(yīng)變狀態(tài)設(shè)置網(wǎng)格為CEP4單元。

以縫高比為0.2的模型對(duì)裂縫周圍進(jìn)行4種數(shù)量網(wǎng)格細(xì)化,分別包括50、100、150和200個(gè),得到的結(jié)果如表1所示,計(jì)算相對(duì)誤差,即解析解KI減去位移外推計(jì)算得到KI的絕對(duì)值與解析解KI的比值,發(fā)現(xiàn)隨網(wǎng)格密度的增加,其相對(duì)誤差越小。根據(jù)解德等[20]實(shí)例中計(jì)算應(yīng)力強(qiáng)度因子時(shí)需要將裂縫尖端奇異點(diǎn)去掉可以提高精確度,以200個(gè)網(wǎng)格為例,繪制Ki-r曲線,根據(jù)公式法計(jì)算出的KI=0.512 MPa·m1/2。圖3中沒有去除奇異點(diǎn)所得到的KI=0.502 MPa·m1/2,相對(duì)誤差為1.95%。發(fā)現(xiàn)裂縫尖端和末端擬合效果較差,因此去除裂尖奇異點(diǎn)所得到的KI=0.511 MPa·m1/2,相對(duì)誤差為0.16%,其決定系數(shù)R2=0.999,擬合直線效果更好。因此使用位移外推法時(shí)可通過進(jìn)行對(duì)模型細(xì)化網(wǎng)格和在結(jié)果處理中去除裂縫尖端和末端的部分?jǐn)?shù)據(jù)來提高結(jié)果的準(zhǔn)確性。

表1 縫高比為0.2的模型不同網(wǎng)格數(shù)量結(jié)果

圖3 Ki-r曲線圖

基于上述要求,計(jì)算不同縫高比下運(yùn)用位移外推法計(jì)算的應(yīng)力強(qiáng)度因子如表2所示,發(fā)現(xiàn)應(yīng)力強(qiáng)度因子隨縫高比的增大而增加,且每組數(shù)據(jù)的相對(duì)誤差均小于1%。

表2 位移外推法下不同縫高比應(yīng)力強(qiáng)度因子結(jié)果

3.1.2 擴(kuò)展有限元法

在Abaqus中建立三維模型,并制定裂縫類型為XFEM裂縫,設(shè)置裂縫不發(fā)生擴(kuò)展并且圍道積分?jǐn)?shù)目設(shè)置為5個(gè),全局采用的網(wǎng)格單元類型為C3D8R。

以縫高比為0.2的模型為例,討論網(wǎng)格密度對(duì)應(yīng)力強(qiáng)度因子的影響:分別進(jìn)行了尺寸為10、5和2.5 mm正方形網(wǎng)格,得到的結(jié)果如表3所示,發(fā)現(xiàn)結(jié)果相對(duì)誤差均在0.5%以內(nèi)且網(wǎng)格細(xì)化可以略提高結(jié)果精度,但是對(duì)于三維模型,過于細(xì)化網(wǎng)格會(huì)導(dǎo)致計(jì)算速度降低,網(wǎng)格尺寸為2.5 mm計(jì)算時(shí)間是5 mm的14倍,因此綜合精度和時(shí)間考慮,在計(jì)算三維模型時(shí)可以對(duì)網(wǎng)格適當(dāng)加密。

表3 縫高比為0.2的三點(diǎn)彎曲梁不同網(wǎng)格尺寸結(jié)果

不同縫高比下運(yùn)用擴(kuò)展有限元法計(jì)算出的應(yīng)力強(qiáng)度因子如表4所示,發(fā)現(xiàn)該方法計(jì)算應(yīng)力強(qiáng)度因子的相對(duì)誤差保持在1%內(nèi)。

表4 XFEM下不同縫高比應(yīng)力強(qiáng)度因子結(jié)果

3.2 矩形板中心裂縫

含有中心預(yù)制裂縫的矩形薄板模型如圖4所示,其中高2H=20 m、寬度2W=10 m、厚度B=0.5 m,預(yù)制裂縫長(zhǎng)度2a=1、2、3、4、5和6 m,即縫寬比α范圍為0.1～0.6,模型上、下表面施加σ=1 MPa拉應(yīng)力,材料屬性為:彈性模量E=30 GPa,泊松比υ=0.167。

圖4 矩形板中心裂縫模型

根據(jù)《應(yīng)力強(qiáng)度因子手冊(cè)》[21],得到該模型下左、右裂尖的計(jì)算公式如下。

(9)

(10)

表5 不同縫寬比下無量綱應(yīng)力強(qiáng)度因子結(jié)果

表6 不同傾角下無量綱應(yīng)力強(qiáng)度因子結(jié)果

3.3 矩形板中心斜裂縫

含有中心預(yù)制裂縫的矩形薄板模型示意圖如圖5所示,其中模型高2H=20 m,寬度2W=10 m,厚度B=0.5 m,為盡量減小邊界對(duì)裂縫的影響,取預(yù)制裂縫長(zhǎng)度2a=1 m,即縫寬比為0.1,裂縫與水平方向傾角β,分別取20°、30°、45°、60°和70°。模型上、下表面施加σ=1 MPa拉應(yīng)力,材料屬性為:彈性模量E=30 GPa,泊松比υ=0.167。

圖5 矩形板中心斜裂縫模型

圖6 矩形板中心界面裂縫模型示意圖

根據(jù)《應(yīng)力強(qiáng)度因子手冊(cè)》[21],得到該模型上、下裂尖的計(jì)算公式如下。

(11)

式(11)中:KI和KII分別為Ⅰ型和Ⅱ型應(yīng)力強(qiáng)度因子,MPa·m1/2。

3.4 矩形板中心界面裂縫

含有中心預(yù)制界面裂縫的矩形薄板模型示意圖如圖11所示,其中模型高2H=20 m、寬度2W=10 m、厚度B=10 m,預(yù)制裂縫長(zhǎng)度2a=1.5和5 m,即縫寬比為0.15和0.5。模型上、下表面施加σ=1 MPa拉應(yīng)力。將模型上下平分給定不同彈性模量E1和E2來形成兩個(gè)材料,其中的E1/E2分別取2、5和10,泊松比υ1=υ2=0.3。

界面材料的應(yīng)力強(qiáng)度因子目前沒有解析解,比例邊界有限元法在計(jì)算應(yīng)力強(qiáng)度因子時(shí)得到廣泛運(yùn)用[22-25]。通過文獻(xiàn)[22]的模型和結(jié)果來與位移外推法和XFEM得到的應(yīng)力強(qiáng)度因子進(jìn)行對(duì)比。

表7 位移外推法下界面應(yīng)力強(qiáng)度因子結(jié)果

使用位移外推法法得到的結(jié)果如表7所示,發(fā)現(xiàn)K1的值遠(yuǎn)大于K2,表明張開型應(yīng)力強(qiáng)度因子占主導(dǎo)因素。因此運(yùn)用擴(kuò)展有限元法計(jì)算應(yīng)力強(qiáng)度因子時(shí)只提取K1。將K1無量綱化處理結(jié)果如表8所示,發(fā)現(xiàn)當(dāng)兩種材料的彈性模量相差較小時(shí),K1與文獻(xiàn)[22]結(jié)果較為接近。但是隨模量比值的增大,其相對(duì)誤差有增大的趨勢(shì)。這可能是由于Abaqus集成的XFEM計(jì)算應(yīng)力強(qiáng)度因子是根據(jù)M積分,而M積分在計(jì)算過程中依賴彈性模量。當(dāng)周圍的彈性模量變化過大可能導(dǎo)致計(jì)算出現(xiàn)較大誤差,從而與文獻(xiàn)和外推方法的結(jié)果相差較大。與位移外推法相比,使用擴(kuò)展有限元法計(jì)算得到的應(yīng)力強(qiáng)度因子略大。

表8 XFEM界面應(yīng)力強(qiáng)度因子結(jié)果

4 結(jié)論

介紹了位移外推法和擴(kuò)展有限元法的基本理論和計(jì)算應(yīng)力強(qiáng)度因子的公式,針對(duì)含預(yù)制裂縫單材料均質(zhì)模型和雙材料界面模型數(shù)值算例,分別采用位移外推法和擴(kuò)展有限元法計(jì)算了裂尖應(yīng)力強(qiáng)度因子。得出如下結(jié)論。

(1)對(duì)于單材料均質(zhì)模型,兩種方法的計(jì)算結(jié)果與解析解均吻合;對(duì)于雙材料界面模型,位移外推法的計(jì)算結(jié)果更加準(zhǔn)確,而擴(kuò)展有限元法在界面兩側(cè)材料彈性模量相差較小時(shí)才表現(xiàn)出良好的精度。

(2)位移外推法可通過去除裂尖的應(yīng)力奇異點(diǎn)和加密裂尖網(wǎng)格密度使計(jì)算結(jié)果更加接近解析解。擴(kuò)展有限元法由于軟件中已經(jīng)集成,計(jì)算簡(jiǎn)單直接,但是目前軟件僅能計(jì)算三維均質(zhì)模型的應(yīng)力強(qiáng)度因子。