李磅磅，鄒志鑒

（200093 上海市 上海理工大學(xué) 機械工程學(xué)院）

0 引言

由于材料自身缺陷或加工過程中造成損傷，結(jié)構(gòu)內(nèi)部不可避免會存在微裂紋，在一定條件下微裂紋會不斷擴展、積聚，最終造成材料斷裂失效[1]。裂紋擴展導(dǎo)致的低應(yīng)力斷裂是造成各種結(jié)構(gòu)、零部件失效及工程中各類重大事故的根源。在斷裂力學(xué)中，應(yīng)力強度因子KI是判斷含裂紋結(jié)構(gòu)的斷裂和計算裂紋擴展速率的重要判據(jù)[2]。應(yīng)力強度因子的計算方法有解析法、有限元法、權(quán)函數(shù)法等。有限元法具有強大的建模能力和計算能力，可以解決多種邊界條件下的復(fù)雜問題，并能獲得較高的精度，已成為確定應(yīng)力強度因子的有效方法[3]。

本文以三點彎曲試樣為研究對象，使用ABAQUS 建立有限元模型，分別采用位移外推法、應(yīng)力外推法、虛擬裂紋閉合法計算I 型裂紋的應(yīng)力強度因子KI，與試樣的解析解進行比較分析，并且研究了載荷P、裂紋長度a及試樣厚度B對I 型裂紋應(yīng)力強度因子的影響。

1 計算方法

1.1 位移外推法

根據(jù)I 型裂紋尖端附近的位移公式[4]，有

用有限元軟件求出位移ui后，代入式（1），即可求出應(yīng)力強度因子KI。令θ=π，代入裂紋張開位移值ν（r，π）得到KI計算公式為

由于式（2）保留了r的奇異項，只在（r→0）時獲得準(zhǔn)確結(jié)果。但r太靠近裂紋時，常規(guī)單元無法正確反映裂紋尖端的奇異性，反而影響有限元的計算結(jié)果。構(gòu)造數(shù)據(jù)對（ri，KIi）繪制曲線，使用外推法擬合數(shù)據(jù)點[5]，擬合直線的截距，即為應(yīng)力強度因子：

1.2 應(yīng)力外推法

根據(jù)I 型裂紋尖端附近的應(yīng)力公式[4]，有

用有限元軟件求出應(yīng)力σij，代入式（4），即可求出應(yīng)力強度因子KI。令θ=0，代入裂紋線上應(yīng)力值σy（r，0）得到KI為

與位移外推法相同，求出不同r處的應(yīng)力，計算得到相應(yīng)的應(yīng)力強度因子KIi，繪制曲線，通過外推法擬合，截距即為應(yīng)力強度因子KI。

1.3 虛擬裂紋閉合法

虛擬裂紋閉合法是根據(jù)Irwin 應(yīng)變能釋放率理論提出的[6-8]，該理論認(rèn)為裂紋擴展中釋放的能量等于閉合裂紋所需要的能量。如圖1 所示，設(shè)裂紋初始長度為a，擴展增量為Δa，根據(jù)Irwin 提出的裂紋閉合積分法，I 型裂紋應(yīng)變能釋放率可表達(dá)為

圖1 虛擬裂紋閉合法示意圖Fig.1 Schematic diagram of virtual crack closing method

式中：B——裂紋體板厚；σy y——沿著閉合裂紋面上的法向應(yīng)力；Δv——閉合裂紋張開時上下裂紋面的相對位移分量。

有限元分析中，根據(jù)虛擬裂紋擴展原理，裂紋向前擴展時應(yīng)力所做總功等于裂紋尖端節(jié)點力在該節(jié)點位移上所做的功，即

式中：Fy1——裂紋擴展前節(jié)點1 上y方向的節(jié)點力；v1,2——裂紋擴展后節(jié)點1、2 間y方向的節(jié)點位移。從而

虛擬裂紋閉合技術(shù)假設(shè)裂紋擴展一個小量Δa后，裂紋尖端的應(yīng)力場、位移場與擴展前相比保持不變。裂紋擴展前，節(jié)點3、4 的位移v3,4近似等于裂紋擴展后的節(jié)點1、2 位移v1,2，且網(wǎng)格擴展增量足夠小，所以式（8）可表示為

應(yīng)變能釋放率與應(yīng)力強度因子的關(guān)系[5]為

2 模型建立與驗證

本文采用的三點彎曲試樣參考GB/T 4161-2007《金屬材料 平面應(yīng)變斷裂韌度KIC試驗方法》中尺寸比例要求。如圖2 所示，寬度W=40 mm，厚度B=20 mm，跨距S=160 mm，試件總長L=200 mm，裂紋長度a=20 mm。材料屬性為：彈性模量E=2×105MPa，泊松比v=0.3；載荷P=200 N。

圖2 三點彎曲試樣Fig.2 Three-point bending specimen

三點彎曲試樣為對稱模型，在ABAQUS 環(huán)境中只需建立右側(cè)1/2 模型進行仿真分析，Part 中建立shell 模型，選擇planar 類型，繪制寬度40 mm，半長100 mm 的試樣；在界面屬性中，添加平面厚度20 mm；采用CPS4R 單元劃分網(wǎng)格，網(wǎng)格長度0.5 mm。建立模型后生成inp 文件，修改inp 文件輸出相關(guān)裂紋節(jié)點的應(yīng)力、位移及節(jié)點力。使用節(jié)點位移外推法、單元應(yīng)力外推法、虛擬裂紋閉合法計算I 型裂紋應(yīng)力強度因子。

根據(jù)GB/T 4161-2007 的規(guī)定，三點彎曲試樣的應(yīng)力強度因子計算表達(dá)式[9]為

經(jīng)仿真與計算，位移外推法、應(yīng)力外推法、虛擬裂紋閉合法與解析法得到的計算結(jié)果分別為17.145 0，16.536 4，16.715 1，16.778 9 MPa·mm1/2，相對誤差分別為2.18%，-1.44%和-0.38%，表明3 種有限元方法都可以準(zhǔn)確計算應(yīng)力強度因子KI值。

3 模型參數(shù)分析

應(yīng)力強度因子的值與試樣的幾何尺寸、裂紋長度及載荷大小有關(guān)，為研究試樣厚度B，裂紋長度a及載荷P對應(yīng)力強度因子的影響，修改inp 文件，計算不同條件下的應(yīng)力強度因子，并進行比較分析。

表1—表4 為不同條件下應(yīng)力強度因子計算結(jié)果對比，圖3 為變化趨勢圖。

表1 不同試樣厚度時KI 的值Tab.1 KI values at different sample thicknesses

由表1 和圖3（a）可知，隨試樣厚度的增加，應(yīng)力強度因子KI越來越小，且減小的趨勢越來越緩慢。3 種有限元法的結(jié)果與解析解一致性較好，誤差分別為-0.38%、2.18%和-1.45%，由式（11）可知，KI與B成反比關(guān)系，曲線變化符合理論結(jié)果。由圖3（b）可知，KI隨著載荷的增大而增大,有限元法與解析解呈線性增長，斜率近乎相同。根據(jù)表2 計算知，VCCT、應(yīng)力外推法和位移外推法的誤差分別為-0.38%、2.17%和-1.44%。

表2 不同載荷時KI 的值Tab.2 KI values at different loads

不同裂紋長度時的應(yīng)力強度因子KI如表3 和圖3（c）所示。隨裂紋長度的增加，應(yīng)力強度因子KI隨之增加，增長速率變快。當(dāng)a/W＜0.6 時，應(yīng)力強度因子KI增長速度較為緩慢。應(yīng)力外推法的結(jié)果整體大于解析解，位移外推法的結(jié)果整體小于解析解，VCCT 的結(jié)果誤差相對更小。

表3 不同裂紋長度時KI 的值Tab.3 KI values at different crack lengths

表4 不同網(wǎng)格長度時KI 的值Tab.4 KI values with different mesh lengths

圖3 不同條件下應(yīng)力強度因子KI 的曲線比較Fig.3 Curve comparison of stress intensity factor KI under different conditions

比較網(wǎng)格尺寸對結(jié)果精度的影響，由表4 和圖3（d）可知，當(dāng)網(wǎng)格尺寸≤1 mm 時，3 種有限元法與解析解高度一致；當(dāng)網(wǎng)格尺寸＞1 mm（Δa/a=0.05）時，隨網(wǎng)格尺寸增大，應(yīng)力外推法和位移外推法的誤差越來越大，而VCCT 的結(jié)果仍然與解析解高度一致。

4 結(jié)論

（1）建立三點彎曲試樣的有限元模型，使用虛擬裂紋閉合法、應(yīng)力外推法、位移外推法來計算應(yīng)力強度因子KI，并與解析解進行比較，誤差都在2.5%以內(nèi)，以上3 種方法在計算I 型裂紋應(yīng)力強度因子方面具有可行性和較高的準(zhǔn)確性。

（2）研究了試樣的幾何尺寸、裂紋長度及載荷大小對I 型裂紋應(yīng)力強度因子KI的影響。KI與載荷P成良好的線性關(guān)系；KI隨試樣厚度B增加而減小，當(dāng)W/B＜2 時，KI受B 影響較大；KI隨裂紋長度a增加而增大，當(dāng)a/W＞0.7 時，KI的增長速度明顯增加。在分析網(wǎng)格尺寸對結(jié)果的影響中，2 種外推法的精度隨網(wǎng)格尺寸增大而降低，虛擬裂紋閉合法始終保持1%以內(nèi)的精度。

（3）分析比較了虛擬裂紋閉合法、應(yīng)力外推法、位移外推法的特點，得到應(yīng)力外推法的誤差始終大于位移外推法，應(yīng)力場需要通過位移場求偏導(dǎo)數(shù)，所以應(yīng)力外推法精度比位移外推法低。而虛擬裂紋閉合法的精度均小于1%，且計算效率高。