综合管廊内天然气泄漏及爆炸过程数值模拟分析

王 奔

(怀化市消防救援支队，湖南 怀化 418000)

0 引言

随着我国城市进程的不断推进，作为现代化城市的“重要生命线工程”地下综合管廊在许多城市的新城区以及新兴二三线城市成为热点建设项目[1]。目前，我国建设的管廊内部环境复杂，包括天然气管线、电力管线、热力管线等多种设施，致使管廊内部的火灾、爆炸等事故发生的危险性进一步增大。火灾、爆炸等事故一旦发生，无论是对管廊内部设施还是对社会经济财产都会造成严重的损失。如2008年2月内蒙古赤峰地下管廊内发生电缆火灾，致使城市大面积断电包括多家重大化工企业，造成巨大的经济损失；2010年7月南京市地下管道发生丙烯泄漏爆炸，直接带来5 000万的经济损失，并导致30余人死亡，100多人受伤；2016年10月德国一家化学公司的地下管廊在维修时，由于工人操作不当，发生乙烯泄漏爆炸，造成10人伤亡。因此，为提高管廊在火灾爆炸发生时的安全保障，有必要对管廊内部的天然气泄漏以及气体爆炸下的冲击波变化规律进行研究。

国内外学者针对管线气体泄漏、气体燃烧爆炸以及管道内的冲击波传播，进行了大量研究。最终气体泄漏高斯扩散模型被提出，以烟团模型表述气体连续泄漏扩散[2]。Lisboa将高斯模型结合计算机模拟对气体扩散规律进行研究，并建立小孔泄漏或管道断裂的泄漏压降模型[3]。杨昭等[4-5]针对管道气体的泄漏速率建立数学模型进行公式推导，并对气体的喷射过程进行分析。李又绿等[6]对典型的气体扩散进行数值模拟，并优化了气体泄漏扩散模型。就气体爆炸方面，Chapman等最早提出C-J爆炸理论模型，详尽地表述了气体的爆炸过程[7]。Fairweather等[8]针对圆形管道内部的甲烷-空气气体，进行气体爆炸传播过程的试验研究；周凯元等[9]则以丙烷-空气混合气体为研究对象，对在直线管道中的爆炸冲击波规律进行试验研究。Chen等[10]建立管廊试验模型，研究了在地震冲击波作用下的破坏规律；毕明树等[7,11]采用数值计算方法研究管道内气体爆炸过程，得到受限空间爆炸压力和升压速率；吴志远等[12]研究了多元混合气体定容、绝热情况下的最大压力简化计算方法。以上学者对管道空间可燃气体泄漏以及爆炸过程进行了大量研究，以数值计算方法居多，主要集中在管道气体的泄漏速率、管道爆炸压力发展、冲击波传播规律等方面。本文通过数学模型计算，探讨管廊无通风条件下天然气管线泄漏模型，研究小孔径条件下的泄漏速度及泄漏量，并计算甲烷-空气混合气体爆炸参数；然后采用数值模拟的方法，构建管廊有限元模型，分析特定燃气泄漏爆炸条件下，管廊舱室内的冲击波传播规律，研究测点的超压变化规律及衰减规律。

1 管廊内天然气泄漏及爆炸计算模型分析

1.1 气体泄漏计算模型

笔者以管廊内天然气管道输送过程的射流形式泄漏为背景，结合气体泄漏模型，并且将气体泄漏的过程视作等熵流动过程，对管廊天然气舱室内部泄漏量和泄漏速度进行分析。如图1所示，图中ρ0、p0、T0、v0的含义分别为管道外气体的密度、压强、温度和比体积，ρ1、p1、T1、v1则分别表示管道内气体的密度、压强、温度以及比体积，A为泄漏面积。

图1 天然气管道气体泄漏示意图

(1)

式中，cf0为气体的泄漏速度；k为气体等熵指数，即比热比，为定压比热与定容比热两者的比值；Rg为气体常数，Rg=R/M0，M0为气体摩尔质量，R为普适理想气体常数，通常取R=8.314 J·mol-1·K-1。

天然气的泄漏量可由式(2)算得：

(3)

式中，α为泄漏射流速度系数，与泄漏形式有关，一般取值在0.6～1之间，缝状口取0.65，圆孔取0.98；Mf为泄漏口的质量流量，kg·s-1。

本文中的管道内部压力根据天然气管道设计现行规范设定，按中级压力等级最低值取0.4 MPa，管道外的压力可按标准大气压取值。根据《城市燃气设计规范》的规定可知，有空气调节的建筑物内取20 ℃，故在管廊内部的天然气管道的输送温度取293 K。

1.2 混合气体爆炸参数计算

本文的爆炸参数计算基于混合气体的合理简化而得，具体假设如下：(1)混合气体为甲烷、空气的均匀混合气体，爆炸模式为稳定的C-J爆轰；(2)假设泄漏的天然气与空气均匀混合，并充满一段管廊。以甲烷的完全反应而定，空气中的氧气体积浓度为21%，即甲烷气体与空气的体积比为9.6%。

甲烷-空气混合气体的C-J爆炸速度计算公式为：

(4)

式中，vDJ为爆炸速度，m·s-1；λ为混合气体比热比，计算得甲烷-空气混合气体的比热比为1.263；Qm为单位质量甲烷的反应生成热，kJ·kg-1。

爆炸波压力的计算公式为：

(5)

式中，PJ为爆炸波波阵面压力，MPa；ρ0为混合气体初始密度，kg·m-3。

1.3 模型建立及算法材料选取

对常见的城市地下管廊结构尺寸进行合理的简化建立综合管廊模型。混凝土、空气、混合气体及覆盖土壤均采用SOLID164实体单元，爆炸气体、空气、管廊结构的网格尺寸为0.2 m×0.2 m，管廊长度总长取100 m，爆炸气体设置在管廊中心，沿纵向10 m填充满管廊的混合气体，管廊模型的横截面如图2所示。在计算中为缩小计算量，设置管廊的1/2模型。管廊周围边界土层厚度取为3 m。空气边界设置为无反射边界条件，土壤表面设置为自由边界，以避免爆炸冲击波在空气边界处反射影响管廊内部空气冲击波的传播规律。

图2 管廊三维模型及横截面示意图(1/2模型)

在LS-DYNA中分别采用*MAT_HIGH_EXPLOSIVE_BURN和*MAT_NULL的关键字，来定义爆炸气体以及空气的材料属性，并采用*EOS_LINEAR_POLINOMIAL来描述其爆炸状态方程[14]，其中混合气体的爆速、爆压、初始内能由上文计算可得，分别为2 105.8 m·s-1、2.42 MPa、4 045 kJ·m-3。土体采用*MAT_SOIL_AND_FOAM模型，混凝土模型采用HJC模型，材料所采用的关键字为*MAT_JOHNSON_HOLMQUIST_CONCRETE，具体参数见表1。

表1 甲烷-空气混合可燃气体与空气材料模型参数

2 管廊天然气爆炸传播过程反应分析

2.1 管廊内冲击波压力云图变化

图3给出了爆炸后70 ms内五个不同时刻管廊内压力分布云图(1/2模型中混合气体的一侧)，展现了管廊内部的爆炸冲击波的传播过程。可以看出，混合气体区域爆炸压强最高，在气体发生爆炸后冲击波以椭圆形的形状向外侧传播，压强最大值是波阵面区域；虽然由于天然气管线的存在，管廊中心截面处冲击波云图与燃气管线所在截面的冲击波云图有所不同，但沿纵向方向上差别并不大。同时，管线的存在造成了冲击波在管廊上方的传播速度要高于下方，从10 ms之后的各个图形都可以看出。在冲击波波阵面迅速传过之后，传过的区域处压强下降甚至形成负压，从而导致部分区域形成负压，在管廊的壁面周围局部压强增大，如图3中20 ms、38 ms、50 ms时刻云图所示。由各个时刻云图中的峰值可以看出，在冲击波传播的过程中压强峰值不断下降，而导致其发生的原因是冲击波在传播过程中撞击上下壁面生成反射波相向而行，可相互叠加抵消，这类反射波抵消过程在冲击波向外传播中重复多次，造成管廊中心的冲击波强度在向外侧传播的过程中逐渐衰减。

图3 管廊中心截面冲击波压力云图(kPa)

2.2 管廊内部超压变化规律分析

图4给出了管廊内超压测点示意图，在截面中心处从混合气体边缘沿管廊纵向方向向外每隔5 m选取一个超压测点。分析各个测点的不同时刻超压变化情况，研究爆炸冲击波在管廊内传播过程中的超压变化规律。

图4 超压测点示意图

如图5所示，分别给出了测点1～8的超压随时间变化曲线图。可以看出各个测点的超压变化总体趋势，在冲击波经过测点时超极速上升，随后在较短时间内下降到一个较低的水平，其中测点1的超压峰值可达415 kPa，测点8的超压峰值可达125 kPa，而超压下降后保持在75 kPa左右，各个测点的衰减速率大小保持一致。由测点1的变化曲线，可以发现超压在下降的初始阶段呈现一定的波动，是由于测点1处于混合爆炸气体边缘而产生的冲击波在横向有着较强的反射波和斜向冲击波导致，其他测点并未出现这样的现象。由图5，测点的超压峰值随着距爆炸气体距离的增加而减小，在超压值大于150 kPa范围，超压衰减速率要大于超压值小于150 kPa的超压衰减速率，如测点6、测点7、测点8的超压衰减近似呈直线衰减。

(a)测点1～4

(b)测点5～8

根据各个测点的超压峰值，应用图形处理软件，对测点超压峰值与时间关系进行拟合，如图6(a)所示。同时，为了分析天然气管线对冲击波传播的影响特设置一无管线管廊的对照组模拟试验，并拟合相同测点的超压峰值与时间关系，如图6(b)所示。可以看出，有天然气管线的管廊中，超压峰值随时间变化拟合曲线方程为二次函数，即超压衰减速率随时间由大减小，最后趋于平缓。而无管线管廊中的超压峰值变化趋势并不相同，衰减的拟合曲线方程为一次函数，呈直线式衰减，衰减速率可视为匀速，说明在管廊内部存在天然气管线对于超压峰值的衰减规律具有一定的影响，会加速超压最大值的衰减。

(a)有管线管廊超压峰值与时间拟合曲线

(b)无管线管廊超压峰值与时间拟合曲线

3 结论

本文对地下管廊中天然气管道舱内的泄漏建立计算模型并计算爆炸参数，且应用LS-DYNA软件对综合管廊内的混合气体爆炸进行数值模拟研究，分析管廊内的冲击波传播过程以及内部测点的超压变化，得到结论为：(1)在管廊内填充的混合气体爆炸后，冲击波以椭圆形的波阵面向外侧传播，截面中心的压强强度要高于管廊其他区域；中心区域冲击波传播速度最快，并且由于天然气管线的存在，在下侧冲击波的传播速度较低。(2)测点最高峰值可达415 kPa，在冲击波传出管廊时，强度削弱了2/3。冲击波强度减弱的原因是横向上反射波多次叠加抵消造成的，冲击波传过区域可形成负压。(3)超压峰值随离混合气体的距离增大而减小，综合管廊内最大超压值随时间变化所拟合曲线函数为二次函数，有天然气管线的管廊中超压衰减速率要大于无管线管廊中的超压衰减速率。