基于有限元法對裂紋尖端應(yīng)力強度因子的計算分析

袁浩，李菁，謝禹鈞，侯汶雨

(1. 遼寧石油化工大學(xué) 機械工程學(xué)院，遼寧 撫順 113001； 2. 黑龍江省北方工具有限公司 機械制造公司，黑龍江 牡丹江 157000； 3. 中國石油天然氣股份有限公司 撫順石化分公司，遼寧 撫順113001)

0 引言

由于加工工藝以及材料本身缺陷等原因，結(jié)構(gòu)構(gòu)件中會不可避免的存在缺陷或裂紋，從而影響結(jié)構(gòu)系統(tǒng)的安全性，甚至?xí)鸾Y(jié)構(gòu)失效。通過對大量破壞事故研究表明，這些裂紋的存在會導(dǎo)致含裂紋構(gòu)件的實際強度遠(yuǎn)低于理論強度[1]。傳統(tǒng)的強度觀點通常把材料視為理想材料，即材料是連續(xù)、均勻、各向同性的，但實際工況中材料很難達(dá)到理想狀態(tài)[2]。為了確保含裂紋構(gòu)件長期穩(wěn)定地安全運行，必須對不可避免存在的裂紋對構(gòu)件的影響進行預(yù)判，從而將發(fā)生損失的風(fēng)險降至最低。在斷裂力學(xué)問題的分析中，應(yīng)力強度因子K是預(yù)判含裂紋構(gòu)件發(fā)生斷裂和裂紋發(fā)生擴展速率的首要判據(jù)[3-4]。獲得應(yīng)力強度因子的方法大致上可分為解析法、數(shù)值法和實驗法。有限元數(shù)值法以計算機為平臺，利用計算機的計算能力和強大的建模能力，可以解決工程中復(fù)雜的幾何條件和邊界條件下的實際問題，而且有限元法不僅局限于線彈性問題，在研究彈塑性斷裂力學(xué)、疲勞和蠕變裂紋擴展速率等問題方面也同樣適用，已經(jīng)成為獲得應(yīng)力強度因子的主要途徑。

本文以薄板中心穿透裂紋為研究對象，利用ANSYS為平臺建立裂紋模型，分別采用J積分法、位移外推法、相互作用積分法等3種方法計算裂紋尖端的應(yīng)力強度因子K，并對3種方法求出的結(jié)果進行比較分析。還研究了受力、裂紋長度、含裂紋薄板的幾何參數(shù)等因素對應(yīng)力強度因子的影響。

1 含裂紋薄板有限元模型的建模方法

在以往的研究中，大多數(shù)裂紋的有限元模型的建立采用的是對稱半裂紋的建模方法，本文采用8節(jié)點四邊形實體PLANE183單元，通過創(chuàng)建位置重合的關(guān)鍵點模擬裂紋尖端的方法建立全裂紋模型。

裂紋尖端附近應(yīng)力場表達(dá)式為[4]：

(1)

圖1 正常單元退化為奇異單元

圖2 裂紋尖端奇異單元排列

需要注意的是，在用相互作用積分法和位移外推法計算K值時，裂紋面要與局部坐標(biāo)x軸方向保持一致，裂紋面法相向與局部坐標(biāo)系y軸保持一致[7]。在運用J積分法計算應(yīng)力強度因子時要保證裂紋面與全局坐標(biāo)系的x軸保持一致。而在裂紋尖端單元的選擇上，除J積分法在裂紋尖端不必采用奇異單元劃分網(wǎng)格外，其余兩種方法均需要在裂紋尖端采用奇異單元[8]。

2 計算應(yīng)力強度因子的幾種方法

由式(1)可知，裂紋尖端應(yīng)力場的強弱主要由K值來描述，它的大小可以說明裂尖附近區(qū)域的安全程度，是預(yù)防含裂紋構(gòu)件的失效和評估結(jié)構(gòu)安全壽命的重要參量。對于應(yīng)力強度因子的計算通常有數(shù)值法和解析法，由于數(shù)值法中的有限元法有計算方便、精度高等特點，如今已被大量地運用在斷裂參量的計算中。根據(jù)有限元思想的基本原理，有相互作用積分法、位移外推法和J積分法等多種方法可被用于應(yīng)力強度因子的計算中。

2.1 相互作用積分法

Wang S.S等[9]首先運用相互作用積分法計算了混合型裂紋的應(yīng)力強度因子。相互作用積分法與J積分法類似，都是通過對積分區(qū)域進行面積分或體積分來獲取應(yīng)力強度因子的方法，其定義式為:

(2)

(3)

由式(3)得，真實場和輔助場共同作用的J積分為：

(4)

減除出真實場和輔助場的J積分項可得附加項式(2)與應(yīng)力強度因子之間的關(guān)系式為：

(5)

(6)

2.2 位移外推法

另一種方法是通過裂紋尖端節(jié)點位移間接推出應(yīng)力強度因子，根據(jù)Ⅰ型裂紋尖端附近的位移表達(dá)式[1]：

(7)

式中，u1=μ,u2=υ，E為彈性模量，r為節(jié)點到裂紋尖端的距離。

k=3-4v平面應(yīng)變狀態(tài)

從式(7)可以看出，應(yīng)用有限元法求出ui后，便可求得應(yīng)力強度因子KⅠ。

(8)

由于節(jié)點的位移在裂紋尖端處最為明顯，所以式(8)只有在裂紋尖端(r→0)時才能獲得較為準(zhǔn)確的結(jié)果，但當(dāng)r→0時，KⅠ→∞，這在數(shù)值模擬中很難完成的。但在r很小的范圍內(nèi)，可近似地把r與KⅠ看成線性關(guān)系。因此，需要沿裂紋面選取多個較小的r值計算節(jié)點位移，從而得到多個相應(yīng)的KⅠ，然后將這些KⅠ值運用曲線擬合或插值的方法最終外推到r→0，便可得到裂紋尖端的KⅠ值。Chen[10]通過計算出圖1的H、D兩1/4分點處的應(yīng)力強度因子進行插值而得到更為精確的應(yīng)力強度因子表達(dá)式：

(9)

式中,uH、uD為H、D兩節(jié)點的位移;L為裂紋面單元的邊長。式(9)較于式(8)，將裂紋尖端附近1/4點處節(jié)點作為一個插值點，與裂紋尖端附近應(yīng)力和應(yīng)變場奇異的特點相符合，因此式(9)的計算結(jié)果更為精確。這種由裂紋尖端附近的節(jié)點位移間接推導(dǎo)出裂紋尖端應(yīng)力強度因子的方法就是位移外推法。

2.3 J積分法

對于二維問題，Rice[4]提出了與路徑無關(guān)的J積分回路積分(圖3)。J積分與應(yīng)力強度因子K一樣具有反映場強的性質(zhì)，其值可以表征裂紋尖端應(yīng)力應(yīng)變場的強弱。張巨偉等[11]也證明了J積分在用ANSYS計算應(yīng)力強度因子中的可行性和準(zhǔn)確性。它的定義表達(dá)式為：

(10)

式中：Γ為起始于裂紋下表面，終止于裂紋下表面的任意一條圍繞裂紋尖端的逆時針回路；Ti為回路Γ上任意一點處的應(yīng)力分量；W為回路Γ上任意一點處的應(yīng)變能密度；n為回路Γ上外法線方向的單位矢量；ui為回路Γ上任意一點處的位移分量；ds為回路Γ上的弧元。經(jīng)過推導(dǎo)證明，在線彈性的情況下，J積分的值與裂紋擴展能量釋放率GⅠ、應(yīng)力強度因子K之間有如下關(guān)系：

圖3 J積分回路示意圖

(11)

運用ANSYS軟件的后處理功能，經(jīng)過映射路徑等一系列ANSYS指令計算出裂紋尖端的J積分，然后將所得的 J積分值代入到式(11)中便得到裂紋尖端的應(yīng)力強度因子。

3 計算實例

如圖4所示一含有中心裂紋長、寬分別為2H=120mm、2W=100mm的薄板，其裂紋長為2a=20mm，材料特性為：彈性模量E=2×1011MPa，泊松比υ=0.3。薄板承受均勻載荷σ=10MPa。

在ANSYS環(huán)境中，建立全裂紋模型，建模中采用PLANE183單元，該單元為8節(jié)點二次實體單元，在裂紋尖端附近則采用退化的PLANE183單元，如圖5、圖6所示。圖中能清楚地看出圍繞在裂紋尖端的是三角形的奇異單元，以及經(jīng)過退化后在1/4分點處的中間節(jié)點，而在遠(yuǎn)離裂紋尖端處則采用正常的四邊形PLANE183單元。通過利用ANSYS求解器以及通用后處理的功能，運用GUI和命令流的操作便可獲得分別用3種方法就算出的KⅠ值。需要說明的是，當(dāng)使用相互作用積分法計算KⅠ時，由于ANSYS14.0暫無相互作用積分相應(yīng)的GUI操作，所以只能采用命令流的方式求解。

圖4 有限寬中心裂紋模型

圖5 全裂紋模型的網(wǎng)格劃分

圖6 裂紋尖端的奇異單元

為了研究載荷和裂紋模型幾何尺寸對上述3種求取KⅠ方法的影響，本文采用控制變量的方法分別以不同的載荷、裂紋半長、薄板長度以及薄板寬度作為變量進行求解，對比相互作用積分法、位移外推法、J積分法的計算結(jié)果，并根據(jù)應(yīng)力強度因子手冊[12]中理論值表達(dá)式：

(12)

計算相對誤差，以得到更為直觀的結(jié)果，如表(1)-表(4)所示，其中KⅠ的單位為MPa·m1/2,表中“()”內(nèi)為與理論值的相對誤差。

表1 不同載荷條件KⅠ的值

表2 不同裂紋長度條件下KⅠ的值

表3 不同薄板寬度條件下KⅠ的值

表4 不同薄板長度條件下KⅠ的值

圖7 不同邊界條件下應(yīng)力強度因子KⅠ曲線比較

4 結(jié)語

1) 通過分析J積分法、位移外推法、相互作用積分法3種計算應(yīng)力強度因子的方法，發(fā)現(xiàn)3種方法與解析法理論值的相對誤差都在合理范圍內(nèi)，證明了以上3種方法在計算應(yīng)力強度因子方面的可行性與準(zhǔn)確性。

2) 3種方法中，J積分法與相互作用積分法較位移外推法求得的應(yīng)力強度因子結(jié)果與理論值的相對誤差較小、精度較高，但由于J積分法需要定義積分路徑，所以相互作用積分法在效率方面要優(yōu)于J積分法。因此，在選擇應(yīng)用有限元法原理計算應(yīng)力強度因子的這3種方法中，相互作用積分法更有優(yōu)勢。

[1] 程靳，趙樹山. 斷裂力學(xué)[M]. 北京：科學(xué)出版社，2006.

[2] 丁遂棟，孫利民. 斷裂力學(xué)[M]. 北京：機械工業(yè)出版社，1997.

[3] Lin X B, Smith R A. Finite Element Modelling of Fatigue Crack Growth of Surface C racked Plates[J]. Engineering Fracture Mechanics , 1999, 63 : 503-522 .

[4] Rice J R. A Path Independent Integral and the Approximate Analysis of Strain Concentration by Notches and Cracks[J]. J ApplMech , 1968, 35 : 379-386 .

[5] 周益春, 鄭學(xué)軍. 材料的宏微觀力學(xué)性能[M]. 北京: 高等教育出版社, 2009 .

[6] Fehl B D, Truman K Z. An Evaluation of Fracture Mechanics Quarter Point Displacement Techniques Used for Computing S tress Intensity F actor s[J]. Engineering Structures, 1999, 21 : 406-415 .

[7] 宮經(jīng)全，張少欽，李禾，等. 基于相互作用積分法的應(yīng)力強度因子計算[J]. 南昌航空大學(xué)學(xué)報(自然科學(xué)版)，2015，29(1)：45-48.

[8] 王偉. Ⅰ- Ⅱ復(fù)合型裂紋 J 積分計算[J]. 石油化工設(shè)備，2011，40(2)：55-57.

[9] Wang S S，Yau J F，Corten H T． A mixed-mode crack analysis of rectilinear anisotropic solids using conservation laws of elasticity[J]. International Journal of Fracture，1980，16( 3) : 247-259.

[10] Chen L S, Kuang J H. A Modified Linear Extrapolation Formula for Determination of Stress Intensity Factor[J]. Int Journ of Fracture, 1992, 54 : 3-8 .

[11] 金曉娜, 張巨偉. J 積分在二維斷裂算例中的計算分析[J]. 機械工程與自動化，2008,149(4):57-59.