王建明，劉偉，呂鶴婷

(1.山東大學(xué)機械工程學(xué)院，山東濟(jì)南250061;2.山東大學(xué) 高效潔凈機械制造教育部重點實驗室，山東 濟(jì)南250061)

構(gòu)件的固有缺陷或加工損傷會引起初始裂紋，初始微裂紋在交變載荷作用下會發(fā)生擴展最終導(dǎo)致結(jié)構(gòu)的斷裂破壞。故研究裂紋擴展規(guī)律尤其是穩(wěn)態(tài)擴展過程中的裂紋擴展行為具有重要意義。

目前國內(nèi)外針對裂紋擴展的數(shù)值研究方法包括有限元法、邊界元法、無網(wǎng)格法、擴展有限元法。而有限元法因其成熟性在線彈性斷裂問題中被廣泛采用，其關(guān)鍵技術(shù)涉及裂紋尖端應(yīng)力奇異場的處理以及適應(yīng)裂紋擴展的網(wǎng)格重構(gòu)，已有眾多學(xué)者做過相關(guān)研究。Bittencourt等［1］基于Franc2D 軟件，針對線彈性材料采用遞歸剖分技術(shù)進(jìn)行網(wǎng)格重構(gòu)。Xu和Needleman［2］提出了單元間裂紋傳播的方法，其裂紋被假設(shè)沿著單元邊界進(jìn)行擴展，雖然不再需要重新劃分網(wǎng)格但裂紋的擴展方向和擴展步長均依賴初始網(wǎng)格的質(zhì)量，計算精度不高。Bouchard等［3］采用節(jié)點釋放技術(shù)模擬了平面裂紋擴展，并對最大周向應(yīng)力(MCSC)、最小應(yīng)變能密度(MSEDC)、最大能量釋放率(MSERRC)3種裂紋擴展準(zhǔn)則進(jìn)行了分析對比，結(jié)果表明MSEDC精確度不如另外2種高，在有限元代碼方面MCSC較為簡單但裂尖網(wǎng)格要求較密，MSERRC最復(fù)雜但所得結(jié)果較為理想。Liu等［4］基于單元子域離散與應(yīng)變光滑技術(shù)提出了光滑有限元法，創(chuàng)建了5節(jié)點三角形奇異單元以模擬裂紋尖端應(yīng)力場，其計算時間較長。Khoei等［5］提出一種在裂紋擴展過程中進(jìn)行誤差估計并優(yōu)化網(wǎng)格尺寸的自適應(yīng)網(wǎng)格劃分方法，盡管提高了計算精度但算法較繁瑣。王慰軍［6］基于ABAQUS進(jìn)行裂紋擴展二次開發(fā)，采用Deluanay三角剖分算法進(jìn)行網(wǎng)格自動劃分，每一步的裂紋擴展增量Δa的大小不對裂紋擴展方向產(chǎn)生影響。王永偉［7］對拉伸載荷作用下的三維表面疲勞裂紋擴展進(jìn)行了模擬，獲得了推導(dǎo)裂紋形狀變化規(guī)律的有效方法。

現(xiàn)有研究裂紋擴展數(shù)值模型大多針對Ⅰ型裂紋，其裂紋擴展方向可事先預(yù)知。本文針對Ⅰ、Ⅱ復(fù)合型裂紋開展裂紋擴展路徑數(shù)值模擬研究?；贏NSYS進(jìn)行二次開發(fā)，采用位移外推法和最大周向應(yīng)力準(zhǔn)則建立了裂紋擴展路徑模擬及疲勞壽命預(yù)測數(shù)值計算模型，取得了理想的仿真結(jié)果。

1 建模相關(guān)理論

1.1 應(yīng)力強度因子計算

應(yīng)力強度因子是描述裂尖區(qū)域位移場和應(yīng)力場的重要參數(shù)，其數(shù)值計算方法主要有J積分法、虛擬裂紋閉合法、虛擬裂紋擴展法、位移外推法等。本文采用節(jié)點位移外推法計算應(yīng)力強度因子。

建立裂尖坐標(biāo)系如圖1所示，根據(jù)線彈性斷裂理論，復(fù)合型裂紋尖端附近位移場可表示為［8］

式中:u、ν分別為平行于裂紋方向(x方向)、垂直于裂紋方向(y方向)的位移分量;r、θ為以裂尖為坐標(biāo)系原點的極坐標(biāo);E為彈性模量，μ為泊松比;KI、KⅡ分別為Ι型、Ⅱ型裂紋的應(yīng)力強度因子;k為膨脹模量，平面應(yīng)力問題取k=3-μ/1+μ，平面應(yīng)變問題取k=3-4μ。

圖1 裂紋尖端坐標(biāo)系Fig.1 The coordinate system at the crack tip

計算KΙ時，令θ=π，由式(2)可得

由式(3)可以看出，KΙ是裂紋尖端r→0( )處的極限值，無法直接由該式得到。在此通過有限元法，利用其節(jié)點位移進(jìn)行外推，結(jié)合最小二乘法數(shù)據(jù)擬合，可形成KΙ的數(shù)值分析方法。在有限元模型中，由于裂尖處應(yīng)力具有奇異性，故圍繞裂紋尖端設(shè)置1/4奇異單元以提高數(shù)值計算精度。

距裂尖ri處，裂紋表面張開位移νi為

式中:νa和νb分別是裂尖后端位于r=ri處的上裂紋面節(jié)點a和對應(yīng)的下裂紋面節(jié)點b的y方向位移如圖2所示。在裂紋面上規(guī)律選取一系列節(jié)點組成數(shù)據(jù)對:

利用最小二乘法對數(shù)據(jù)對進(jìn)行擬合處理，假定KΙi和ri之間用線性關(guān)系近似:

當(dāng)r=0時，。根據(jù)最小二乘法原理，其最佳擬合應(yīng)滿足:

類似地，

圖2 裂尖后端裂紋表面張開位移示意圖Fig.2 The deformation of crack surfaces behind the crack tip

1.2 裂紋擴展方向判定

由于裂紋方位不對稱或者載荷分布不對稱，裂紋并非簡單地沿著直線向前擴展，而是根據(jù)KΙ與KⅡ的關(guān)系表現(xiàn)為曲線軌跡。目前裂紋擴展方向的判斷主要有3種方法:1)最大周向應(yīng)力準(zhǔn)則;2)最大能量釋放率準(zhǔn)則;3)最小應(yīng)變能密度準(zhǔn)則。本文采用最大周向應(yīng)力準(zhǔn)則計算裂紋擴展方向角θ，即認(rèn)為裂紋沿周向正應(yīng)力σθθ最大的方向擴展。在裂尖極坐標(biāo)系中，裂尖附近的應(yīng)力場可表示為［8］

令 ?σθθ/?θ=0，τrθ=0 可得

求解上述方程可得裂紋擴展角:

Keq與材料斷裂韌性Kc比較可作為裂紋是否發(fā)生失穩(wěn)擴展的判據(jù)。

1.3 Paris準(zhǔn)則

Paris用應(yīng)力強度因子幅定量描述穩(wěn)態(tài)區(qū)疲勞裂紋擴展速率，其表達(dá)式為

式中:C和m是實驗確定的材料常數(shù);ΔKeq是等效應(yīng)力強度因子幅值，可用 ΔKΙ、ΔKⅡ代替方程(13)中的求得。利用Paris公式模擬疲勞裂紋擴展有2種處理方法:

1)在各裂紋擴展步中設(shè)定載荷循環(huán)次數(shù)ΔN，計算相應(yīng)的裂紋擴展增量Δa:

2)在各裂紋擴展步中設(shè)定裂紋擴展增量Δa，計算相應(yīng)的載荷循環(huán)次數(shù)ΔN:

本文采用方法2)，即在各裂紋擴展步中給定裂紋擴展增量Δa，計算相應(yīng)的載荷循環(huán)次數(shù)ΔN。

2 算法流程

本文利用ANSYS平臺，基于上述建模理論，通過APDL編程形成裂紋擴展路徑模擬和疲勞壽命預(yù)測專用分析模塊，其計算流程可分為以下4個步驟:

1)建立含裂紋的參數(shù)化幾何模型;

2)網(wǎng)格自動劃分，其中圍繞裂紋尖端采用1/4奇異單元，其余部分采用二次三角形單元，參數(shù)化施加載荷及位移邊界條件;

3)進(jìn)行有限元計算，得到裂尖后端上下表面對應(yīng)節(jié)點位移，分別利用式 ( 7)、式 ( 8)計算應(yīng)力強度因子KΙ、KΙΙ，利用式(13)計算等效應(yīng)力強度因子Keq;

4)通過設(shè)置多個裂紋擴展步以跟蹤裂紋擴展路徑。在每個裂紋擴展步中設(shè)定裂紋擴展增量Δa，計算相應(yīng)的Keq(或ΔKeq)，利用式 16( )計算該擴展步對應(yīng)的載荷循環(huán)次數(shù)ΔN，當(dāng)Keq≥Kc時停止計算，視為裂紋穩(wěn)態(tài)擴展階段結(jié)束;否則根據(jù)裂紋擴展方向角θ及裂紋擴展增量Δa形成新的裂紋擴展步。如此循環(huán)計算直至滿足Keq≥Kc。裂紋擴展增量△a取值應(yīng)與成反比，Δa越小路徑模擬結(jié)果越精確，但計算時間加長，一般可取初始裂紋長度的 10%～20%［9］。

圖3給出該裂紋擴展路徑模擬和疲勞壽命預(yù)測專用分析模塊的分析流程，其中在計算裂尖位置并形成新的裂紋擴展步模塊中，首先根據(jù)式(12)確定在新擴展步中的裂紋方向，再根據(jù)給定的擴展增量確定新的裂尖位置，根據(jù)新的裂尖位置調(diào)整計算網(wǎng)格，完成相應(yīng)計算。

圖3 裂紋擴展模擬流程圖Fig.3 Flow chart of crack growth path prediction

3 數(shù)值算例

本節(jié)將從應(yīng)力強度因子計算、復(fù)合裂紋擴展路徑模擬和疲勞壽命預(yù)測三方面通過典型算例驗證前文提出的算法及程序設(shè)計的正確性。

3.1 單邊裂紋平板應(yīng)力強度因子計算

單邊裂紋平板幾何模型如圖4所示，裂紋長度a=0.4 mm，平板寬度b=1 mm，高度h=2 mm，軸向載荷σ=1 MPa。材料參數(shù):彈性模量E=30 MPa，泊松比μ=0.3，按平面應(yīng)變問題計算，其Ι型裂紋應(yīng)力強度因子解析解為

式中:修正系數(shù)Cf取:

圖5為Von Mises等效應(yīng)力云圖。裂尖后端5點對應(yīng)力強度因子數(shù)值結(jié)果如表1所示。

將表1數(shù)據(jù)代入式(7)，得到位移外推法應(yīng)力強度因子數(shù)值解為2.370 MPa·mm1/2。由式(17)可得其解析解為2.358 MPa·mm1/2，兩者相對誤差僅為0.5%，說明本文所給出的應(yīng)力強度因子數(shù)值計算方法具有較高的精度。

圖4 單邊裂紋平板幾何模型Fig.4 A plate with single-edge crack under tension

圖5 Von Mises等效應(yīng)力云圖Fig.5 The contours of Von Mises stress distribution

表1 裂尖點對KΙi數(shù)值解Table 1 The KIinumerical solutions at crack tip points

3.2 有機玻璃試樣裂紋擴展路徑模擬

本算例相關(guān)模型參數(shù)均來源于文獻(xiàn)［1］，以便于將仿真結(jié)果與文獻(xiàn)［1］中的實驗結(jié)果進(jìn)行對照分析。選用有機玻璃PMMA試樣模擬Ⅰ、Ⅱ型復(fù)合裂紋的擴展過程，同時考慮孔對裂紋擴展的影響。試樣長度L=20 mm，寬度W=8 mm，3個孔的半徑r=0.25 mm，位置如圖6所示，將試樣在底邊兩點處固定并在頂邊中點施加集中載荷P=1 N。材料參數(shù)為:彈性模量E=3×104MPa，泊松比μ=0.3。裂紋初始長度a=1 mm、距試樣左端距離b=4 mm。

圖7給出了工況Ⅰ下裂紋擴展第2步、第6步時的擴展路徑，分析中取裂紋擴展增量Δa=0.5 mm。由模擬結(jié)果可知初始裂紋經(jīng)過8個裂紋擴展步后最終從左邊趨向于中間孔，模擬路徑與文獻(xiàn)［1］實驗結(jié)果對比如圖8所示，二者非常接近，同時與文獻(xiàn)［9-10］中的數(shù)值結(jié)果吻合較好。圖9為工況ⅠVon Mises等效應(yīng)力云圖。圖10、圖11為裂紋擴展過程中Ⅰ、Ⅱ型應(yīng)力強度因子。

圖6PMMA試樣幾何模型Fig.6 Geometric model of the PMMA specimen

圖7PMMA試樣裂紋擴展路徑Fig.7 Crack growth path of the PMMA specimen

圖8 裂紋擴展路徑模擬結(jié)果與實驗結(jié)果對比Fig.8 Comparison of the crack growth path between numerical results and experimental results

圖9 Von Mises等效應(yīng)力云圖Fig.9 The contours of Von Mises stress distribution

圖10 Ⅰ型應(yīng)力強度因子Fig.10 Type I stress intensity factor

圖11 Ⅱ型應(yīng)力強度因子Fig.11 TypeⅡstress intensity factor

3.3 單邊斜裂紋疲勞壽命預(yù)測

單邊斜裂紋幾何模型如圖12所示，長度b=100 mm，寬度h=200 mm，裂紋傾斜角度θ=40°，初始長度a0=20 mm，將模型視為平面應(yīng)變問題，試樣承受脈沖循環(huán)載荷，Δσ=40 MPa(σmax=40 MPa，σmin=0)，材料參數(shù)為:E=74 000 MPa，泊松比μ=0.3。Paris模型參數(shù)C=2.087 136 × 10-13，m=3.32，平面應(yīng)變斷裂韌性Kc=1 897 MPa·mm1/2。

圖13所示為應(yīng)力強度因子隨裂紋長度的變化曲線圖，各步裂紋擴展增量取Δa=1 mm，從圖13可以看出應(yīng)力強度因子KⅡ經(jīng)過微幅上升之后迅速變?yōu)?，之后裂紋由Ι-Ⅱ復(fù)合型變?yōu)榧儲┬汀Ａ鸭y擴展軌跡如圖14所示，裂紋開裂時其方向發(fā)生突變，隨后保持直線擴展直至試樣失穩(wěn)斷裂。

圖12 單邊斜裂紋平板承受循環(huán)載荷Fig.12 A plate with edge-angled crack under a cyclic tension

圖13 應(yīng)力強度因子隨裂紋長度變化曲線Fig.13 Stress intensity facors vs.crack length

圖14 單邊斜裂紋擴展路徑Fig.14 Crack growth path of the edge-angled crack problem

圖15為裂紋擴展的a-N曲線，當(dāng)裂紋等效應(yīng)力強度因子達(dá)到臨界值Kc時，其裂紋長度為61.4 mm，載荷循環(huán)次數(shù)為130 016次，與文獻(xiàn)［11］采用無網(wǎng)格法模擬得到的結(jié)果:裂紋長度61.3 mm，載荷循環(huán)次數(shù)131 411次比較，兩者基本吻合，驗證了本文所得方法的正確性與有效性。

圖15 單邊斜裂紋疲勞壽命曲線Fig.15 Fatigue life curves of the edge-angled crack problem

4 結(jié)束語

本文基于ANSYS編制了裂紋擴展路徑模擬和疲勞壽命預(yù)測專用分析模塊，采用位移外推法數(shù)值計算應(yīng)力強度因子，并結(jié)合最大周向應(yīng)力準(zhǔn)則和Paris準(zhǔn)則形成模擬復(fù)合型裂紋擴展路徑及疲勞壽命預(yù)計的有效算法。該分析模塊由APDL參數(shù)化編程，在裂紋擴展的每個計算步中采用全自動網(wǎng)格劃分，其分析過程簡便高效，為裂紋擴展模擬和疲勞壽命預(yù)測提供有利工具。

［1］BITTENCOURT T N，WAWRZYNEK P A，INGRAFFEA A R，et al.Quasi-automatic simulation of crack propagation for 2D LEFM problems［J］.Engineering Fracture Mechanics，1996，52(2):321-334.

［2］XU X P，NEEDLEMAN A.Numerical simulations of fast crack growth in brittle solids［J］.Journal of the Mechanics and Physics of Solids，1994，42(9):1397-1434.

［3］BOUCHARD P O，BAY F，CHASTEL Y.Numerical modelling of crack propagation:automatic remeshing and comparison of different criteria［J］.Computer Methods in Applied Mechanics and Engineering，2003，192(35/36):3887-3908.

［4］LIU G R，NOURBAKHSHNIA N，CHEN L，et al.A novel general formulation for singular stress field using the ESFEM method for the analysis of mixed-mode cracks［J］.International Journal of Computational Methods，2010，7(1):191-214.

［5］KHOEI A R，AZADI H，MOSLEMI H.Modeling of crack propagation via an automatic adaptive mesh refinement based on modified superconvergent patch recovery technique［J］.Engineering Fracture Mechanics，2008，75(10):2921-2945.

［6］王慰軍.基于ABAQUS的裂紋擴展仿真軟件及應(yīng)用［D］.杭州:浙江大學(xué)，2006:17-44.WANG Weijun.The program and its applications for simulation of crack propagation based on ABAQUS［D］.Hangzhou:Zhejiang University，2006:17-44.

［7］王永偉.結(jié)構(gòu)疲勞裂紋擴展的數(shù)值模擬［D］.大連:大連理工大學(xué)，2005:26-42.WANG Yongwei.Numerical simulations of fatigue crack growth of structure［D］.Dalian:Dalian University of Technology，2005:26-42.

［8］程靳，趙樹山.斷裂力學(xué)［M］.北京:科學(xué)出版社，2006:21-45.

［9］NOURBAKHSHNIA N，LIU G R.A quasi-static crack growth simulation based on the singular ES-FEM［J］.International Journal for Numerical Methods in Engineering，2011，88(5):473-492.

［10］XUAN H N，LIU G R，BORDAS S，et al.An adaptive singular ES-FEM for mechanics problems with singular field of arbitrary order［J］.Computer Methods in Applied Mechanics and Engineering，2013，253:252-273.

［11］DUFLOT M，DANG H N.A meshless method with enriched weight functions for fatigue crack growth［J］.International Journal for Numerical Methods in Engineering，2004，59(14):1945-1961.