浅埋三舱管廊甲烷爆炸的地面响应规律

王桂林，欧阳啸天，翟 俊，孙 帆

（1. 重庆大学土木工程学院，重庆 400045；2. 库区环境地质灾害防治国家地方联合工程研究中心，重庆 400045；3. 重庆大学三峡库区生态环境教育部重点实验室，重庆 400045）

城市地下综合管廊含有燃气管道，管道燃气泄漏爆炸事故时有发生，如2014 年中国台湾高雄燃气管廊爆炸[1]（图1），2016 年德国综合管廊燃气爆炸等[2]。由于管廊常位于人员密集区，地下管廊一旦发生爆炸，将会给地面人员的生命和财产造成巨大损失，因此，研究地下管廊燃气爆炸地面响应对安全生产和生活具有重要意义。

由于实验成本高且安全风险大，针对地下空间爆炸问题，学者们一般采用有限元数值模拟方法[3]开展研究。曲树盛等[4]利用AUTODYN 有限元软件模拟了地铁车站内爆炸波的传播过程和衰减规律，廖维张等[5-6]利用ANSYS/LS-DYNA 有限元软件模拟了爆炸波在二层地铁站的传播规律和中柱的动力响应，田力等[7]利用ABAQUS 和LS-DYNA 软件分析了地下空间爆炸产生的荷载对建筑物的影响。然而在解决由爆炸引起的岩土体大变形问题方面，有限元方法常采用输入应力波的方式模拟爆炸应力，这与气体爆炸冲击岩土体的实际情况存在较大差异，此外，岩土体大变形将造成有限元网格畸变。物质点法采用拉格朗日质点和欧拉网格双重描述，将连续体离散成一组质点，每个质点代表一块材料区域并携带该材料区域的所有物理信息，计算网格仅用于动量方程的求解和空间导数的计算[8]，同时不需要对物质交界面进行特殊处理，自动满足无滑移的接触条件[9]，能够直接模拟物体之间的冲击现象。这些特点使得物质点法在处理爆炸冲击[10-12]和岩土大变形问题上[13-14]较有限元方法更有优势。

图1 2014 年高雄地下丙烯管道爆炸现场[1]Fig. 1 Explosion scene of the underground propylene pipeline in Kaohsiung in 2014[1]

本研究依托重庆市某地下综合管廊试点工程，采用点火增长模型，运用物质点法模拟浅埋管廊甲烷气体泄漏后爆炸冲击管廊本体结构和围岩的过程，以获得地面压强与位移的响应特性。

1 模型设置

1.1 接触算法

标准格式的物质点法可以自动处理物体之间的无滑移接触[8]，但为考虑物体之间的分离和滑动，Bardenhagen 等[15]将接触算法引入物质点法中，通过在网格节点上引入多重速度场的方式，模拟物体之间由于速度不同而产生的分离和滑动。该算法中，接触判据表达式为

如果物体B 与C 之间发生接触和穿透，则需要施加接触力以防止穿透发生，物体C 对B 产生的接触力计算公式为

通过在单步计算中遍历对单个网格节点上产生影响的多个物体间的两两接触，将得到的接触力进行矢量叠加，从而实现处理多个物体之间接触问题。计算示意图如图2 所示。

图2 多物体接触计算示意图Fig. 2 Schematic diagram of multi-object contact algorithm

1.2 甲烷爆轰模型设置

甲烷的燃烧有定压燃烧、爆燃、爆轰等形式，其中爆轰最剧烈，且对周围环境的破坏性最强。本研究采用点火增长模型模拟管廊中泄漏甲烷气体的爆轰过程，并将两种常见的定常爆轰反应率函数结合使用。

1.2.1 反应率函数设定

Wilkins 函数表达式为

式中：ps为反应完全时爆炸产物压强函数，E0为气体的单位体积内能。

1.2.2 人工体积黏性

在爆轰过程中冲击波前后的物质点参数变化巨大，在应力项上施加人工体积黏性项q[16]，将爆轰波的强间断面化成一个急剧但连续变化的波阵面，起到光滑间断面的作用。

人工体积黏性计算公式表示为

式中：c0、c1为输入常数，分别取 3.0、0.3；c 为声速； ρ为密度；为体积应变率。

1.3 模型及参数

重庆市某地下综合管廊试点工程的综合管廊标准断面高3.3 m，宽9.0 m，管廊包括燃气舱、电力舱和水信舱3 个舱，如图3 所示，顶板位于途经道路设计路面标高以下约3.4 m。建立的模型如图4 所示。

图3 某地下综合管廊现场照片Fig. 3 Picture of an underground comprehensive pipe gallery

图4 数值模拟模型内部示意图Fig. 4 Internal schematic diagram of numerical simulation model

模型尺寸为40.0 m × 30.0 m × 17.0 m，质点均匀分布，质点间距0.5 m，正方体网格边长1.0 m，质点总数为157 632 个。燃气舱、电力舱和水信舱的截面尺寸分别为2.0 m × 3.0 m、4.0 m × 3.0 m、2.0 m × 3.0 m，管廊结构厚度均为0.5 m。由于实际情况中管廊上覆岩土体与管廊之间的黏结力小，故模拟中不考虑。参考Zhu 等[17]开展的混合气体爆炸模拟实验，设置发生泄漏的甲烷-空气混合气体充满管廊中部，长度为4.0 m 的区域，起爆点位于混合气体正中心。模型四周采用黏弹性边界模拟侧向无限岩土体，底部采用对称边界模拟基岩面，顶部采用自由边界模拟地面。

模型围岩和管廊本体结构模拟采用DP 强度准则，甲烷-空气混合气体模拟采用高能燃烧模型。参考工程地质勘察报告和文献[18]，岩土物理力学参数取值如表1 所示，混合气体参数取值如表2 所示。其中：E 为弹性模量， ν为泊松比， σt为 岩石抗拉强度， qφ、 Kφ为由摩尔库伦参数换算得到的材料参数，ψ为剪胀角，Vf为体积分数。

表1 围岩和管廊本体结构力学参数Table 1 Mechanical parameters of surrounding rock and pipe gallery structure

表2 甲烷-空气混合气体高能燃烧模型计算参数Table 2 Calculation parameters of high-energy combustion model of methane-air mixture

2 模型验证

参照 Zhu 等[17]、Xu 等[19]和 Qu 等[20]的实验，在一个长900.0 m、一端封闭、截面面积为7.2 m2的正方形隧道中，充满区域长度约为14 m、体积为100 m3、体积分数为10%的甲烷-空气混合气体。点火点在封闭隧道的一端。图5 为物质点法模型示意图，其中质点间距0.3 m，网格间距0.6 m，质点总数363 001。表3 为数值模拟结果及其与文献实验结果的对比，其中： δ为相对误差。

图5 模型截面示意图Fig. 5 Model section schematic

表3 数值模拟验证结果Table 3 Numerical simulation verification results

分析表3 可以发现，物质点法数值模拟结果与实验结果的偏差约为9%，较能反映实际情况[21]，产生误差的原因主要有以下3 点：（1）由于物质点法的算法特性，单个测量质点不能代表具体的物质形状，相较于实验中测量元件具有明确的尺寸，物质点法中测量质点得出的结果与实际情况存在出入；（2）本研究中质点间距为0.3 m，网格边长为0.6 m，对于精确模拟来说精度不够，然而随着物质点数量增加，网格划分密度增大，模拟结果会更加贴近实际情况；（3）实验中隧道与高速气体存在一定的摩擦力，这种摩擦力会影响爆炸压强，而具体区域的摩擦力大小受气体的流速和隧道壁的粗糙程度影响，致使数值模拟结果大于实际情况。

3 爆炸作用地面响应特性

3.1 地面压强响应特性

3.1.1 地面响应压强平面分布特性

图6 为混合气体引爆0.02、0.20、0.50 s 时对应的地面压强分布情况。由图6 可以看出，起爆后0.02 s，地面上以起爆点竖直方向为中心的椭圆形区域出现方向向上的响应压强，并向四周辐射，这是由于混合气体引爆以后，受到爆轰压力的管廊衬砌产生向上的形变，从而冲击地面造成的。起爆后0.20 s，爆炸应力波向四周传播，管廊上覆地面响应压强方向整体向下，这是由于地表为自由边界，应力波传播至地表后发生反射导致向下传播；部分区域出现应力集中现象，是由于应力波在管廊水信舱及电力舱空洞中发生了反射；紧邻管廊上覆地表的岩土体出现应力集中现象，是由于应力波在通过管廊和岩土接触截面时发生了折射。观察地面上整体响应压强分布情况，可以发现燃气舱一侧地面响应压强分布集中，数值较大；电力舱一侧响应压强分布范围更广，但数值较小。起爆后0.50 s，地面响应压强分布逐渐稳定，在模型四周出现的响应压强是由于吸收边界未能完全吸收反射波。地面中部出现的应力集中现象是由于管廊未考虑阻尼，在管廊内部反射的应力波未完全衰减。

图6 爆炸后不同时刻地面压强分布Fig. 6 Distribution of ground pressure at different times after explosion

3.1.2 管廊横向压强分布特性

图7 为管廊横向距起爆点水平距离-4.0、0、4.0、8.0 m 处质点压强随时间的变化曲线。由图7可知，从响应时间来看，地面质点的响应时间随距起爆点水平距离的增大而延长。管廊上覆地面质点（0、4.0、8.0 m）在经历初始响应压强峰值后均会出现方向向下的负压强，而围岩上覆地面（-4.0、-8.0 m）则不会出现这种情况。此外，除了起爆点正上方质点外，其他质点在爆炸过程中出现了由于接触界面反射作用形成的次生响应压强波峰，且随着距离的增大，次生波峰峰值逐渐增大。根据主波和次生波产生的阶段，可以将爆炸响应分成3 个阶段：爆炸响应阶段、次生响应阶段和渐趋稳定阶段。

图8 为管廊横向质点压强峰值随起爆点距离的变化曲线。由图8 可以看出，爆炸引起管廊横向的地面峰值响应压强受地下情况影响明显。其中，响应压强峰值的最大值偏离起爆点正上方，管廊上覆地面距起爆点2.0 m 处峰值压强最大，为37.16 kPa。从响应压强峰值分布来看，围岩上覆地面响应峰值的最大值小于管廊上覆地面，但随距起爆点距离的增大，其减小幅度小于管廊上覆地面。响应压强最大值峰值点的偏移，是由于浅埋管廊与岩土体接触面和舱体临空面的存在，使得爆炸应力波传输到管廊本体结构时发生折射，导致响应压强峰值出现偏移现象。

图7 管廊横向质点压强随时间变化曲线Fig. 7 Variation curves of particle pressure in the horizontal direction of pipe gallery with time

图8 管廊横向质点压强峰值随起爆点距离的变化Fig. 8 Peak particle pressure in the horizontal direction of pipe gallery with distance of initiation point

3.1.3 管廊纵向压强分布特性

由于模型在纵向对称，取纵向距起爆点水平距离0、2.0、4.0、6.0 m 处质点的响应压强随时间的变化曲线，如图9 所示。由图9 可知，地面质点的响应时间随距起爆点水平距离的增大而延长，但与横向的响应特性不同，纵向各质点在受爆炸作用后产生的次生波振幅较小，且随距离的变化不大；此外，处于管廊正上方纵向各质点的压强均竖直向下。

图10 为管廊纵向质点压强峰值随起爆点距离的变化曲线。从图10 中可以看出，沿管廊纵向的地面响应峰值压强随与起爆点水平距离的增大而经历先增大后减小的过程，且随着距起爆点距离的增大，压强峰值变化速率逐渐减小。最大值出现在距管廊起爆点2.0 m 处，峰值压强为26.29 kPa。出现这种现象的原因是：混合气体在横向尺寸较大，气体发生化学反应的先后不同导致应力波叠加，出现峰值压强最大值偏离起爆点；距离起爆点越远，混合气体中甲烷消耗得越多，释放的压强越小，因而对距离起爆点较远的地面质点影响更小。

图9 管廊纵向质点压强随时间变化曲线Fig. 9 Variation curves of particle pressure in the longitudinal direction of the pipe gallery with time

图10 管廊纵向质点压强峰值随起爆点距离变化Fig. 10 Peak particle pressure of the pipe gallery with the distance of the initiation point in the longitudinal direction

3.2 地面竖向位移响应特性

3.2.1 地面响应竖向位移分布特性

图11 为截取地表各质点爆炸作用下的位移情况。（1）从绝对位移来看，模型地面整体出现了竖直向下的位移，这是由于管廊可燃性气体爆炸时，应力波向四周围岩传播，管廊底部围岩受力产生压缩变形，使得地面发生沉降。（2）从地面隆起形态来看，在隆起范围方面，由于管廊燃气舱中甲烷爆炸作用，管廊燃气舱上覆地面立即相对隆起，而其他区域均下沉；在隆起幅度方面，靠近起爆点的管廊上覆土层隆起剧烈，而远离起爆点的管廊上覆土层隆起幅度较小。

图11 不同时刻地面竖向位移响应分布情况Fig. 11 Distribution of ground vertical displacement at different times

3.2.2 管廊纵向地表竖向位移分布特性

由于浅埋管廊爆炸中地表竖向位移远大于水平位移，因此着重分析地表竖向位移的响应。图12为管廊纵向距起爆点水平距离0、2.0、4.0、6.0、8.0 m 处各质点竖向位移随时间的变化曲线。由图12可知，管廊纵向质点竖向位移随距起爆点水平位移距离的增大而下降明显。结合3.1 节的压强分析，爆炸过程中管廊纵向地面各质点位移响应经历了3 个阶段：0～0.05 s，地面质点响应压强快速增加，而竖向位移增加缓慢，其中0 和2.0 m 处质点位移向上，4.0、6.0、8.0 m 处质点位移向下；0.05～0.39 s，地面质点响应压强度过峰值阶段，而质点竖向位移增加迅速，靠近起爆点正上方的质点由于模型整体沉降作用而呈现先上升后下降的过程；0.39～0.50 s，地面响应压强逐渐稳定至较小值，质点位移增加速率减小。

图13 为管廊纵向以地面上距起爆点最远处质点为参考对象，质点相对竖向位移随起爆点距离的变化曲线。由图13 可知，管廊纵向相对竖向位移随距起爆点距离呈较明显的阶段性特征，其中0～2.5 m 处地面向上剧烈隆起，管廊由于甲烷爆炸作用发生裂损，气体涌出直接冲击上覆岩土层；2.5～7.5 m 处为过渡阶段，地面隆起程度随距离的过渡性变化较为明显；7.5～12.5 m处为平缓阶段，地面隆起平缓，这是由于爆炸作用使得管廊整体发生振动引起地面隆起。整体上，管廊纵向相对竖向位移随距起爆点水平距离的增大而减小，相对位移峰值出现在起爆点正上方，相对隆起高度为1.3 m。

图12 管廊纵向质点竖向位移随时间变化曲线Fig. 12 Variation curves of vertical displacement of particles in the longitudinal direction of pipe gallery with time

图13 管廊纵向质点相对竖向位移随起爆点距离的变化Fig. 13 Relative vertical displacement of particles in the longitudinal direction of pipe gallery with distance of initiation point

3.2.3 管廊横向地表竖向位移分布特性

图14 为沿管廊横向距起爆点水平距离-8.0、-4.0、0、4.0、8.0 m 处质点相对竖向位移随时间的变化曲线。与纵向相似，除了起爆点正上方外，各质点竖直位移呈现3 个阶段：缓慢增加阶段、加速增加阶段和减速稳定阶段。

图15 为管廊横向以地面上距起爆点最远处质点为参考对象，质点相对竖向位移随距起爆点距离的变化曲线。管廊横向相对竖向位移最大值出现在起爆点斜上方，这是由于管廊上覆地面岩土体与管廊黏结力较小，而岩土体之间的黏结力使得横向位移最大值偏移起爆点正上方；各质点相对竖向位移随距起爆点的水平距离增大而减小，其中管廊上覆地面的竖向位移要大于岩层上覆地面，且位移减少速率更小；与各质点响应峰值压强呈现明显分段性不同，上述散点分布的分段性不明显，这是由于材料的连续性；随着质点距起爆点距离的增大，相对竖向位移减少速度逐渐减小。

图14 管廊横向质点相对竖向位移随时间变化曲线Fig. 14 Variation curves of relative vertical displacement of particles in the horizontal direction of pipe gallery with time

图15 管廊横向质点相对竖向位移随距起爆点距离变化Fig. 15 Rrelative vertical displacement of particles in the horizontal direction of pipe gallery with distance of initiation point

4 结 论

采用点火增长模型模拟管廊内甲烷气体泄漏爆炸冲击管廊本体结构和围岩的过程，研究了地面压强和位移响应规律，得出以下主要结论。

（1）在响应压强方面：爆炸过程中，管廊及围岩中出现由接触面反射和折射产生的次生应力波，这种次生波在管廊横向方向，振幅随距起爆点水平距离的增大而增大；而在管廊纵向上产生的次生波振幅较小，且随距离增大变化较小；管廊横纵方向上响应压强峰值发生点均偏离起爆点；管廊上覆地面在爆炸作用后出现由地面反射波造成的负向压强。

（2）在响应位移方面：爆炸作用造成整体地面沉降，但在起爆点中心附近地面隆起，这种隆起由管廊衬砌破裂、气体直接冲击岩土体形成的剧烈隆起，以及管廊整体振动形成的轻微隆起两部分组成，其中围岩上覆地面位移小于管廊上覆地面；管廊横向的竖向位移偏移于起爆点正上方。