复合材料由于具有比强度高、比模量大、耐疲劳及抗破损等特点，在航空航天、交通工程和电气工程等领域得到广泛应用.但是，复合材料在复杂应力状态下，内部极易产生裂纹并扩展，最终导致材料的断裂而引发事故.因此，对复合材料板中的裂纹缺陷问题进行分析具有重要的意义.
摘要：为提高求解断裂参数的精度和效率，将光滑有限元法与虚拟裂纹闭合法相结合，提出光滑有限元虚拟裂纹闭合法.对含不同长度和角度的倾斜裂纹复合材料圆板的断裂参数进行了求解，并与有限元虚拟裂纹闭合法计算结果进行了对比.数值算例结果验证了该方法具有高精度，裂纹尖端处单元不需特殊处理，对网格尺寸要求低等优点，是分析断裂问题简洁高效的数值计算方法.
关键词：职称论文发表平台投稿,光滑有限元法,正交各向异性,虚拟裂纹闭合法,应变能释放率,应力强度因子
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.
