基于擴展有限元的多裂紋擴展分析

?

基于擴展有限元的多裂紋擴展分析

束一秀，李亞智，姜薇，賈雨軒

(西北工業大學航空學院，陜西西安710072)

摘要:研究基于擴展有限元的多裂紋擴展分析方法以及ABAQUS環境下的程序實現。在位移函數中增加擴充項以描述裂紋周圍的不連續位移場。初始裂紋以點集或方程的形式給出，并使用水平集函數將多裂紋信息離散到單元節點上，水平集函數還用來追蹤裂紋擴展路徑。使用局部水平集更新方法減小了計算規模，改進了裂尖單元的判斷準則，借助商業軟件Tecplot軟件實現了裂紋擴展的動態顯示功能。使用交互積分法計算混合模式下的應力強度因子，用最大周向拉應力準則判斷裂紋擴展方向。探討了網格密度和積分域尺寸對方法的影響，數值算例表明擴展有限元方法能夠模擬任意形狀的裂紋，并且能夠反映多裂紋擴展的規律。

關鍵詞:ABAQUS;裂尖;疲勞裂紋擴展;有限元;流程圖;矩陣代數;剛度矩陣;應力強度因子;交互積分;水平集法;擴展有限元法

研究疲勞載荷下的裂紋擴展規律對于結構安全性評價非常重要。國內外學者已經提出了很多的方法來計算裂紋擴展問題，比如邊界元法［1］、無網格方法［2］、有限元網格更新法［3-4］等，這些方法在求解裂紋擴展問題時都有良好的表現。基于有限元方法的網格更新法雖然適應性強，但是需要隨著裂紋的擴展不斷更新裂紋周圍的網格，工作量大，計算效率較低。擴展有限元方法［5-6］基于傳統有限元方法，并且克服了網格更新方法的缺點，通過擴充位移項來描述不連續的位移場，使裂紋獨立于網格存在，從而能夠模擬任意形狀的裂紋。此外，由于擴展有限元方法僅僅在有限元法的基礎上擴充了位移項，可以通過對現有商業有限元軟件進行二次開發實現其功能，一定程度上減少了工作量。

基于商業軟件的擴展有限元功能開發是普遍采用的一種方法，其優勢是可以省去編寫求解平衡方程程序的過程。ABAQUS從6．9版本開始內置了擴展有限元方法，但是只考慮了單元被裂紋完全貫穿的情況，而且在使用上存在諸多限制［7］。方修君等［8］在ABAQUS環境下開發了擴展有限元功能，并對三點彎曲梁的開裂過程進行了模擬，但是文章只研究了Ⅰ型開裂問題，且沒有考慮裂尖位于單元內部的情況。謝海等［9］解決了ABAQUS環境下求解混合型裂紋的擴展有限元方法，但是沒有實現裂紋擴展問題的模擬。Giner等［10］使用用戶子程序功將XFEM方法嵌入到ABAQUS平臺中，實現了對多裂紋擴展的模擬。

本文推導并編程實現了基于XFEM的多裂紋擴展模擬，改進了水平集更新以及裂尖單元判斷算法，解決了可視化后處理問題。

1　模擬裂紋擴展的擴展有限元方法

1．1含裂紋區域的位移場形式

在XFEM中，對于包含裂紋的區域，通過添加擴充項描述包含間斷的位移場［6］

式中，ui為普通節點自由度，n為所有節點，nJ是屬于被單元完全貫穿的單元的節點集，稱為階躍擴充節點，nT是屬于被裂尖刺中的單元的節點集，稱為裂尖擴充節點，擴充節點和裂紋的位置關系如圖1所示。ai是階躍擴充節點的附加自由度(擴充函數H(x)在裂紋一側取+ 1，另一側取－1)是裂尖擴充節點的擴充自由度，其擴充函數Ψα(x)的形式為

式中，r和θ為裂尖局部坐標參數。

圖1　單元節點擴充類型

1．2剛度矩陣和幾何矩陣

根據(1)式定義的位移場，單元剛度矩陣K的形式為

式中，B為單元幾何矩陣，對于普通節點

對于階躍擴充節點，幾何矩陣附加項Bi，jump=(H(x) －H(xi) ) Bi，c，對于裂尖擴充節點，幾何矩陣附加項為

1．3裂紋擴展的水平集表征

使用擴展有限元方法模擬斷裂問題，不需要直接對裂紋進行建模，通常使用水平集方法［11］將單元內部的裂紋信息以符號距離的形式離散到節點上，然后在后處理中通過形函數插值還原裂紋信息。

區域內的一條裂紋的幾何形狀由切向水平集函數φ和法向水平集函數ψ確定，如圖2所示裂紋，它的形狀可以描述為:

對于內嵌裂紋(2個端點)，需要使用2個法向水平集函數ψ1和ψ2來確定裂尖的位置。

圖2　裂紋的水平集函數表征

對于多裂紋情況，為了節省計算空間，為區域內的每條裂紋分配一個影響區域，對處在區域內的單元的水平集函數值進行初始化。如果裂紋位置比較靠近，其影響區域相互疊加，則該區域內的單元節點需要分別計算相對于兩條或更多條裂紋的水平集函數值，如圖3所示。

圖3　多裂紋情況的水平集函數分布

對于裂紋經過的單元，計算單元節點的水平集函數值φi和ψi，則通過形函數插值可以得到單元內部水平集函數的表達式:

式中，n為單元節點數。裂紋擴展后，對處于裂紋影響區域內的單元的水平集函數進行更新，即可通過插值還原擴展后的裂紋幾何形狀。水平集的更新可以通過求解(8)式得到［7］

式中，v為裂尖擴展速率，Δt為時間增量，Ωupdate為切向水平集函數的更新域，其范圍如圖4所示，取為裂紋擴展方向的一個矩形區域，其垂直于裂紋增量方向的邊長為裂尖單元邊長的4倍，沿著裂紋增量方向的邊長為裂紋增量的2～3倍，這樣處理可以減少更新切向水平集函數φ的計算量。

圖4　切向水平集更新區域

1．4間斷單元處理

使用單元節點的水平集值可以方便地判斷間斷單元的類型:對于完全貫穿單元，滿足φminφmax＞0，ψmax＜0;對于裂尖單元，滿足ψminψmax＜0，ψminψmax＜0。其中，φmin、φmax、ψmin、ψmax分別為單元節點的最小和最大水平集函數值。

實際問題中，簡單的使用上述方法并不能準確判斷單元類型，注意圖5所示的情況，圖5中的灰色單元滿足裂尖單元判據，但不是裂尖單元，單元類型判斷的錯誤會賦予節點錯誤的附加自由度，從而影響模擬的精度，需要增加額外的判斷條件來提高方法的穩定性。

圖5　單元擴充類型判斷

對于初步判定為裂尖擴充形式的單元，判斷裂尖點與單元的位置關系，若包含裂尖，則可以判定為裂尖單元。所以裂尖單元的判定標準可以定義為(10)式。

式中，xtip為裂尖位置，Ωele為單元區域。

由于包含裂紋信息的單元內部不連續，使用Gauss積分時，無法得到精確解，一般采用劃分積分子區域的方法進行積分。單元與裂紋可能存在4種交匯關系，針對不同的交匯關系對單元進行分割，然后分別對單元的各子區域積分，如圖6所示。裂紋與單元的交點可以由(11)式得到

式中，n為單元節點數，ξi和ηi為單元節點的等參坐標。

圖6　四節點單元分割方案

1．5應力強度因子計算

基于J積分的交互積分方法［12］被廣泛用于混合型裂紋的應力強度因子的求解。交互積分在真實場的基礎上引入輔助場，輔助場可以是任意滿足變形方程和本構方程的應力-變形場，在求解斷裂問題時取為裂尖奇異場的形式。則交互積分的表達式為

式中，A為積分區域，ui和分別為真實位移場和輔助位移場，σij和分別為真實和輔助應力場，W ==εij。q為積分區域的權重函數，在積分域內取1，在積分域外取0。

交互積分和應力強度因子的關系如(13)式所示

對于復合型裂紋，裂紋擴展方向可以由最大周向拉應力準則計算

1．6程序結構

有別于傳統的在ABAQUS中嵌入XFEM功能的方案，本文以模塊化的形式編寫裂紋擴展分析程序:使用FORTRAN分別編寫了模型創建、裂紋初始化、水平集更新、斷裂參數計算以及可視化處理(借助商業軟件Tecplot，實現了裂紋擴展過程的動態顯示)等模塊并留出接口，方便不同裂紋擴展問題的創建以及后續功能的開發。圖7為程序的詳細流程圖。

圖7　裂紋擴展分析流程圖

2　方法精度

為驗證上述方法的計算精度，計算中心斜裂紋平板在不同裂紋傾角下的應力強度因子。模型構形如圖8所示。

圖8　中心斜裂紋平板

受垂直方向均勻拉伸作用，裂紋長度與板寬滿足a/W = 0．8，板長H = 2W。裂紋與載荷夾角為θ，取其值為0°～90°，計算裂尖的無量綱應力強度因子并與應力強度因子手冊結果［13］進行比較。

2．1網格密度的影響

為了研究網格密度對結果的影響，考慮疏密程度不同的網格劃分，網格類型為均勻的四邊形網格，單元數目分別為1 250(粗網格)、2 450(普通網格) 和5 000(細網格)。圖9為不同單元數目下的無量綱應力強度因子(K/σ)與角度θ的關系曲線，由于對稱性，本文只給出一個裂尖的結果。

圖9　不同網格密度下的無量綱應力強度因子

從結果可以看出，當網格較粗時(網格數為1 250)，角度80°～85°時KⅠ值計算誤差較大，這是因為此時裂尖離模型邊緣很近，從裂尖到模型邊緣之間只存在1～2層單元，計算交互積分時積分域內單元數過少，導致了計算精度不夠，可以通過加密模型邊緣網格的方法解決。除此之外，應力強度因子的計算解與理論解吻合較好，表明方法受網格密度影響較小。

2．2積分域尺寸的影響

為了考察交互積分的積分域半徑對計算精度的影響，對網格數5 000的模型分別取積分半徑為2～4倍的特征長度(裂尖單元面積的平方根)，以取4倍特征長度計算的結果為基準，計算它們的相對誤差，結果如圖10所示，其中R為積分域半徑，H為單元特征長度。

圖10　不同積分半徑下的無量綱應力強度因子誤差

從結果可以看出，取3種積分半徑計算應力強度因子的結果相對誤差非常小，同時可以說明，取2倍的特征長度已經可以取得不錯的計算精度。

3　裂紋擴展模擬

3．1帶孔板邊緣裂紋擴展模擬

對圖11所示的邊緣裂紋模型進行模擬，圖中長度單位均為mm，預制初始裂紋為略高于中心線的邊緣裂紋，其尺寸a0= 10 mm。模型材料為7075-T6，彈性模量E = 7．17×104N/mm2，泊松比υ = 0．33，集中載荷P = 20 kN，用于對比的數值和試驗結果取自文獻［16］。

圖11　帶孔平板裂紋擴展

采用控制裂紋增量的方式計算裂紋擴展，取Δa =3 mm，當裂紋穿入圓孔時停止計算，計算結果如圖11b)所示，圖11c)為后處理云圖結果。我們將計算得到的裂紋路徑點按照坐標繪制出來，并與文獻中的計算路徑作了對比，兩者結果與試驗吻合較好，本文結果略優于文獻中的模擬結果。

3．2多裂紋擴展模擬

考慮圖12所示的雙邊裂紋擴展模型，兩孔關于模型中心對稱分布，初始裂紋對于模型中心對稱分布于模型兩側，長度為1 mm。模擬裂紋擴展路徑并與更新網格方法(Bouchard等［6］和Khoei［7］)進行對比。矩形板的彈性模量E= 2×105N/mm2，泊松比v =0．3。對模型下緣進行固支，對上緣施加等位移加載，同時約束上緣節點水平方向的自由度，使它們在水平方向不能收縮。我們在程序中可以使用裂紋增量或者應力循環數作為計算裂紋的控制參數，對于此算例，由于對稱性，使用控制裂紋增量的形式模擬裂紋擴展。取Δa=1 mm，擴展路徑如圖13所示。

圖12　孔邊裂紋擴展模型

圖13　裂紋擴展路徑

從裂紋擴展路徑可以看出，我們計算的結果與Khoei等人計算的結果吻合較好并且2條裂紋的擴展路徑對稱分布，而Bouchard等人計算的結果右側裂紋出現了較大的拐點，三者裂紋擴展的走勢基本相符。裂紋剛開始朝孔的方向擴展，然后由于2條裂紋間的相互影響，裂紋開始偏離圓孔并相向擴展，最終裂紋在模型中部交錯和匯合。

為了研究裂紋與孔的距離對裂紋擴展路徑的影響，我們將初始裂紋的位置向孔靠近(距離模型端面分別為3．15 mm、3．45 mm)，對這2種情況進行模擬并與原來的情況進行對比，由于使用XFEM模擬裂紋不依賴與網格，因此相比于更新網格法，可以方便地改變初始裂紋的位置。模擬結果如圖14所示。

圖14　裂紋與模型邊緣距離不同時的擴展路徑

從圖中可以看出，當裂紋離孔較近時會直接向孔靠近并最終刺入孔中;當裂紋離孔較遠時，會呈現圖13中討論的裂紋路徑形態;當裂紋與孔的距離為某一值時(圖14中實線的裂紋)，裂紋在擴展到孔邊緣時，會受對面裂紋的影響而偏離孔，并最終形成同一水平面上相向擴展的裂紋。以上關于初始裂紋位置對擴展路徑影響的模擬和討論雖然沒有試驗對比或文獻可資參考，但是其擴展規律的合理性是明顯的。

4　結論

本文基于XFEM基本理論，利用商業軟件ABAQUS和Tecplot的程序接口，編寫了用于多裂紋擴展模擬的Fortran程序，在傳統XFEM方法上進行了部分的優化: 1．使用局部水平集更新算法，減小了計算量; 2．改進了裂尖單元的判斷準則，提高了方法的穩定性。3．編寫了可視化處理模塊，彌補了ABAQUS不能顯示用戶單元的缺陷，方便觀測含裂紋體的應力分布狀態。數值算例證明了本文的程序計算結果精度較高且適應性強，與更新網格方法的對比證明使用XFEM可以靈活地調整裂紋位置，方便研究不同位置的裂紋擴展規律。

參考文獻:

［1］Portela A，Aliabadi M H，Rooke D P．The Dual Boundary Element Method: Effective Implementation for Crack Problems［J］．International Journal for Numerical Methods in Engineering，1992，33(6) : 1269-1287

［2］Belytschko T，Gu L，Lu Y Y．Fracture and Crack Growth by Element Free Galerkin Methods［J］．Modelling and Simulation in Materials Science and Engineering，1994，2(3A) : 519-534

［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) : 3887-3908

［4］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

［5］Belytschko T，Black T．Elastic Crack Growth in Finite Elements with Minimal Remeshing［J］．International Journal for Numerical Methods in Engineering，1999，45(5) : 601-620

［6］Dolbow J，Belytschko T．A Finite Element Method for Crack Growth without Remeshing［J］．Int J Numer Meth Engng，1999，46: 131-150

［7］Pais M J．Variable Amplitude Fatigue Analysis Using Surrogate Models and Exact XFEM Reanalysis［D］．University of Florida，2011

［8］方修君，金峰．基于ABAQUS平臺的擴展有限元法［J］．工程力學，2007，24(7) : 6-10 Fang Xiujun，Jin Feng．Extended Finite Element Method Based on ABAQUS［J］．Engineering Mechanics，2007，24(7) : 6-10 (in Chinese)

［9］謝海，馮淼林．擴展有限元的ABAQUS用戶子程序實現［J］．上海交通大學學報，2009 (10) : 1644-1648 Xie Hai，Feng Miaolin．Implementation of Extended Finite Element Method with ABAQUS User Subroutines［J］．Journal of Shanghai Jiaotong University，2009 (10) : 1644-1648 (in Chinese)

［10］Giner E，Sukumar N，Tarancón J E，et al．An Abaqus Implementation of the Extended Finite Element Method［J］．Engineering Fracture Mechanics，2009，76(3) : 347-368

［11］Stolarska M，Chopp D L，Mo?s N，et al．Modelling Crack Growth by Level Sets in the Extended Finite Element Method［J］．International Journal for Numerical Methods in Engineering，2001，51(8) : 943-960

［12］Rao B N，Rahman S．An Interaction Integral Method for Analysis of Cracks in Orthotropic Functionally Graded Materials［J］．Computational Mechanics，2003，32(1/2) : 40-51

［13］Sih G C．Handbook of Stress-Intensity Factors［M］．Lehigh University，Institute of Fracture and Solid Mechanics，1973

Analyzing Multiple Crack Propagation Using Extended Finite Element Method (X-FEM)

Shu Yixiu，Li Yazhi，Jiang Wei，Jia Yuxuan

(College of Aeronautics，Northwestern Polytechnical University，Xi＇an 710072，China)

Abstract:A numerical method and its implementation in ABAQUS using X-FEM to analyze multiple crack propagation was studied．The displacement discontinuity was approximated by appending enriched items to standard field．We give the initial cracks in the form of point-set or curve functions and dispersed them to the element nodes using Level Set Method (LSM)．LSM was also used to track crack propagation paths．The local level set update technique was used to reduce the computing scale and the Judgmental Approach has been improved．The display of the process of crack propagation was realized using the commercial software Tecplot．Interaction integral technique was used to calculate mixed mode stress intensity factors．The maximal circumferential stress criterion was used to calculate the kinking angles of the propagating cracks．Influence of mesh and integral zone on simulation results were taken into consideration．Numerical examples were presented to demonstrate the benefits of the proposed implementation．

Key words:ABAQUS，crack tips，fatigue crack propagation，finite element method，flowcharting，matrix algebra，stiffness matrix，stress intensity factors; interaction integral，level set method，X-FEM (Extended Finite Element Method)

作者簡介:束一秀(1988—)，西北工業大學博士研究生，主要從事金屬疲勞斷裂數值算法研究。

收稿日期:2014-09-28

文章編號:1000-2758(2015) 02-0197-07

文獻標志碼:A

中圖分類號:O346．1