0 引　言

由于复杂载荷、环境以及内部缺陷（裂纹）等因素协同作用下，船舶与海洋工程等重要机械结构在服役过程中往往会发生突然断裂事故，造成不可挽回的“灾难性”后果^[1]。针对断裂问题的研究，参数应力强度因子（stress intensity factor, SIF）十分重要，常用于确定裂纹扩展状态和稳定性^[2]。因此，SIF的准确获取对含缺陷材料和结构的裂纹扩展状态以及寿命预测至关重要。

SIF的获取方法主要有实验、公式解以及数值模拟法。王强胜等^[3]通过位错法求解了SIF，该方法成本高且实验周期长，对SIF的求解具有一定的局限性。解析法往往针对简单问题进行求解，Chen等^[4]，Eshraghi等^[5]分别采用权函数法、复变函数法对SIF的求解进行了研究。即使针对均一材料，公式中的积分和矩阵运算仍有一定的难度^[6]。因此通过数值模拟法准确快速的获取SIF已被广泛采用，在数值模拟法中，主要的计算方法有围道积分法和扩展有限元法。为提高计算精度，Wang等^[7]在扩展有限元法中引入高阶富集函数来精确计算SIFS。针对复杂问题，扩展有限元法能够保证足够的可行性与计算效率。同时，它还可以用来获得弯曲裂纹的SIF，对于这一点其他方法很难应用或不适用。

综上所述，很多学者在裂纹SIF的获取方法上做了大量研究。然而，在使用扩展有限元法获取SIF的计算上仍然缺乏相关的文献对其阐述。因此，本文基于扩展有限元法研究三维裂纹前缘SIF的求解，并与围道积分法和公式解析解的计算结果进行对比，分别从收敛性、精度以及网格密度等对结果的影响角度分析，研究基于扩展有限元法所获取的三维裂纹前缘SIF的规律。通过该方法的研究，使得复杂任意裂纹前缘SIF的获取变得简单、快捷，为含裂纹缺陷的结构完整性分析提供了一种新的技术手段。

1 扩展有限元简介 1.1 扩展有限元法的基本原理

在传统有限元的基础上，以单位分解理论为基础的扩展有限元法降低了对网格的依赖度，并引入富集函数，通过富集裂纹附近的少量单元以反映非连续位移场^[8]。因此，应用扩展有限元法分析断裂问题时，可直接划分网格，用能反映裂纹尖端位移场奇异性的附加函数以及反映裂纹面不连续性的阶跃函数来描述单元内部的局部属性^[9]。位移场近似的扩展有限元法公式可定义为：

$ u=\sum _{i=1}^{N}{N}_{I}\left(x\right)\left[{u}_{I}+H\left(x\right){a}_{I}+\sum _{a=1}^{4}{F}_{a}\left(x\right){b}_{I}^{a}\right]。$

(1)

其中：N_I(x)为节点位移型函数； $ {u}_{I} $ 为节点位移； $ H\left(x\right){a}_{I} $ 用来描述裂纹上下表面特性， $ H\left(x\right) $ 为Heaviside函数，是一种间断跳跃函数（见式（2））； $ {a}_{I} $ 为裂纹两边增强节点自由度向量； $ {F}_{a}\left(x\right){b}_{I}^{a} $ 用来描述裂尖的奇异性； $ {F}_{a}\left(x\right) $ 为裂尖的应力渐进函数； $ {b}_{I}^{a} $ 为裂尖增强节点自由度向量。

Heaviside增强函数在裂纹面的一侧取值1，另一侧取值为−1，函数表达式为：

$ H\left(x\right)=\left\{\begin{array}{c}1,裂纹面之上，\\ -1,裂纹面之下。\end{array}\right. $

(2)

裂尖渐进函数 $ {F}_{a}\left(x\right) $ 的极坐标表达式为：

$ {F}_{a}\left(x\right)=\left[\sqrt{r}\mathit{\rm sin}\frac{\theta }{2},\sqrt{r}\mathit{\rm cos}\frac{\theta }{2},\sqrt{r}\mathit{\rm sin}\frac{\theta }{2}\mathit{\rm sin}\theta ,\sqrt{r}\mathit{\rm sin}\theta \mathit{\rm cos}\frac{\theta }{2}\right]。$

(3)

1.2 相互作用积分法

工程结构内部不可避免的含有许多细微裂纹，在服役过程中裂纹有一个扩展的过程。而裂纹前缘的SIF与裂纹的启裂及扩展息息相关，通过对SIF的分析来表征裂纹尖端的位移、应力和应变场^[11]。扩展有限元法中与线弹性断裂力学相关的常用方法是通过J积分及其变体计算SIF。相互作用积分法作为J积分的扩展，是在真实场J积分的基础上引入一个辅助场，将真实场和辅助场叠加后代入J积分表达式，获取两者因相互作用而产生的J积分分量，将这个分量称为相互作用积分^[12]。定义如下^[13]：

$ {M}^{\left(\mathrm{1,2}\right)}={\int }_{\Gamma }^{}\left({W}^{\left(\mathrm{1,2}\right)}{\delta }_{1j}-{\sigma }_{ij}^{\left(1\right)}\frac{{\partial u}_{i}^{\left(2\right)}}{{\partial x}_{1}}-{\sigma }_{ij}^{\left(2\right)}\frac{{\partial u}_{i}^{\left(1\right)}}{{\partial x}_{1}}\right){n}_{j}{\rm{d}}\Gamma。$

(4)

式中：(1)代表真实场；(2)代表辅助场； $ {W}^{\left(\mathrm{1,2}\right)} $ 为应变能密度；Γ为围绕裂纹尖端的任一逆时针回路； $ {n}_{j} $ 为Γ的单位外法线向量。

2 三维裂纹有限元模型

选取恒载荷1T-CT试样作为计算的标准试样，1T-CT试样几何尺寸如图1所示。扩展有限元法模型对网格依赖度低，无需对裂纹前缘网格做特殊处理，全局网格采用C3D8单元，如图1(b)所示。分析网格尺寸参数对裂纹前缘SIF的影响程度，最终确定采用扩展有限元法计算裂纹前缘SIF的优越性。通过有限元软件Abaqus模拟加载过程和计算过程，试样材料选择316L不锈钢，弹性模量 $ E $ 为206000 MPa，泊松比 $ v $ 为0.28。模型左端固定，对2个施加载荷的参考点施加Y方向的载荷约束，同时约束其在X，Z方向的自由度，限制其绕X轴或Y轴的转动，试样两端承受垂直裂纹面方向的均匀拉伸载荷作用，使模型符合实际工况约束条件。

图中，净宽度W=50，厚度B=0.5 W，a=0.45−0.55 W，孔直径D=0.25 W。可以得到CT试样应力强度因子的理论解如下式：

$ \begin{split}{K}_{I}\frac{B{W}^{\frac{1}{2}}}{P}=&29.6\left(\frac{a}{W}\right)^{\frac{1}{2}}-185.5\left(\frac{a}{W}\right)^{\frac{3}{2}}+\\ &655.7\left(\frac{a}{W}\right)^{\frac{5}{2}}-1\;017\left(\frac{a}{W}\right)^{\frac{7}{2}}+638.9\left(\frac{a}{W}\right)^{\frac{9}{2}}。\\[-15pt]\end{split} $

(5)

其中，P为所施加在CT试样上的恒载荷。

此公式在 $0.2\leqslant \dfrac{a}{W}\leqslant 1.0$ 范围内使用，其准确度为 $\pm 0.5\text{%}，\dfrac{H}{W}=0.6$ 。

3 结果对比与分析 3.1 扩展有限元法计算收敛性的分析

如图2所示，通过对比不同网格尺寸下2种算法可知，当采用围道积分法进行SIF计算，网格尺寸为0.8和1.5 mm，围线数达到3圈时，计算结果逐渐收敛且趋于稳定；当网格尺寸为2.5和3.0 mm时，计算值收敛性明显降低且随着网格变粗，误差亦随之增大；当全局网格尺寸为4.0和5.0 mm时，计算结果发生明显的不收敛。然而，当采用扩展有限元法计算SIF时，随着网格变粗，单元数减少，计算效率提高，且计算结果仍然收敛。因此，在同等网格尺寸作用下，扩展有限元法具有良好的收敛性。

3.2 扩展有限元法计算精度的分析

分别对比裂纹长度a为8，9，10 mm时，不同网格尺寸(size)下其精度的计算值，如图3所示。由图可知，2种计算方法具有良好的一致性，相对于围道积分法，扩展有限元法的收敛性以及精度更能达到理想的程度。沿厚度方向，SIF大小以两边小中间大的规律对称分布，从而证明裂纹总是从中间开裂，与实际情况相符。当a=9 mm，size=0.8 mm时，围道积分法和扩展有限元法的计算精度误差最小，为0.18%；当a=8 mm，size=1.5 mm时，围道积分法和扩展有限元法的计算精度误差最大，为3.06%。因此，扩展有限元法在计算SIF时可以保证足够的精度。

3.3 扩展有限元法计算精度与网格尺寸的关系

如图4所示，使用扩展有限元法计算SIF时，计算结果具有良好的收敛性，且对网格的依赖度较低。通过和理论值对比发现可知，使用扩展有限元法求解应力强度因子时，并非网格越细求解越精确。由表1所示，通过对比分析可知，当a=8 mm，size为3 mm时，计算值与理论值的误差为4.2%；当a=9 mm，size为4 mm时，计算值与理论值的误差为1.4%；当a=10 mm，size为3 mm时，计算值与理论值的误差为0.27%；当a=11 mm，size为4 mm时，计算值与理论值的误差为2.3%，均在允许的误差范围内。因此该方法可以准确且高效的获取裂纹前缘SIF。

4 结　语

本文以Abaqus软件为平台，通过采用紧凑拉伸试样，研究了基于扩展有限元法获取裂纹前缘SIF的数值解，并分别与围道积分法和公式解析解进行了对比验证分析，主要结论如下：

1）通过对不同裂纹长度以及不同网格尺寸下SIF的计算值可知，扩展有限元法与围道积分法具有良好的一致性，相对于围道积分法，扩展有限元法的收敛性以及精度更能达到理想的程度。当a=9 mm，size=0.8时，围道积分法和扩展有限元法的计算精度误差最小，为0.18%；当a=8 mm，size=1.5时，围道积分法和扩展有限元法的计算精度误差最大，为3.06%。因此，扩展有限元法在计算SIF时可以保证足够的精度。

2）沿厚度方向SIF先递增至最大值而后逐渐递减，且以两边小中间大的规律对称分布，即裂纹总是从中间开裂，与实际情况相符。

3）使用扩展有限元法计算SIF时，计算结果具有良好的收敛性，且对网格的依赖程度较低，随着网格尺寸的增加，扩展有限元法可以保证足够的收敛性与稳定性，提高计算效率。