基于擴展有限元法的裂紋擴展數值模擬及程序實現

崔 韋，熊 欣，喻濺鑒，王建軍

(1.中國直升機設計研究所，江西 景德鎮 333001；2 北京航空航天大學 能源與動力工程學院，北京 100191)

基于擴展有限元法的裂紋擴展數值模擬及程序實現

崔 韋1，熊 欣1，喻濺鑒1，王建軍2

(1.中國直升機設計研究所，江西 景德鎮 333001；2 北京航空航天大學 能源與動力工程學院，北京 100191)

常規有限元法在解決裂紋擴展問題時需要重構網格而使其應用受到局限。擴展有限元法(XFEM)使用擴充形函數描述不連續的位移場，具有在模擬裂紋擴展時無需更新網格模型的顯著優點。基于Matlab語言實現了平面問題裂紋擴展模擬程序的開發，通過解析解算例和對某鈦合金裂紋擴展試驗的模擬驗證了程序的精度。

擴展有限元法；裂紋擴展；數值模擬；Matlab

0 引言

疲勞引起的裂紋萌生和擴展是直升機結構失效的主要原因，近年來損傷容限設計廣泛應用于直升機結構設計中[1]。常規有限元法(FEM)在模擬裂紋擴展、移動界面等不連續問題時，需要進行計算網格的重構，自動化的網格重構是困難而耗費計算資源的：一方面,全自動的網格重構無法保證全部劃分結構化網格，目前計算力學界的主要策略為子模型內局部采用非結構化網格[2-4]；另一方面,網格重構將耗費較多計算資源。為了達到在模擬裂紋擴展時不進行網格重構，使用擴充形函數描述不連續位移場的擴展有限元法(eXtended Finite Element Method，簡稱XFEM)應運而生[5,6]。本文對擴展有限元法解決裂紋擴展問題的基本原理進行了闡述，基于Matlab語言實現了平面問題裂紋擴展模擬程序的開發，通過算例驗證了程序求解應力強度因子和模擬裂紋擴展路徑的精度，最后將程序用于某國產鈦合金裂紋擴展試驗的數值模擬。

1 擴展有限元法的原理

結構內任意形狀的單一裂紋如圖1所示，有限元網格的所有節點的集合記為S，所有圍繞裂紋尖端的單元的節點記為SC，所有被裂紋貫穿的單元的節點(非SC)記為SH，其余節點為常規有限元節點。單元相應分為四類：(1)裂尖單元，單元節點全部由SC組成的單元；(2)貫穿單元，單元節點全部由SH組成的單元；(3)混合單元，單元節點由SC和/或SC與常規有限元節點組成；(4)常規有限元單元。其中，涉及到SC和SC節點的單元需要使用擴充形函數，以下分別介紹各類單元的擴充形函數。

圖1 二維裂紋問題中的XFEM擴充節點

由常規有限元形函數、階躍擴充形函數和漸進位移場擴充形函數，可以將包含二維平面不連續位移場表示為如下形式：

式中，I表示所有節點的集合，J為被貫穿單元的節點集合，K為裂尖單元所有節點的集合；NI，NJ和NK分別代表以上三種節點的形函數；uI代表常規有限元位移自由度，aJ和bK為節點附加自由度，分別代表對貫穿單元和裂尖單元位移場的擴充，成為Heaviside自由度和裂尖自由度；H(x)為階躍函數，Φα(x)為漸進位移場函數。

2 平面問題裂紋擴展的程序實現

作者以商用數學應用程序Matlab為平臺，使用m語言編寫了二維問題的XFEM程序。該程序可以實現應力強度因子計算和裂紋擴展模擬的功能。在不引入擴充形函數時，該程序可以退化為常規有限元程序進行求解。程序進行裂紋擴展模擬的運行流程如圖2所示，通過斷裂韌度和門檻值作為循環結束判斷標準的過程，在XFEM程序中進行實現。程序除主函數外，各子函數按功能可以分為初始化、前處理、求解器和后處理四類。在程序網格劃分方面，由于XFEM描述裂紋時網格不依賴裂紋邊界，而且研究的矩形計算域拓撲結構簡單，因此，程序使用了結構化網格進行劃分。結構化網格具有網格生成速度快、質量高和數據結構簡單的優點。

2.1 算例：斜裂紋矩形板受均勻單向拉伸

為了驗證自編程序對I、II復合型裂紋應力強度因子的計算精度，本算例模擬了斜裂紋矩形板受均勻單向拉伸問題。圖3(a)所示矩形平板長H=5m，寬W=4m，裂紋長度a=1m，裂紋與載荷方向夾角度θ=67.5°，在上、下邊施加均布拉力σ=1Pa。彈性模量E200GPa，泊松比v=0.3。程序網格如圖3 (b)所示，約束右上角節點UX自由度，約束右下角節點UX、UY自由度，模擬上下邊可轉動單軸拉伸邊界條件。計算了不同網格密度下2個裂紋尖端的應力強度因子KⅠ、KⅡ，計算結果如表1所示。

表1 斜裂紋矩形平板(θ=67.5°)單軸拉伸問題計算結果(SIF單位

圖2 XFEM程序流程圖

由表1可見，網格密度在160×200時計算結果具有較高精度，KⅠ、KⅡ與理論解的相對誤差僅為1%左右，程序適合平面復合型問題的求解。將真實變形量放大1e10倍后繪制變形后網格如圖 4(a)、(b)所示。

2.2 算例：三點彎曲裂紋擴展

程序對裂紋擴展模擬的有效性通過簡支裂紋梁的三點彎曲算例進行驗證，如圖 5示意。簡支梁長L=200mm,高h=40mm，寬度b=2mm，梁跨度中點布置一長a=10mm的初始裂紋。使用平面應力單元離散，網格規模為201×40。圖示集中力載荷F=50N。Paris公式系數C=1.5e-11，指數m=3，裂紋擴展判據使用最大周向應力準則。在裂紋擴展計算中，假設斷裂韌度KIC無限大，設定裂紋增量Δa=2mm，裂紋從10mm擴展至35mm分為5個裂紋擴展增量步進行計算。

計算結果表明，裂紋開裂角度θc與純I型三點彎曲的理論值0°差距較小，均在1°左右。裂紋擴展路徑如圖 6所示，可見程序較好地模擬了純I型問題的裂紋擴展問題。

圖3 斜裂紋矩形板受均勻拉伸

圖4 變形后的網格(圖示為真實變形放大1e10倍)：

圖5 三點彎曲裂紋擴展問題

3 某國產鈦合金裂紋擴展性能試驗數值模擬

為進一步驗證本文擴展有限元程序的有效性與精度，對預研課題[9]某型國產鈦合金裂紋擴展性能試驗進行了數值模擬。鈦合金M(T)試件及試驗現場照片如圖7所示。

取試驗中兩件應力比R=0.5試件的結果，對Paris公式進行擬合,如表2所示。

圖6 三點彎曲梁的裂紋擴展路徑

Paris公式RCnda/dN=C(ΔK)n0.57.516×10-82.906

圖7 某國產鈦合金M(T)試樣及試驗現場照片

使用擴展有限元程序對該裂紋擴展過程進行了模擬，目的在于將試驗獲得的Paris公式參數作為輸入，驗證程序對裂紋擴展長度的計算精度。有限元的數值模擬中，計算模型尺寸取試驗件平均值，寬W=75.12mm，厚B=5.16mm,計算中初始半裂紋長取a=5.28mm，計算載荷加載均布應力σ=103.2MPa，計算網格為340*150。計算循環數與試驗數據點對應，計算裂紋長度-循環數曲線如圖8所示。

圖8 鈦合金試件裂紋擴展過程XFEM與

結果表明XFEM與試驗裂紋長度最大差異在4%以內。可見擴展有限元程序在模擬裂紋擴展方面具有較好的精度。

4 結語

本文介紹了擴展有限元(XFEM)求解包含不連續位移場斷裂力學問題的理論基礎，基于Matlab編寫了求解二維問題的XFEM程序。該程序可進行平面問題的裂紋擴展模擬。通過算例驗證了XFEM方法求解二維問題I型和I、II復合型應力強度因子具有較好的精度；以三點彎曲算例驗證了XFEM模擬裂紋擴展路徑具有較高精度。通過對某國產鈦合金裂紋擴展試驗的數值模擬，驗證了XFEM程序計算裂紋擴展具有較好的精度。XFEM在模擬裂紋擴展時不需重構網格，該獨特優勢使其具有較好的工程應用前景。

[1] 穆志韜, 曾本銀. 直升機結構疲勞[M]. 北京：國防工業出版社, 2009.

[2] 賈學明, 王啟智. 三維斷裂分析軟件FRANC3D[J]. 計算力學學報, 2004, 21(6): 764-768.

[3] Hou J, Goldstraw M, Maan S, et al. An evaluation of 3D crack growth using ZENCRACK[R]. Defence Science and Technology Organisation Victoria (Australia) Aeronautical and Maritime Research Lab, 2011.

[4] Richard H A, Fulland M, Sch?llmann M, et al. Simulation of fatigue crack growth using ADAPCRACK3D[C]. Fatigue. 2002: 1405-1412.

[5] Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing [J]. International Journal for Numerical Methods in Engineering. 1999(86): 1141-1151.

[6] Belytschko T, Gracie R, Ventura G. A review of extended/generalized finite element methods for material modeling[J]. Modelling and Simulation in Materials Science and Engineering, 2009, 17(4): 043001.

[7] sher S, Fedkiw R P. Level Set Methods and Dynamic Implicit Surfaces [M]. Springer, 2003.

[8] 仇仲翼. 應力強度因子手冊[M]. 北京: 科學出版社, 1993.

[9] 鈦合金XX材料損傷容限性能研究[S]. 直升機設計研究所,2015.

Crack Propagation Numerical Simulation and Implementation by eXtended Finite Element Method

CUI Wei1, XIONG Xin1, YU Jianjian1, WANG Jianjun2

(1.China Helicopter Research and Development Institute, Jingdezhen 333001, China;2.School of Energy and Power Engineering, Beihang University, Beijing 100191, China)

Due to the demands of remeshing tasks, standard Finite element method is confined in crack growth simulations. Based upon research on enriched shape function of discontinues describing, a framework of crack propagation simulation was proposed based upon extended finite element method (XFEM). The theories of crack propagation simulation in XFEM were explained. Then a two-dimensional XFEM implementation within Matlab was presented. The accuracy of XFEM implementation was verified by static fracture mechanics examples and experiment of titanium alloy crack growth test.

eXtended Finite Element Method； crack growth; numerical simulation; Matlab

2016-02-19 作者簡介：崔 韋(1985-)，男，天津人，博士，工程師，主要研究方向：直升機結構疲勞與損傷容限。

1673-1220(2016)03-013-05

V215.5+2；TB115.1

A