?

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

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

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

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

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

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

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

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

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

1．1含裂紋區(qū)域的位移場形式

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

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

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

圖1　單元節(jié)點擴充類型

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

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

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

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

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

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

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

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

圖2　裂紋的水平集函數(shù)表征

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

圖3　多裂紋情況的水平集函數(shù)分布

對于裂紋經(jīng)過的單元，計算單元節(jié)點的水平集函數(shù)值φi和ψi，則通過形函數(shù)插值可以得到單元內(nèi)部水平集函數(shù)的表達式:

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

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

圖4　切向水平集更新區(qū)域

1．4間斷單元處理

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

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

圖5　單元擴充類型判斷

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

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

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

式中，n為單元節(jié)點數(shù)，ξi和ηi為單元節(jié)點的等參坐標。

圖6　四節(jié)點單元分割方案

1．5應力強度因子計算

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

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

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

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

1．6程序結構

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

圖7　裂紋擴展分析流程圖

2　方法精度

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

圖8　中心斜裂紋平板

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

2．1網(wǎng)格密度的影響

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

圖9　不同網(wǎng)格密度下的無量綱應力強度因子

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

2．2積分域尺寸的影響

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

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

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

3　裂紋擴展模擬

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

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

圖11　帶孔平板裂紋擴展

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

3．2多裂紋擴展模擬

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

圖12　孔邊裂紋擴展模型

圖13　裂紋擴展路徑

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

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

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

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

4　結論

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

參考文獻:

［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用戶子程序?qū)崿F(xiàn)［J］．上海交通大學學報，2009 (10) : 1644-1648 Xie Hai，F(xiàn)eng 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—)，西北工業(yè)大學博士研究生，主要從事金屬疲勞斷裂數(shù)值算法研究。

收稿日期:2014-09-28

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

文獻標志碼:A

中圖分類號:O346．1