王建明，劉偉，呂鶴婷

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

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

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

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

1 建模相關理論

1.1 應力強度因子計算

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

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

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

圖1 裂紋尖端坐標系Fig.1 The coordinate system at the crack tip

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

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

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

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

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

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

類似地，

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

1.2 裂紋擴展方向判定

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

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

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

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

1.3 Paris準則

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

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

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

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

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

2 算法流程

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

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

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

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

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

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

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

3 數(shù)值算例

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

3.1 單邊裂紋平板應力強度因子計算

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

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

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

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

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

圖5 Von Mises等效應力云圖Fig.5 The contours of Von Mises stress distribution

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

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

本算例相關模型參數(shù)均來源于文獻［1］，以便于將仿真結果與文獻［1］中的實驗結果進行對照分析。選用有機玻璃PMMA試樣模擬Ⅰ、Ⅱ型復合裂紋的擴展過程，同時考慮孔對裂紋擴展的影響。試樣長度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。由模擬結果可知初始裂紋經(jīng)過8個裂紋擴展步后最終從左邊趨向于中間孔，模擬路徑與文獻［1］實驗結果對比如圖8所示，二者非常接近，同時與文獻［9-10］中的數(shù)值結果吻合較好。圖9為工況ⅠVon Mises等效應力云圖。圖10、圖11為裂紋擴展過程中Ⅰ、Ⅱ型應力強度因子。

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

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

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

圖9 Von Mises等效應力云圖Fig.9 The contours of Von Mises stress distribution

圖10 Ⅰ型應力強度因子Fig.10 Type I stress intensity factor

圖11 Ⅱ型應力強度因子Fig.11 TypeⅡstress intensity factor

3.3 單邊斜裂紋疲勞壽命預測

單邊斜裂紋幾何模型如圖12所示，長度b=100 mm，寬度h=200 mm，裂紋傾斜角度θ=40°，初始長度a0=20 mm，將模型視為平面應變問題，試樣承受脈沖循環(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，平面應變斷裂韌性Kc=1 897 MPa·mm1/2。

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

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

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

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

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

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

4 結束語

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

［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的裂紋擴展仿真軟件及應用［D］.杭州:浙江大學，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］王永偉.結構疲勞裂紋擴展的數(shù)值模擬［D］.大連:大連理工大學，2005:26-42.WANG Yongwei.Numerical simulations of fatigue crack growth of structure［D］.Dalian:Dalian University of Technology，2005:26-42.

［8］程靳，趙樹山.斷裂力學［M］.北京:科學出版社，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.