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

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

　　關鍵詞：職稱論文發表平臺投稿,光滑有限元法,正交各向異性,虛擬裂紋閉合法,應變能釋放率,應力強度因子

　　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

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

　　本文基于SFEM并結合VCCT，提出了SFEMVCCT法，對含傾斜裂紋復合材料圓板的斷裂參數進行了數值分析，并與FEMVCCT計算結果進行了對比.

　　1Cellbased光滑有限元法 　　均勻正交各向異性彈性力學平面問題的光滑Galerkin弱形式[5]可表示為：

　　∫ΩδT()()dΩ-∫ΩδTdΩ-∫ΓδTdΓ=0.(1)

　　式中：Ω為求解域;δ為變分符號;T為矩陣的轉置;為應變矩陣;為彈性矩陣(與柔度矩陣互逆);為廣義位移;為體力;為力邊界Γ上的面力.

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

　　廣義位移場為：

　　=uvT=∑npi=1Nii.(2)

　　式中：i=uivi為廣義節點位移;Ni為形函數對角矩陣;np=4.

　　光滑應變為：

　　xc=∫ΩcxΦx-xcdΩ.(3)

　　式中：Φ為光滑函數，取

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

　　0，xΩc.(4)

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

　　將式(4)代入式(3)，由分部積分得：

　　xc=12Ac∫Γcuinj+ujnidΓ.(5)

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

　　將式(2)代入式(5)，可得：

　　xc=∑nci=1ixcqi.(6)

　　式中：nc為光滑單元個數.

　　ixc=1Ac∑nbb=1NixGbnx0

　　0NixGbny

　　NixGbnyNixGbnx lcb.(7)

　　式中：xGb和lcb分別為光滑邊界Γcb的中點(高斯點)和長度;nb為每個光滑元的邊界總數.

　　FEM通過對單元形函數矩陣求導得到單元應變矩陣，通常采用高斯數值積分計算單元域積分.由式(7)可見，Cellbased光滑有限元計算光滑應變矩陣時無需確定形函數在光滑域內解析函數式及其導數，只需利用光滑域邊界各高斯點處的形函數，將形函數導數的域內積分轉化為形函數的邊界線上的積分，提高了數值計算的精度和收斂性.

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

　　=F.(8)

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

　　ij=∑Nsk=1TiDjAk;(9)

　　F為力向量.

　　F=∫ΩNTdΩ-∫ΓNTdΓ.(10)

　　由上式可見，Cellbased光滑有限元法的形函數選取簡單，計算應變矩陣時只需用形函數本身，對網格質量要求低，編程簡單，容易實現.

　　2虛擬裂紋閉合法

　　如圖2所示，長度為a的主裂紋前端虛擬擴展了長度為Δa的微小子裂紋，在此過程中裂紋虛擬擴展Δa時釋放的能量等于裂紋從a+Δa閉合到初始實際裂紋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為原始裂紋尖端處的應力分量;Δv2為裂紋虛擬擴展到a+Δa時裂紋面上的張開位移;Δu2為裂紋虛擬擴展到a+Δa時裂紋面上的相對滑動位移;B為裂紋體厚度;GI和GⅡ分別為Ⅰ型和Ⅱ型裂紋的應變能釋放率分量.

　　如圖3所示，基于光滑有限元網格，虛擬裂紋線上節點力在節點位移上做的功等于應力所做的功，即

　　F1y1v21，1′=∫Δa0σ1yyxΔv2xdx.(13)

　　式中：v21，1′=v1-v′1，為節點1和1′之間的垂直位移變化量;Fy1為節點1上豎直方向的節點力;上標(1)和(2)分別為初始實際裂紋和虛擬擴展裂紋;應力σ1yyx=Ax，其中，A為常數;對于線性四邊形單元，位移為Δv2x=1-x/Δav21，1′.

　　經整理得：

　　F1y1=∫Δa0Ax1-xΔadx=43AΔa.(14)

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

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

　　應力分布為：

　　σ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)的近似表達為：

　　GI=12BΔaFy1Δv3，4. (19)

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

　　GⅡ=12BΔaFx1Δu3，4.(20)

　　對于二維平面內傾斜裂紋的虛擬裂紋閉合法可采用斷裂單元[14].當裂紋方向與各向異性材料某一對稱軸重合時，能量釋放率與應力強度因子的關系為：

　　GⅠ=K2IS11S22212S22S1112+2S12+S662S1112，(21)

　　GⅡ=K2ⅡS112S22S1112+2S12+S662S1112. (22)

　　式中：Sij為柔度系數.

　　3數值算例

　　為驗證SFEMVCCT的正確性與有效性，采用文獻[15]的算例，含中心斜裂紋復合材料圓板受集中載荷作用，幾何構型和加載方式如圖4所示，裂紋長度為2a， SymbolaA@ 為裂紋傾斜角，板厚B=1.0，材料參數E11=0.1，E22=1.0，G12=0.5，v12=0.03. 　　圖5僅給出了2a=2，α=0o時，取四邊形單元數分別為4 138和1 500時單元分布情況，裂紋尖端單元正常離散.表1給出了當α=0o，2a=2，2a=4，2a=6和2a=8時，采用光滑有限元虛擬裂紋閉合法(SFEMVCCT)和有限元虛擬裂紋閉合法(FEMVCCT)的單元個數及KI值.SFEMVCCT相對FEMVCCT也不需要對裂尖單元特殊處理，與 FEMVCCT所得結果基本一致，當單元數為4 138和1 500時，KI值分別為22.106和21.997，與FEMVCCT計算結果的相對誤差僅為2%和2.5%，可見，該方法完全繼承了VCCT的優點，不需要對裂尖單元特殊處理，對網格尺寸要求低，精度高.

　　圖6和圖7分別給出當2a=2，α=0°，α=15°，α=30°和α=45°，對應的單元個數分別為3 814，4 167，4 045和4 184時，采用SFEMVCCT和FEMVCCT得到的GⅠ和GⅡ值，所得結果基本一致，可見，SFEMVCCT法是正確有效的.

　　4結論

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

　　1)SFEMVCCT計算時形函數簡單，對網格質量要求低，形函數導數的域內積分轉化為形函數的邊界線上的積分，編程簡單，容易實現.

　　2)SFEMVCCT不需要對裂尖單元特殊處理，單元數為4 138和1 500時，KⅠ值分別為22.106和21.997，與FEMVCCT計算結果的相對誤差僅為2%和2.5%，完全繼承了VCCT的優點.

　　參考文獻

　　[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法的動態斷裂力學問題[J]. 湖南大學學報：自然科學版， 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.