基于二次插值重构有限元法的动态裂纹扩展模拟

基于二次插值重构有限元法的动态裂纹扩展模拟
作者：王玉光李忠文吴圣川臧晓蕾王燕
来源：《计算机辅助工程》2017年第03期
摘要：
基于二次插值重构有限元法（Twiceinterpolation Finite Element Method， TFEM）分析动态断裂力学问题并进行数值实验，考察TFEM在裂纹动态扩展模拟中的准确性和可靠性.由于TFEM保证节点处梯度场的连续性，因此裂尖附近的应力场可以得到较好的逼近.把该算法成功移植到自主开发的三维裂纹扩展仿真软件（ZonCrack）中.利用ZonCrack进行的裂纹扩展，分析结果表明：TFEM得到裂尖应力强度因子（Stress Intensity Factor， SIF）与解析解基本一致；裂纹扩展的模拟结果与实验值吻合良好.
关键词：
插值重构有限元；断裂力学；疲劳裂纹扩展；裂纹尖端场；损伤容限； ZonCrack软件
中图分类号： TB115；O242.21
文献标志码： B
0 引言
动态疲劳裂纹的模拟分析对工程结构服役安全性和剩余寿命评价有重要的工程指导意义.[1]为模拟裂纹扩展，一般采用节点松弛方法，即在计算中把节点沿着裂缝路径一分为二.[2]由于这种技术适用性差，学者们便又提出一些新型计算方法，如无网格[3]、边界单元法[4]，以及
有效的单元分割策略[5].这些方法虽然可在一定程度上提高计算精度，但计算效率却降低，显然不利于工程应用.
为此，ZHENG等[6]提出一种二次插值重构有限元法（Twiceinterpolation Finite Elements Method， TFEM）求解线弹性断裂力学问题.为提高缺陷体的分析精度，通过2次连续插值得
到形函数，其主要特点是把标准有限元法（Finite Elements Method， FEM）插值函数作为第二阶段的插值参数，或者说把传统FEM不连续的节点梯度（节点导数）作为插值参数，这样，新的形函数便具有连续的节点梯度场.结果发现，TFEM比经典FEM有更好的计算精度和收敛性.
该算法植入新一代基于自适应加密的扩展有限元（eXtend FEM，XFEM）的断裂仿真软件ZonCrack.ZonCrack基于概率疲劳断裂力学理论，对含裂纹工程的结构进行裂纹扩展仿真，获
取无损探伤数据，为工程师提供科学的基础数据和判据.该算法能够预测静载荷或疲劳载荷作用下的开裂行为，确定工程结构损伤容限，不仅可以进行传统有限元强度分析，而且提供疲劳断裂及扩展寿命的分析功能.
1 基于平面三角形的TFEM
选取对复杂几何形状适应性好的三角形单元构建插值函数以求解断裂力学问题.
1.1 一次插值函数
给定一个典型的常应变的单纯形三角形，单元由节点i，j，m逆时针排列构成.插值重构法在二维单纯形单元中的示意见图1.
如图1所示，Si，Sj和Sm为有公共节点i，j和k的关系单元集合，组成支撑域.对u（x）有贡献的节点集（无单元法中一般称为支撑点）正是Si，Sj和Sm所涵盖的全部节点.设qs为支撑节点的位移向量，
1.2 二次插值函数
根据经典FEM得到节点导数后，再进入二次插值获取一个新的形函数
1.3 关于C0节点
从式（16）看出，二次插值得到的试函数在节点处具有C1连续性.这一特点具有提高计算精度的优势，但位于材料界面和位移边界上的节点给求解带来困难.为此，设节点i具有C0连续性，则对于插值点x位于单元e内的节点，可令
分析发现，当问题域内所有节点进行C0连续性处理后，
TFEM型函数将退化为经典FEM型函数，即在某些情况下2种算法可以耦合.也就是说，对于含缺陷体需要更高插值精度的局部区域可采用TFEM，而在远离裂纹前缘区则可采用一般精度FEM，从而实现计算精度和效率的平衡考虑.
2 裂纹扩展计算关键技术
发展TFEM方法用于求解断裂力学问题，并植入ZonCrack软件中.专业断裂仿真软件ZonCrack是一个面向动态断裂过程的仿真实验平台，其以成熟的疲劳断裂力学为理论基础，对含缺陷工程部件进行仿真模拟.[7]为减少对经典FEM源代码的改动，引入节点松弛技术模拟材料的劈开过程.
2.1 应力强度因子
采用区域互作用积分计算混合型裂纹的应力强度因子（Stress Intensity Factor， SIF）解[8]，以最大周向应力判断裂纹方向[9]，则极坐标形式的裂尖应力场分量可有I型和II型裂纹应力强度因子表示
2.2 节点投影线技术
研究发现，网格质量对插值重构有限元法的计算精度影响较小（网格敏感性低），也就是说调整裂纹面附近的网格对SIF的影响也较小.
当前求解断裂力学问题的通用处理方法是确保裂纹面/线与单元边/节点的重合.本文采用节点投影方法准确模拟裂纹的扩展路径，具体操作过程为：网格整体上不调整，仅对裂纹附近节点投影到裂纹面/线上（见图2）.一般来说，裂纹与三角元有3种关系：第一，裂纹与2个边相交，把相交边的公共节点移至裂纹；第二，裂纹与1个边相交，把离裂纹近的节点移至裂纹；第三，裂尖位于单元内，则移动单
元上距裂纹尖端最近的节点.
采用上述处理后，问题域内所有单元边/节点都会位于裂纹面/线上.自此，采用经典节点分离法人为把裂纹面/线上的节点一分为二，并赋予不同节点号.通过这样处理，裂纹扩展路径上就出现一个强不连续性位移，即实现裂纹张开模拟.
3 算例验证
为考察新型二次插值有限元的求解精度和稳定性，定义2范数位移误差和能量误差为
式中：下标e表示精确解或者理论解；下标n表示数值方法参考解或者TFEM解.在TFEM 实现中，单元刚度矩阵积分采用4点Hammer积分方法.
在ZonCrack中调用TFEM求解器分析一些典型的断裂力学问题，以测试TFEM法准确性与合理性.文中所用变量均采用国际单位制.
3.1 单边直裂纹
一个单边直裂纹的无限大平板，两端受均匀拉力，假设平面应变状态.模型的几何及加载见图3，且有H=L=P=1.力学性能参数E=1010，ν=0.25.解析解为
应力强度因子解对比曲线见图4.由此可知，对于不同扩展长度的裂纹，TFEM数值解与理论解和经典FEM解吻合良好，且TFEM法的计算精度高于传统FEM法.TFEM与传统FEM的裂纹尖端应力场分布下见图5.由此可见，TFEM结果与理论解给出的裂纹简单应力场吻合更好，而传统有限元法却存在着明显的应力跳跃现象（即不连续）.
3.2 中心斜裂纹
为进一步考察TFEM求解SIF的准确性与可靠性，以中心斜裂纹板为例研究，见图6.在无限大平板x方向和y方向分别搜到σ1和σ2大小的拉应力，β为斜裂纹与x轴的夹角，显然，KI和KII为关于裂纹倾角β的函数，理论解为
有限元计算选取的计算域为边长L=10的方板，所含裂纹的长度a=0.25，假设轴向拉力
P=1.考虑到裂纹长度远小于矩形板的长度，所以可以视为无限大平板的中心斜裂纹问题.
3.3 4点剪切含裂纹矩形板
4点剪切受载问题的裂纹扩展研究较多，几何模型及受载情况见图7，L=4，b=1，P=1.材料参数为E=200 MPa，ν=0.25，假设平面应变条件.第12扩展步后的裂纹路径见图8.由此可见，插值重构法的计算结果与实验基本吻合.
4 结论
将新型TFEM植入高级三维裂纹扩展仿真与分析系统ZonCrack中用于求解动态断裂问题.通过标准算例分析，得出以下结论.
（1）在ZonCrack中调用TFEM求解器，能够获得与理论解可比的I型、II型及混合型SIF解，且在同一单元网格下，计算精度高于传统FEM.
（2）与传统FEM解相比较，TFEM可在裂尖区域得到更加准确和光滑连续的应力分布.
（3）TFEM法预测的裂纹路径与实测解基本一致，表明TFEM算法的可靠性.
必须指出，与传统有限元比较，虽然插值重构有限元法使用经典FEM的导数进行二次插值，但并未因此显著增加问题的总求解自由度.
参考文献：
[1] 杨小彬，庄茁，庄传晶. 富气输送管道裂纹动态扩展的数值模拟[J]. 清华大学学报（自然科学版）， 2008， 48（8）： 13551358. DOI： 10.3321/j.issn：10000054.2008.08.032.
YANG X B， ZHUANG Z， ZHUANG C J. Numerical study of crack propagation in rich gas transporting pipelines[J]. Journal of Tsinghua University （Science and Technology）， 2008， 48（8）： 13551358. DOI： 10.3321/j.issn：10000054.2008.08.032.
[2] 解德，钱勤，李长安. 断裂力学中的数值计算方法及工程应用[M]. 北京：科学出版社， 2009.
[3] BELYTSCHKO T， LU Y Y， GU L. Crack propagation by elementfree Galerkin methods[J]. Engineering Fracture Mechanics， 1995， 51（2）： 295315. DOI：
10.1016/00137944（94）001539.
[4] 李俊，冯伟哲，高效伟. 一种基于直接计算高阶奇异积分的断裂力学双边界积分方程分析法[J]. 力学学报， 2016， 48（2）： 387398.DOI： 10.6052/0459187915342.
LI J， FENG W Z， GAO X W. A dual boundary integral equation method based on direct evaluation of higher order singular integral for crack problems[J]. Chinese Journal of Theoretical and Applied Mechanics， 2006， 48（2）： 387398.DOI： 10.6052/0459187915342.
[5] MELENK J M， BABUS K. The partition of unity finite element method： basic theory and applications[J]. Computer Methods in Applied Mechanics and Engineering， 1996， 139（14）：289314. DOI： 10.1016/S0*******（96）010870.
[6] ZHENG C， WU S C， TANG X H， et al. A novel twice interpolation finite element method for solid mechanics problems[J]. Acta Mechanica Sinica， 2010， 26（2）： 265278. DOI： 10.1007/s1040900902653.
[7] 吴圣川，吴玉程. 断裂分析及CAE软件的现状与发展[J]. 计算机辅助工程， 2011， 20（1）： 12.
WU S C， WU Y C. Current status and development of fracture mechanics based CAE software[J]. Computer Aided Engineering， 2011， 20（1）： 12.
[8] BUDIANSKY B， RICE J R. Conservation laws and energyrelease rates[J]. Journal of Applied Mechanics， 1973， 40（2）： 201203. DOI： 10.1115/1.3422926.
[9] FLEMING M， CHU Y A， MORAN B， et al. Enriched elementfree Galerkin methods for cracktip fields[J]. International Journal for Numerical Methods in Engineering， 1997， 40（8）：14831504. DOI： 10.1002/（SICI）10970207（19970430）40：83.0.CO；26.
（编辑武晓英）。