随着社会生产的发展，新工艺、新技术、高强度材料得到广泛应用，工程部件和结构在服役期间或由于制造缺陷会产生微裂纹，之后由于应力集中、疲劳、腐蚀等因素的影响，微裂纹不断成核、扩展、生长为宏观裂纹，最终导致材料断裂失效。断裂准则是含裂纹材料在各种应力作用下断裂和破坏的重要判据，解释了裂纹何时起裂和判断裂纹扩展方向的问题，因此在预测裂纹扩展路径方面得到广泛应用。传统线弹性断裂力学中常用的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]可以计算裂纹尖端的应力强度因子K_I、K_II，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型裂纹扩展实验结果建立了最大周向应力准则，即在脆性材料中，裂纹将沿着裂尖周围的最大周向应力方向θ₀开裂。CHAO等 ^[8-9]提出了经过T应力修正的最大周向应力准则，修正后的准则考虑了K和常数项T应力的影响。GMTS(Generalized maximum tangential stress)准则可以表示如下：

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

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

由式(3)可知，如果${{{T_{\text{c}}}} \mathord{\left/ {\vphantom {{{T_{\text{c}}}} {{K_{\text{c}}}}}} \right. } {{K_{\text{c}}}}} < \dfrac{3}{8}\dfrac{1}{{\sqrt {2{\text{π }}r} }}$，则I型断裂条件下的裂纹将沿${\theta _0} = 0$方向直线扩展；如果${{{T_{\text{c}}}} \mathord{\left/ {\vphantom {{{T_{\text{c}}}} {{K_{\text{c}}}}}} \right. } {{K_{\text{c}}}}} > \dfrac{3}{8}\dfrac{1}{{\sqrt {2{\text{π }}r} }}$，则裂纹将发生偏转，裂纹起裂角度由式(3)的第二个解确定。综上，I型断裂条件下的裂纹弯曲将发生在T应力超过某个正值时。

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

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

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

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

式中：${{{\text{d}}W} \mathord{\left/ {\vphantom {{{\text{d}}W} {{\text{d}}V}}} \right. } {{\text{d}}V}}$为单位面积的应变能；r_c为临界半径；ν、E分别为材料的泊松比和弹性模量；σ_x、σ_y、σ_z为空间直角坐标系中的正应力分量；τ_xy、τ_xz、τ_yz为空间直角坐标系中的剪切应力分量。将包含常数项T应力的裂纹尖端应力场表达式代入式(5)~式(7)中，得到GSED(Generalized Minimum Strain Energy Density)准则，经求导化简后其表达式为：

其中：

式中：$\alpha = \sqrt {\dfrac{{2{r_{\rm c}}}}{a}}$，$B = \dfrac{{T\sqrt {{\text{π }}a} }}{{{K_{{\text{eff}}}}}}$，a为裂纹长度。平面应变问题：$\kappa = 3 - 4\nu$；平面应力问题：$\kappa = ( {3 - \nu } ) / ({1 + \nu } )$。在K_Ⅰ、K_Ⅱ和T应力已知的情况下，求解上述式(8)即可得到I型或Ⅰ-Ⅱ复合型裂纹扩展角θ_m。

2   裂纹扩展模拟方法

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

2.1   裂纹自动扩展程序流程

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

图  1  程序流程图

Figure  1.  Program flow chart

下载: 全尺寸图片 幻灯片

裂纹自动扩展程序功能主要包含以下几个部分：应力分析及裂纹扩展检测、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个粒子，用v_i为粒子的速度，x_i为粒子的位置，速度和位置更新公式如下：

速度更新公式为：

位置更新公式为：

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

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

图  2  粒子群算法流程图

Figure  2.  Flow chart of particle swarm algorithm

下载: 全尺寸图片 幻灯片

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

2.3   程序编制

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

#MTS+T 计算裂纹起裂角度

K_c=T/K₁

rad=acos((1+sqrt(1+(1024*pi*r_c*(K_c**2))/9))/((512*pi*r_c*(K_c**2))/9))

Cpd=rad*180/pi

#SED+T 计算裂纹起裂角度(粒子数N=20，迭代次数iter_N=100)

α=sqrt(2*0.86/a)

B=T*sqrt(pi*arf)/K₁

pso=PSO(N=20, iter_N=100,v=v,B=B,$\text{α}$=$\text{α}$)

def fitness(self, x):

return (4−8*self.v)*np.cos(x)−np.cos(2*x)+self.B *self.α*((6−16*self.v)*np.cos(x/2)+2* np.cos (5* x/2))

pso.implement()

Cpd=pso.gbest*180/pi

#SED+T 计算裂纹扩展角度(粒子数N=20，迭代次数iter_N=100)

pso=PSO(N=20, iter_N=100, v=v, K₁=K₁, K₂= K₂)

def fitness(self, x):

return ((self.$\text{κ}$−1)*np.cos(x)−0.5*np.cos(2*x))* (self.K₁**2)+2*((1−self.$\text{κ}$)*np.sin(x)+ np.sin (2*x))*self. K₁*self.K₂+((1−self.$\text{κ}$)* np.cos(x) +1.5*np.cos(2*x))*(self.K₂**2)

pso.implement()

Cpd=pso.gbest*180/pi

3   应用验证

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

图  3  试样几何形状

Figure  3.  Specimen geometry

下载: 全尺寸图片 幻灯片

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

表  1  模型相关信息^[13]

Table  1.  Model related information^[13]

试样类型	模型长度W/mm	模型宽度h/mm	裂纹长度a/mm	断裂载荷Fc/N
DCB1	90	30	45	141
DCB2	150	30	75	87

下载: 导出CSV 

| 显示表格

由于试样的厚度(10 mm)相对于其他尺寸比例较大，应采用平面应变条件进行建模。模型圆孔中心采用集中力形式进行约束，在ABAQUS中实现的具体操作：在圆孔中心设置参考点，在相互作用模块中将参考点和圆孔内边界进行耦合(耦合形式为Kinematic Coupling)，在加载模块中在参考点处施加集中力荷载以及边界条件(其中在上圆孔的参考点处约束U1方向的位移/转角，在下圆孔的参考点处约束U1、U2方向的位移/转角)。

裂尖网格采用Isoparametric 8-node quadrilateral elements(CPE8R)单元，并将裂尖周围单元设置为奇异性单元，以满足裂纹尖端应力与应变的奇异性要求^[13]，模型有限元网格如图4所示。最后进行了收敛性研究，以确认在有限元建模中使用了适当数量的元素。

图  4  模型有限元网格

Figure  4.  Model Finite Element Mesh

下载: 全尺寸图片 幻灯片

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

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

图  5  Mises应力云图(DCB2)

Figure  5.  Mises stress nephogram(DCB2)

下载: 全尺寸图片 幻灯片

图  6  DCB1裂纹扩展路径对比图

Figure  6.  DCB1 Comparison of crack propagation paths

下载: 全尺寸图片 幻灯片

图  7  DCB2裂纹扩展路径对比图

Figure  7.  DCB2 Comparison of crack propagation paths

下载: 全尺寸图片 幻灯片

表  2  理论结果与实验结果对比

Table  2.  Comparison of theoretical and experimental results

试样类型	Bα	断裂起始角θ0 (文献[13])	起裂角θ0 (GSED)	起裂角θ0 (GMTS)
DCB1	0.22	−9	−7	0
DCB2	0.25	−56	−49	0

下载: 导出CSV 

| 显示表格

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

图  8  断裂起始角θ₀相对于Bα的变化^[[28-29]]

Figure  8.  Variation of fracture initiation angle θ₀ 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时才能预测裂纹的弯曲现象。