（海軍航空工程學(xué)院 青島校區(qū)，山東 青島 266041）

基于Seam型裂紋的復(fù)合型裂紋起裂角預(yù)測(cè)的有限元法

王躍，穆志韜，劉濤，牛勇

（海軍航空工程學(xué)院 青島校區(qū)，山東 青島 266041）

目的對(duì)處于多軸應(yīng)力狀態(tài)下的復(fù)合型裂紋起裂角進(jìn)行預(yù)測(cè)。方法考慮材料的彈性性能，利用商業(yè)有限元軟件ABAQUS中的Seam裂紋模擬穿透型裂紋損傷，建立雙軸載荷作用下中心帶孔邊裂紋板的三維有限元模型，通過(guò)位移外推法得到裂紋尖端的應(yīng)力強(qiáng)度因子（Stress Intensity Factor，SIF），并利用求得的SIF基于最大周向應(yīng)力斷裂準(zhǔn)則預(yù)測(cè)孔邊裂紋的起裂角。結(jié)果通過(guò)與相關(guān)文獻(xiàn)結(jié)果進(jìn)行對(duì)比，發(fā)現(xiàn)通過(guò)該有限元方法計(jì)算得到的應(yīng)力強(qiáng)度因子與文獻(xiàn)結(jié)果的最大誤差在2%以內(nèi)，預(yù)測(cè)的裂紋起裂角與文獻(xiàn)結(jié)果的最大誤差在3%以內(nèi)。結(jié)論該有限元方法計(jì)算的裂紋起裂角與文獻(xiàn)結(jié)果一致，因此文中求解復(fù)合型裂紋起裂角所使用的有限元方法是可靠有效的。

復(fù)合型裂紋；起裂角；Seam型裂紋；應(yīng)力強(qiáng)度因子；最大周向應(yīng)力準(zhǔn)則

機(jī)身結(jié)構(gòu)需要承受由飛行條件引起的機(jī)身剪力和彎曲力，開口損傷逐漸形成裂紋。這不僅對(duì)機(jī)身結(jié)構(gòu)的使用功能和結(jié)構(gòu)完整性產(chǎn)生消極的影響，甚至導(dǎo)致整個(gè)結(jié)構(gòu)系統(tǒng)失效，釀成飛行事故。形成的裂紋依然處于多向載荷的共同作用下，裂紋的擴(kuò)展往往表現(xiàn)為復(fù)合型裂紋的特點(diǎn)[1—4]。復(fù)合型裂紋的擴(kuò)展不同于單純張開型裂紋，具體表現(xiàn)在其擴(kuò)展軌跡不再沿原裂紋的軌跡方向，而是與原裂紋方向成某一角度[1—12]。斷裂力學(xué)是研究此類問(wèn)題的有效手段之一，它是利用合適的斷裂準(zhǔn)則和應(yīng)力強(qiáng)度因子（SIF）確定裂紋起裂角和裂紋擴(kuò)展時(shí)的臨界狀態(tài)。Erdogan[5]等提出最大周向應(yīng)力準(zhǔn)則，假定裂紋沿最大軸向力擴(kuò)展。Sih[6]提出最小應(yīng)變能密度準(zhǔn)則，對(duì)復(fù)合型裂紋的起裂進(jìn)行了預(yù)測(cè)。Theocaris[7]等改進(jìn)最小應(yīng)變能密度準(zhǔn)則，提出了裂紋沿著最大彈性應(yīng)變能密度開裂的T準(zhǔn)則。Khan[8]等對(duì)裂紋尖端的塑性區(qū)進(jìn)行研究，發(fā)現(xiàn)裂紋的擴(kuò)展路徑與從裂紋尖端到彈塑性邊界最短路徑一致。由于建模能力強(qiáng)大，并可模擬復(fù)雜幾何邊界條件和本構(gòu)關(guān)系等優(yōu)點(diǎn)，有限元方法越來(lái)越成為工程人員所青睞的最有效的SIF求解方法。Henshell[9]等提出了四分之一結(jié)點(diǎn)奇異等參元，利用裂紋尖端附近的應(yīng)力場(chǎng)（或位移場(chǎng)），進(jìn)行回歸分析、外推計(jì)算等方法間接獲得裂紋尖端的SIF。楊綠峰[10]等對(duì)Williams級(jí)數(shù)進(jìn)行了改進(jìn)，建立了奇異域子單元的局部位移場(chǎng)，通過(guò)解Ⅰ-Ⅱ混合型Williams單元?jiǎng)偠确匠糖蠼鈶?yīng)力強(qiáng)度因子。

有限元分析軟件ABAQUS中有三種定義裂紋的方式[11—12]。Seam型裂紋可以求解裂紋尖端的應(yīng)力應(yīng)變場(chǎng)，但無(wú)法輸出裂紋尖端的應(yīng)力強(qiáng)度因子。使用Contour定義裂紋求解裂紋尖端SIF時(shí)，需要事先對(duì)裂紋擴(kuò)展方向的相關(guān)參數(shù)進(jìn)行設(shè)定，而裂紋的擴(kuò)張方向?qū)IF的影響很大。使用XFEM定義擴(kuò)展型裂紋預(yù)測(cè)擴(kuò)展路徑時(shí)，需要的材料斷裂參數(shù)繁多，較難獲得。文中利用Seam型裂紋受力后裂紋尖端的應(yīng)變，采用位移外推法得到裂紋尖端附近的SIF，進(jìn)行回歸分析得到裂紋尖端處的SIF，并與相關(guān)文獻(xiàn)進(jìn)行對(duì)比。采用Erdogan等提出最大周向應(yīng)力準(zhǔn)則對(duì)裂紋的起裂角進(jìn)行預(yù)測(cè)。

1 應(yīng)力強(qiáng)度因子的求解

線彈性斷裂力學(xué)中，SIF是分析裂紋擴(kuò)展和失效的一個(gè)重要指標(biāo)。文中利用裂紋尖端附近的位移場(chǎng)，使用外推的方法對(duì)裂紋尖端的SIF進(jìn)行求解。機(jī)身結(jié)構(gòu)件中的裂紋可以使用三維模型來(lái)研究，裂紋的三維模型中裂紋前沿任意點(diǎn)處局部坐標(biāo)系和位移場(chǎng)如圖1所示[13]。

圖1 三維模型中裂紋前沿任意點(diǎn)處局部坐標(biāo)系和位移場(chǎng)

由圖1可以看出，三維裂紋模型中裂紋前沿任意點(diǎn)處的局部直角坐標(biāo)系由x，y和z三坐標(biāo)軸構(gòu)成，其中x-y坐標(biāo)平面構(gòu)成了裂紋前沿任意點(diǎn)處的法平面，y-z坐標(biāo)平面則構(gòu)成了裂紋前沿任意點(diǎn)處的切平面。模型中裂紋前沿任意點(diǎn)處的局部柱坐標(biāo)系的三坐標(biāo)軸為r，θ和z，柱坐標(biāo)系下裂紋前沿任意點(diǎn)處的法平面由r和θ所構(gòu)成。使用三維裂紋模型中的局部柱坐標(biāo)系來(lái)表示局部直角坐標(biāo)系下裂紋前沿任意點(diǎn)處的位移場(chǎng)，可得[13]：

式中：μ，ν和ω分別對(duì)應(yīng)于裂紋尖端局部坐標(biāo)的位移；r和θ為計(jì)算點(diǎn)在局部柱坐標(biāo)的坐標(biāo)值；G為剪切模量；λ為泊松比，對(duì)于平面應(yīng)力，k=3-λ/（1+λ），而對(duì)于平面應(yīng)變，k=3-4λ；KⅠ，KⅡ和KⅢ分別為Ⅰ型Ⅱ型和Ⅲ型應(yīng)力強(qiáng)度因子；()ro為高階無(wú)窮小量。

根據(jù)式（1），若裂紋表面的位移已知，可以導(dǎo)出對(duì)稱裂紋的應(yīng)力強(qiáng)度因子計(jì)算公式[14—15]：

式中：Δμ，Δν和Δω為一裂紋面相對(duì)于另一裂紋面的位移。通過(guò)建立三維線彈性有限元模型，可以求解結(jié)構(gòu)在單元節(jié)點(diǎn)處的位移及應(yīng)力，將得到的位移代入式（2）就可求得該節(jié)點(diǎn)處的應(yīng)力強(qiáng)度因子。通過(guò)對(duì)不同節(jié)點(diǎn)處的應(yīng)力強(qiáng)度因子進(jìn)行插值，就可以得到r→0時(shí)，KⅠ，KⅡ和KⅢ的值。

2 裂紋開裂角預(yù)測(cè)的斷裂準(zhǔn)則

利用Erdogan最大周向力斷裂準(zhǔn)則對(duì)復(fù)合型裂紋的初始起裂角進(jìn)行預(yù)測(cè)。關(guān)于最大周向應(yīng)力斷裂準(zhǔn)則的兩個(gè)基本假設(shè)為[1,13]：

1）裂紋沿最大周向應(yīng)力的方向起裂；

2）當(dāng)此方向的最大周向力達(dá)到了臨界值，裂紋失穩(wěn)擴(kuò)展。

由彈性力學(xué)和斷裂力學(xué)的知識(shí)，可以得到Ⅰ-Ⅱ型復(fù)合型裂紋尖端的周向應(yīng)力σθ表達(dá)式為：

由式（4）解得起裂角為：

將位移外推法求得的裂紋尖端應(yīng)力強(qiáng)度因子KⅠ和KⅡ代入式（5），即可求得裂紋的開裂角，結(jié)合有限元模型的裂紋尖端的應(yīng)力分布特點(diǎn)，就可預(yù)測(cè)裂紋開裂的起始方向。

3 雙軸載荷狀態(tài)下裂紋起裂角分析

為了驗(yàn)證文中所描述的裂紋起裂角預(yù)測(cè)方法的有效性和準(zhǔn)確性，對(duì)含孔邊裂紋中心帶孔板建立雙軸載荷狀態(tài)下的有限元分析模型，計(jì)算裂紋尖端的SIF和裂紋的起裂角，并將計(jì)算結(jié)果與文獻(xiàn)[2]的計(jì)算結(jié)果進(jìn)行對(duì)比。雙軸載荷狀態(tài)下含孔邊裂紋中心帶孔板結(jié)構(gòu)模型如圖2所示。

圖2 含孔邊裂紋中心帶孔板的結(jié)構(gòu)模型

選用雙軸載荷下的300 mm（寬W）×400 mm（長(zhǎng)L）×3 mm（厚tp）的中心帶孔板模型，中心孔半徑R=20 mm，孔邊含有一初始裂紋長(zhǎng)度a0=30 mm的平直裂紋。在模型的長(zhǎng)度方向的一端施加σ=100 MPa拉伸載荷和τ=100 MPa剪切載荷，另一端固定，板寬度方向的兩邊自由。板的材料為2024-T3鋁合金，其彈性模量E= 73 084 MPa，泊松比λ= 0.33。

采用ABAQUS 6.10建立中心帶孔板的三維線彈性有限元模型。對(duì)板的建模采用適合彈塑性分析的C3D8R三維實(shí)體單元，這是由于該單元具有位移求解結(jié)果精確、不容易體積自鎖等特點(diǎn)[12]。使用Seam型裂紋模擬孔邊平直穿透型裂紋損傷，全局種子的尺寸設(shè)定為5 mm，在裂紋尖端處進(jìn)行網(wǎng)格密度細(xì)化，設(shè)置局部種子的大小為0.25 mm，得到含孔邊裂紋中心帶孔板的有限元模型，如圖3所示。

圖3 含孔邊裂紋中心帶孔板有限元模型

選取裂紋尖端附近的單元節(jié)點(diǎn)，如圖4所示。節(jié)點(diǎn)1為裂紋尖端，利用節(jié)點(diǎn)2和節(jié)點(diǎn)6、節(jié)點(diǎn)3和節(jié)點(diǎn)7、節(jié)點(diǎn)4和節(jié)點(diǎn)8、以及節(jié)點(diǎn)5和節(jié)點(diǎn)9四組點(diǎn)在x軸和y軸方向的位移差。通過(guò)式（2）可以求得沿裂紋方向距裂紋尖端不同位移處的KⅠ和KⅡ的值，并對(duì)其進(jìn)行直線擬合，如圖5所示，擬合后的方程為：

圖4 裂紋尖端的節(jié)點(diǎn)

圖5 沿裂紋方向距裂紋尖端不同位移處的應(yīng)力強(qiáng)度因子擬合曲線

當(dāng)r=0時(shí)，代入式（6）求得裂紋尖端的應(yīng)力強(qiáng)度因子KⅠ=1394 MPa·mm1/2，KⅡ=865 MPa·mm1/2。為方便與文獻(xiàn)[2]的結(jié)果進(jìn)行對(duì)比，將應(yīng)力強(qiáng)度因子無(wú)量，與文獻(xiàn)[2]中FRANC2D模型的計(jì)算結(jié)果相比，相差1.4%和1.5%。說(shuō)明文中使用的應(yīng)力強(qiáng)度因子求解方法是有效的。

由于板的尺寸tp/W=0.01，屬于薄板的平面應(yīng)力狀態(tài)，將計(jì)算的KⅠ和KⅡ的值代入式（6）可以求得裂紋的起裂角為44.7°。與文獻(xiàn)[2]中FRANC2D模型的計(jì)算結(jié)果相比，相差2.82%。圖6為有限元模型受力后的應(yīng)力分布情況，可以看出，裂紋尖端處應(yīng)力集中，應(yīng)力分布呈8字形，裂紋將沿圖中所標(biāo)的方向開裂，這與文獻(xiàn)[2]中的FRANC2D模型中的裂紋初始起裂方向一致。說(shuō)明本文使用的裂紋起裂角的預(yù)測(cè)方法是有效的。

圖6 有限元模型的應(yīng)力分布

4 結(jié)語(yǔ)

建立了雙軸載荷下含孔邊裂紋中心帶孔板的線彈性三維有限元模型，利用位移外推法求解了裂紋尖端的應(yīng)力強(qiáng)度因子，基于最大周向力斷裂準(zhǔn)則，對(duì)有限元模型復(fù)合型裂紋的起裂方向進(jìn)行了預(yù)測(cè)，預(yù)測(cè)結(jié)果與文獻(xiàn)結(jié)果基本一致，文中的計(jì)算結(jié)果是可靠有效的。

[1] 詹福宇, 楊偉, 王生楠, 等. 基于Abaqus二次開發(fā)的機(jī)身開口損傷容限分析[J]. 機(jī)械強(qiáng)度, 2014, 36(1): 116—120.

[2] 張志賢, 王生楠, 詹福宇. 基于有限元的雙軸載荷狀態(tài)下自動(dòng)裂紋擴(kuò)展分析[J]. 航空計(jì)算技術(shù), 2012, 42(2): 65—68.

[3] 馬家升. 三維復(fù)合型裂紋應(yīng)力強(qiáng)度因子有限元分析[D].哈爾濱: 哈爾濱工程大學(xué), 2008.

[4] 翁劍成, 謝煌生. 復(fù)合型裂紋應(yīng)力強(qiáng)度因子的有限元分析[J]. 龍巖學(xué)院學(xué)報(bào), 2011, 29(2): 44—47.

[5] ERDOGAN F, SIH G C. On the Crack Extension in Plates under Plane Loading and Transverse Shear[J]. Journal of Basic Engineering, 1963, 85(4): 519—527.

[6] SHI G C. Strain-energy-density Factor Applied to Mixed Mode Crack Problems[J]. International Journal of Fracture, 1974, 10(3): 305—321.

[7] THEOCARIS P S, ANDRIANOPOULOS N P. The Mises Elastic-plastic Boundary as the Core Region in Fracture Criteria[J]. Engineering Fracture Mechanics, 1982, 16(2): 425—432.

[8] KHAN S M A, KHRAISHEH M K. A New Criterion for Mixed Mode Fracture Initiation Based on the Crack Tip Plastic Core Region[J]. International Journal of Plasticity, 2004, 20(1) :55—48.

[9] HENLLEN T K, BLACKBURN W S. The Calculation of Stress Intensity Factors for Combined Tensile and Shear Loading[J]. International Journal of Fracture, 1975, 11(4): 605—617.

[10] 楊綠峰, 徐華, 佘振平, 等. Willams單元分析Ⅰ-Ⅱ混合型裂紋應(yīng)力強(qiáng)度因子[J]. 船舶力學(xué), 2014, 18(1/2): 116—123.

[11] 徐慧, 伍曉贊, 程仕平, 等. 復(fù)合裂紋的應(yīng)力強(qiáng)度因子有限元分析[J]. 中南大學(xué)學(xué)報(bào)(工程科學(xué)版), 2007, 38(1): 79—83.

[12] 石亦平, 周玉蓉. ABAQUS有限元分析實(shí)例詳解[M].北京: 機(jī)械工業(yè)出版社, 2006: 151—169.

[13] 李慶芬, 胡勝海, 朱世范. 斷裂力學(xué)及其工程應(yīng)用[M].哈爾濱: 哈爾濱工程大學(xué)出版社, 2005.

[14] 王濤金, 車有金. Ⅰ-Ⅱ型復(fù)合裂紋初始開裂角的三角函數(shù)表達(dá)式及其相應(yīng)斷裂判據(jù)分析[J]. 中國(guó)鍋爐壓力容器安全, 2002, 19(4): 27—30.

[15] 王超, 朱哲明, 劉宏杰. 裂紋起裂角預(yù)測(cè)方法研究[J].四川大學(xué)學(xué)報(bào)(工程科學(xué)版), 2015, 47(2): 28—35.

Finite Element Method for Predicting Crack Initiation Angle of Mixed Mode Crack Based on Seam Crack Model

WANG Yue, MU Zhi-tao, LIU Tao, NIU Yong
(Qingdao Branch of Naval Aeronautical Academy, Qingdao 266041, China)

ObjectiveTo predict the initiation angle of mixed-mode crack under multiaxial stress.MethodsThe Seam crack in commercial finite element software ABAQUS was used to simulate penetrating crack damage according to elastic property of materials. A 3-D finite element model of center-hole with edge cracking plate under multiaxial stress was established. Displacement and extrapolation method and linear fitting method were used to calculate the SIF of crack tip. Initiation angles of cracks were predicted by maximum circumferential stress criterion, which were compared with those obtained by relative papers.ResultsThe maximum error of SIF between the result obtained by finite element models and that of document was within 2%, while the maximum error of initiation angles was within 3%.ConclusionThe results obtained by the finite element method in this paper are in line with document; therefore, the finite element method used for calculating the initiation angle of mixed-mode crack is reliable.

mixed-mode crack; crack initiation angle; seam crack; SIF; maximum circumferential stress criterion

10.7643/ issn.1672-9242.2017.02.021

TJ07

A

1672-9242(2017)02-0105-04

2016-08-24；

2016-09-04

王躍（1989—），男，山東人，博士研究生，主要研究方向?yàn)轱w機(jī)結(jié)構(gòu)腐蝕疲勞。