基于光滑有限元的含裂紋復(fù)合材料的虛擬裂紋閉合法＊

周立明，孟廣偉，王 暉，李 鋒?，郭學(xué)東

（1．吉林大學(xué) 機械科學(xué)與工程學(xué)院，吉林 長春 130025；2．吉林大學(xué) 交通學(xué)院，吉林 長春 130025）

復(fù)合材料由于具有比強度高、比模量大、耐疲勞及抗破損等特點，在航空航天、交通工程和電氣工程等領(lǐng)域得到廣泛應(yīng)用．但是，復(fù)合材料在復(fù)雜應(yīng)力狀態(tài)下，內(nèi)部極易產(chǎn)生裂紋并擴展，最終導(dǎo)致材料的斷裂而引發(fā)事故．因此，對復(fù)合材料板中的裂紋缺陷問題進行分析具有重要的意義．

計算斷裂參數(shù)是進行斷裂分析的第一步，許多數(shù)值計算方法比如有限元法（Finite Element Meth－od，F(xiàn)EM）、擴展有限元法、有限差分法、邊界元法和無網(wǎng)格法等［1－2］都被嘗試用來計算斷裂參數(shù)，其中有限元法已經(jīng)成為求解斷裂參數(shù)的有效方法．由于采用位移有限元法理論得到的位移解偏小，文獻［3］提出了將形函數(shù)導(dǎo)數(shù)的域內(nèi)積分轉(zhuǎn)化為形函數(shù)的邊界線上的積分、網(wǎng)格劃分要求低和位移解更加準確的光滑有限元法（Smoothed Finite Element Method，SFEM）．虛擬裂紋閉合法（Virtual Crack Closure Technique，VCCT）［4］具有裂尖單元不需特殊處理和對網(wǎng)格尺寸要求低的優(yōu)點．SFEM 是Liu等［5］將光滑應(yīng)變措施引入有限元法，改進有限元法剛度結(jié)構(gòu)的一種方法，具有形函數(shù)簡單、對網(wǎng)格質(zhì)量要求低、計算精度高等優(yōu)點，現(xiàn)已廣泛應(yīng)用于各個領(lǐng)域［6－8］．VCCT 由Rybicki和Kanninen［9］于1977年提出的．Xie等［10－13］對VCCT 做了大量研究工作．VCCT 比外推法、等效積分區(qū)域積分法以及全局或局部虛擬裂紋擴展法求解斷裂參數(shù)具有明顯優(yōu)勢，它僅利用節(jié)點力與節(jié)點位移來計算應(yīng)變能釋放率，且只需要一步數(shù)值分析，最大程度地簡化了問題，具有高精度、高效率、裂尖單元不需特殊處理和對網(wǎng)格尺寸要求低等優(yōu)點．

本文基于SFEM 并結(jié)合VCCT，提出了SFEMVCCT 法，對含傾斜裂紋復(fù)合材料圓板的斷裂參數(shù)進行了數(shù)值分析，并與FEM－VCCT 計算結(jié)果進行了對比．

1 Cell－based光滑有限元法

均勻正交各向異性彈性力學(xué)平面問題的光滑Galerkin弱形式［5］可表示為：

式中：Ω為求解域；δ為變分符號；T 為矩陣的轉(zhuǎn)置；為應(yīng)變矩陣；為彈性矩陣（與柔度矩陣互逆）；為廣義位移；為體力；為力邊界Γ上的面力．

將求解域Ω離散為Ne個四邊形單元，節(jié)點個數(shù)為為空集，再將劃分為Ns＝4個光滑區(qū)域，如圖1所示，●為節(jié)點，□為光滑節(jié)點，○為高斯點，（N1N2N3N4）為該點處的位移形函數(shù)值．

廣義位移場為：

圖1 光滑域的劃分Fig．1 Division of an element into smoothing cells

光滑應(yīng)變?yōu)椋?/p>

式中：Φ為光滑函數(shù)，取

式中：Ac為第c光滑區(qū)域的面積，

將式（4）代入式（3），由分部積分得：

式中：Γc為光滑域Ωc的邊界；ni和nj分別為積分段外法向向量的分量．

將式（2）代入式（5），可得：

式中：nc為光滑單元個數(shù)．

FEM 通過對單元形函數(shù)矩陣求導(dǎo)得到單元應(yīng)變矩陣，通常采用高斯數(shù)值積分計算單元域積分．由式（7）可見，Cell－based光滑有限元計算光滑應(yīng)變矩陣時無需確定形函數(shù)在光滑域內(nèi)解析函數(shù)式及其導(dǎo)數(shù)，只需利用光滑域邊界各高斯點處的形函數(shù)，將形函數(shù)導(dǎo)數(shù)的域內(nèi)積分轉(zhuǎn)化為形函數(shù)的邊界線上的積分，提高了數(shù)值計算的精度和收斂性．

將式（6）和式（2）代入式（1），可得離散方程為：

式中：為整體光滑剛度矩陣，可由光滑單元剛度矩陣組裝得到．

F為力向量．

由上式可見，Cell－based光滑有限元法的形函數(shù)選取簡單，計算應(yīng)變矩陣時只需用形函數(shù)本身，對網(wǎng)格質(zhì)量要求低，編程簡單，容易實現(xiàn)．

2 虛擬裂紋閉合法

如圖2所示，長度為a的主裂紋前端虛擬擴展了長度為Δa的微小子裂紋，在此過程中裂紋虛擬擴展Δa時釋放的能量等于裂紋從a＋Δa閉合到初始實際裂紋a所需做的功．Irwin 的裂紋閉合積分為：

圖2 VCCT 示意圖Fig．2 The virtual crack closure technique

如圖3所示，基于光滑有限元網(wǎng)格，虛擬裂紋線上節(jié)點力在節(jié)點位移上做的功等于應(yīng)力所做的功，即

圖3 虛擬裂紋擴展法計算應(yīng)變能釋放率Fig．3 Virtural crack extension technique for strain energy release rate

經(jīng)整理得：

由于虛擬擴展裂紋尖端后面的張開位移和初始實際裂紋尖端后面的張開位移近似相等，式（11）可改寫為：

應(yīng)力分布為：

位移分布為：

將式（15）整理得：

式（18）的近似表達為：

類似地，Ⅱ型裂紋的計算公式為：

對于二維平面內(nèi)傾斜裂紋的虛擬裂紋閉合法可采用斷裂單元［14］．當(dāng)裂紋方向與各向異性材料某一對稱軸重合時，能量釋放率與應(yīng)力強度因子的關(guān)系為：

式中：Sij為柔度系數(shù)．

3 數(shù)值算例

為驗證SFEM－VCCT 的正確性與有效性，采用文獻［15］的算例，含中心斜裂紋復(fù)合材料圓板受集中載荷作用，幾何構(gòu)型和加載方式如圖4所示，裂紋長度為2a，α為裂紋傾斜角，板厚B＝1．0，材料參數(shù)E11＝0．1，E22＝1．0，G12＝0．5，v12＝0．03．

圖4 含中心斜裂紋復(fù)合材料圓板受集中載荷作用Fig．4 Geometry and boundary conditions for an orthotropic disk with inclined central crack subjected to point loads

圖5僅給出了2a＝2，α＝0o時，取四邊形單元數(shù)分別為4 138和1 500時單元分布情況，裂紋尖端單元正常離散．表1給出了當(dāng)α＝0o，2a＝2，2a＝4，2a＝6和2a＝8時，采用光滑有限元－虛擬裂紋閉合法（SFEM－VCCT）和有限元－虛擬裂紋閉合法（FEMVCCT）的單元個數(shù)及KI值．SFEM－VCCT相對FEM－VCCT 也不需要對裂尖單元特殊處理，與FEM－VCCT 所得結(jié)果基本一致，當(dāng)單元數(shù)為4 138和1 500 時，KI值分別為22．106 和21．997，與FEM－VCCT 計算結(jié)果的相對誤差僅為2% 和2．5%，可見，該方法完全繼承了VCCT 的優(yōu)點，不需要對裂尖單元特殊處理，對網(wǎng)格尺寸要求低，精度高．

圖6和圖7分別給出當(dāng)2a＝2，α＝0°，α＝15°，α＝30°和α＝45°，對應(yīng)的單元個數(shù)分別為3 814，4 167，4 045 和4 184 時，采 用SFEM－VCCT 和FEM－VCCT 得到的GⅠ和GⅡ值，所得結(jié)果基本一致，可見，SFEM－VCCT 法是正確有效的．

圖5 含中心斜裂紋復(fù)合材料圓板離散方式（α＝0°）Fig．5 Discretization of the orthoropic disk with an inclined central crack

表1 應(yīng)力強度因子、單元數(shù)與裂紋長度的關(guān)系（α＝0°）Tab．1 Stress intensity factors and elements corresponding to different inclined cracks in an orthotropic disk subjected to point loads（α＝0°）

圖6 GI與斜裂紋角度α的關(guān)系Fig．6 GIvalues corresponding to different central crack anglesαin the orthotropic disk

圖7 GⅡ與斜裂紋角度α的關(guān)系Fig．7 GⅡvalues corresponding to different central crack anglesαin the orthotropic disk

4 結(jié) 論

本文提出光滑有限元－虛擬裂紋閉合法，對含不同長度和角度的傾斜裂紋復(fù)合材料圓板的斷裂參數(shù)進行了模擬，并與有限元－虛擬裂紋閉合法計算結(jié)果進行了對比，得到如下結(jié)論：

1）SFEM－VCCT 計算時形函數(shù)簡單，對網(wǎng)格質(zhì)量要求低，形函數(shù)導(dǎo)數(shù)的域內(nèi)積分轉(zhuǎn)化為形函數(shù)的邊界線上的積分，編程簡單，容易實現(xiàn)．

2）SFEM－VCCT 不需要對裂尖單元特殊處理，單元數(shù)為4 138和1 500時，KⅠ值分別為22．106和21．997，與FEM－VCCT 計算結(jié)果的相對誤差僅為2%和2．5%，完全繼承了VCCT 的優(yōu)點．

［1］MOTAMEDI D，MOHAMMADI S．Dynamic analysis of fixed cracks in composites by the extended finite element method［J］．Engng Fract Mech，2010，77（17）：3373－3393．

［2］龍述堯，張國虎．基于MLPG 法的動態(tài)斷裂力學(xué)問題［J］．湖南大學(xué)學(xué)報：自然科學(xué)版，2012，39（11）：41－45．

LONG Shu－yao，ZHANG Guo－h(huán)u．An analysis of the dynamic fracture problem by the meshless local Petrov－Galerkin method［J］．Journal of Hunan University：Natural Sciences，2012，39（11）：41－45．（In Chinese）

［3］LIU G R，NGUYEN T T，DAI K Y，etal．Theoretical aspects of the smoothed finite element method（SFEM）［J］．Int J Numer Meth Eng，2007，71：902－930．

［4］RAJU I S．Calculation of strain－energy release rates with highorder and singular finite－elements［J］．Engineering Fracture Mechanics，1987，28：251－274．

［5］LIU G R，DAI K Y，NGUYEN T T．A smoothed finite element method for mechanics problems［J］．Computation Mechanics，2007，39（6）：859－877．

［6］CHEN J S，WU C T，YOON S．A stabilized conforming nodal integration for Galerkin meshfree method［J］．Int J Numer Meth Eng，2001，50：435－466．

［7］姚凌云，于德介，藏獻國．二維聲學(xué)數(shù)值計算的光滑有限元法［J］．機械工程學(xué)報，2010，46（18）：115－120．

YAO Ling－yun，YU De－jie，ZANG Xian－guo．Smoothed finite element method for two－dimensional acoustic numerical computation［J］．Journal of Mechanical Engineering，2010，46（18）：115－120．（In Chinese）

［8］周立明，孟廣偉，周振平，等．含孔復(fù)合材料層合板Cell－Based光滑有限元法［J］．東北大學(xué)學(xué)報：自然科學(xué)版，2013，34（Sup2）：90－93．

ZHOU Li－ming，MENG Guang－wei，ZHOU Zhen－ping，etal．Cell－Based smoothed finite element method for laminated composited plates with circular opening［J］．Journal of Northeastern University：Natural Science，2013，34（Sup2）：90－93．（In Chinese）

［9］RYBICKI E F，KANNINEN M F．A finite element calculation of stress intensity factors by a modified crack closure integral［J］．Engineering Fracture Mechanics，1977，9：931－938．

［10］XIE D，BIGGERSJR S B．Progressive crack growth using interface element based on the virtual crack closure technique［J］．Finite Elements in Analysis and Design，2006，42（11）：977－984．

［11］XIE D，SHERRILL B，BIGGERS J．Strain energy release rate calculation for a moving delamination front of arbitrary shape based on the virtual crack closure technique．part I：formulation and validation［J］．Engng Fracture Mech，2006，73（6）：771－785．

［12］XIE D，SHERRILL B，BIGGERS J．Strain energy release rate calculation for a moving delamination front of arbitrary shape based on the virtual crack closure technique．partⅡ：sensitivity study on modeling details［J］．Engng Fracture Mech，2006，73（6）：786－801．

［13］GHORASHI S S，MOHAMMADI S，SABBAGH－YAZDI S．Orthotropic enriched element free Galerkin method for fracture analysis of composites［J］．Engineering Fracture Mechanics，2011，78（9）：1906－1927．

［14］孟令兵，陳普會．層壓復(fù)合材料分層擴展分析的虛擬裂紋閉合技術(shù)及其應(yīng)用［J］．復(fù)合材料學(xué)報，2010，27（1）：190－195．

MENG Ling－bing，CHEN Pu－h(huán)ui．Virtual crack closure technique for delamination growth analysis of laminated composites and its application［J］．Acta Materiae Compositae Sinica，2010，27（1）：190－195．（In Chinese）

［15］ASADPOURE A，MOHAMMADI S．Developing new enrichment functions for crack simulation in orthotropic media by the extended finite element method［J］．Int J Numer Meth Eng，2007，69：2150－2172．