基于ABAQUS 平台考虑T 应力的I 型裂纹扩展模拟开发

杨立云，韦 鹏，王青成，陈美霞，杨登辉

(中国矿业大学(北京)，力学与建筑工程学院，北京 100083)

随着社会生产的发展，新工艺、新技术、高强度材料得到广泛应用，工程部件和结构在服役期间或由于制造缺陷会产生微裂纹，之后由于应力集中、疲劳、腐蚀等因素的影响，微裂纹不断成核、扩展、生长为宏观裂纹，最终导致材料断裂失效。断裂准则是含裂纹材料在各种应力作用下断裂和破坏的重要判据，解释了裂纹何时起裂和判断裂纹扩展方向的问题，因此在预测裂纹扩展路径方面得到广泛应用。传统线弹性断裂力学中常用的MTS[1]和SED[2]准则仅采用应力强度因子K表征裂纹尖端应力集中的程度[3]，但经过研究发现[4]，当裂纹受到高度约束时，使用单个参数(即临界应力强度因子)来表征线弹性材料的脆性断裂韧度是有缺陷的。SUN 等[5]研究发现边界约束能显著改变裂纹尖端应力场，尤其是由奇异应力场(K-主导区)主导的区域大小，随着K-主导区尺寸的减小，奇异应力可能不足以描述裂纹的起裂，需要引入一个新的参数来表征裂纹尖端附近非奇异应力场对断裂的影响。

国外研究者LIM 等[6]和KHAN 等[7]在分析线弹性材料的脆性断裂时发现，不同几何构型的试件通过断裂实验得到的材料断裂韧度差别很大，而且应用传统的断裂准则计算的裂纹起裂角和扩展角与试验结果有一定的偏差，传统断裂力学理论无法对这一现象作出解释。CHAO 等[8]指出I 型裂纹扩展发生偏转是由于裂纹尖端约束效应的作用，即试件厚度、平面几何形状和加载方式等与材料的断裂行为的相关性，并用裂纹尖端Williams级数展开式的第2 项-平行于裂纹方向的常数项即T应力来表征约束效应，发现其对裂纹起裂和扩展的影响显著。研究结果表明[8-9]：经过T应力修正后的断裂准则能更准确地预测裂纹的偏转角和扩展路径。

II 型裂纹与I-II 复合型裂纹的扩展相对于I 型裂纹更加复杂，需要确定裂纹起裂角度以及方向。对于II 型裂纹，T应力对裂纹的断裂韧度有较大影响，SMITH 等[10]基于理论与实验研究了T应力对线弹性脆性材料II 型断裂韧度的影响，结果表明：较大的负T应力使II 型断裂韧度增加，T应力和断裂韧性之间几乎呈线性反比关系。

对于I-II 复合型裂纹，随着T应力的增大，裂纹扩展角逐渐变大。高文和王生楠[11]将理论预测结果和试验数据进行比较得出结论：负T应力可以提高裂纹扩展阻力，正T应力则可以降低材料的裂纹扩展阻力。

ABAQUS 是目前国际上通用的非线性有限元分析软件之一，它提供了基于常规有限元法或扩展有限元法的断裂力学参数计算方法，采用ABAQUS软件内置的交互积分方法[12]可以计算裂纹尖端的应力强度因子KI、KII，J积分和T应力，在模拟裂纹起裂和扩展问题时具有可靠性。ABAQUS 中模拟裂纹扩展路径的方法主要有：1) 增量裂纹扩展法(The incremental crack growth)[13]是在利用传统有限元法求解断裂力学问题时，将裂纹尖端邻域网格单元进行细化，并设置为奇异性单元，模拟裂纹扩展时还需要对网格进行重划分，耗时且效率低；2)扩展有限元法 (Extended finite element method, XFEM)[14]消除了裂纹对网格的依赖，可以自动模拟裂纹扩展路径，但不能准确得到应力强度因子等断裂力学参数；3)虚拟裂纹闭合法(Virtual crack closure technique, VCCT)[15]基于线弹性断裂力学，具有对有限元网格尺寸不敏感和无需对裂纹尖端应力单元特殊处理等优点，分析相对稳定，但同样无法求取准确的裂尖应力参数；4)离散元法[16](Discrete element method, DEM)作为分析不连续介质问题的数值方法，通过将颗粒相互黏结形成的不连续黏结体来模拟各种材料，可以较为真实地模拟非线性变形和裂纹扩展特征；5)近场动力学法[17](Peridynamics, PD)是一种基于非局部作用思想并以积分形式为运动方程的数值方法，解决了传统连续介质力学在处理不连续问题时的奇异性问题；6)相场法[18](Phase field method,PDM)相对于描述离散裂纹的计算裂纹方法，无需额外的裂纹起裂准则，可以描述从无裂纹到裂纹起裂扩展的全过程。

ABAQUS 为用户提供了两种二次开发接口：1)使用Fortran 语言二次开发本构模型或单元[19-20]；2)使用Python 语言二次开发ABAQUS 脚本接口。很多国内外学者利用Python 语言对ABAQUS 进行二次开发，实现裂纹自动扩展数值模拟，SONG等[21]采用子模型技术二次开发ABAQUS 软件，成功模拟了混凝土裂缝扩展过程；李强等[22]通过ABAQUS 脚本接口进行二次开发，对闭口斜裂纹的分支裂缝的扩展过程进行了数值模拟；WANG等[23]基于Python 程序进行ABAQUS 预处理二次开发，建立了中尺度轻骨料混凝土有限元模型，模拟结构中离散裂纹的起裂和扩展行为；王海鹏和吴永东[24]利用Python 语言编写了裂纹自动扩展的程序模块；竹宇波[25]利用围道积分及网格重剖分技术，实现了裂纹自由扩展数值模拟功能，并对SCB 模型进行了应力强度因子及裂纹扩展路径数值仿真。关于ABAQUS 裂纹扩展模块的开发主要是采用最大周向应力准则(MTS)或最小应变能密度因子准则(SED)判断裂纹的扩展方向，却没有考虑T应力的影响。

基于此，本文采用Python 语言编写基于修正最大周向应力和最小应变能密度因子准则的程序模块，将T应力引入裂纹起裂角的计算中，研究T应力对起裂角度的影响，并利用编写的程序进行裂纹扩展路径模拟，与相关文献中含I 型裂纹的实验结果进行对比，验证程序的有效性。

1 理论分析

1.1 修正的最大周向应力准则

ERDOGAN 等[1]根据纯II 型裂纹扩展实验结果建立了最大周向应力准则，即在脆性材料中，裂纹将沿着裂尖周围的最大周向应力方向θ0开裂。CHAO 等[8-9]提出了经过T应力修正的最大周向应力准则，修正后的准则考虑了K和常数项T应力的影响。GMTS(Generalized maximum tangential stress)准则可以表示如下：

求解式(1)和式(2)，两个可能的结果可以写为：

式中：Tc、Kc分别为断裂时的临界T应力和应力强度因子；rc为裂纹尖端的临界距离或断裂过程区域的半径。

裂纹一旦发生弯曲，裂纹模式将由I 型转变为I-II 复合型，此时应力强度因子KI、KII为应力场的主项，因此不考虑高阶项T应力的影响。根据MTS 准则并且仅考虑Williams 级数展开的第一奇异项时，裂纹扩展角θ0可以通过以下公式求得：

1.2 修正的最小应变能密度准则

SIH[2]于1973 年提出了基于应变能密度场的断裂概念，得到了应变能密度因子准则，该准则在预测含裂纹结构的疲劳和断裂方面得到广泛应用。SED 准则认为当以裂尖为圆心、临界距离rc为半径的小圆周上的应变能密度因子达到临界值Sc时裂纹失稳扩展，裂纹沿最小的应变能密度因子方向扩展，记为θm。rc和Sc均为材料常数。SED 准则可表示为：

式中，S是应变能密度因子，表达式为：

式中： dW/dV为单位面积的应变能；rc为临界半径；ν、E分别为材料的泊松比和弹性模量；σx、σy、σz为空间直角坐标系中的正应力分量；τxy、τxz、τyz为空间直角坐标系中的剪切应力分量。将包含常数项T应力的裂纹尖端应力场表达式代入式(5)～式(7)中，得到GSED(Generalized Minimum Strain Energy Density)准则，经求导化简后其表达式为：

其中：

2 裂纹扩展模拟方法

ABAQUS 作为常用的非线性有限元分析软件之一，其内置程序为二次开发用户提供了丰富的数据库，通过Python 语言来调用这些数据库可以实现从建模到后处理等过程的全程控制，还可以开发用户自定义界面以实现交互式操作等功能。

2.1 裂纹自动扩展程序流程

利用Python 语言对ABAQUS 进行二次开发，编写裂纹自动扩展程序。设置裂纹扩展增量，自动提取新裂尖坐标，删除上一步模型的网格，根据ABAQUS 软件内置的交互积分方法[12]得到的裂尖应力强度因子和T应力计算裂纹起裂和扩展角度，自动生成新的模型以及裂尖网格细化，重复上述步骤，直至裂纹扩展到模型边界，完成自动扩展程序模块的编写。程序运行流程图如图1 所示，图中Cpd 为内置程序计算的裂纹起裂角和扩展角，Seam 为指定裂缝，Remeshing 为重新划分网格。

裂纹自动扩展程序功能主要包含以下几个部分：应力分析及裂纹扩展检测、ODB 文件数据收集、更新模块、网格重划分、将结果写入指定文件。建立模型，定义材料，施加荷载和边界条件等裂纹体前处理操作均由ABAQUS/CAE 和ABAQUS/Standard 完成。运行脚本程序并读取模型信息，自动获取裂纹及其尖端位置，建立分析步、划分网格后提交计算，根据计算结果和断裂准则判断裂纹是否扩展，如果不扩展则终止程序；如果扩展则根据断裂准则模块计算的角度和给定的扩展增量得到新裂纹尖端的坐标，返回定义裂纹的步骤并重新划分网格，往复循环，直至裂纹扩展到模型边界或止裂时程序终止。通过导入CAE 文件建立初始模型进行分析，是为了拓宽程序模块的使用范围，例如对含中心斜裂纹、多条裂纹等试件的模拟，只需提前在ABAQUS 中建立前处理模型，再运行脚本程序即可实现对其裂纹扩展路径的分析。

2.2 粒子群算法

考虑T应力对裂尖应力场的影响，对最大MTS或最小SED 准则进行修正，需要利用Python 语言将修正后的准则编入裂纹自动扩展程序脚本，判断裂纹是否起裂。这涉及利用Python 语言求解三角函数方程的问题，由于式(8)所示的方程没有解析解，Python 中的科学计算库Sympy 模块无法解决，本文通过Python 调用Matlab 中的函数解决方程问题，再将计算结果返回原Python 程序，此种数据转换方法的运算效率非常低。

1997 年KENNEDY[26]根据对鸟群觅食行为的研究和模拟，建立了基于群体智能的全局随机搜索算法——粒子群算法(PSO)，主要用于解决最优化问题，PSO 算法可以描述为n个粒子在N 维搜索空间里，通过不断更新速度和位置来寻找最优解的过程。在N 维空间中有n个粒子，用vi为粒子的速度，xi为粒子的位置，速度和位置更新公式如下：

速度更新公式为：

位置更新公式为：

式中：w为惯性权重(有速度就有运动惯性)；pbest为本粒子历史上最好的位置；gbest 为种群中所有粒子当前最好的位置；c1和c2为学习因子，分别向本粒子历史最好位置和种群当前最好位置进行学习；r1和r2为两个随机函数，取值范围[0, 1]，以增加更新速度的随机性。

PSO 算法具有较强的函数极值求解能力，其算法流程如图2 所示，首先定义粒子速度和位置的取值范围并随机初始化每个粒子，根据适应度函数(需要求解极值的函数方程)计算粒子适应度值，通过粒子间的相互搜索、学习得到个体极值和群体极值；利用式(10)～式(11)更新粒子的速度和位置，并相应地更新个体极值和群体极值，反复迭代，直至输出全局最优结果并结束程序。

图2 粒子群算法流程图Fig.2 Flow chart of particle swarm algorithm

为了求解式(8)所示的三角函数方程，本文适应度函数应为其微分后得到的原函数，适应度值为其函数值，种群粒子数为20，每个粒子的维数为2，算法迭代进化次数为100，输入相关参数后利用Python 编写的PSO 算法进行计算，得到的关于原函数的最优解(极值点)即为三角函数方程的数值解，经过换算可得裂纹的起裂角或扩展角。

2.3 程序编制

利用Python 将修正后的MTS 或SED 准则编入裂纹自动扩展程序脚本，基于智能群体全局随机搜索的粒子群算法，计算裂纹起裂角度和扩展角度，部分程序如下：

3 应用验证

AYATOLLAHI 等[13]在室温静态荷载下利用含I 型裂纹的双悬臂梁(Double cantilever beam, DCB)试件进行了断裂试验并给出了扩展路径，试件模型如图3 所示，其材质为聚甲基丙烯酸甲酯(PMMA 或有机玻璃)，该试件模型的断裂韧度为1.3 MPa·m0.5。本节将利用ABAQUS 模拟相同条件下的裂纹扩展路径，使用增量裂纹扩展法自动迭代处理，并将模拟结果与文献中两种尺寸的DCB 试样实验结果进行对比，以此验证程序的有效性。

图3 试样几何形状Fig.3 Specimen geometry

在文献[13]中初始裂纹长度与模型长度比值a/W设定为0.5，弹性模量E=2.9 GPa，泊松比ν=0.35，临界半径rc=0.1 mm，模型其他参数见表1。研究发现[27]裂纹偏转具有随机性，可能会向上或者向下偏转，如图3 所示。由于参考的文献实验结果是向下扩展，所以本文模拟选取向下的偏转角度。

表1 模型相关信息[13]Table 1 Model related information[13]

由于试样的厚度(10 mm)相对于其他尺寸比例较大，应采用平面应变条件进行建模。模型圆孔中心采用集中力形式进行约束，在ABAQUS中实现的具体操作：在圆孔中心设置参考点，在相互作用模块中将参考点和圆孔内边界进行耦合(耦合形式为Kinematic Coupling)，在加载模块中在参考点处施加集中力荷载以及边界条件(其中在上圆孔的参考点处约束U1 方向的位移/转角，在下圆孔的参考点处约束U1、U2 方向的位移/转角)。裂尖网格采用Isoparametric 8-node quadrilateral elements(CPE8R)单元，并将裂尖周围单元设置为奇异性单元，以满足裂纹尖端应力与应变的奇异性要求[13]，模型有限元网格如图4 所示。最后进行了收敛性研究，以确认在有限元建模中使用了适当数量的元素。

图4 模型有限元网格Fig.4 Model Finite Element Mesh

裂纹以设置的1 mm 的增量长度进行扩展，为了判断裂纹的扩展方向，采用上述修正准则计算每一增量步的偏转角度：① 对于GMTS 准则，将KI、KII和T应力的数值代入式(3)和式(4)即可求得下一增量步的裂纹扩展角；② 对于GSED 准则，将数值结果代入式(8)即可计算裂纹扩展角。

裂纹在不同扩展阶段下(选取扩展步1、6、11、15)的Mises 应力云图如图5 所示(以DCB2 为例)。图6～图7 分别给出了两种模型的GMTS 和GSED 准则预测的理论裂纹扩展路径与实验结果之间的对比，具体结果见表2。由对比结果可知，两种模型的对比结果相似，其中GMTS 准则未能预测实验中发生的裂纹偏转，而是沿初始直线路径扩展，与文献[13]中的计算结果吻合。基于GSED准则的实验结果与理论预测结果之间可以观察到良好的一致性，裂纹的整个扩展过程与试验数据基本符合，证明了程序的有效性。

表2 理论结果与实验结果对比Table 2 Comparison of theoretical and experimental results

图6 DCB1 裂纹扩展路径对比图Fig.6 DCB1 Comparison of crack propagation paths

图7 DCB2 裂纹扩展路径对比图Fig.7 DCB2 Comparison of crack propagation paths

进一步研究发现上述不同准则预测裂纹扩展路径存在差异的主要原因是SED 准则考虑了裂纹尖端附近6 个应力分量的作用，计算裂纹尖端附近局部的应变能密度，建立裂纹失稳判据。GSED准则相对于GMTS 准则，还考虑泊松比的影响。AYATOLLAHI 等[28]基于式(8)给出了断裂起始角θ0相对于Bα 的变化规律(ν=0.35)，如图8 所示。对于GMTS 准则，SMITH 等[29]在考虑Bα 对裂纹起始角θ0的影响时得出了类似的结论并绘制于图8，当Bα＞0.375 时(Bα 见式(8))，I 型裂纹弯曲；当0.22＜Bα＜0.375 时，GMTS 准则无法预测裂纹的偏转；而GSED 准则可以预测在此Bα 范围内裂纹的扩展路径。

图8 断裂起始角θ0 相对于Bα 的变化[[28-29]]Fig.8 Variation of fracture initiation angle θ0 relative to Bα[[28-29]]

4 结论

本文采用Python 语言对ABAQUS 进行二次开发，将考虑T应力的GMTS、GSED 准则引入裂纹起裂角的计算中并研究T应力对起裂角度的影响，将模拟结果与相关文献中含I 型裂纹的实验结果进行对比，得到以下结论：

(1)通过ABAQUS 二次开发接口，利用Python语言编写裂纹自动扩展程序，实现了模拟裂纹扩展过程中的自动化功能，并将考虑T应力的GSED和GMTS 准则编入程序脚本中，实现多种断裂准则对裂纹扩展路径的模拟。

(2)分别利用GSED 和GMTS 断裂准则研究了I 型载荷下裂纹扩展路径的差异，并将数值模拟结果与相关文献结果进行对比，其中GSED 准则能较好地预测该模式下的裂纹扩展路径，表明了程序的有效性。对于GMTS 准则，当Bα＞0.375 时才能预测裂纹的弯曲现象。