0 引言

伴随“一带一路”和“油气管网”等政策的逐步实施，作为能源输送主要方式的管道将迎来蓬勃发展期。由于路由限制的影响，双管或者多管并行敷设方式不可避免。然而，管道不可避免地会受到腐蚀、第三方破坏^[1-3]及自然灾害^[4-6]等导致的多种损伤，容易引发泄漏，从而导致气云爆炸，与之相邻的并行管道极有可能受到威胁从而造成冲击破坏，引发二次事故。因此，开展爆炸载荷作用下埋地并行输气管道的动力响应研究就显得尤为重要。

当前，国外对爆炸载荷作用于并行邻管的研究相对较少，多以弹坑模型^[7-8]界定并行管道的安全间距。国内则开展了大量研究，笔者所在团队^[9-12]前期针对并行管道泄漏及物理爆炸对邻管的冲击破坏开展了一系列研究。此外，部分国内学者^[13-18]也开展了相关的研究工作。以上研究对于模型精度的验证大都基于经验描述，缺乏基于试验的精确验证，即使部分涉及试验验证，数值模拟的精度也有待进一步提高，其次在ALE数值分析中不能很好地描述土体变形与管道瞬态响应的内在联系。

为此，笔者首先基于土中爆炸试验提高ALE算法的精度；其次选用TNT当量模型估算气云爆炸能量，建立了管土-炸药耦合三维模型，分析了泄漏量和管间距对管道动力响应规律的影响；最后基于应力和应变失效准则评估了管道的安全状况。所得结论可为并行输气管道的安全评估提供参考。

1 小孔泄漏模型与TNT当量法 1.1 小孔泄漏模型

天然气管道的泄漏通常为小孔泄漏。图 1中管输流量为Q₁，1点为管段起始点，距起点L处的点3处发生泄漏，泄漏点与点2垂直于管道轴线。假设理想气体在管内做一维稳态绝热流动，且气体在孔口处做等熵流动。

依据理想气体状态方程和欧拉运动微分方程及临界压力比，得到临界流的小孔泄漏模型计算公式，即有：

(1)

式中：Q为气体泄漏率，kg/s; C_D为泄漏系数，取值范围0.61~1.00，圆形孔口取1.00；A为泄漏孔口面积，m²；φ为流速系数，表示实际流速与理论流速之比，由于实际流速小于理论流速，一般取0.97~0.98；p₁为泄漏起点的压力，Pa；k为气体绝热系数，燃气取1.3；M为气体摩尔质量，燃气取0.016 kg/mol；R为摩尔气体常数，取8.314 J/(mol·K)；T₁为泄漏点温度，K。

1.2 TNT当量法

定量计算蒸汽云爆炸能量的方法主要有TNT当量法、数值模拟、球形数值模型及自相似等方法。TNT当量法依据能量相当法则，将气云爆炸产生的威力等效为凝聚相TNT爆炸产生的威力，可直接通过数学转换方式进行转换。TNT当量计算公式为：

(2)

式中：W_TNT为TNT当量，kg；α为蒸汽云爆炸当量系数，取值范围0.01~0.10，燃气通常取0.04；W_f为参与爆炸的蒸汽云质量，kg；Q_f为可燃物燃烧值，取54×10⁶ J/kg；Q_TNT为TNT爆炸热，取4.52×10⁶ J/kg。

2 材料本构与管体失效准则 2.1 材料本构

土体选用MAT_SOIL_AND_FOAM模型，本构方程为：

(3)

式中：S_ij为Cauchy偏应力张量；δ_ij为Kronecker符号；a₀为与土体摩擦角相关的参数，取3.4×10⁹ Pa²；a₁为与土体黏聚力相关的参数，取7.033×10⁴ Pa；a₂为与土体动载效应相关的参数，取0.3；p为土体剪压力，GPa。

此外，建模中土体密度为1.8 g/cm³，剪切模量为63.8 MPa，体积模量为30 GPa。

埋地管线为西二线常用的X80管线钢，输送压力12 MPa，规格为ø1 219 mm×26.4 mm。本构方程采用随动硬化双线性弹塑模型，该模型遵从Von Mises屈服准则，表达式如下：

(4)

式中：σ为管材应力，MPa；E_s为弹性模量，取210 GPa；E_tan为切线模量，取1 500 GPa；σ_y为屈服应力，取590 MPa；ε、ε_e分别为应变和弹性极限应变。

建模时管材密度为7 900 kg/m³，泊松比为0.3。

建模时炸药选用MAT_HIGH_EXPLOSIVE_BURN本构模型与JWL状态方程，方程表达式为：

(5)

式中：p₁为爆轰压力，MPa；E₀为单位体积爆轰产物的内能，取0.06 J/m³；V为爆轰产物的相对体积；A、B、R₁、R₂和ω为常数，分别取值373.77 GPa、3.747 1 GPa、4.15、0.9和0.35。

建模中炸药密度为1.63 g/cm³，爆速为6 930 m/s，爆压为21 GPa。

2.2 管体失效准则

《输气管道工程设计规范》(GB 50251—2015)规定：工作载荷条件下，管道可能出现的最大应力应小于其屈服应力。第四强度理论综合考虑了管道的环向应力、轴向应力和径向应力的作用，表达式为：

(6)

式中：σ_φ为管道的环向应力，σ_r为径向薄膜应力，σ_θ为轴向薄膜应力，[σ_t]为管体许用应力，单位均为MPa。

管道变形包含椭圆化(度)变形、拉伸变形和压缩变形。评判管体因变形引发的失效常用椭圆化失效准则，表达式为^[19]：

(7)

式中：δ_θ为椭圆度，无量纲常数；δ_θ^crit为临界椭圆度，通常取3%^[20]；D_max为管道外径最大值，mm；D_min为管道外径最小值，mm。

3 计算模型 3.1 模型精度验证

加拿大防务研究所对100 g炸药在土中进行了一系列爆炸试验^[21]，试验装置及装置简图分别如图 2a和图 2b所示。试验主要包含一个装满土体的铁桶、土中不同埋深下的3个压力传感器以及大气中设置的5个压力传感器(见图 2b黑色方块)。研究仅针对地下部分，因此建立如图 2c所示的三维实体模型，模型顶面为自由表面，对称面上施加对称约束，其余面上施加非反射边界约束。

图 3为土中3个压力传感器的实测值与模拟结果对比图。

由图 3可知，模拟所得结果与试验值存在一定偏差，但整体接近，偏差分别为1.12%、7.00%和8.86%，处于工程可接受范围，因此该试验定量验证了数值模拟的可行性。

3.2 有限元模型

若处于同一高度上的3条输气管道等间距同沟并行敷设，中间管道发生小孔泄漏，高速气流与孔壁间产生摩擦生热，存在爆炸安全隐患。选取TNT当量法作为气云泄漏爆炸的估算模型，利用对称性建立整体有限元计算模型，如图 4所示。其中药包置于爆管形心处，药包形心至并行邻管形心的距离为R(爆心距)。管体有限元模型及网格划分如图 5所示。

整体有限元模型顶面为自由表面，在XOY及YOZ平面内施加如图 6a所示的对称约束。为了模拟无限域环境，在其余面上施加如图 6b所示的非反射边界条件，管道内壁施加如图 6c所示的内压载荷，载荷大小为12 MPa。

4 结果分析 4.1 压力波在土中的传播规律

假定燃气泄漏一定时间后，气云爆炸的等效TNT当量为15 kg，炸药为直径16 cm、高46 cm的柱形药包，建立爆心距离R为220 cm的有限元模型进行数值计算。图 7为药包起爆后压力波在土中的传播规律。由图 7可知：压力波在土中以球面波的形式向外传播，随着时间t的推移波前区域面积逐渐增大；压力波在t=1.5 ms时刻抵达管道，由于管道阻碍了其横向传播，因而发生扰流现象，向上透过地表耗散于大气中，向下朝土体下方扩散进而耗尽。此外，药包起爆后土中形成一个爆炸空腔，随着时间的推移腔体逐渐增大，地表产生土丘状隆起，腔体形状由球形向倒梨形过渡。

4.2 动力响应分析

考虑到气云泄漏爆炸的复杂性，基于能量相当法则得到不同泄漏时长下装药量G为36、52、68及84 kg，简化药包为半径12 cm，高度为49、70、92和114 cm的圆柱形，建立爆心距R为205 cm的数值模型进行仿真。

4.2.1 应力响应分析

由于爆炸载荷下管道的迎爆面为受力最敏感的部位，这里提取了不同炸药量下管体迎爆点的应力时程曲线(见图 8)。由图 8可知，药包起爆后，4种炸药量下管体迎爆点的应力时程曲线趋势基本一致，均呈现先减小后陡然增大，然后再急剧减小后小幅度增长，最后缓慢减小并逐渐趋于稳定的趋势。

这是因为管道起始时刻在纯内压的作用下处于应力拉伸状态，4种工况下的初始应力为256.5 MPa，药包起爆后，冲击波率先抵达管道，抵消了部分由于内压产生的拉应力，4种工况下的应力瞬间降至极小值。此后，冲击波得到源源不断的补充，管道受力又开始反向急剧增大，达到极大值后由于冲击波得不到补充，且管道变形部分恢复，应力又开始急剧减小达到第二处极小值，各工况下的应力极值如表 1所示。此刻之后由于爆腔继续扩大，对管道的作用力增大，管道应力又开始缓慢增大，最终随着能量的耗尽与爆腔扩张速率的减慢管道应力又开始减小。

将表 1中不同药量下管体迎爆点的最大应力进行函数拟合，拟合曲线近似呈线性关系。拟合公式如下：

(8)

式中：y为不同药量下管体迎爆点的最大应力，MPa；x为药包的质量，kg。

为了确保管道安全平稳运行，需满足应力校核准则，即管道的最大应力应小于其抗拉强度719 MPa。由公式(8)得到管道的极限抗爆能力(TNT当量)为58.4 kg。

4.2.2 位移响应分析

(1) X向位移响应分析。迎爆点X方向位移时程曲线如图 9所示。正对爆心管道截面上迎爆点、背爆点、管顶及管底4个观测点X向位移随时间的推移均呈现先增大后减小的趋势。以图 9为例，测点位移存在一极大值点，该点之前位移随着时间的变化而增大，之后随着时间的变化而减小。不同药量下各测点不同时刻对应的位移极值如表 2所示。由表 2可以看出：药量越大管体位移响应越明显；管底和管顶的位移响应基本一致，而迎爆点、管顶(管底)和背爆点的位移响应依次滞后。

(2) Y向位移响应分析。Y方向的位移时程曲线与测点位置相关，其中迎爆点与背爆点的位移随时间的变化呈现先增大后减小的趋势，仅是响应时刻及响应值存在微小差别。因此，笔者只给出了迎爆点的时程曲线，如图 10所示。

图 11为管底Y方向位移时程曲线。

由图 11可知，管底位移呈现由负向位移(先增大后减小)向正向位移(先增大后减小到再增大后减小)过渡的现象。这是因为球形冲击波及爆炸空腔的挤压作用，管道一开始产生负向位移，之后爆腔向上发展速度较其他方向快，爆腔失去球形对称性，迫使管道被抬升，产生向上的位移。图 12为管顶Y方向位移时程曲线。由图 12可以看出，管顶位移呈现先增大后减小的二次脉动现象。这是因为管道首先受到爆炸波的冲击应力急剧增大，之后冲击波得不到持续补充，管道受力开始急剧衰减，爆炸空腔继续扩张，对管道挤压作用增强，管道受力又开始增大，最终随着能量耗散，管道变形部分恢复，因而受力减小。各测点不同时刻Y向位移极值如表 3所示。

4.3 定量失效风险评估

管道作为国家能源输送的命脉，确保其安全运行具有重要意义，这里分别基于应变和应力失效判据评估不同装药量及不同爆心距下管体的运行状态，评估结果分别如表 4和表 5所示。

由表 4可以看出，同一爆心距下，管道的失效风险均随着装药量的增大而增加，基于应力的失效判据较应变失效判据相对偏保守。考虑到工程实际中管体本身的缺陷，研究工况下建议采用椭圆度失效准则。

不同爆心距下，管体各测点的参数时间历程与不同装药量下的规律大体相同，在此不做赘述，仅评估其安全状态。由表 5可知，管道的失效风险随着爆心距的增大而减小。这是因为爆心距越大，管道间的覆土越多，冲击波及土壤的塑性挤压作用越弱。为确保管道留有一定的安全余量，研究工况下建议管道的并行安全间距不小于240 cm。

5 结论

(1) 基于加拿大土中爆炸试验验证了ALE法的精度，ALE法可以较好地描述爆炸波在土中的传播规律，以及爆炸空腔及地表土丘的宏观演化过程。

(2) 同一爆心距下，应力和位移均随着装药量的增大而增大；应力随着时间的变化先增大后减小；X向位移随着时间的变化先增大后减小；Y向位移与时间的关系取决于测点的位置，迎爆点与背爆点的位移随着时间的变化先增大后减小，管顶位移随着时间的变化呈现先增大后减小的二次脉冲现象，管底位移随着时间的变化由负向位移向正向位移过渡。

(3) 基于应力的管道安全评定较应变失效准则偏保守，评定准则的选取应视具体工况而定。