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

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

收稿日期：20131119

基金項(xiàng)目：國(guó)家重大科學(xué)儀器設(shè)備開(kāi)發(fā)專(zhuān)項(xiàng)項(xiàng)目（2012YQ030075）；國(guó)家自然科學(xué)基金資助項(xiàng)目（51305157）；吉林省科技廳基金資助項(xiàng)目（20130305006GX，201215048，201205011）

作者簡(jiǎn)介：周立明（1982-），男，吉林白山人，吉林大學(xué)講師，博士，在站博士后

通訊聯(lián)系人，Email：fengli@jlu.edu.cn

摘要：為提高求解斷裂參數(shù)的精度和效率，將光滑有限元法與虛擬裂紋閉合法相結(jié)合，提出光滑有限元虛擬裂紋閉合法.對(duì)含不同長(zhǎng)度和角度的傾斜裂紋復(fù)合材料圓板的斷裂參數(shù)進(jìn)行了求解，并與有限元虛擬裂紋閉合法計(jì)算結(jié)果進(jìn)行了對(duì)比.數(shù)值算例結(jié)果驗(yàn)證了該方法具有高精度，裂紋尖端處單元不需特殊處理，對(duì)網(wǎng)格尺寸要求低等優(yōu)點(diǎn)，是分析斷裂問(wèn)題簡(jiǎn)潔高效的數(shù)值計(jì)算方法.

關(guān)鍵詞：光滑有限元法；正交各向異性；虛擬裂紋閉合法；應(yīng)變能釋放率；應(yīng)力強(qiáng)度因子

中圖分類(lèi)號(hào)：TB115 文獻(xiàn)標(biāo)識(shí)碼：A

Virtual Crack Closure Technique Based on Smoothed Finite

Element Method for Composite Materials with Cracks

ZHOU Liming1，MENG Guangwei1，WANG Hui1，LI Feng1，GUO Xuedong2

（1.School of Mechanical Science and Engineering， Jilin Univ， Changchun， Jilin130025， China；

2.College of Traffic， Jilin Univ， Changchun， Jilin130025， China）

Abstract：To improve the solution efficiency and accuracy of the fracture parameters， the smoothed finite element methodvirtual crack close method was proposed based on the combination of these two methods. The fracture parameters of the composite material circular plate with different length and angle oblique cracks were solved， and the result was compared with that of the finite element methodvirtual crack closure method. The calculation result confirms that this method has advantages of high accuracy， no special treatment on the elements at crack tip and less mesh resolution requirement. This method is simple but efficient for the calculation of fracture problems.

Key words：smoothed finite element method；orthotropic；Virtual Crack Closure Technique（VCCT）；strain energy release rate；stress intensity factor

復(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ù)合材料板中的裂紋缺陷問(wèn)題進(jìn)行分析具有重要的意義.

計(jì)算斷裂參數(shù)是進(jìn)行斷裂分析的第一步，許多數(shù)值計(jì)算方法比如有限元法（Finite Element Method， FEM）、擴(kuò)展有限元法、有限差分法、邊界元法和無(wú)網(wǎng)格法等[1-2]都被嘗試用來(lái)計(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ù)簡(jiǎn)單、對(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)勢(shì)，它僅利用節(jié)點(diǎn)力與節(jié)點(diǎn)位移來(lái)計(jì)算應(yīng)變能釋放率，且只需要一步數(shù)值分析，最大程度地簡(jiǎn)化了問(wèn)題，具有高精度、高效率、裂尖單元不需特殊處理和對(duì)網(wǎng)格尺寸要求低等優(yōu)點(diǎn).

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

1Cellbased光滑有限元法

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

∫ΩδT（）（）dΩ-∫ΩδTdΩ-∫ΓδTdΓ=0.（1）

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

將求解域Ω離散為Ne個(gè)四邊形單元，節(jié)點(diǎn)個(gè)數(shù)為Nd，Ω=∪Nei=1Ωei，Ωei∩Ωej=，i≠j，為空集，再將Ωei劃分為Ns=4個(gè)光滑區(qū)域，如圖1所示，●為節(jié)點(diǎn)，□為光滑節(jié)點(diǎn)，○為高斯點(diǎn)，（N1 N2N3N4）為該點(diǎn)處的位移形函數(shù)值.

廣義位移場(chǎng)為：

=uvT=∑npi=1Nii.（2）

式中：i=uivi為廣義節(jié)點(diǎn)位移；Ni為形函數(shù)對(duì)角矩陣；np=4.

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

xc=∫ΩcxΦx-xcdΩ.（3）

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

Φx-xc=1/Ac，x∈Ωc，

0，xΩc.（4）

式中：Ac為第c光滑區(qū)域的面積，Ac=∫ΩcdΩ.

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

xc=12Ac∫Γcuinj+ujnidΓ.（5）

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

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

xc=∑nci=1ixcqi.（6）

式中：nc為光滑單元個(gè)數(shù).

ixc=1Ac∑nbb=1NixGbnx0

0NixGbny

NixGbnyNixGbnx lcb.（7）

式中：xGb和lcb分別為光滑邊界Γcb的中點(diǎn)（高斯點(diǎn)）和長(zhǎng)度；nb為每個(gè)光滑元的邊界總數(shù).

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

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

=F.（8）

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

ij=∑Nsk=1TiDjAk；（9）

F為力向量.

F=∫ΩNTdΩ-∫ΓNTdΓ.（10）

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

2虛擬裂紋閉合法

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

GⅠ=lim Δa→012BΔa∫Δa0σ1yyΔa-r，0Δv2r，πdr，（11）

GⅡ=lim Δa→012BΔa∫Δa0τ1yyΔa-r，0Δu2r，πdr.（12）

式中：σ1yy和τ1yy為原始裂紋尖端處的應(yīng)力分量；Δv2為裂紋虛擬擴(kuò)展到a+Δa時(shí)裂紋面上的張開(kāi)位移；Δu2為裂紋虛擬擴(kuò)展到a+Δa時(shí)裂紋面上的相對(duì)滑動(dòng)位移；B為裂紋體厚度；GI和GⅡ分別為Ⅰ型和Ⅱ型裂紋的應(yīng)變能釋放率分量.

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

F1y1v21，1′=∫Δa0σ1yyxΔv2xdx.（13）

式中：v21，1′=v1-v′1，為節(jié)點(diǎn)1和1′之間的垂直位移變化量；Fy1為節(jié)點(diǎn)1上豎直方向的節(jié)點(diǎn)力；上標(biāo)（1）和（2）分別為初始實(shí)際裂紋和虛擬擴(kuò)展裂紋；應(yīng)力σ1yyx=Ax，其中，A為常數(shù)；對(duì)于線性四邊形單元，位移為Δv2x=1-x/Δav21，1′.

經(jīng)整理得：

F1y1=∫Δa0Ax1-xΔadx=43AΔa.（14）

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

GI=lim Δa→012BΔa∫Δa0σ1yyΔa-r，0Δv1r，πdr.（15）

應(yīng)力分布為：

σ1yyΔa-r=AΔa-r.（16）

位移分布為：

Δv1r，π=rΔaΔv13，4.（17）

將式（15）整理得：

GI=lim Δa→012BΔa∫Δa0AΔa-rrΔaΔv13，4dr=

lim Δa→012BΔaF1y1Δv13，4.（18）

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

GI=12BΔaFy1Δv3，4. （19）

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

GⅡ=12BΔaFx1Δu3，4.（20）

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

GⅠ=K2IS11S22212S22S1112+2S12+S662S1112，（21）

GⅡ=K2ⅡS112S22S1112+2S12+S662S1112. （22）

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

3數(shù)值算例

為驗(yàn)證SFEMVCCT的正確性與有效性，采用文獻(xiàn)[15]的算例，含中心斜裂紋復(fù)合材料圓板受集中載荷作用，幾何構(gòu)型和加載方式如圖4所示，裂紋長(zhǎng)度為2a， SymbolaA@ 為裂紋傾斜角，板厚B=1.0，材料參數(shù)E11=0.1，E22=1.0，G12=0.5，v12=0.03.

圖5僅給出了2a=2，α=0o時(shí)，取四邊形單元數(shù)分別為4 138和1 500時(shí)單元分布情況，裂紋尖端單元正常離散.表1給出了當(dāng)α=0o，2a=2，2a=4，2a=6和2a=8時(shí)，采用光滑有限元虛擬裂紋閉合法（SFEMVCCT）和有限元虛擬裂紋閉合法（FEMVCCT）的單元個(gè)數(shù)及KI值.SFEMVCCT相對(duì)FEMVCCT也不需要對(duì)裂尖單元特殊處理，與FEMVCCT所得結(jié)果基本一致，當(dāng)單元數(shù)為4 138和1 500時(shí)，KI值分別為22.106和21.997，與FEMVCCT計(jì)算結(jié)果的相對(duì)誤差僅為2%和2.5%，可見(jiàn)，該方法完全繼承了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í)，采用SFEMVCCT和FEMVCCT得到的GⅠ和GⅡ值，所得結(jié)果基本一致，可見(jiàn)，SFEMVCCT法是正確有效的.

4結(jié)論

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

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

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

參考文獻(xià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]龍述堯， 張國(guó)虎. 基于MLPG法的動(dòng)態(tài)斷裂力學(xué)問(wèn)題[J]. 湖南大學(xué)學(xué)報(bào)：自然科學(xué)版， 2012， 39（11）：41-45.

LONG Shuyao， ZHANG Guohu. An analysis of the dynamic fracture problem by the meshless local PetrovGalerkin 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， et al. 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 strainenergy release rates with highorder and singular finiteelements[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)國(guó). 二維聲學(xué)數(shù)值計(jì)算的光滑有限元法[J]. 機(jī)械工程學(xué)報(bào)，2010，46（18）：115-120.

YAO Lingyun， YU Dejie， ZANG Xianguo. Smoothed finite element method for twodimensional acoustic numerical computation[J]. Journal of Mechanical Engineering， 2010，46（18）：115-120.（In Chinese）

[8]周立明，孟廣偉，周振平，等. 含孔復(fù)合材料層合板CellBased 光滑有限元法[J]. 東北大學(xué)學(xué)報(bào)：自然科學(xué)版，2013，34（Sup2）：90-93.

ZHOU Liming， MENG Guangwei， ZHOU Zhenping， et al. CellBased 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， SABBAGHYAZDI 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 Lingbing，CHEN Puhui. 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.

[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， SABBAGHYAZDI 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 Lingbing，CHEN Puhui. 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.

[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， SABBAGHYAZDI 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 Lingbing，CHEN Puhui. 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.