ABAQUS User Subroutine Development for Energy Releasing Rate of the Functionally Graded Plate with Cracks
近年来, 人们对含裂纹功能材料体在碰撞事件中受到冲击载荷作用时的响应日益关注.为了精确评价动态载荷作用下结构的断裂力学问题, 提出了动态断裂参数的概念, 如动态应力强度因子、动态能量释放率等.对于简单的几何模型、理想的材料模型和特殊加载模型的动态断裂参数可由解析法求得.但往往由于结构或边界条件比较复杂, 解析法并不可行, 实验测试又费时费钱, 数值计算就成为研究该类问题的有力工具.
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) |
式中:M(z)代表功能梯度材料物参; Mh代表功能梯度材料结构底部的物参; Mb代表功能梯度材料结构顶部的物参; φc为底部材料的体积分数; n为分布规律的形状因子.图 1给出了不同形状因子下, 无量纲厚度z/h和体积分数的关系, 当n=0时, 该材料为均匀材料, n>0时, 该材料为功能梯度材料，z为沿z轴坐标的长度, h为板宽.
2 非均匀有限元法
非均匀有限元方程为
|
(3) |
式中:K为刚度矩阵, 由单元刚度矩阵Ke(e=1, 2, …, n, n是离散单元的个数)集成; F为载荷列阵, 由单元载荷列集成; U为位移列阵.
单元刚度矩阵Ke的表达式为
|
(4) |
式中:Ωe为单元域; Be为单元应变矩阵; De(x)为与坐标x=[x, y, z]有关的单元材料弹性矩阵.
式(4) 采用高斯积分求解, 二维问题可写成
|
(5) |
式中:Jij为雅可比行列式; wi和wj为二维高斯积分的权系数.De(x)可由插值求得, 其弹性模量E和泊松比ν可表达为
|
(6) |
|
(7) |
该部分在ABAQUS中由UMAT编写完成.
3 虚拟裂纹闭合法
计算平面功能梯度材料内任意线状裂纹的应变能释放率, 可采用哑节点断裂单元, 如图 2所示.断裂单元中裂尖的节点1和节点2坐标重合; 节点3和节点4坐标重合; 节点5在裂尖前面, Δa为虚拟裂纹增量, Δc为裂尖后单元尺寸, 取Δc=Δa.
节点位移矢量为
|
(8) |
虚拟裂纹的增量Δa为
|
(9) |
式中, (x1, y1)和(x5, y5)分别为整体坐标(x, y)下节点1和节点5的坐标.
x和x′轴的夹角为
|
(10) |
|
(11) |
在节点1和节点2处放置特殊刚度的弹簧, 位移矢量如图 3所示, 裂尖处节点力和张开位移在局部坐标系(x′, y′)下为
|
(12) |
|
(13) |
|
(14) |
|
(15) |
式中:
|
(16) |
|
(17) |
|
(18) |
|
(19) |
式中, Kx, Ky分别为x和y方向的弹簧刚度.
应变能释放率为
|
(20) |
|
(21) |
式中, B为裂纹体的厚度.
动态载荷作用下的固定裂纹问题的大量实验结果表明无需在以上求解公式中增加惯性项或在裂尖处施加特殊质量的质点[1].静态的虚拟裂纹闭合法计算公式完全适用于计算动态载荷作用下的固定裂纹问题, 计算过程中注意各个物理单位的匹配.
动态应力强度因子和动态能量释放率的关系为
|
(22) |
|
(23) |
式中:KⅠdyn(t)和KⅡdyn(t)分别为Ⅰ型和Ⅱ型动态应力强度因子; GⅠdyn(t)和GⅡdyn(t)分别为Ⅰ型和Ⅱ型动态应变能释放率.平面应力状态, E=Etip; 平面应变状态, E=Etip (1－(νtip)2), Etip和νtip分别为裂尖处材料的弹性模量和泊松比.该部分在ABAQUS中由UEL编写完成.
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求解动态能量释放率后转化成动态应力强度因子, 并进行无量纲化处理, 即
|
(24) |
图 6给出了形状因子n=0.0, 0.2, 0.5, 1.0和5.0时, 采用本文方法计算求解所得的含中心裂纹功能梯度结构的无量纲动态应力强度因子KⅠ.当n=0时, 材料梯度无变化, 梯度板的材料仅为材料Y-TZP, 与有限差分法(FDM)和文献[7]的计算结果作对比, 从图中可以看出, n=0时, 本文方法与文献[1]的计算结果基本相同, 与FDM计算结果也吻合很好, 从而证明了本文方法的正确性.可见本文方法具有求解精度高、计算方便、裂尖单元无需特殊处理等优点.当形状因子n由0增加到5.0, KⅠ值变化不大, 是由于功能梯度材料的组分材料Y-TZP和γ-TiAl弹性模量相差不大和形状因子变化范围小导致的.图 7给出了不同方法、不同时间步时中心裂纹功能梯度板的变形图, 从而验证了本文方法的正确性与可靠性.
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, 无阻尼.
将结构均匀划分成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给出了不同方法、不同时间步时倾斜裂纹功能梯度板的变形图, 可见本文方法是正确有效的.
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 |
2017, Vol. 38
