位移外推法和擴展有限元法計算應力強度因子的比較

張亞洲, 鐘紅, 王立強

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

混凝土作為非均質材料,其表面和內部存在大量微裂縫甚至宏觀缺陷[1]。混凝土壩中的裂縫可能損害大壩的整體性,改變其受力狀態,降低其耐腐蝕性,從而縮短其使用壽命[2-4]。對混凝土抗裂性能的研究有利于延緩、減輕混凝土的開裂以及對已開裂混凝土采取止裂措施[5],從而提高混凝土結構的耐久性和安全性[6]。其中混凝土材料中存在的微裂縫可能影響裂尖應力強度因子。作為判斷結構失效和裂縫的擴展,高效準確地計算裂尖應力強度因子尤為重要。

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

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

1 位移外推法計算應力強度因子

外推法最先由Chan等[17]利用應力場和位移場中奇異項進行展開,基于最小二乘法繪制K與r曲線擬合所得應力強度因子,公式如下。

(1)

(2)

2 擴展有限元法計算應力強度因子

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

圖1 擴展有限元法的各單元示意圖

(3)

(4)

(5)

3 實例計算

3.1 三點彎曲梁模型

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

圖2 三點彎曲梁模型

應力強度因子的計算公式[19]如下。

(6)

(7)

(8)

3.1.1 位移外推法

在Abaqus中建立二維模型用于位移外推法計算應力強度因子,裂尖網格加密處理,考慮平面應變狀態設置網格為CEP4單元。

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

表1 縫高比為0.2的模型不同網格數量結果

圖3 Ki-r曲線圖

基于上述要求,計算不同縫高比下運用位移外推法計算的應力強度因子如表2所示,發現應力強度因子隨縫高比的增大而增加,且每組數據的相對誤差均小于1%。

表2 位移外推法下不同縫高比應力強度因子結果

3.1.2 擴展有限元法

在Abaqus中建立三維模型,并制定裂縫類型為XFEM裂縫,設置裂縫不發生擴展并且圍道積分數目設置為5個,全局采用的網格單元類型為C3D8R。

以縫高比為0.2的模型為例,討論網格密度對應力強度因子的影響:分別進行了尺寸為10、5和2.5 mm正方形網格,得到的結果如表3所示,發現結果相對誤差均在0.5%以內且網格細化可以略提高結果精度,但是對于三維模型,過于細化網格會導致計算速度降低,網格尺寸為2.5 mm計算時間是5 mm的14倍,因此綜合精度和時間考慮,在計算三維模型時可以對網格適當加密。

表3 縫高比為0.2的三點彎曲梁不同網格尺寸結果

不同縫高比下運用擴展有限元法計算出的應力強度因子如表4所示,發現該方法計算應力強度因子的相對誤差保持在1%內。

表4 XFEM下不同縫高比應力強度因子結果

3.2 矩形板中心裂縫

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

圖4 矩形板中心裂縫模型

根據《應力強度因子手冊》[21],得到該模型下左、右裂尖的計算公式如下。

(9)

(10)

表5 不同縫寬比下無量綱應力強度因子結果

表6 不同傾角下無量綱應力強度因子結果

3.3 矩形板中心斜裂縫

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

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

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

根據《應力強度因子手冊》[21],得到該模型上、下裂尖的計算公式如下。

(11)

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

3.4 矩形板中心界面裂縫

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

界面材料的應力強度因子目前沒有解析解,比例邊界有限元法在計算應力強度因子時得到廣泛運用[22-25]。通過文獻[22]的模型和結果來與位移外推法和XFEM得到的應力強度因子進行對比。

表7 位移外推法下界面應力強度因子結果

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

表8 XFEM界面應力強度因子結果

4 結論

介紹了位移外推法和擴展有限元法的基本理論和計算應力強度因子的公式,針對含預制裂縫單材料均質模型和雙材料界面模型數值算例,分別采用位移外推法和擴展有限元法計算了裂尖應力強度因子。得出如下結論。

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

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