研究疲劳载荷下的裂纹扩展规律对于结构安全性评价非常重要。国内外学者已经提出了很多的方法来计算裂纹扩展问题，比如边界元法[1]、无网格方法[2]、有限元网格更新法[3, 4]等，这些方法在求解裂纹扩展问题时都有良好的表现。基于有限元方法的网格更新法虽然适应性强，但是需要随着裂纹的扩展不断更新裂纹周围的网格，工作量大，计算效率较低。扩展有限元方法[5, 6]基于传统有限元方法，并且克服了网格更新方法的缺点，通过扩充位移项来描述不连续的位移场，使裂纹独立于网格存在，从而能够模拟任意形状的裂纹。此外，由于扩展有限元方法仅仅在有限元法的基础上扩充了位移项，可以通过对现有商业有限元软件进行二次开发实现其功能，一定程度上减少了工作量。
基于商业软件的扩展有限元功能开发是普遍采用的一种方法，其优势是可以省去编写求解平衡方程程序的过程。ABAQUS从6.9版本开始内置了扩展有限元方法，但是只考虑了单元被裂纹完全贯穿的情况，而且在使用上存在诸多限制[7]。方修君等[8]在ABAQUS环境下开发了扩展有限元功能，并对三点弯曲梁的开裂过程进行了模拟，但是文章只研究了Ⅰ型开裂问题，且没有考虑裂尖位于单元内部的情况。谢海等[9]解决了ABAQUS环境下求解混合型裂纹的扩展有限元方法，但是没有实现裂纹扩展问题的模拟。Giner等[10]使用用户子程序功将XFEM方法嵌入到ABAQUS平台中，实现了对多裂纹扩展的模拟。
本文推导并编程实现了基于XFEM的多裂纹扩展模拟，改进了水平集更新以及裂尖单元判断算法，解决了可视化后处理问题。
1 模拟裂纹扩展的扩展有限元方法 1.1 含裂纹区域的位移场形式在XFEM中，对于包含裂纹的区域，通过添加扩充项描述包含间断的位移场[6]
式中，ui为普通节点自由度，n为所有节点，nJ是属于被单元完全贯穿的单元的节点集，称为阶跃扩充节点，nT是属于被裂尖刺中的单元的节点集，称为裂尖扩充节点，扩充节点和裂纹的位置关系如图 1所示。ai是阶跃扩充节点的附加自由度(扩充函数H(x)在裂纹一侧取+1，另一侧取-1)，biα是裂尖扩充节点的扩充自由度，其扩充函数Ψα(x)的形式为
 |
| 图 1 单元节点扩充类型 |
根据(1)式定义的位移场，单元刚度矩阵K的形式为
式中，B为单元几何矩阵，对于普通节点 对于阶跃扩充节点，几何矩阵附加项Bi,jump=(H(x)-H(xi))Bi,c，对于裂尖扩充节点，几何矩阵附加项为 1.3 裂纹扩展的水平集表征
使用扩展有限元方法模拟断裂问题，不需要直接对裂纹进行建模，通常使用水平集方法[11]将单元内部的裂纹信息以符号距离的形式离散到节点上，然后在后处理中通过形函数插值还原裂纹信息。
区域内的一条裂纹的几何形状由切向水平集函数φ和法向水平集函数ψ确定，如图 2所示裂纹，它的形状可以描述为：
 |
| 图 2 裂纹的水平集函数表征 |
对于内嵌裂纹(2个端点)，需要使用2个法向水平集函数ψ1和ψ2来确定裂尖的位置。
对于多裂纹情况，为了节省计算空间，为区域内的每条裂纹分配一个影响区域，对处在区域内的单元的水平集函数值进行初始化。如果裂纹位置比较靠近，其影响区域相互叠加，则该区域内的单元节点需要分别计算相对于两条或更多条裂纹的水平集函数值，如图 3所示。
 |
| 图 3 多裂纹情况的水平集函数分布 |
对于裂纹经过的单元，计算单元节点的水平集函数值φi和ψi，则通过形函数插值可以得到单元内部水平集函数的表达式：
式中，n为单元节点数。裂纹扩展后，对处于裂纹影响区域内的单元的水平集函数进行更新，即可通过插值还原扩展后的裂纹几何形状。水平集的更新可以通过求解(8)式得到[7] 式中，v为裂尖扩展速率，Δt为时间增量，Ωupdate为切向水平集函数的更新域，其范围如图 4所示，取为裂纹扩展方向的一个矩形区域，其垂直于裂纹增量方向的边长为裂尖单元边长的4倍，沿着裂纹增量方向的边长为裂纹增量的2~3倍，这样处理可以减少更新切向水平集函数φ的计算量。
 |
| 图 4 切向水平集更新区域 |
使用单元节点的水平集值可以方便地判断间断单元的类型：对于完全贯穿单元，满足φminφmax>0,ψmax < 0；对于裂尖单元，满足ψminψmax < 0,ψminψmax < 0 。其中，φmin、φmax、ψmin、ψmax分别为单元节点的最小和最大水平集函数值。
实际问题中，简单的使用上述方法并不能准确判断单元类型，注意图 5所示的情况，图 5中的灰色单元满足裂尖单元判据，但不是裂尖单元，单元类型判断的错误会赋予节点错误的附加自由度，从而影响模拟的精度，需要增加额外的判断条件来提高方法的稳定性。
 |
| 图 5 单元扩充类型判断 |
对于初步判定为裂尖扩充形式的单元，判断裂尖点与单元的位置关系，若包含裂尖，则可以判定为裂尖单元。所以裂尖单元的判定标准可以定义为(10)式。
式中，xtip为裂尖位置，Ωele为单元区域。
由于包含裂纹信息的单元内部不连续，使用Gauss积分时，无法得到精确解，一般采用划分积分子区域的方法进行积分。单元与裂纹可能存在4种交汇关系，针对不同的交汇关系对单元进行分割，然后分别对单元的各子区域积分，如图 6所示。裂纹与单元的交点可以由(11)式得到
 |
| 图 6 四节点单元分割方案 |
基于J积分的交互积分方法[12]被广泛用于混合型裂纹的应力强度因子的求解。交互积分在真实场的基础上引入辅助场，辅助场可以是任意满足变形方程和本构方程的应力-变形场，在求解断裂问题时取为裂尖奇异场的形式。则交互积分的表达式为
式中，A为积分区域，ui和uiaux分别为真实位移场和辅助位移场，σij和σijaux分别为真实和辅助应力场，W=σijεijaux=σijauxεij。q为积分区域的权重函数，在积分域内取1，在积分域外取0。
交互积分和应力强度因子的关系如(13)式所示
分别取KⅠaux和KⅡaux为特定值即可得到应力强度因子的值:
对于复合型裂纹，裂纹扩展方向可以由最大周向拉应力准则计算
1.6 程序结构
有别于传统的在ABAQUS中嵌入XFEM功能的方案，本文以模块化的形式编写裂纹扩展分析程序：使用FORTRAN分别编写了模型创建、裂纹初始化、水平集更新、断裂参数计算以及可视化处理(借助商业软件Tecplot，实现了裂纹扩展过程的动态显示)等模块并留出接口，方便不同裂纹扩展问题的创建以及后续功能的开发。图 7为程序的详细流程图。
 |
| 图 7 裂纹扩展分析流程图 |
为验证上述方法的计算精度，计算中心斜裂纹平板在不同裂纹倾角下的应力强度因子。模型构形如图 8所示。
 |
| 图 8 中心斜裂纹平板 |
受垂直方向均匀拉伸作用，裂纹长度与板宽满足a/W=0.8，板长H=2W。裂纹与载荷夹角为θ，取其值为0°~90°，计算裂尖的无量纲应力强度因子并与应力强度因子手册结果[13]进行比较。
2.1 网格密度的影响为了研究网格密度对结果的影响，考虑疏密程度不同的网格划分，网格类型为均匀的四边形网格，单元数目分别为1 250(粗网格)、2 450(普通网格)和5 000(细网格)。图 9为不同单元数目下的无量纲应力强度因子(
)与角度θ的关系曲线，由于对称性，本文只给出一个裂尖的结果。
 |
| 图 9 不同网格密度下的无量纲应力强度因子 |
从结果可以看出，当网格较粗时(网格数为1 250)，角度80°~85°时KⅠ值计算误差较大，这是因为此时裂尖离模型边缘很近，从裂尖到模型边缘之间只存在1~2层单元，计算交互积分时积分域内单元数过少，导致了计算精度不够，可以通过加密模型边缘网格的方法解决。除此之外，应力强度因子的计算解与理论解吻合较好，表明方法受网格密度影响较小。
2.2 积分域尺寸的影响为了考察交互积分的积分域半径对计算精度的影响，对网格数5 000的模型分别取积分半径为2~4倍的特征长度(裂尖单元面积的平方根)，以取4倍特征长度计算的结果为基准，计算它们的相对误差，结果如图 10所示，其中R为积分域半径，H为单元特征长度。
 |
| 图 10 不同积分半径下的无量纲应力强度因子误差 |
从结果可以看出，取3种积分半径计算应力强度因子的结果相对误差非常小，同时可以说明，取2倍的特征长度已经可以取得不错的计算精度。
3 裂纹扩展模拟 3.1 带孔板边缘裂纹扩展模拟对图 11所示的边缘裂纹模型进行模拟，图中长度单位均为mm，预制初始裂纹为略高于中心线的边缘裂纹，其尺寸a0=10mm。模型材料为7075-T6,弹性模量E=7.17×104N/mm2，泊松比υ=0.33,集中载荷P=20 kN，用于对比的数值和试验结果取自文献[16]。
 |
| 图 11 带孔平板裂纹扩展 |
采用控制裂纹增量的方式计算裂纹扩展，取Δa=3 mm，当裂纹穿入圆孔时停止计算，计算结果如图 11b)所示，图 11c)为后处理云图结果。我们将计算得到的裂纹路径点按照坐标绘制出来，并与文献中的计算路径作了对比，两者结果与试验吻合较好，本文结果略优于文献中的模拟结果。
3.2 多裂纹扩展模拟考虑图 12所示的双边裂纹扩展模型，两孔关于模型中心对称分布，初始裂纹对于模型中心对称分布于模型两侧，长度为1 mm。模拟裂纹扩展路径并与更新网格方法(Bouchard等[6]和Khoei[7])进行对比。矩形板的弹性模量E=2×105 N/mm2，泊松比v=0.3。对模型下缘进行固支，对上缘施加等位移加载，同时约束上缘节点水平方向的自由度，使它们在水平方向不能收缩。我们在程序中可以使用裂纹增量或者应力循环数作为计算裂纹的控制参数，对于此算例，由于对称性，使用控制裂纹增量的形式模拟裂纹扩展。取Δa=1 mm，扩展路径如图 13所示。
 |
| 图 12 孔边裂纹扩展模型 |
 |
| 图 13 裂纹扩展路径 |
从裂纹扩展路径可以看出，我们计算的结果与Khoei等人计算的结果吻合较好并且2条裂纹的扩展路径对称分布，而Bouchard等人计算的结果右侧裂纹出现了较大的拐点，三者裂纹扩展的走势基本相符。裂纹刚开始朝孔的方向扩展，然后由于2条裂纹间的相互影响，裂纹开始偏离圆孔并相向扩展，最终裂纹在模型中部交错和汇合。
为了研究裂纹与孔的距离对裂纹扩展路径的影响，我们将初始裂纹的位置向孔靠近(距离模型端面分别为3.15 mm、3.45 mm)，对这2种情况进行模拟并与原来的情况进行对比，由于使用XFEM模拟裂纹不依赖与网格，因此相比于更新网格法，可以方便地改变初始裂纹的位置。模拟结果如图 14所示。
 |
| 图 14 裂纹与模型边缘距离不同时的扩展路径 |
从图中可以看出，当裂纹离孔较近时会直接向孔靠近并最终刺入孔中；当裂纹离孔较远时，会呈现图 13中讨论的裂纹路径形态；当裂纹与孔的距离为
某一值时(图 14中实线的裂纹)，裂纹在扩展到孔边缘时，会受对面裂纹的影响而偏离孔，并最终形成同一水平面上相向扩展的裂纹。以上关于初始裂纹位置对扩展路径影响的模拟和讨论虽然没有试验对比或文献可资参考，但是其扩展规律的合理性是明显的。
4 结 论本文基于XFEM基本理论，利用商业软件ABAQUS和 Tecplot的程序接口，编写了用于多裂纹扩展模拟的Fortran程序，在传统XFEM方法上进行了部分的优化：1. 使用局部水平集更新算法，减小了计算量；2. 改进了裂尖单元的判断准则，提高了方法的稳定性。3. 编写了可视化处理模块，弥补了ABAQUS不能显示用户单元的缺陷，方便观测含裂纹体的应力分布状态。数值算例证明了本文的程序计算结果精度较高且适应性强，与更新网格方法的对比证明使用XFEM可以灵活地调整裂纹位置，方便研究不同位置的裂纹扩展规律。
| [1] | Portela A, Aliabadi M H, Rooke D P. The Dual Boundary Element Method: Effective Implementation for Crack Problems[J]. International Journal for Numerical Methods in Engineering, 1992, 33(6): 1269-1287 |
| Click to display the text | |
| [2] | Belytschko T, Gu L, Lu Y Y. Fracture and Crack Growth by Element Free Galerkin Methods[J]. Modelling and Simulation in Materials Science and Engineering, 1994, 2(3A): 519-534 |
| Click to display the text | |
| [3] | Bouchard P O, Bay F, Chastel Y. Numerical Modelling of Crack Propagation:Automatic Remeshing and Comparison of Different Criteria[J]. Computer Methods in Applied Mechanics and Engineering, 2003, 192(35): 3887-3908 |
| Click to display the text | |
| [4] | Khoei A R, Azadi H, Moslemi H. Modeling of Crack Propagation via an Automatic Adaptive Mesh Refinement Based on Modified Superconvergent Patch Recovery Technique[J]. Engineering Fracture Mechanics, 2008, 75(10): 2921-2945 |
| Click to display the text | |
| [5] | Belytschko T, Black T. Elastic Crack Growth in Finite Elements with Minimal Remeshing[J]. International Journal for Numerical Methods in Engineering, 1999, 45(5): 601-620 |
| Click to display the text | |
| [6] | Dolbow J, Belytschko T. A Finite Element Method for Crack Growth without Remeshing[J]. Int J Numer Meth Engng, 1999, 46: 131-150 |
| Click to display the text | |
| [7] | Pais M J. Variable Amplitude Fatigue Analysis Using Surrogate Models and Exact XFEM Reanalysis[D]. University of Florida, 2011 |
| Click to display the text | |
| [8] | 方修君, 金峰. 基于ABAQUS平台的扩展有限元法[J]. 工程力学, 2007, 24(7): 6-10 Fang Xiujun, Jin Feng. Extended Finite Element Method Based on ABAQUS[J]. Engineering Mechanics, 2007, 24(7): 6-10 (in Chinese) |
| Cited By in Cnki (105) | Click to display the text | |
| [9] | 谢海, 冯淼林. 扩展有限元的 ABAQUS 用户子程序实现[J]. 上海交通大学学报, 2009 (10): 1644-1648 Xie Hai, Feng Miaolin. Implementation of Extended Finite Element Method with ABAQUS User Subroutines[J]. Journal of Shanghai Jiaotong University, 2009 (10): 1644-1648 (in Chinese) |
| Cited By in Cnki (11) | Click to display the text | |
| [10] | Giner E, Sukumar N, Tarancón J E, et al. An Abaqus Implementation of the Extended Finite Element Method[J]. Engineering Fracture Mechanics, 2009, 76(3): 347-368 |
| Click to display the text | |
| [11] | Stolarska M, Chopp D L, Mos N, et al. Modelling Crack Growth by Level Sets in the Extended Finite Element Method[J]. International Journal for Numerical Methods in Engineering, 2001, 51(8): 943-960 |
| Click to display the text | |
| [12] | Rao B N, Rahman S. An Interaction Integral Method for Analysis of Cracks in Orthotropic Functionally Graded Materials[J]. Computational Mechanics, 2003, 32(1/2): 40-51 |
| Click to display the text | |
| [13] | Sih G C. Handbook of Stress-Intensity Factors[M]. Lehigh University, Institute of Fracture and Solid Mechanics, 1973 |