基于有限元法對裂紋尖端應力強度因子的計算

楊 巍，張 寧，許 良,3

(1.中航工業(yè)沈陽飛機工業(yè)(集團)有限公司 制造工程部，沈陽 110136； 2.沈陽航空航天大學 機電工程學院，沈陽 110136；3.沈陽航空航天大學 航空制造工藝數(shù)字化國防重點實驗室，沈陽 110136)

基于有限元法對裂紋尖端應力強度因子的計算

楊 巍1，張 寧2，許 良2,3

(1.中航工業(yè)沈陽飛機工業(yè)(集團)有限公司 制造工程部，沈陽 110136； 2.沈陽航空航天大學 機電工程學院，沈陽 110136；3.沈陽航空航天大學 航空制造工藝數(shù)字化國防重點實驗室，沈陽 110136)

基于ANSYS有限元軟件通過相互作用積分法建立了求解三維穿透裂紋應力強度因子的有限元模型，將有限元法和解析法求得的應力強度因子進行比較驗證了模型的準確性。研究了載荷、裂紋長度、試樣寬度、厚度對裂紋尖端應力強度因子的影響，在對比結果的基礎上分析了裂紋尖端應力強度因子的三維效應。結果表明：在不同條件下有限元模型都可以很好的模擬出應力強度因子的值，二維狀態(tài)時應力強度因子的分布規(guī)律與三維狀態(tài)時的分布規(guī)律有較大差異，出于安全的考慮不應忽略應力強度因子的三維效應，對三維應力強度因子的有限元求解有一定的指導意義。

相互作用積分方法；三維穿透裂紋；應力強度因子；ANSYS

斷裂是工程構件最危險的一種失效方式，尤其是脆性斷裂，它是突然發(fā)生破壞，斷裂前沒有明顯的征兆，這就常常引起災難性的事故。大量斷裂事故分析中發(fā)現(xiàn)，斷裂皆與結構中存在缺陷或裂紋有關。由缺陷或裂紋所引起的機械、結構的斷裂失效,是工程中最重要、最常見的失效模式[1]。傳統(tǒng)強度理論是建立在假設材料無缺陷的基礎上的,但工程實際中很多按傳統(tǒng)強度理論設計的結構由于裂紋的產生和擴展,出現(xiàn)大量斷裂失效。為保證含裂紋構件的安全性和可靠性,必須預測裂紋的擴展速率和構件的斷裂強度,在斷裂力學的工程應用中,應力強度因子K是判斷含裂紋結構的斷裂和計算裂紋擴展速率的重要參數(shù)。應力強度因子的計算方法有解析法、有限元法、邊界元法、權函數(shù)法等[2]。解析法可以得到理論上的精確解,但不適應工程實際當中復雜的條件(如結構形狀、載荷情況、裂紋形式等)。而有限元法具有強大的建模能力和可以充分利用計算機的計算能力,可以解決復雜幾何狀況,各種邊界條件下的平面和空間問題以及各向異性、熱應力和非線性問題,并能獲得較高的精度,已成為確定應力強度因子的最有效方法。

本文以平板類穿透性裂紋為研究對象，利用有限元軟件ANSYS建立了含裂紋平板的有限元模型，采用相互作用積分法計算了裂紋尖端的應力強度因子KI，與解析法的計算結果進行了對比，并且研究了載荷、裂紋長度、試樣寬度、厚度對裂紋尖端應力強度因子的影響。

1 相互作用積分法

1.1 相互作用積分的定義

如式(1)所示，RICE提出了J積分[3]，

(1)

(2)

整理得到：

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

(3)

(4)

(5)

1.2 應力強度因子的提取

對于線彈性情況，J積分、應變能釋放率和應力強度因子存在如下關系[4]：

(6)

其中,E和υ為裂紋尖端處的彈性模量和泊松比。類似地，相互作用積分

(7)

2 應力強度因子的有限元計算

在ANSYS中求解斷裂力學問題,首先要進行彈性分析或彈塑性靜力分析,然后再用特殊的后處理命令，或宏命令計算所需的斷裂參數(shù)。

2.1 有限元模型的建立

本文模型的建立是基于含中心穿透裂紋平板試樣(M(T)試樣)。如圖1所示，寬度為2W，厚度為2B，穿透裂紋長度為2a。由于裂紋尖端存在較高的應力梯度且尖端處具有奇異性,Barsoum、Henshell和Shaw[5]分別發(fā)現(xiàn)通過把裂紋尖端附近二階單元的中間節(jié)點沿裂紋尖端方向移至靠近裂尖1/4分點處(如圖2中B點)就可以使裂尖附近應力場達到奇異性。因此，裂紋尖端附近采用奇異單元劃分網格，遠離裂紋的余下部分采用常規(guī)單元進行劃分[6]。圖3為有限元模型的網格劃分，尖端網格尺寸為a/12，圍繞裂紋尖端一圈的網格數(shù)量為15，每層網格尺寸的增長比例為1∶1。建模過程中平面單元采用plane183，實體單元采用solid186，沿板厚度方向單元層數(shù)取5。裂紋所在延長線上分別由裂紋尖端至兩端網格寬度逐漸增大，比率為10∶1。模型的底端固定，板遠端作用有均布應力，彈性模量和泊松比分別為E=195 GPa和μ=0.3。通過ANSYS求解出裂紋尖端的應力、應變和位移場(圖4為等效應力云圖)并將其與相互作用積分的公式結合，再通過相互作用積分與應力強度因子的關系導出KI值。

2.2 有限元模型的驗證

GB/T6398-2000中給出的求解斷裂力學中有限大平板的應力強度因子公式為：

(8)

本文中的解析解均按照式(8)計算得出。

圖1 中心穿透裂紋平板試樣示意圖

圖2 裂紋尖端奇異單元

圖3 模型的網格劃分

圖4 裂紋尖端等效應力云圖

表1算例參數(shù)

參數(shù)符號2W/mmt/mmH/mma/mmE/GPaμσ/MPa參數(shù)值502200101950.3200

3 應力強度因子的有限元分析

隨著施加載荷、裂紋長度和試樣尺寸的改變，應力強度因子的求解結果也會隨之改變，為了研究這些因素對有限元方法計算KI值的影響，本文對比了不同條件下的有限元法和解析法的計算結果。

表2至表4為不同條件下應力強度因子計算結果對比，圖5為變化趨勢圖。

3.1 二維應力強度因子的影響因素

由表2和圖5a所示，有限元法和解析法在不同載荷的情況下具有相同的變化趨勢，均隨著載荷的增加呈線性增長。且各應力水平下的相對誤差均低于2%。模型在相同應力水平σ=200 MPa的前提下不同裂紋長度所對應的應力強度因子KI如表3和圖5(b)所示。由表3可知不同裂紋長度下KI值的相對誤差低于2.5%，說明該模型可以很好的仿真出不同裂紋長度下的KI值；在相同的裂紋長度增幅下所求得的KI值前后相差量逐漸增大即隨著裂紋長度的增加應力強度因子KI的增長速率變快。當只改變平板寬度W，取W分別為50 mm、75 mm、100 mm、130 mm、170 mm，即a/W=0.4、0.267、0.2、0.154、0.118，

表2 不同載荷時KI計算結果

表3 不同裂紋長度時KI計算結果

表4 不同平板寬度時KI計算結果

圖5 不同條件下應力強度因子KI曲線的比較

由表4和圖5c可知，隨著a/W的減小應力強度因子KI值也在減小，當a/W<0.2時KI值逐漸趨于一定值，符合斷裂力學中“無限大”平板的假設。

3.2 三維應力強度因子的影響因素

隨著斷裂力學學科的發(fā)展，許多研究者不滿足于將裂尖附近的三維問題等效為二維問題來看待，進而開展了許多關于三維應力強度因子的研究[7-14]，這些研究都表明裂紋尖端的應力場有著強烈的三維效應。圖5d為裂紋長度a/W=0.4不同厚度(t=2,4,6,10 mm)含有中心穿透裂紋板應力強度因子沿厚度變化規(guī)律。顯然，應力強度因子沿厚度的方向分布是不均勻的，但從中不難發(fā)現(xiàn)一些規(guī)律：無論式樣的厚度如何變化，在裂紋前緣的中心處(z/t=0.5)附近應力強度因子值最大，從中心分別向兩表面處(z/t=0,1)延伸應力強度因子值逐漸減小且減小幅度越來越大，在沿試樣厚度方向的兩表面處達到最小值；從整體上看隨著試樣厚度的增加裂紋前緣各處的應力強度因子的值均有所增大。

相互作用積分類似于J積分都是基于能量平衡理論產生的方法。金屬材料抵抗裂紋擴展的能力分為兩部分，分別是形成裂紋新表面所需的表面能和裂紋擴展所需的塑性變形能。試樣的中心面由于受厚度方向的約束較大，所以塑性區(qū)尺寸較表面處的塑性區(qū)要小即需要克服的塑性變形能減少，而隨著試樣厚度的增加厚度方向的約束愈來愈大，導致在其他條件相同的情況下試樣越厚其塑性區(qū)尺寸越小所需克服的塑性變形能越小[15]，這也與圖4d中不同厚度含有中心穿透裂紋板應力強度因子沿厚度變化規(guī)律相符合。

由以上分析可以看出，二維狀態(tài)時應力強度因子的分布規(guī)律與三維狀態(tài)時的分布規(guī)律有較大差異，在工程應用中人們通常會忽略這些差異從而使構件裂紋擴展速率的計算以及剩余壽命的評估不夠準確，容易留下隱患。

4 結論

(1)基于相互作用積分法建立了求解應力強度因子的三維有限元模型，并將有限元數(shù)值解與解析解進行了對比誤差僅為1.7%，證明了有限元模型求解的準確性。

(2)分別通過有限元法與解析法計算了不同載荷、裂紋長度、試樣寬度情況下的應力強度因子，結果表明有限元法可以很好的模擬不同條件下的應力強度因子，且當a/W<0.2時有限元求得的解趨于定值，此時的平板可以近似看作“無限大”板。

(3)根據(jù)所建立的模型對不同厚度試樣的應力強度因子進行求解，結果表明：應力強度因子沿厚度的分布不是均勻的，有限厚度板材板厚度中心應力強度因子最大且沿厚度方向逐漸減?。辉谒炞C的厚度范圍內試樣越厚其應力強度因子越大。因此，對有限厚度裂紋問題僅考慮二維狀態(tài)是不安全的。

[1]陳傳堯.疲勞與斷裂[M].武漢:華中科技大學出版社，2002：1-14.

[2]程勒.斷裂力學[M].北京:科學出版社，2006：8-12.

[3]Rice J R.A path independent integral and the approximate analysis of strain concentration by notches and cracks[J].J.Appl.Mech,1968(35)：379-386.

[4]于紅軍，果立成，國峰楠,等.材料屬性連續(xù)性對邊裂紋應力強度因子的影響[J].哈爾濱工業(yè)大學學報，2011,43(S1)：16-20.

[5]Fehl B D，Truman K Z.An Evaluation of Fracture Mechanics Quarter Point Displacement Techniques Used for Computing Stress Intensity Factors[J].Engineering Structures，1999，21(3)：406-415.

[6]王峰.三維應力強度因子分析及干涉預應力影響研究[D].西安：西北工業(yè)大學，2007：25-30.

[7]湯英，張曉晶，吳學仁.單邊缺口拉伸試樣噴丸強化殘余應力及其三維應力強度因子分析[J].航空學報，2011，33(7)：1265-1274.

[8]謝偉，黃其青.塑性區(qū)阻止效應對三維裂紋應力強度因子的影響研究[J].西北工業(yè)大學學報，2005，23(5)：585-588.

[9]楊政，霍春勇,郭萬林,等.有限厚度板穿透裂紋前緣附近三維彈性應力場分析[J].應用力學學報，2003，20(3)：74-78.

[10]何宇廷，傅祥炯.鋁合金 CCT 試樣裂紋擴展三維尺寸效應金相分析[J].空軍工程大學學報，2002，3(2)：5-8.

[11]Guo W.Three-dimensional analyses of plastic constraint for through-thickness cracked bodies[J].Eng Fract Mech，1999，62：383-407.

[12]Sternberg E，Sadowsky M A.Three-dimensional solution for the stress concentration around a circular hole in a plate of arbitrary thickness[J].J Appl Mech，1949，16：27-36.

[13]Hartranft R J，Sih G C.An approximate three-dimensional theory of plates with application to crack problems[J].Int J Eng Sci，1970(87)：11-29.

[14]Guo W.Elastoplastic three dimensional crack border field-I singular structure of the field[J].Eng Fract Mech，1993(46)：93-104.

[15]張詩捷,朱亦鋼,黃新躍,等.試件厚度對鋁合金疲勞裂紋擴展的影響[J].航空學報，1994，6(6)：757-760.

(責任編輯：宋麗萍 英文審校：劉紅江)

CalculationofstressintensityfactorsofcracktipsbasedonFEM

YANG Wei1,ZHANG Ning2,XU Liang2,3

(1.The Manufacturing Engineering Department,Aviation Industry Corporation of China,Shenyang 110136,China; 2.College of Mechanical and Electrical Engineering,Shenyang Aerospace University,Shenyang 110136,China; 3.Key Laboratory of Fundamental Science for National Defense of Aeronautical Digital Manufacturing Process,Shenyang 110136,China)

Based on ANSYS software,the calculation model of stress intensity factors(SIF)was built utilizing the interaction integral method to predict the SIF for 3D straight through cracks.The theoretical calculation results are used to compare with the results calculated by using the FEM method.The results of the comparison testify the accuracy of the model.The influences of load,crack length,width and thickness on crack-tip stress intensity factor are studied.Based on the comparison,the 3D effects of the crack tip stress factors are analyzed.According to the results,the calculation model of SIF can simulate the value of SIF accurately in different situations.There are huge differences between the 2D and the 3D on the distribution of SIF.The method is conductive for calculating crack-tip stress intensity factor.

interaction integral method;3D straight through crack;SIF;ANSYS

2014-04-11

楊巍(1982-)，男，遼寧沈陽人，工程師，主要研究方向：結構強度及完整性性評定，E-mail:421852615@qq.com。

2095-1248(2014)03-0019-05

O346.1

A

10.3969/j.issn.2095-1248.2014.03.004