基于多相耦合的含缺陷埋地管道外爆破坏影响因素研究

丁宇奇，戴子威，芦 烨，杨 明，张 强，叶碧涛

(东北石油大学 机械科学与工程学院，黑龙江大庆 163318)

0 引言

埋地管道在制造生产、焊接和安装过程中，不可避免地会在表面产生各种缺陷。当遭受临近管道爆炸、工程爆破等偶然外爆载荷作用，缺陷处可能会被破坏，在管内流体作用下缺陷可能开始扩展并造成流体泄漏和管道大面积破坏，因此研究爆炸等冲击性载荷作用下含缺陷埋地管道的破坏问题对预防和控制此类事故具有重要意义[1-7]。

含缺陷埋地管道在外爆载荷作用下的破坏问题涉及到炸药-土体-管道-管内流体的多相耦合作用和爆炸载荷对管道的破坏问题。对于外爆多相耦合方面的研究，宋娟等[8]利用共节点、接触和ALE等算法探讨了土中爆炸数值计算方法的优劣，得出ALE算法更适合模拟土中爆炸问题的结论；梁博等[9]用SPH-FEM耦合算法研究埋地管道近场爆炸的冲击响应，实现了土体与管道之间不同算法的耦合，详细描述了最大冲击应力的截面分布及其在不同时刻的出现位置；Parviz等[10]对含水埋地管道进行了外爆冲击模拟，得到了管内水压力变化和振动情况。在含缺陷管道破坏问题方面的研究中，辛健[11]通过试验研究得到了不同形状的贯穿缺陷在内爆载荷作用下的破坏形状，并用数值模拟进行了对比验证；陈良路[12]通过数值模拟方法研究了炸药直接接触对管道的冲击，得到了管道局部的花瓣状破口，且破口翻唇程度与壁厚有关；田晓建等[13]通过SPH-FEM耦合法研究了含缺陷输气管道在外爆载荷下的冲击破坏情况，得出缺陷尺寸与位置相比缺陷深度和爆心距对管道的扰动程度小的结论；Jayasinghe等[14]探究了爆炸载荷作用下含缺陷圆柱桩基的失效问题，得到了桩基断裂时的极限爆炸载荷；杜洋[15]利用数值模拟方法研究了管内可燃气体爆炸后含缺陷管道的动态断裂过程，并提出了管道爆炸后果的评估方法；Aursand等[16]通过自编程序实现了CO2气体压力致使管道缺陷破坏并伴随气体泄漏减压的过程，提出了预防CO2管道破坏的方法。

通过上述文献对比分析可以看出，现有对埋地管道在外爆载荷作用下的破坏研究主要集中在管道受到冲击后的变形和应力分析，且部分研究对象为不含缺陷的管道。而对于管道的破坏问题，材料的损伤是一个动态过程，对由冲击引起的管道厚度方向损伤积累过程导致破坏的研究目前还很少。因此，通过建立土体-含缺陷管道-流体的多相耦合动力学模型，分析含缺陷埋地管道在外爆载荷作用下缺陷沿管道厚度方向的损伤程度，并探讨改变缺陷深度与形状、炸药量与爆心距对含缺陷管道破坏的影响，可为管道的破坏预测提供参考。

1 含缺陷埋地管道多相耦合模型和管道破坏判定准则

1.1 含缺陷埋地管道多相耦合模型

建立土体-含缺陷管道-流体的外爆多相耦合动力学模型，管道长度5 m，外径1 219 mm，壁厚16.3 mm，并在管道外表面正中间预置一个长100 mm、宽度20 mm的菱形缺陷，该缺陷深度8 mm，约等于壁厚的一半。土体长为5 000 mm，宽高都为3 219 mm，并在土体上表面放置正方体型TNT炸药，边长600 mm(352 kg)，管道缺陷正对着炸药。所有模型均采用实体单元进行离散，管道与土体采用拉格朗日网格描述，炸药和流体采用欧拉网格描述。土体模型用六面体单元划分，单元尺寸为50 mm，长度方向约100个单元，宽高方向约60个单元，可保证土体准确传播爆炸波。管道缺陷处用四面体单元划分，单元尺寸为3 mm，可准确描述缺陷的形状并保证厚度方向有多层单元；管道其余部分用六面体单元划分，单元尺寸为20 mm，在减少单元数量的同时，确保管道受力的精确性，混合单元的过渡部分为渐进过渡，避免应力梯度的突变。欧拉网格要小于拉格朗日网格，单元尺寸为20 mm，确保流固耦合作用的正常传递，网格尺寸经过网格无关性验证，精度满足要求。土体除上表面为自由边界外，其他面均施加无反射边界条件，流体设定12 MPa输送压力，起爆点设在炸药正方体中心，整体模型见图1。

炸药采用TNT材料模型，并采用JWL状态方程来描述爆炸过程，其状态方程表达式为：

(1)

式中P——TNT爆炸压力,Pa；

A，B，R1，R2，w——TNT材料常数；

η——无量纲量，η=ρ/ρ0;

ρ——密度，kg/m3;

ρ0——参考密度,kg/m3。

各参数具体数值见表1。

(a)埋地管道外爆多相耦合整体模型 (b)管道与预置缺陷模型

图1 含缺陷埋地管道外爆多相耦合动力学模型

表1 TNT材料参数

在高速剧烈载荷下，管道材料应表现出应变率强化效应，故采用Johnson-Cook强化模型来描述管道的塑性性能，该模型广泛用于描述材料在高速载荷冲击下的屈服与强化[17-18]，其表达式为：

(2)

式中Y——屈服强度,Pa；

A，B，C，m，n——管道材料常数；

εp——有效塑性应变,s-1；

TH——温度,K;

本文不考虑温度影响，管材选取X90钢的参数[18]，具体数值见表2。

表2 X90钢材料参数

土体材料采用Drucker-Prager材料模型模拟土壤，它考虑了材料的剪切膨胀特性，并可以描述土壤的长期塑性变形[19]，表3为土体材料参数。

表3 土壤材料参数

管内流体为天然气，采用理想气体状态方程来描述：

p=(γ-1)ρe

(3)

式中p——天然气压力,Pa；

γ——绝热系数，γ=1.4；

ρ——天然气密度，kg/m3,ρ=0.717 kg/m3。

1.2 土体-管道-流体载荷传递方法

管土之间采用接触描述，将管道与土体轴向摩擦系数视为无穷大，即土体表面与管道外表面之间没有切向滑动位移，同时管土界面法向接触采用罚函数[20]的方式进行计算；管内流体与管道内壁的相互作用采用流固耦合算法[21]，模型整体的多相耦合作用如图2所示。

图2 炸药-土体-管道-流体多相耦合作用示意

爆炸波冲击土体使土体向管道面运动，管道与土体之间会产生一个接触力FN，土体接触面穿透管道目标面时产生一段穿透距离xp，而罚函数算法会在接触面与目标面节点之间加上一个弹簧，该弹簧具有罚刚度kN，此时接触力FN可表示为：

FN=kNxp

(4)

当弹簧罚刚度kN很大时，穿透距离xp很小，以至于xp接近为零，此时接触力为准确值，管道与土体之间则通过这种方式传递接触力。

管内壁与流体之间的耦合界面上需要满足运动学条件：

Df=Ds

(5)

式中Df,Ds——流体、管道位移向量。

流体的力通过下式传递到管道结构上：

τf=Hfsτs

(6)

式中τf——管道内流体压力向量；

Hfs——管道内流体域到管道结构域的传递矩阵；

τs——管道结构的压力向量。

管道结构接收到流体的力会产生变形，并通过下式反馈给流体：

df=Hsfds

(7)

式中df——管道结构的位移向量；

Hsf——管道结构域到管道内流体域的传递矩阵；

ds——管道结构域变形后流体边界处的位移向量。

1.3 管道破坏判定准则

预置缺陷穿透则视为管道被破坏，而缺陷的穿透本质是材料的破坏。材料的破坏不是瞬间完成的，而是有一个损伤积累过程。管道受冲击作用时，当应力超过材料的屈服强度，该处会出现不可恢复的塑性变形；随着冲击作用的加深，塑性应变会积累，而材料本身的性质决定了在外力作用下断裂应变的大小。随着爆炸波的传播，管道外壁塑性应变会逐渐增大，塑性应变积累到一定程度，达到了断裂应变，材料即被破坏；当内壁塑性应变也达到了断裂应变，管道即被穿透。取塑性应变和断裂应变的比值作为损伤值，用来衡量材料破坏的程度。管道受到爆炸载荷这种高速剧烈的冲击载荷引起的破坏，断裂应变必须要考虑到应变率的影响。由于管道材料使用了Johnson-Cook强化模型来描述管道的应变率效应，故相应地使用Johnson-Cook失效准则来判定材料的破坏与否。Johnson-Cook失效准则的表达式为：

(8)

式中εf——管道单元有效断裂应变；

σ*——管道单元的应力，Pa；

T*——温度,K；

D1～D5——材料本身的参数，通过试验来测出，这里使用X90动态断裂试验拟合出来数据[20]，具体参数见表4；

Δε——管道有效塑性应变增量；

D——损伤积累比率。

表4 X90钢的Johnson-Cook失效准则参数

式(8)中3个括号部分分别描述压力效应、应变率效应和温度效应，综合三者计算出爆炸波冲击管道所产生的断裂应变。

管道的变形超过弹性阶段后，产生的有效塑性应变与断裂应变的比值为该时刻管道材料的损伤比率，通过求和的方式描述在一段时间内管道的损伤积累程度即为D值。D值表示为比率形式，取值范围为0～1，当达到1时，管道单元被视作失效并删除。当壁厚方向的单元损伤都达到了1，管道即被穿透。

为了监测管道的损伤积累比率D值，如图1(b)所示设置两个监测点和一条路径，A点为缺陷表面中心，B点为缺陷在管道厚度方向内壁位置；路径L由点A指向点B，目的是为了监测厚度方向的D值变化。

2 埋地含缺陷管道在外爆载荷下破坏的计算

对于管道结构，在前5 ms仅有管内流体作用，从第5 ms起爆炸波开始在土中传播，第6.75 ms时爆炸波刚刚传播到管道缺陷处，第10 ms时爆炸波已穿过整个管道截面，计算总时间为15 ms，其中爆炸波的传播时间共计10 ms。爆炸波的传播过程与相应时间管道缺陷处的损伤积累比率见图3。

图3 爆炸波的传播过程与相应时间管道缺陷处的损伤积累比率图

由图3(a)可以看出，第5 ms时炸药爆炸，爆炸波压迫土体使之变形形成凹坑，但爆炸波未传播到管道表面。此时管道仅受流体压力作用，管道整体的损伤集中在缺陷处。缺陷处的损伤集中在尖端，最大值为0.049，说明损伤很小。图3(b)所示第6.75 ms时，土体的凹坑变大，爆炸波此时经过土体刚刚到达管道表面，管道迎爆面由于受到爆炸波的压力，且与流体压力方向相反，损伤分布表现出了中间大、两边小的趋势。缺陷处A点处于正中心，此处的损伤值最大，达到了0.83。图3(c)所示第10 ms时，土体凹坑变得更大，爆炸波前缘经过管道整个截面，管道缺陷位置处的土体向下挤压管道，此时管道已被挤压出一个凹坑，凹坑处的损伤值大于其他部位。缺陷处出现了大面积的损伤，A点处损伤值仍然最大，达到了0.92。图3(d)所示第15 ms时，爆炸波穿越了整个土体，管道的迎爆面变形达到最大，此时缺陷处损伤加剧，A点处损伤积累比率达到1，单元被删除，管道被破坏。

由图3可以看出，贯穿位置发生在A点处，为了分析缺陷的贯穿情况，取A点、B点的损伤积累比率D随时间变化的曲线如图4所示，并取路径L在不同时刻的D值分布如图5所示。

图4 监测点处损伤积累比率随时间变化曲线

由图4可以看出，5 ms之前只有流体压力作用，故两个点处的D值很小。第5 ms炸药爆炸，大约第6.75 ms时爆炸波到达管道表面，两个点处D值有明显的升高。点A在缺陷外表面，比B点更早、更剧烈地受到爆炸波的冲击，因此A点处的D值要比B点大。从6.75～10 ms之间，由于爆炸波接触到管道表面后，爆炸波压力一开始并不很大，管内流体会抵消一部分爆炸波的能量，D值升高的速度变慢；随着爆炸波的传播，爆炸波压力上升，流体对爆炸波的抵抗能力越来越小，D值升高速度变快，因此导致D值在上升过程中，曲线呈阶梯性上升而不是线性上升。10 ms后爆炸波前缘已经过缺陷处，穿过整个管道截面，此时管道的受力又以流体为主导，所以D值趋于平稳。两点的D值在12 ms后都达到了1，为极限D值，此时管内外单元完全失效被删除。

图5 不同时刻路径L上的损伤积累比率分布

图5提取了4个时间点的路径L的D值，路径L平分为50个点，代表缺陷厚度方向的50份数据。第6.75,10,15 ms时间点为图3中爆炸波作用于管道的3个关键时刻，第11.25 ms为A点刚刚破坏而B点尚未破坏的时刻。图中4条不同时刻曲线都是呈线性趋势，A点比B点的D值大，说明D值在壁厚方向由外向内越来越小，即缺陷在爆炸波作用下的损伤在壁厚方向由外向内越来越小。随着时间的增长，曲线的斜率越来越小，说明两点之间的差值越来越小，表5给出了4个时刻下A点、B点处的D值。

从表5可以看出，第6.75 ms时A点D值比B点大得多，约为B点的2.5倍。随着时间的增长，B点D值增大，第10 ms时A点D值只比B点大6.9%。结合图4和表5，11.25 ms时A点已达到极限D值，B点D值为0.98，缺陷接近贯穿；第12 msA点和B点的D值已经为1，说明第12 ms时缺陷已贯穿，因此12 ms之后的曲线呈水平状态。

表5 不同时刻监测点D值汇总

本节通过计算爆心距为1 m的352 kg炸药，对含深度为8 mm的菱形缺陷管道的破坏情况分析，得到了缺陷处A,B两点及其厚度方向的损伤积累比率D的变化规律。而不同缺陷深度、形状和炸药量、爆心距等因素会影响缺陷的破坏，因此以下在该模型的基础上做不同工况的对比，以研究不同因素对缺陷破坏的影响。

3 不同预置缺陷对管道破坏的影响

3.1 初始缺陷深度对破坏情况的影响

分别设定缺陷深度为4,8,12 mm，缺陷形状为菱形，炸药量与爆心距不变，得到3种缺陷深度下路径L的D值曲线如图6所示，表6为3个时刻下各监测点处D值汇总。

图6 不同缺陷深度下第10 ms时路径L的D值分布

当爆炸波完全传播穿过管道截面的时刻(第10 ms)，从图6和表6可以看出，随着缺陷深度的增加，在路径L上D值曲线斜率逐渐降低，损伤比率也逐渐降低，说明在管道缺陷的厚度方向上，损伤积累随深度的增加而增加。这是由于缺陷深度越深，管道壁厚越薄，对爆炸波的抵抗能力越弱。对于4 mm缺陷深度，全程没有被贯穿，8 mm缺陷深度在15 ms时已经贯穿，12 mm缺陷深度在10 ms时已经贯穿。

表6 不同缺陷深度下各监测点处不同时刻的D值汇总

由表6可看出，对于4 mm深度缺陷，6.75 ms时路径L损伤比率高至318%，说明缺陷处两点D值的差值大，厚度方向损伤积累较小；10 ms时降至56.8%，损伤积累增大；15 ms时A点最终的损伤值为0.716 4，未达到极限值，说明4 mm深度缺陷未破坏。对于8 mm深度缺陷，6.75 ms时路径L损伤比率为251%，相比318%变小了，说明管道壁厚减薄，损伤积累增加了；10 ms时降至6.9%，下降程度和速度都很大，此时A点损伤值达到0.927 6，接近极限值；11.25 ms时A点达到了极限值而B点未到极限，该时刻缺陷即将贯穿，12 ms时就已经贯穿了。对于12 mm深度缺陷，6.75 ms时路径L损伤比率为176%，8.9 ms时就降到了12.1%，并且A点达到了极限值1；10 ms 时两点都达到极限，缺陷贯穿，此时4 mm和8 mm深度缺陷都未贯穿，说明12 mm深度缺陷不仅不能抵抗爆炸波的冲击，而且受到冲击后的损伤积累也很快，更加表明缺陷深度增加的危险性。

3.2 初始缺陷形状对破坏情况的影响

分别设定缺陷形状为菱形、椭圆形、圆形，缺陷深度为8 mm，尺寸如图7所示，得到3种缺陷形状下路径L的D值曲线如图8所示。表7为4个时刻下各监测点处D值汇总。

图7 不同缺陷形状尺寸

图8 不同缺陷形状下第10 ms时路径L的D值分布

表7 不同缺陷形状下各监测点处不同时刻的D值汇总

从图8和表7可以看出，第10 ms时椭圆形缺陷的路径L曲线D值要远远小于菱形和椭圆形的曲线；椭圆形和菱形缺陷的A点损伤值相差0.05，圆形和菱形缺陷的A点损伤值相差为0.7；但在厚度方向，圆形缺陷的路径L损伤比率为94.5%，远远大于其他两种缺陷。这是由于菱形缺陷有尖角，会在尖端产生损伤值集中，导致缺陷处的损伤整体偏大，椭圆形缺陷也会在长轴端产生一定程度的损伤集中。圆形缺陷将损伤近似平均地分布在了边缘处，使中心点A处的损伤值降低，圆面面积大比较均匀地吸收了爆炸波的破坏作用，因此在厚度方向损伤值下降得比较快。

由表7可以看出，菱形缺陷在第6.75 ms时A点损伤值为0.433 1，第10 ms时为0.927 6，提高了1.14倍；椭圆形缺陷在第6.75 ms时A点损伤值提高了1.16倍；圆形缺陷在第6.75 ms时A点损伤值仅提高了0.25倍。这说明菱形和椭圆形缺陷在爆炸波的作用过程中受到的损伤相近，圆形缺陷受到的损伤相比之下较小。第15 ms时菱形缺陷已经贯穿；椭圆形缺陷A点损伤值达到了0.908 2，路径L损伤比率为2.7%，近乎贯穿；圆形缺陷此时的A点损伤值仅为0.230 7，路径L损伤比率为82.8%，远达不到贯穿效果。

4 炸药量和爆心距对缺陷破坏的影响

4.1 不同炸药量作用下管道缺陷贯穿的情况

分别设定炸药量为104,203,352 kg(分别对应炸药立方体模型边长400,500,600 mm)，爆心距为1 m，缺陷形状为菱形，深度为8 mm，得到3种炸药量下路径L的D值曲线如图9所示。表8为4个时刻下各监测点处D值汇总。

图9 不同炸药量下第10 ms时路径L的D值分布

表8 不同炸药量下各监测点处不同时刻的D值汇总

由图9和表8可以看出，第10 ms时，3种炸药量的路径L曲线斜率随炸药量的增加而减少，D值随炸药量增加而增加。第6.75～10 ms时间段内104 kg炸药作用下管道的路径L损伤比率从231%降至207%，平均每毫秒下降7.4%；203 kg炸药作用下平均每毫秒下降53.5%；352 kg炸药作用下平均每毫秒下降75%。下降速度越快，表明B点D值接近A点的速度越快，这说明在该时间段内，炸药量越大，管道的破坏速度越快。第10～第15 ms，352 kg炸药作用下平均每毫秒下降1.4%，203 kg作用下平均每毫秒下降0.9%，说明该段时间内管道的破坏速度变慢；而104 kg炸药作用下平均每毫秒下降27.6%，说明该段时间内管道破坏速度相较352 kg和203 kg炸药量的作用要快。这是由于损伤积累比率是塑性应变和断裂应变的比值，大炸药量下塑性应变积累速度快，D值升高速度变快，因此6.75～10 ms之间炸药量越大，损伤比率降速越快；而10～15 ms之间，由于塑性应变不能无限增加，故大炸药量作用下剩余塑性变形很少，导致损伤比率降速很小。并且由式(8)可以看出，断裂应变与管道的应力和应变率有关，但并不是二者的线性函数，因此炸药量的线性增加并不能表现出管道损伤的线性增加。

4.2 不同爆心距下管道缺陷贯穿的情况

分别设定爆心距为1,2,3 m，炸药质量为352 kg，得到3种爆心距下路径L的D值曲线见图10。表9为4个时刻下各监测点处D值汇总。

图10 不同爆心距下第10 ms时路径L的D值分布

表9 不同爆心距下各监测点处不同时刻的D值汇总

由图10、表9可看出，第10 ms时曲线D值随爆心距的增加而减少，2 m和3 m爆心距的路径L斜率都小于1 m的曲线，1 m爆心距的A点损伤值为0.927 6，比2 m的损伤值0.089 7高出967%，比3 m的损伤值0.041 4高出2 140%，说明爆心距的增加会极大地减少A点处的损伤，这是由于2 m和3 m爆心距工况下的炸药爆炸后的冲击波在传播过程中能量耗散得很严重，导致到达缺陷时的能量大大减小，因此曲线D值整体都很小。

由表9可看出，第6.75 ms时，2 m和3 m爆心距下，由于炸药距离变远，爆炸波还未达到管道表面，此时的损伤值由流体产生，且损伤值相同。第6.75～10 ms之间，2 m爆心距下路径L损伤比率从95.1%减少到15.5%，而3 m爆心距下从95.1%增加到132%，说明此段时间内2 m爆心距下的损伤累积速度变快，3 m爆心距下的损伤累积速度已经变慢。第10～15 ms之间，2 m爆心距下A点D值增加了18.4%，B点仅增加了4%，路径L损伤比率从15.5%增加到了31.4%；3 m爆心距下A点增加20%，B点仅增加2.8%，路径L损伤比率从132%增加到了161%，损伤累计速度都已经变慢。与不同炸药量作用工况同理，爆心距的线性变化未表现出损伤积累的线性变化。

以上通过不同炸药质量和爆心距的研究得知，炸药量变化与爆心距变化相比，每增加100 kg左右的炸药量后管道损伤值仅为增加前的1～2倍，说明改变爆心距对缺陷破坏的影响要比改变炸药量更敏感。

5 结论

(1)考虑了含缺陷管道与管内流体的流固耦合作用、土体与管道间的接触作用，建立了土体-含缺陷管道-流体的外爆多相耦合模型。整体模型采用实体单元进行离散，采用考虑应变率效应的Johnson-Cook强化准则描述管道性能，并使用基于损伤积累的Johnson-Cook失效准则来计算管道的损伤。通过计算管道缺陷处在外爆载荷作用下的损伤积累比率，分析了管道破坏的损伤积累过程，并以该失效准则为标准判定管道的破坏程度。

(2)以含8 mm深度菱形缺陷的管道为例，对管道在352 kg炸药量、1 m爆心距的外爆载荷下缺陷破坏的过程进行计算，监测A点、B点的损伤积累比率值D，得到了两点处和路径L的D值变化规律。分析3个关键时刻的D值，发现缺陷损伤在壁厚方向外侧大于内侧，以路径L损伤比率作为壁厚破坏的程度，增幅越小、破坏越大。随着时间的推移，监测A，B两点之间的损伤差值越来越小，路径L损伤比率由251%降为6.9%，最终降为0，缺陷贯穿。

(3)通过对不同缺陷深度的计算结果表明，缺陷越深，管道越容易破坏，且破坏的时间越早；通过对不同缺陷形状的计算结果表明，菱形、椭圆形、圆形3种缺陷在爆炸波刚开始作用时的损伤相差不大，随着爆炸波的传递，菱形缺陷被贯穿，椭圆形缺陷所受最大损伤达到了0.9，圆形缺陷所受损伤最小，说明圆形缺陷对爆炸波的抵抗能力最强。同时，缺陷深度的变化比形状变化在外爆载荷作用下破坏影响更大。

(4)探讨了不同炸药量和爆心距对缺陷破坏的影响，结果表明：炸药量越大，管道损伤值越大，且管道损伤速度随炸药量增加而增加；不同爆心距下由于炸药的能量耗散严重，损伤值随爆心距的增加而减小。同时，1 m爆心距增加到3 m后的管道损伤值比1 m爆心距高出2 140%，而每增加100 kg左右的炸药量后管道损伤值仅为增加前的1～2倍，说明爆心距变化比炸药量变化对缺陷的破坏更敏感。