楊 巍，張 寧，許 良,3

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

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

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

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

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

相互作用積分方法；三維穿透裂紋；應(yīng)力強(qiáng)度因子；ANSYS

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

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

1 相互作用積分法

1.1 相互作用積分的定義

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

(1)

(2)

整理得到：

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

(3)

(4)

(5)

1.2 應(yīng)力強(qiáng)度因子的提取

對于線彈性情況，J積分、應(yīng)變能釋放率和應(yīng)力強(qiáng)度因子存在如下關(guān)系[4]：

(6)

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

(7)

2 應(yīng)力強(qiáng)度因子的有限元計(jì)算

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

2.1 有限元模型的建立

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

2.2 有限元模型的驗(yàn)證

GB/T6398-2000中給出的求解斷裂力學(xué)中有限大平板的應(yīng)力強(qiáng)度因子公式為：

(8)

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

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

圖2 裂紋尖端奇異單元

圖3 模型的網(wǎng)格劃分

圖4 裂紋尖端等效應(yīng)力云圖

表1算例參數(shù)

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

3 應(yīng)力強(qiáng)度因子的有限元分析

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

表2至表4為不同條件下應(yīng)力強(qiáng)度因子計(jì)算結(jié)果對比，圖5為變化趨勢圖。

3.1 二維應(yīng)力強(qiáng)度因子的影響因素

由表2和圖5a所示，有限元法和解析法在不同載荷的情況下具有相同的變化趨勢，均隨著載荷的增加呈線性增長。且各應(yīng)力水平下的相對誤差均低于2%。模型在相同應(yīng)力水平σ=200 MPa的前提下不同裂紋長度所對應(yīng)的應(yīng)力強(qiáng)度因子KI如表3和圖5(b)所示。由表3可知不同裂紋長度下KI值的相對誤差低于2.5%，說明該模型可以很好的仿真出不同裂紋長度下的KI值；在相同的裂紋長度增幅下所求得的KI值前后相差量逐漸增大即隨著裂紋長度的增加應(yīng)力強(qiáng)度因子KI的增長速率變快。當(dāng)只改變平板寬度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 不同載荷時(shí)KI計(jì)算結(jié)果

表3 不同裂紋長度時(shí)KI計(jì)算結(jié)果

表4 不同平板寬度時(shí)KI計(jì)算結(jié)果

圖5 不同條件下應(yīng)力強(qiáng)度因子KI曲線的比較

由表4和圖5c可知，隨著a/W的減小應(yīng)力強(qiáng)度因子KI值也在減小，當(dāng)a/W<0.2時(shí)KI值逐漸趨于一定值，符合斷裂力學(xué)中“無限大”平板的假設(shè)。

3.2 三維應(yīng)力強(qiáng)度因子的影響因素

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

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

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

4 結(jié)論

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

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

(3)根據(jù)所建立的模型對不同厚度試樣的應(yīng)力強(qiáng)度因子進(jìn)行求解，結(jié)果表明：應(yīng)力強(qiáng)度因子沿厚度的分布不是均勻的，有限厚度板材板厚度中心應(yīng)力強(qiáng)度因子最大且沿厚度方向逐漸減小；在所驗(yàn)證的厚度范圍內(nèi)試樣越厚其應(yīng)力強(qiáng)度因子越大。因此，對有限厚度裂紋問題僅考慮二維狀態(tài)是不安全的。

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

[2]程勒.斷裂力學(xué)[M].北京:科學(xué)出版社，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ù)性對邊裂紋應(yīng)力強(qiáng)度因子的影響[J].哈爾濱工業(yè)大學(xué)學(xué)報(bào)，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]王峰.三維應(yīng)力強(qiáng)度因子分析及干涉預(yù)應(yīng)力影響研究[D].西安：西北工業(yè)大學(xué)，2007：25-30.

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

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

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

[10]何宇廷，傅祥炯.鋁合金 CCT 試樣裂紋擴(kuò)展三維尺寸效應(yīng)金相分析[J].空軍工程大學(xué)學(xué)報(bào)，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]張?jiān)娊?朱亦鋼,黃新躍,等.試件厚度對鋁合金疲勞裂紋擴(kuò)展的影響[J].航空學(xué)報(bào)，1994，6(6)：757-760.

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

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-)，男，遼寧沈陽人，工程師，主要研究方向：結(jié)構(gòu)強(qiáng)度及完整性性評定，E-mail:421852615@qq.com。

2095-1248(2014)03-0019-05

O346.1

A

10.3969/j.issn.2095-1248.2014.03.004