含裂纹管道结构的三维应力强度因子数值仿真

Numerical Simulation of Three-Dimensional Stress Intensity Factor of Cracked Pipe Structure

摘要: 含裂纹管道结构的应力强度因子对结构的安全评估至关重要，对于非穿透型裂纹来说，很难得到其理论解，必须采用数值仿真得到近似解。本文基于ANSYS软件给出了管道结构三维应力强度因子的数值建模方法，并利用含穿透型裂纹直管结构验证了本文方法的精度和可行性。进一步针对含非穿透型裂纹的弯管结构的三维应力强度因子进行了参数影响研究，分析了轴向及环向裂纹深度和长度对裂尖应力强度因子的影响规律。本文所建立的三维应力强度因子数值仿真方法可推广用于其它工程结构应力强度因子的计算评估。

Abstract: The stress intensity factor is very important to the safety evaluation of cracked pipe structures, and it is very hard to obtain the theoretical solution for part-through cracks, so numerical simulation is always employed to obtain the approximate solution. The present paper shows the numerical modeling approach for the three-dimensional stress intensity factor of pipe structure with ANSYS software, and the cracked straight pipe is employed to verify the accuracy and feasibility of the proposed approach. The influence of parameters on the three-dimensional stress intensity factor of elbow pipe structure with a part-through crack was investigated. The influence of depth and length of axial and radial part-through cracks on stress intensity factor was analyzed. The numerical simulation approach of the three-dimensional stress intensity factor proposed in this paper could be used for calculating and evaluating the stress intensity factor of other engineering structures.

1. 引言

管道结构广泛应用于石油天然气、化工、航空航天等各个工业领域，其安全性受到了学术界和工业界的普遍关注 [1]。在实际的应用过程中，管道结构中不可避免地会出现各种各样的缺陷，比如腐蚀性凹坑或者各种裂纹。针对含裂纹的管道结构，可以基于断裂力学进行计算分析 [2]。在断裂力学中，用应力强度因子K来反映裂纹尖端处的应力场的奇异性效果，其既反映了平均应力的效果，同时还包括了裂纹特征尺寸的影响，是表征裂纹尺寸是否扩展的一个参量。应力强度因子一般是工作应力、裂纹尺寸和裂纹几何的函数。

临界状态时的应力强度因子一般用 $K_{\text{C}}$ 表示，临界状态是指裂纹从缓慢扩展到失稳扩展的转折点(失稳扩展为载荷突然下降的裂纹的快速扩展)。 $K_{\text{C}}$ 也称为断裂韧性，它是材料常数。在实际工程中，I型裂纹最常见，也最危险，一般把其它类型的裂纹转化为I型裂纹来处理，统一采用I型裂纹结构的断裂判据：

$K_{\text{I}}=K_{\text{IC}}$ (1)

对于含裂纹的管道结构，可以计算出裂纹尖端的应力强度因子 $K_{\text{I}}$，并根据式(1)的裂纹扩展判据，为结构设计、选材和检验提供依据，可以对管道进行断裂安全评估。材料的断裂韧性 $K_{\text{IC}}$ 可查阅相关手册或进行测量得到。

对于复杂的工程结构，精确的应力强度因子往往很难获得，因此，需要采用有限元方法来求得数值解。王永伟等研究了有限元方法在求解裂纹尖端应力强度因子时的网格划分问题，并基于ANSYS软件建立了含表面裂纹的平板模型 [3]。伊军利用有限元分析软件ANSYS建立了内压作用下弯管环向穿透裂纹模型 [4]。喻健良等人基于ANSYS建立了含外表面轴向裂纹的直管模型，并对内压下裂纹尖端的应力强度因子进行了参数影响研究 [5]。鞠峰等人针对含轴向穿透裂纹的三通管道应力强度因子，进行了有限元分析 [6]。王冉等人基于ANSYS对内压作用下含纵向半穿透裂纹斜接弯管的应力强度因子，进行了三维数值仿真 [7]。于桂杰等人基于ABAQUS软件建立了直管外表面斜裂纹的有限元模型 [8]。俞树荣等人采用SOLID186单元建立了基于ANSYS的表面裂纹的三维模型及应力强度因子计算方法 [9]。贾旭等人提出了复杂载荷下偏心穿透裂纹应力强度因子的计算方法，给出了通用的权函数系数 [10]。Subbaiah利用二维轴对称有限元模型，针对圆柱形压力容器表面倾斜裂纹的应力强度因子进行了数值仿真 [11]。

ANSYS中有两种方法来计算应力强度因子：交互积分法(Interaction Integral Method)和位移外推法(Displacement Extrapolation Method)。交互积分法在求解时进行应力强度因子的计算，并将其存储到结果文件中用于后处理。由于围道积分在远离裂纹间断的点上进行计算，因此，交互积分法适用于混合型应力强度因子的高精度计算，同时也可用于非线性材料。位移外推法在后处理过程中进行应力强度因子的计算，只适用于裂尖区域是各向同性线弹性应力状态的问题。本文采用交互积分法进行应力强度因子的计算。

2. ANSYS计算三维应力强度因子的步骤

在采用ANSYS软件进行结构三维应力强度因子计算时，其建模、分网、施加载荷和边界条件与常规的有限元分析一样，值得注意的是，划分网格时需要利用KSCON命令对裂纹尖端的网格进行控制，使得网格呈发射状外扩。其使用方法为：KSCON,NPT,DELR,KCTIP,NTHET,RRAT。这里，NPT是网格密集处的关键点编号，也可以是一个关键点组的名称。DELR是以关键点为圆心的第一行单元的半径。KCTIP是裂纹尖端奇异性参数，0表示单元中节点位置不变，1表示第一行单元的中节点移动到1/4边长处，用于模拟裂尖的奇异性。NTHET是环向单元的份数(默认大概45度一个单元，如果KCTIP = 1，那么30度一个单元)。RRAT是第二行单元尺寸与第一行单元的比值(默认是0.75，如果KCTIP = 1，那么取为0.5)。

在进入求解模块/SOLU后，与应力强度因子相关的操作如下所示。

1) 创建一个新的应力强度因子计算任务并设定求解类型。

CINT,NEW,CRACK_NUM

CINT,TYPE,SIFS

这里，CRACK_NUM是裂纹编号，SIFS表示求解的是应力强度因子。

2) 定义裂纹尖端信息。

CINT,CTNC,CRACKTIP,Par2

应力强度因子的计算需要定义裂尖节点组CRACKTIP和裂纹扩展方向。裂纹扩展方向通过定义裂纹面的法向来实现。裂尖节点组需要包括裂纹尖端的所有节点，且这些节点是连接在一起的，ANSYS基于这些节点形成的裂纹尖端进行围道积分的单元。这里，Par2是位于裂纹张开侧的任意一个节点。

3) 设定围道个数。

CINT,NCON,COUNTER_NUM

这里，COUNTER_NUM是围道个数。

4) 设定对称性。

CINT,SYMM,ON

一般来说，对于混合型裂纹不存在对称性。如果设置了对称性，那么K_II和K_III均强制为0，计算得到的应力强度因子乘以2即为K_I。

5) 设定裂纹面的法向。

CINT,NORM,Par1,Par2

这里，Par1是定义的局部坐标系编号，Par2是该局部坐标系的坐标轴，默认值为2，对应于局部坐标系的y轴。

在后处理模块/POST1，通过对应力计算结果进行后处理即可得到裂纹的应力强度因子。主要的控制命令有：

*GET,KI_Cal,CINT,CRACK_NUM,CTIP,NODE_NUM,,COUNTER_NUM,,K1

这里，KI_Cal是存储在参数KI的变量名，CINT表明计算的是应力强度因子，CTIP是固定标识符。NODE_NUM是需要计算应力强度因子的节点编号。K1表明需要计算的是KI。注意，对于弯管内裂纹，需要选择裂纹尖端不同位置的节点进行应力强度因子的计算，反复使用本命令可以把不同节点处的应力强度因子赋给不同的变量名KI_Cal。

3. 数值算例

基于上述数值仿真方法，本部分首先针对含环向和轴向穿透裂纹的直管道结构，将应力强度因子的理论解与基于本文方法的数值仿真结果进行对比，验证本文方法的精度与可靠性。进而针对含非穿透型裂纹的典型弯管结构，计算不同几何尺寸轴向及环向裂纹的裂纹应力强度因子。

3.1. 含穿透裂纹直管道的应力强度因子

对于直管道来说，可以近似认为壁厚上的应力沿厚度方向是均匀分布的，与平板的受力状态相同。因此，对于直管道上的裂纹，其应力强度因子近似按平板相应裂纹的应力强度因子公式计算，在此基础上，考虑到曲率与膨胀效应的修正，可得到与实际情况基本相符的结果。按照裂纹长度方向不同，管道中的裂纹分为轴向裂纹和环向裂纹。轴向裂纹的裂纹长度方向与管道轴向同向，环向裂纹的长度方向沿着管道的周向。

对于直管道上的穿透型裂纹，其应力强度因子为 [12]：

$K_{\text{I}}=M\sigma \sqrt{\pi a}$ (2)

这里，M为鼓胀因子，a为裂纹半长。均匀拉伸应力 $\sigma$ 在轴向裂纹时取环向应力 ${\sigma }_1$，在环向裂纹时取轴向应力 ${\sigma }_{\text{2}}$，且有：

${\sigma }_1=\frac{pD}{2t},{\sigma }_2=\frac{pD}{4t}$ (3)

其中，p为管道内压，D为管道外径，t为管道壁厚。值得注意的是，式(3)仅适用于薄壁管道。

轴向和环向穿透裂纹对应的鼓胀因子 $M_1$ 和 $M_2$ 分别为：

$M_1=\sqrt{1+3.21\frac{a^2}{Dt}},M_2=\sqrt{1+0.64\frac{a^2}{Dt}}$ (4)

直管道外径D = 200 mm，厚度t = 10 mm，裂纹长度2a = 20 mm，内压p = 2 MPa。管道材料的弹性模量为210 GPa，泊松比为0.3。采用带中间节点的二次六面体实体单元SOLID186。基于式(2)得到该管道轴向和环向穿透裂纹的应力强度因子分别为 $\text{3}\text{.82MPa}\cdot \sqrt{\text{m}}$ 和 $\text{1}\text{.80MPa}\cdot \sqrt{\text{m}}$。在采用前节所述数值仿真方法分别对含轴向和环向穿透型裂纹结构进行仿真，得到的裂纹尖端的等效应力场分别如图1和图2所示，进一步计算得到的应力强度因子分别为 $\text{3}\text{.56MPa}\cdot \sqrt{\text{m}}$ 和 $\text{1}\text{.83MPa}\cdot \sqrt{\text{m}}$。

可以看出，基于本文数值仿真方法得到的结合和经验解吻合得很好，出现的误差在于经验解基于管道厚度上应力是均匀分布的假设，而实际上对于有一定厚度的管道来说，其厚度方向上的应力分布并不完全是均匀的，由此也导致了沿着裂纹尖端不同位置得到的应力强度因子是不同的。

Figure 1. Equivalent stress field of the crack tip of the axial penetrated crack

图1. 轴向穿透裂纹尖端的等效应力场云图

Figure 2. Equivalent stress field of the crack tip of the radial penetrated crack

图2. 环向穿透裂纹尖端的等效应力场云图

3.2. 含非穿透裂纹弯管结构的应力强度因子

本节研究的弯管管道外径为406 mm，壁厚10 mm，弯曲部分的曲率半径为2 m，对应的圆心角为20度，其几何模型如图3所示。管道内压8.5 MPa，其下部端面固定，上部端面施加30 MPa的轴向拉力来模拟管道其它部分对模型的作用。

Figure 3. Geometry of the pipe

图3. 弯管几何模型

3.2.1. 轴向裂纹

对于轴向的非穿透裂纹，这里假设裂纹位于管道内壁的内侧。裂纹为椭圆形，如图4所示，这里B为管道壁厚。沿管道轴向的裂纹长度为2a = 20 mm，其深度b = 5 mm。

Figure 4. Geometry of the part-through crack

图4. 非穿透型裂纹几何

基于上述几何模型、载荷条件和材料属性得到的管道结构裂纹尖端的等效应力场如图5所示。经过后处理得到裂纹顶端的应力强度因子为 $\text{17}\text{.62MPa}\cdot \sqrt{\text{m}}$。

Figure 5. Equivalent stress field of the crack tip of the axial part-through crack

图5. 非穿透轴向裂纹尖端的等效应力场云图

保持轴向裂纹深度5 mm不变，不同轴向裂纹长度时裂纹尖端应力强度因子的变化趋势如图6所示。保持裂纹长度20 mm不变，不同轴向裂纹深度时裂纹尖端应力强度因子的变化趋势如图7所示。

Figure 6. Variations of stress intensity factor with the length of the axial part-through crack tip

图6. 非穿透轴向裂纹尖端应力强度因子随裂纹长度的变化趋势

Figure 7. Variations of stress intensity factor with the depth of the axial part-through crack tip

图7. 非穿透轴向裂纹尖端应力强度因子随裂纹深度的变化趋势

从图6和图7可以看出，随着轴向裂纹长度的增加，裂尖应力强度因子基本上呈线性增加趋势，而随着裂纹深度的增加，在裂纹深度不大的情况下，裂尖应力强度因子随其深度的增加，基本上呈线性增加趋势，但当裂纹深度比较大时，裂尖应力强度因子增加趋势变缓。

3.2.2. 环向裂纹

对于环向的非穿透裂纹，这里也假设裂纹位于管道内壁的内侧。裂纹几何形状与轴向裂纹相同均为椭圆形。沿管道环向的裂纹长度为2a = 20 mm，其深度b = 6 mm。基于上述几何模型、载荷条件和材料属性得到的管道结构裂纹尖端的等效应力场如图8所示。经过后处理得到裂纹顶端的应力强度因子为 $\text{6}\text{.37MPa}\cdot \sqrt{\text{m}}$。

Figure 8. Equivalent stress field of the crack tip of the radial part-through crack

图8. 非穿透环向裂纹尖端的等效应力场云图

保持环向裂纹深度6 mm不变，不同环向裂纹长度时裂纹尖端应力强度因子的变化趋势如图9所示。保持裂纹长度20 mm不变，不同环向裂纹深度时裂纹尖端应力强度因子的变化趋势如图10所示。

Figure 9. Variations of stress intensity factor with the length of the radial part-through crack tip

图9. 非穿透环向裂纹尖端应力强度因子随裂纹长度的变化趋势

Figure 10. Variations of stress intensity factor with the depth of the radial part-through crack tip

图10. 非穿透环向裂纹尖端应力强度因子随裂纹深度的变化趋势

从图9和图10可以看出，随着环轴向裂纹长度和深度的增加，裂尖应力强度因子基本上呈线性增加趋势。当裂纹深度或长度比较大时，裂尖应力强度因子的增加有变缓的趋势。

4. 结论

本文给出了利用ANSYS进行三维应力强度因子数值仿真的关键步骤及实施过程，同时针对含穿透型裂纹的直管道的裂尖应力强度因子进行了计算，并与文献中的工程近似解进行了对比，二者结果吻合得很好，从而验证了本文数值仿真方法的精度和可靠性。本文进一步针对含非穿透型裂纹弯管结构的应力强度因子进行了数值仿真，并分别针对轴向和环向裂纹的长度及深度对其应力强度因子的影响进行了参数分析。数值结果表明，随着非穿透型裂纹长度和深度的增加，裂纹尖端应力强度因子会逐渐增加，但是当裂纹深度或长度比较大时，此增加的趋势会变缓。本文所建立的三维裂纹应力强度因子数值仿真方法可进一步应用于其它类型工程结构的安全性评估。