基于光滑有限元的含裂纹复合材料的虚拟裂纹闭合法＊

周立明，孟广伟，王 晖，李 锋†，郭学东

（1．吉林大学 机械科学与工程学院，吉林 长春 130025；2．吉林大学 交通学院，吉林 长春 130025）

复合材料由于具有比强度高、比模量大、耐疲劳及抗破损等特点，在航空航天、交通工程和电气工程等领域得到广泛应用．但是，复合材料在复杂应力状态下，内部极易产生裂纹并扩展，最终导致材料的断裂而引发事故．因此，对复合材料板中的裂纹缺陷问题进行分析具有重要的意义．

计算断裂参数是进行断裂分析的第一步，许多数值计算方法比如有限元法（Finite Element Meth－od，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 法，对含倾斜裂纹复合材料圆板的断裂参数进行了数值分析，并与FEM－VCCT 计算结果进行了对比．

1 Cell－based光滑有限元法

均匀正交各向异性弹性力学平面问题的光滑Galerkin弱形式［5］可表示为：

式中：Ω为求解域；δ为变分符号；T 为矩阵的转置；为应变矩阵；为弹性矩阵（与柔度矩阵互逆）；为广义位移；为体力；为力边界Γ上的面力．

将求解域Ω离散为Ne个四边形单元，节点个数为为空集，再将划分为Ns＝4个光滑区域，如图1所示，●为节点，□为光滑节点，○为高斯点，（N1N2N3N4）为该点处的位移形函数值．

广义位移场为：

图1 光滑域的划分Fig．1 Division of an element into smoothing cells

光滑应变为：

式中：Φ为光滑函数，取

式中：Ac为第c光滑区域的面积，

将式（4）代入式（3），由分部积分得：

式中：Γc为光滑域Ωc的边界；ni和nj分别为积分段外法向向量的分量．

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

式中：nc为光滑单元个数．

FEM 通过对单元形函数矩阵求导得到单元应变矩阵，通常采用高斯数值积分计算单元域积分．由式（7）可见，Cell－based光滑有限元计算光滑应变矩阵时无需确定形函数在光滑域内解析函数式及其导数，只需利用光滑域边界各高斯点处的形函数，将形函数导数的域内积分转化为形函数的边界线上的积分，提高了数值计算的精度和收敛性．

将式（6）和式（2）代入式（1），可得离散方程为：

式中：为整体光滑刚度矩阵，可由光滑单元刚度矩阵组装得到．

F为力向量．

由上式可见，Cell－based光滑有限元法的形函数选取简单，计算应变矩阵时只需用形函数本身，对网格质量要求低，编程简单，容易实现．

2 虚拟裂纹闭合法

如图2所示，长度为a的主裂纹前端虚拟扩展了长度为Δa的微小子裂纹，在此过程中裂纹虚拟扩展Δa时释放的能量等于裂纹从a＋Δa闭合到初始实际裂纹a所需做的功．Irwin 的裂纹闭合积分为：

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

如图3所示，基于光滑有限元网格，虚拟裂纹线上节点力在节点位移上做的功等于应力所做的功，即

图3 虚拟裂纹扩展法计算应变能释放率Fig．3 Virtural crack extension technique for strain energy release rate

经整理得：

由于虚拟扩展裂纹尖端后面的张开位移和初始实际裂纹尖端后面的张开位移近似相等，式（11）可改写为：

应力分布为：

位移分布为：

将式（15）整理得：

式（18）的近似表达为：

类似地，Ⅱ型裂纹的计算公式为：

对于二维平面内倾斜裂纹的虚拟裂纹闭合法可采用断裂单元［14］．当裂纹方向与各向异性材料某一对称轴重合时，能量释放率与应力强度因子的关系为：

式中：Sij为柔度系数．

3 数值算例

为验证SFEM－VCCT 的正确性与有效性，采用文献［15］的算例，含中心斜裂纹复合材料圆板受集中载荷作用，几何构型和加载方式如图4所示，裂纹长度为2a，α为裂纹倾斜角，板厚B＝1．0，材料参数E11＝0．1，E22＝1．0，G12＝0．5，v12＝0．03．

图4 含中心斜裂纹复合材料圆板受集中载荷作用Fig．4 Geometry and boundary conditions for an orthotropic disk with inclined central crack subjected to point loads

图5仅给出了2a＝2，α＝0o时，取四边形单元数分别为4 138和1 500时单元分布情况，裂纹尖端单元正常离散．表1给出了当α＝0o，2a＝2，2a＝4，2a＝6和2a＝8时，采用光滑有限元－虚拟裂纹闭合法（SFEM－VCCT）和有限元－虚拟裂纹闭合法（FEMVCCT）的单元个数及KI值．SFEM－VCCT相对FEM－VCCT 也不需要对裂尖单元特殊处理，与FEM－VCCT 所得结果基本一致，当单元数为4 138和1 500 时，KI值分别为22．106 和21．997，与FEM－VCCT 计算结果的相对误差仅为2% 和2．5%，可见，该方法完全继承了VCCT 的优点，不需要对裂尖单元特殊处理，对网格尺寸要求低，精度高．

图6和图7分别给出当2a＝2，α＝0°，α＝15°，α＝30°和α＝45°，对应的单元个数分别为3 814，4 167，4 045 和4 184 时，采 用SFEM－VCCT 和FEM－VCCT 得到的GⅠ和GⅡ值，所得结果基本一致，可见，SFEM－VCCT 法是正确有效的．

图5 含中心斜裂纹复合材料圆板离散方式（α＝0°）Fig．5 Discretization of the orthoropic disk with an inclined central crack

表1 应力强度因子、单元数与裂纹长度的关系（α＝0°）Tab．1 Stress intensity factors and elements corresponding to different inclined cracks in an orthotropic disk subjected to point loads（α＝0°）

图6 GI与斜裂纹角度α的关系Fig．6 GIvalues corresponding to different central crack anglesαin the orthotropic disk

图7 GⅡ与斜裂纹角度α的关系Fig．7 GⅡvalues corresponding to different central crack anglesαin the orthotropic disk

4 结 论

本文提出光滑有限元－虚拟裂纹闭合法，对含不同长度和角度的倾斜裂纹复合材料圆板的断裂参数进行了模拟，并与有限元－虚拟裂纹闭合法计算结果进行了对比，得到如下结论：

1）SFEM－VCCT 计算时形函数简单，对网格质量要求低，形函数导数的域内积分转化为形函数的边界线上的积分，编程简单，容易实现．

2）SFEM－VCCT 不需要对裂尖单元特殊处理，单元数为4 138和1 500时，KⅠ值分别为22．106和21．997，与FEM－VCCT 计算结果的相对误差仅为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 Shu－yao，ZHANG Guo－hu．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］姚凌云，于德介，藏献国．二维声学数值计算的光滑有限元法［J］．机械工程学报，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］周立明，孟广伟，周振平，等．含孔复合材料层合板Cell－Based光滑有限元法［J］．东北大学学报：自然科学版，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］孟令兵，陈普会．层压复合材料分层扩展分析的虚拟裂纹闭合技术及其应用［J］．复合材料学报，2010，27（1）：190－195．

MENG Ling－bing，CHEN Pu－hui．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．