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

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

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

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

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

1 Cell－based光滑有限元法

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

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

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

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

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

式中：為整體光滑剛度矩陣，可由光滑單元?jiǎng)偠染仃嚱M裝得到．

F為力向量．

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

2 虛擬裂紋閉合法

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

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

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

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

經(jīng)整理得：

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

應(yīng)力分布為：

位移分布為：

將式（15）整理得：

式（18）的近似表達(dá)為：

類似地，Ⅱ型裂紋的計(jì)算公式為：

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

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

3 數(shù)值算例

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

圖6和圖7分別給出當(dāng)2a＝2，α＝0°，α＝15°，α＝30°和α＝45°，對(duì)應(yīng)的單元個(gè)數(shù)分別為3 814，4 167，4 045 和4 184 時(shí)，采 用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)力強(qiá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é) 論

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

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

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

［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 法的動(dòng)態(tài)斷裂力學(xué)問題［J］．湖南大學(xué)學(xué)報(bào)：自然科學(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］姚凌云，于德介，藏獻(xiàn)國．二維聲學(xué)數(shù)值計(jì)算的光滑有限元法［J］．機(jī)械工程學(xué)報(bào)，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é)報(bào)：自然科學(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］孟令兵，陳普會(huì)．層壓復(fù)合材料分層擴(kuò)展分析的虛擬裂紋閉合技術(shù)及其應(yīng)用［J］．復(fù)合材料學(xué)報(bào)，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．