基于Abaqus计算界面裂纹应力强度因子的可行性分析

复合材料是由两种或多种不同性质的材料用物理和化学方法在宏观尺度上组成具有新性能的材料，在复合材料中往往存在多种形式的材料界面。因此，在评估复合材料的断裂性能时，需要计算材料界面间裂纹尖端的断裂控制参量，如线弹性断裂力学中的应力强度因子等。

目前，有限元软件Abaqus已经被广泛应用于计算复杂结构中裂纹尖端的应力强度因子。Abaqus采用围线积分法（contour integral）计算应力强度因子，该方法是基于J积分法建立的，而传统的J积分法仅适用于无体积力的单一均质材料，且裂纹面不存在面力的情况。虽然目前已经建立了多种能够准确计算复杂材料界面情况下裂纹尖端断裂参量的交互积分法，但Abaqus内置的围线积分算法能否准确提取界面裂纹尖端的应力强度因子尚缺乏验证。为此，本文以含界面裂纹的双材料均质板和非均质板为例，采用有限元软件Abaqus内置的围线积分法，分别计算界面裂纹尖端的应力强度因子，并与现有文献给出的参考解进行对比，以考察基于Abaqus计算界面裂纹应力强度因子的可行性。

考虑如图1所示的中心裂纹板，平板的长为2L，宽为2W，中心裂纹长度为2a，平板的上、下端部受到均匀拉应力σ₀的作用。平板由两种材料构成，其弹性模量分别为E₁和E₂，泊松比分别为ν₁和ν₂。计算示例中用到的数据为：2W=100 mm，2L=200 mm，归一化裂纹长度a/W=0.4，E₁=2.058×10⁵ MPa，E₁/E₂=1~100，ν₁=ν₂=0.3，σ₀=9.8 MPa。假设平板处于平面应力状态，下面采用Abaqus内置的围线积分法计算界面裂纹尖端的应力强度因子，所使用的Abaqus版本为Abaqus 2020。

 

图1 含界面裂纹的双材料均质板示意图

该计算示例在建模方面并不存在较大难度，仅需在普通中心裂纹板有限元模型的基础上分别在不同区域赋予相应的材料属性即可，因此本文不再赘述相应的建模过程，关于裂纹建模和应力强度因子计算的相关内容可参考《二维中心裂纹板的应力强度因子分析》。

通过计算，图2给出了当上、下材料的弹性模量比E₁/E₂为100时，含界面裂纹均质板的米塞斯应力分布。

 

图2 含界面裂纹的均质板的米塞斯应力云图（E₁/E₂=100）

从图中可以看到，由于平板上、下材料的刚度差异明显，因此上、下裂纹面的裂纹张开位移并非是对称分布的，并且位于裂纹面上、下两侧的应力分布也出现了明显的不连续现象，这与均质材料中的中心裂纹有很大区别。

基于围线积分法，图3给出了裂纹尖端应力强度因子随围线的变化规律。这里，不同围线代表了不同的积分区域，围线编号越大则围线距离裂纹尖端越远，积分区域越大，从而能够避免过于靠近裂纹尖端而引起的数值误差。

 

(a)K₁

 

(b)K₂

图3 均质板间界面裂纹应力强度因子随围线变化

从图中可以看到，当围线编号大于3后，K₁和K₂均趋于稳定值，这证明在本算例中围线积分的守恒性得到了很好地满足。注意，由于界面裂纹中剪切和拉伸断裂模式无法分离，因此这里采用K₁和K₂来区分传统的I型和II型应力强度因子。这里，我们可以对均质材料间界面裂纹的守恒性做一个简单证明。考虑如图4所示具有界面裂纹的双材料。

 

图4 双材料界面裂纹尖端的J积分路径

已知J积分的定义为：

式中：W为应变能密度，u_i为位移矢量的分量，ds为积分路径上的微段。

分别研究裂纹尖端附近的三条路径：

由于路径Γ₁和Γ₂没有围绕裂纹尖端，并且处于单一的均质材料中，因此由J积分的守恒性有：

注意到

因此有

在材料界面CD上，满足dy=0，并且由于Y方向的应力和位移都是连续的，因此有：

代入可得

可以看到，对于双材料均质板，界面裂纹J积分的守恒性可以得到满足，这也是为什么通过Abaqus计算得到的应力强度因子值可以趋于稳定值。当然，这种守恒性仅当裂纹平行于材料界面时才能够成立，当裂纹不平行时，由于上式中的J积分不再等于零，因此守恒性不成立。

需要注意的是，围线积分的守恒性成立并不代表Abaqus给出的应力强度因子是准确的，因此仍需要将Abaqus计算得到的应力强度因子与精确值进行对比验证。表1给出了在不同弹性模量比E₁/E₂下，含界面裂纹的双材料均质板的归一化应力强度因子与参考解的对比结果。表中Abaqus给出的应力强度因子是通过剔除前三个围线上的应力强度因子并将余下的应力强度因子取平均值得到的，而归一化应力强度因子的参考解来自于文献1，其中K₀的定义为：

表1 含界面裂纹的双材料均质板的归一化应力强度因子（平面应力状态，a/W=0.4）

从表中可以看到，对于均质材料（E₁/E₂=1），Abaqus给出的结果与参考解吻合良好；而当两种材料存在刚度差异时，对于拉伸模式引起的应力强度因子K₁，Abaqus给出的计算结果与参考解吻合较好；但对于剪切模式引起的应力强度因子K₂，两者差别较大，并且误差随着弹性模量比E₁/E₂的增大而增大。这说明Abaqus内置的围线积分法并不能准确计算界面裂纹的应力强度因子。当两种材料的刚度差异较小时，可以采用Abaqus做近似计算。

考虑如图5所示的含界面裂纹的双材料非均质板，平板长为2L，宽为2W，中心裂纹长度为2a，平板的下端固定，上端受到均匀拉应力σ₀的作用。平板由两种材料构成，其弹性模量E₁和E₂分别可表示为：

式中：E₀为参考弹性模量，δ的定义为

两种材料的泊松比分别为ν₁和ν₂。

图5 含界面裂纹的双材料非均质板示意图

计算示例中用到的数据为：L=200 mm，W=50 mm，归一化裂纹长度a/W=0.4，E₀=1000 MPa，E₁(W)/E₂(-W)=10~100，ν₁=ν₂=0.3，σ₀=1 MPa。假设平板处于平面应变状态，下面需要计算界面裂纹尖端的应力强度因子。

由于两种材料的弹性模量在平板宽度方向并不是均匀变化的，因此在建立有限元模型时必须建立全模型。建模中比较难以处理的部分为材料的弹性模量与材料点的坐标存在关联，必须使用Abaqus用户子程序USDFLD（子程序代码参见附录）以自定义场的形式来给定弹性模量。因此，在定义材料属性时，需要在General中选择User Defined Field，并在Elastic中定义场变量（field variable），通过将场变量Field 1与材料的弹性模量和泊松比关联在一起，即可通过用户子程序USDFLD定义场变量值来实现预定的材料属性值，如图6所示。同时，本文还通过在Depvar中定义一个状态变量STATEV来实现材料弹性模量的绘制，从而检验用户子程序定义的弹性模量是否正确，如图7所示。

图6 用户自定义场设置

图7 自定义状态变量

为了保证用户自定义场变量和状态变量能够在后处理中以云图的形式呈现，还需要在场输出中勾选SDV和FV场变量，如图8所示。

图8 场变量输出设置

同时，还需要在inp文件中添加如下所示的关键字：

这里，VARIABLE表示场变量的数量。

基于上述设置进行计算，计算中取E₁(W)/E₂(-W)=10。计算完成后，得到采用状态变量SDV1绘制的弹性模量分布如图9所示。从图中可以看到，材料弹性模量的分布与图5中给出的弹性模量分布基本一致，证明子程序定义的弹性模量场无误。

图9 非均质板中弹性模量分布

基于上图所示的弹性模量场，图10给出了含界面裂纹的非均质板的米塞斯应力分布。可以看到，平板呈现S形的弯曲变形，这与均质板的变形模式有很大区别。

 

图10 含界面裂纹的非均质板的米塞斯应力分布

表2给出了采用Abaqus内置的围线积分法计算得到的归一化应力强度因子与参考解的对比结果，其中应力强度因子的参考解来源于文献2。结果显示，Abaqus给出的应力强度因子与参考解的误差很大，这说明对于非均质板，Abaqus内置的围线积分法同样不能准确计算界面裂纹尖端的应力强度因子。

表2 含界面裂纹的双材料非均质板的归一化应力强度因子（平面应变状态，a/W=0.4）

图11给出了当E(W)/E(-W)=10时，由Abaqus给出的应力强度因子随围线编号的变化。可以看到对于拉伸断裂模式引起的应力强度因子K₁，随着围线编号的增大应力强度因子呈现收敛的趋势；而对于剪切断裂模式引起的应力强度因子K₂，应力强度因子并没有呈现收敛的趋势，说明对于非均质材料的界面裂纹，围线积分的守恒性也无法得到满足。

 

(a)K₁

 

(b)K₂

图11 非均质板间界面裂纹应力强度因子随围线变化

本文以含界面裂纹的双材料均质板和非均质板为例，采用Abaqus 2020内置的围线积分法计算了裂纹尖端的应力强度因子，并与现有文献给出的参考解进行了对比。结果表明，对于含界面裂纹的双材料均质板，围线积分的守恒性能够得到很好地满足，但两种材料的刚度差异较大时，计算结果有较大误差；对于含界面裂纹的双材料非均质板，围线积分的守恒性无法得到满足，并且计算结果误差很大，说明当前版本的Abaqus内置的围线积分法尚不具备准确计算界面裂纹应力强度因子的功能。

[1] T. Nagashima, Y. Omoto and S. Tani. Stress intensity factor analysis of interface cracks using X-FEM. International Journal for Numerical Methods in Engineering. 2003, 56: 1151-1173.

[2] 于红军. 含复杂界面非均匀材料断裂力学研究[D]. 哈尔滨工业大学, 2010.