吉林大学 机械科学与工程学院, 吉林 长春 130025
收稿日期：2016-03-15
基金项目：国家自然科学基金资助项目(11502092)；吉林省科技厅基金资助项目(20160520064JH, 20170101043JC)；中央高校基本科研业务费专项资金资助项目(451170306066)。
作者简介：周立明(1982-), 男, 吉林白山人, 吉林大学副教授;
孟广伟(1959-), 男, 吉林长春人, 吉林大学教授。
摘要：为求解含裂纹功能梯度板的动态能量释放率, 提出了基于非均匀有限元法的虚拟裂纹闭合法.基于有限元软件ABAQUS平台, 开发了动态载荷作用下功能梯度材质的哑节点断裂单元, 编写了用户自定义子程序UMAT和UEL, 求解了动态载荷作用下的含裂纹功能梯度板能量释放率分量.对不同形状因子下的中心和倾斜裂纹功能梯度板的动态能量释放率分量进行了求解, 并转化为动态应力强度因子, 与差分法作比较.数值算例结果表明, 动态载荷作用下功能梯度材质的哑节点断裂单元具有精度高、简捷方便、裂尖无需使用奇异元或特殊单元和子程序通用性强等优点, 为工程师求解动态断裂参数提供了新途径.
关键词：裂纹动态能量释放率虚拟裂纹闭合法动态应力强度因子非均匀有限元法
ABAQUS User Subroutine Development for Energy Releasing Rate of the Functionally Graded Plate with Cracks
ZHOU Li-ming, REN Shu-hui, MENG Guang-wei, LI Rong-jia
School of Mechanical Science and Engineering, Jilin University, Changchun 130025, China
Corresponding author: LI Rong-jia, E-mail: lmzhou@jlu.edu.cn
Abstract: To obtain the dynamic energy releasing rate of the functionally graded plate with cracks, a virtual crack closure technique was proposed based on the inhomogeneous finite element method. A fracture dummy node element of the functionally graded material under dynamic loading was developed, user-defined subroutines named UMAT and UEL were programmed, and the component of energy releasing rate for the functionally graded plate caused by dynamic loading was calculated based on the finite element software ABAQUS platform. For the different shape factors such as the center-crack and inclined-crack inside the functionally graded plates, the component of the dynamic energy releasing rate was solved, and then it was transformed into the dynamic stress intensity factor, which was used for comparing the results calculated by the finite difference method. The numerical calculating results showed that the dummy node element of the functionally graded material fracture caused by dynamic loading is characterized by high precision, simplicity and convenience without using the singular element or special element for the crack tip and the versatility of its subroutines, which may provide a new approach for the engineers to get the dynamic fracture parameters.
Key Words: crackdynamic energy releasing ratevirtual crack closure techniquedynamic stress intensity factorinhomogeneous finite element method
近年来, 人们对含裂纹功能材料体在碰撞事件中受到冲击载荷作用时的响应日益关注.为了精确评价动态载荷作用下结构的断裂力学问题, 提出了动态断裂参数的概念, 如动态应力强度因子、动态能量释放率等.对于简单的几何模型、理想的材料模型和特殊加载模型的动态断裂参数可由解析法求得.但往往由于结构或边界条件比较复杂, 解析法并不可行, 实验测试又费时费钱, 数值计算就成为研究该类问题的有力工具.
Xie等[1-2]基于ABAQUS软件提出了界面单元, 求解了冲击载荷作用下均匀材料板的动态能量释放率, 拓展了虚拟裂纹闭合法的应用范围.Wimmer等[3-5]对层压复合材料组件的强度和断裂问题进行了数值模拟.Msekh等[6]采用ABAQUS建立了脆性断裂的相场模型.Wang等[7]采用实验观察和有限元法对氧化铝陶瓷的动态疲劳和断裂问题进行了分析.Monfared等[8-9]对含多裂纹的功能梯度涂层的动态断裂参数进行了求解.Manolis等[10]基于边界元法对含裂纹光滑非均匀材料板的动态断裂问题进行了研究.Lee等[11]解决了正交异性功能梯度材料的动态裂纹扩展问题.Abdollahifar等[12]采用无网格局部彼得洛夫伽辽金法求解了功能梯度板的动态应力强度因子.Rouzegar等[13]建立了Kirchhoff板壳动态裂纹扩展的扩展有限元模型.Bui等[14]采用NURBS法对含裂纹压电材料的动态裂纹问题进行了研究.Zhao等[15]对冲击载荷作用下反平面功能梯度材料无限板的动态应力强度因子进行了求解.Bhardwaj等[16]提出了基于XIGA的NURBS法, 对不同载荷和边界条件下的功能梯度裂纹板的动态问题进行了模拟.C, eribas, s等[17]对超椭圆均匀功能梯度材料薄板的动态问题进行了模拟.可见, 目前新兴的计算方法对研究问题的针对性较强, 不具有通用性, 推广起来相对比较困难; 基于ANSYS，ABAQUS等软件的二次开发分析动态裂纹问题主要集中在均匀材料上, 对功能梯度材料的动态裂纹问题研究还有待进一步拓展.
本文基于商业有限元软件ABAQUS, 结合非均匀有限元法, 将虚拟裂纹闭合法(virtual crack closure technique, VCCT)扩展到非均匀材质中, 开发了功能梯度材料的哑节点断裂单元, 编写了用户自定义子程序UMAT和UEL, 求解了动态载荷作用下的含裂纹功能梯度板能量释放率分量.
1 功能梯度材料功能梯度材料物参(如弹性模量、泊松比等)可遵循如下梯度分布规律:
 | (1) |
 | (2) |
图 1(Fig. 1)
 | 图 1 无量纲厚度z/h和体积分数的关系Fig.1 Variation of the volume fraction versus the non-dimensional thickness z/h |
2 非均匀有限元法非均匀有限元方程为
 | (3) |
单元刚度矩阵Ke的表达式为
 | (4) |
式(4) 采用高斯积分求解, 二维问题可写成
 | (5) |
 | (6) |
 | (7) |
3 虚拟裂纹闭合法计算平面功能梯度材料内任意线状裂纹的应变能释放率, 可采用哑节点断裂单元, 如图 2所示.断裂单元中裂尖的节点1和节点2坐标重合; 节点3和节点4坐标重合; 节点5在裂尖前面, Δa为虚拟裂纹增量, Δc为裂尖后单元尺寸, 取Δc=Δa.
图 2(Fig. 2)
 | 图 2 哑节点断裂单元Fig.2 Fracture element of the dummy node |
节点位移矢量为
 | (8) |
 | (9) |
x和x′轴的夹角为
 | (10) |
 | (11) |
图 3(Fig. 3)
 | 图 3 哑节点断裂单元的位移矢量Fig.3 Displacement vector for fracture element of the dummy node |
 | (12) |
 | (13) |
 | (14) |
 | (15) |
 | (16) |
 | (17) |
 | (18) |
 | (19) |
应变能释放率为
 | (20) |
 | (21) |
动态载荷作用下的固定裂纹问题的大量实验结果表明无需在以上求解公式中增加惯性项或在裂尖处施加特殊质量的质点[1].静态的虚拟裂纹闭合法计算公式完全适用于计算动态载荷作用下的固定裂纹问题, 计算过程中注意各个物理单位的匹配.
动态应力强度因子和动态能量释放率的关系为
 | (22) |
 | (23) |
4 数值算例4.1 算例1图 4为含中心裂纹的功能梯度板受单向冲击载荷作用图, 板长L=40 mm, 板宽W=20 mm, 板厚B=0.1 mm, 裂纹长度2a=4.8 mm.功能梯度材料由Y-TZP和γ-TiAl两种材料构成, 物参服从式(1), Y-TZP的弹性模量Mh=200 GPa, γ-TiAl的弹性模量Mb=186 GPa, 两种材料的泊松比ν均为0.3, 密度ρ=5×103 kg/m3; 载荷为Heaviside阶跃函数的形式, 如图 5所示, 载荷幅值为σ0=0.4 GPa, 无阻尼.采用VCCT求解动态能量释放率后转化成动态应力强度因子, 并进行无量纲化处理, 即
图 4(Fig. 4)
 | 图 4 受冲击载荷作用的中心裂纹功能梯度板Fig.4 Center-crack functionally graded plate under impact loading> |
图 5(Fig. 5)
 | 图 5 阶跃函数的载荷形式Fig.5 Loading type of Heaviside step function |
 | (24) |
图 6(Fig. 6)
 | 图 6 不同形状因子下中心裂纹功能梯度板的动态应力强度因子KⅠFig.6 Dynamic stress intensity factor KⅠ of the center-crack functionally graded plate with different shape factors |
图 7(Fig. 7)
 | 图 7 不同时间步时中心裂纹功能梯度板的变形图Fig.7 Deformation figures of the center-crack functionally graded plate under different time steps |
4.2 算例2图 8是受单向冲击载荷作用的含中心倾斜裂纹板, 板长L=60 mm, 板宽W=30 mm, 单位板厚, 裂纹长度a=14.14 mm, 倾斜角α=45°.功能梯度材料由Y-TZP和γ-TiAl两种材料构成, 物参服从式(1).Y-TZP的弹性模量Mh=200 GPa, γ-TiAl的弹性模量Mb=186 GPa, 两种材料的泊松比ν均为0.3, 密度ρ=5×103 kg/m3, 载荷形式如图 5所示, 幅值为σ0=0.4 GPa, 无阻尼.
图 8(Fig. 8)
 | 图 8 受冲击载荷作用的倾斜裂纹功能梯度板Fig.8 Inclined-crack functionally graded plate under impact loading |
将结构均匀划分成200×100网格, ABAQUS中采用CPE4单元, 时间步长Δt=0.2 μs, 弹簧刚度设置为Mh的1 000倍.
图 9和图 10分别给出了形状因子n=0.0, 0.2, 0.5, 1.0和5.0时, 采用本文方法计算求解所得的含倾斜裂纹功能梯度板的无量纲动态应力强度因子KⅠ和KⅡ.当n=0时, 材料梯度无变化, 梯度板的材料仅为材料Y-TZP, 与有限差分法(FDM)和文献[1]的计算结果作对比, 从图中可以看出, n=0时, 本文方法与文献[1]的计算结果基本相同, 与FDM计算结果吻合得也很好, 进一步证明了本文方法的正确性, 可见本文方法具有精度高和求解简单的优点.同时通过借助ABAQUS计算结果提取相关信息来计算动态应力强度因子, 也使得程序具有了较强的通用性, 较大的拓展空间.图 11给出了不同方法、不同时间步时倾斜裂纹功能梯度板的变形图, 可见本文方法是正确有效的.
图 9(Fig. 9)
 | 图 9 不同形状因子下倾斜裂纹功能梯度板的动态应力强度因子KⅠFig.9 Dynamic stress intensity factor KⅠ of the inclined-crack functionally graded plate with different shape factors |
图 10(Fig. 10)
 | 图 10 不同形状因子下倾斜裂纹功能梯度板的动态应力强度因子KⅡFig.10 Dynamic stress intensity factor KⅡ of inclined-crack functionally graded plate with different shape factors |
图 11(Fig. 11)
 | 图 11 不同时间步时倾斜裂纹功能梯度板的变形图Fig.11 Deformation figures of inclined-crack functionally graded plate under different time steps |
5 结论1) 本文方法求解动态能量释放率时具有精度高、计算方便、裂尖单元无需特殊处理等优点.
2) 本文编写了基于非均质材料的UMAT和UEL, 通过借助ABAQUS计算结果提取相关信息来计算动态应力强度因子, 使得程序具有了较强的通用性、较大的拓展空间.
3) 当功能梯度材料的组分材料弹性模量值和形状因子变化不大时, 动态应力强度因子的变化也较小.
参考文献
| [1] | Xie D, Biggers S B. Calculation of transient strain energy release rates under impact loading based on the virtual crack closure technique[J].International Journal of Impact Engineering, 2007, 34(6): 1047–1060.DOI:10.1016/j.ijimpeng.2006.02.007 |
| [2] | Xie D, Biggers S B. Strain energy release rate calculation for a moving delamination front of arbitrary shape based on the virtual crack closure technique.part Ⅰ:formulation and validation[J].Engineering Fracture Mechanics, 2006, 73(6): 771–785.DOI:10.1016/j.engfracmech.2005.07.013 |
| [3] | Wimmer G, Schuecker C, Pettermann H E. Numerical simulation of delamination in laminated composite components—a combination of a strength criterion and fracture mechanics[J].Composites Part B:Engineering, 2009, 40(2): 158–165.DOI:10.1016/j.compositesb.2008.10.006 |
| [4] | Sundaram B M, Tippur H V. Dynamic crack growth normal to an interface in bi-layered materials:an experimental study using digital gradient sensing technique[J].Experimental Mechanics, 2016, 56(1): 37–57.DOI:10.1007/s11340-015-0029-x |
| [5] | Ahmed A, Sluys L J. A computational model for the simulation of dynamic fracture in laminated composite plates[J].Journal of Composite Materials, 2015, 49(14): 1717–1738.DOI:10.1177/0021998314539777 |
| [6] | Msekh M A, Sargado J M, Jamshidian M, et al. Abaqus implementation of phase-field model for brittle fracture[J].Computational Materials Science, 2015, 96(Part B): 472–484. |
| [7] | Wang Z, Li P. Dynamic failure and fracture mechanism in alumina ceramics:experimental observations and finite element modelling[J].Ceramics International, 2015, 41(10): 12763–12772.DOI:10.1016/j.ceramint.2015.06.110 |
| [8] | Bagheri R, Fard K M, Monfared M M. Dynamic fracture analysis in an orthotropic half-plane with FGM coating containing several moving cracks[J].Acta Mechanica, 2015, 226(6): 1725–1736.DOI:10.1007/s00707-014-1285-z |
| [9] | Monfared M M, Ayatollahi M. Dynamic stress intensity factors of multiple cracks in an orthotropic strip with FGM coating[J].Engineering Fracture Mechanics, 2013, 109: 45–57.DOI:10.1016/j.engfracmech.2013.07.002 |
| [10] | Manolis G D, Dineva P S, Rangelov T V. Dynamic fracture analysis of a smoothly inhomogeneous plane containing defects by BEM[J].Engineering Analysis with Boundary Elements, 2012, 36(5): 727–737.DOI:10.1016/j.enganabound.2011.11.010 |
| [11] | Lee K H. Analysis of a propagating crack tip in orthotropic functionally graded materials[J].Composites Part B:Engineering, 2016, 84(2): 83–97. |
| [12] | Abdollahifar A, Nami M R. Determination of dynamic stress intensity factor in FGM plates by MLPG method[J].IJST, Transactions of Mechanical Engineering, 2014, 38(M1): 181–194. |
| [13] | Rouzegar S J, Mirzaei M. Modeling dynamic fracture in Kirchhoff plates and shells using the extended finite element method[J].Scientia Iranica, 2013, 20(1): 120–130. |
| [14] | Bui T Q. Extended isogeometric dynamic and static fracture analysis for cracks in piezoelectric materials using NURBS[J].Computer Methods in Applied Mechanics and Engineering, 2015, 295: 470–509.DOI:10.1016/j.cma.2015.07.005 |
| [15] | Zhao W, Hu Z, Zhang X, et al. The dynamic stress intensity factor around the anti-plane crack in an infinite strip functionally graded material under impact loading[J].Theoretical and Applied Fracture Mechanics, 2014, 74: 1–6.DOI:10.1016/j.tafmec.2014.05.002 |
| [16] | Bhardwaj G, Singh I V, Mishra B K, et al. Numerical simulation of functionally graded cracked plates using NURBS based XIGA under different loads and boundary conditions[J].Composite Structures, 2015, 126: 347–359.DOI:10.1016/j.compstruct.2015.02.066 |
| [17] | ?eriba?s S. Static and dynamic analyses of thin uniformly loaded super elliptical FGM plates[J].Mechanics of Advanced Materials and Structures, 2012, 19(5): 323–335.DOI:10.1080/15376494.2010.528160 |
