地下资源的开发利用可以有效解决城市空间的扩大与土地资源紧缺之间的矛盾，而综合管廊作为城市地下空间利用的重要组成部分也取得了长足发展。综合管廊容纳多种市政管线，但由于燃气管道的特殊性，在管廊内必须设置独立的腔室^[1]。若燃气管道发生破坏致使燃气介质泄入燃气仓，又由于处于较密闭空间，气体容易积聚而形成爆炸环境。管廊一般处于市区建筑人口密集地段，一旦发生爆炸会导致电力信息中断、人员伤亡、地上建筑倒塌等一系列事故，对人们生活和社会秩序产生不利影响。对于密闭、半密闭空间内的气体爆炸，诸多学者在该领域开展了较多研究并取得一些成果。陈国华等^[2]通过模拟地下密闭空间，探究了密闭空间内燃气爆炸中冲击波变化规律，并得到了相应的爆炸波传播衰减公式。路长^[3]、秦波涛^[4]等研究了不同浓度的可燃气体在管道内的传播规律，卢捷^[5]、王成^[6]等研究了管道内混合气体爆炸火焰传播加速现象和加速机理。左哲等^[7]模拟了天然气在受限空间内的爆炸，并研究了空间尺寸、初始压力等因素对爆炸冲击波的影响规律。周宁等^[8]研究了不同点火能对管道内可燃气体爆炸的影响规律，杜扬^[9]、蒋新生^[10]等研究了分支结构对管道内可燃气体爆炸强度的影响。

与凝聚炸药相比，气体爆炸受环境影响较大，不同的密闭结构内可能有完全不同的爆炸形式^[11]。本研究以地下管廊为研究对象，利用LS-DYNA软件，基于流固耦合计算方法对其进行模拟分析，探讨管廊结构动力响应规律以及结构内冲击波传播规律和特点，为地下管廊结构抗爆设计提供部分参考。

1.   地下管廊爆炸计算模型建立

1.1   工程概况

平潭综合试验区环岛路管线工程中山大道至新湖路段中E型综合管廊，管廊截面宽度10700 mm，高度为3500 mm，混凝土结构外墙厚400 mm，内墙厚为300 mm。工程为盾构开挖，土层性质未改变，管廊底部至路面，为3层不同性质的土体，从上往下分别为素填土、全风化花岗岩以及砂土状强风化花岗岩，土壤参数见表 1, 混凝土材料参数见表 2。管廊共有3个仓室，燃气仓由于其特殊性独立成仓，位于管廊最右边。

表  1  土壤材料参数

Table  1.  Parameters of soil material

Material	Thickness/ cm	Density/ (g·cm-3)	Cohesion/ kPa	Internal frictionangle /(°)	Elastic modulus/MPa	Poisson'sratio
Plain fill	130	1.8	6	10	4.2	0.30
Fully weathered granite	210	1.9	20	25	20	0.24
Sandy weathered granite	620	2.0	30	32	54000	0.21

下载: 导出CSV 

| 显示表格

表  2  混凝土材料参数

Table  2.  Parameters of concrete

Density/(g·cm-3)	Elastic modulus/GPa	Poisson's ratio	Yield strength/MPa	Shear modulus/GPa
2.5	30	0.22	33.8	12.5

下载: 导出CSV 

| 显示表格

1.2   有限元模型及参数

有限元模型共有4部分组成，分别为爆炸气体、空气、土体以及混凝土。采用SOLID164实体单元建模，混凝土结构用Lagrange算法，爆炸混合气体、土体以及空气采用欧拉网格，欧拉网格包含Lagrange网格，单元使用多物质ALE(Arbitrary Lagrange Euler)算法，采用LS-DYNA有限元程序提供的CONSTRAINED_LAGRANGE_IN_SOLID耦合方式进行流固耦合计算，以模拟爆炸冲击波与结构的相互作用。

建模时采用cm-g-μs单位制，在保证合适的计算精度，同时考虑合适的计算机成本要求，本模拟采用20 cm的控制尺寸划分网格，由于混凝土结构为主要的受力构件，且容易产生应力集中，为进行有效合理的模拟，对网格进行加密。考虑对称性问题，以起爆点的xOy面作为对称面建立1/2模型，模型对称面采用对称约束，其余界面设置为无反射边界以模拟无限域的环境条件，有限元模型及相关尺寸如图 1所示。

图  1  有限元模型

Figure  1.  Finite element model

下载: 全尺寸图片 幻灯片

天然气的主要成分为甲烷(CH₄)，天然气的爆炸即是甲烷管廊内甲烷-空气(CH₄-Air)混合气体在一定点火能量下发生爆轰。假设爆炸性气体达到最严重的爆炸气体浓度，选取甲烷体积浓度为9.5%，燃气位于燃气仓中间位置，为160 cm×300 cm×500 cm长方体。甲烷和空气混合气体在爆炸反应前后，其原始反应物和爆炸产物分别遵循理想气体的状态变化特性，因此，可以通过定义理想气体的γ律方程来描述爆炸后气体产物的状态方程。空气采用*MAT_NULL空物质材料模型定义，爆炸气体采用*MAT_HIGH_EXPLOSIVE_BURN爆炸燃烧材料模型，二者均采用线性多项式状态方程来描述爆炸气体产物的状态变化过程，线性多项式状态方程为

式中：ρ为当前密度，ρ₀为初始密度，E₀为爆炸气体内能，C₀~C₆为状态方程参数。线性多项式状态方程的相关参数见表 3。

表  3  线性多项式状态方程参数^[12]

Table  3.  Parameters for linear polynomial equation of state^[12]

Material	ρ/(mg·cm-3)	C0/MPa	C1	C2	C3	C4	C5	C6	E0/(MJ·m-3)	V0
Air	1.234	-0.1	0	0	0	0.400	0.400	0	0.250	1.0
CH4-Air	1.293	0	0	0	0	0.274	0.274	0	3.408	1.0

下载: 导出CSV 

| 显示表格

混凝土采用*MAT_PLASTIC_KINEMATIC随动强化模型，该模型是带有失效的材料模型，可以通过定义单元的失效来模拟混凝土受力破坏的过程。

土壤采用*MAT_FHWA_SOIL材料模型，这是一种各向同性损伤材料，该材料适用于固体单元且允许失效^[13]。该模型有一个修正的Mohr-Coulomb面用以决定依赖峰值的剪切强度，修正后Mohr-Coulomb屈服表面表达式为

式中：σ_v为屈服面，p为压力，φ为内摩擦角，J₂为应力偏张量第二不变量，K(q)为张量平面角的函数，A为D-P准则系数，c为黏聚力。

*MAT_FHWA_SOIL材料模型采用的修正后的Mohr-Coulomb屈服准则，该材料可以同时考虑塑性硬化、塑性软化以及应变速率效应等，相对于LS-DYNA提供的其他土壤模型，该模型计算结果更加准确。

2.   地下管廊爆炸计算结果分析

2.1   冲击波传播规律

图 2为不同时刻下燃气仓yOz面爆炸冲击波超压云图。从图 2(a)中可以看出，在混合气体爆炸开始阶段，波阵面以立方体形状向外扩散，爆炸中心压力最大，由里向外压力逐渐减小。从图 2(b)和图 2(c)可以看出，随着爆炸波向周围不断扩散，冲击波的压力峰值从爆炸中心转移到波阵面上，且在转移过程中压力峰值逐渐衰减，这是由于爆炸发生后，冲击波迅速向周围传播，一部分沿管廊轴向传播，另一部分冲击波遇到四周管壁发生反射，其中燃气仓内上管壁和下管壁的反射波向两个相对方向传播，相互叠加后导致能量耗散，致使冲击波能量逐渐降低。同理左右管壁的反射波也由于向两个相对方向传播，相互叠加后导致能量耗散。从图 2(b)、图 2(c)、图 2(d)和图 2(e)中可以看出，在冲击波波阵面向远处传播的过程中，压力峰值从中心向两侧转移，爆源附近压力小于波阵面压力，在内部产生负压。产生负压是由于爆炸气体起爆后，爆炸产物以极高的速度向周围扩散，压力峰值从中心向周围转移，形成初始冲击波。随着冲击波的向前传播，压力降至周围气体压力为扰动时的初始压力，随着冲击波向前传播，此时爆炸产物的平均压力低于介质未经扰动时的初始压力，则出现负压区。在图 2(f)中，t=100 ms时，管廊各位置处的爆炸压力相差不大，廊内压力基本相同。

图  2  管廊内冲击波超压云图

Figure  2.  Overpressure contour of shock wave in gallery at different simulation times

下载: 全尺寸图片 幻灯片

图 3给出了不同时刻下管廊xOy截面的压力云图。从图 3中可以看出，由于波阵面最初以立方体形状向外扩散，所以燃气仓内四面管壁受压状态大致相同，但随着爆炸的不断发生，管廊燃气仓左侧混凝土结构受压逐渐增大，在t=5.5 ms时，部分混凝土发生破坏。

图  3  管廊结构所受压力云图

Figure  3.  Pressure contour of gallery structure at different simulation times

下载: 全尺寸图片 幻灯片

2.2   管廊结构动力响应分析

为进一步探究混凝土受力过程，在xOy面上选取部分测点记录爆炸过程中不同位置的压力变化，测点位置如图 4所示，不同时刻下各个测点的超压时程曲线见图 5。

图  4  xOy面数据采集点位置示意图

Figure  4.  Schematic diagram of data acquisition points in the xOy plane

下载: 全尺寸图片 幻灯片

由图 5(a)可以看出，测点L₁出现强烈的正负压振荡，而测点L₁₁、L₁₂和R₁均未出现压力波强烈振荡。测点L₁在t=7.8 ms时，超压峰值达到最大，为18.65 MPa，而测点R₁压力上升速度较为缓慢，在t=10 ms左右超压值达到最大，为8.99 MPa，随后超压值逐渐趋于零。与L₁和R₁相比，L₁₁和L₁₂的压力上升速度比较慢，超压峰值较低。从图 5(b)和图 5(c)中可以看出，T₁和B₁压力变化大致相同，都在短时间内达到压力峰值，分别为1.78 MPa和1.87 MPa，随后出现较大的负压，经过几次正负压振荡后逐渐趋于零。T₁₁、T₁₂、B₁₁以及B₁₂节点的超压现象不太明显。

图  5  各节点超压时程曲线变化图

Figure  5.  Time history diagram of node overpressure

下载: 全尺寸图片 幻灯片

由图 5可以看出，与其他仓测点相比，燃气仓内各测点的受压较大，所以在爆炸荷载下燃气仓结构可能最先发生破坏。故沿燃气仓纵向选取测点记录燃气仓内爆炸过程中不同时刻的压力变化，测点位置见图 6。各测点的超压时程曲线见图 7。

图  6  数据采集点位置示意图

Figure  6.  Position diagram of data acquisition points

下载: 全尺寸图片 幻灯片

由图 7可以看出：燃气仓内墙上测点L₁的超压峰值最高，测点L₄的超压峰值最小为7.5 MPa，内墙上的各测点正负压变化比较明显，且持续时间较长；燃气仓右侧结构上的各测点与内墙测点L₁~L₄相比，超压峰值不高，且负压峰值也较低，正负压振荡持续时间较短，压力消散较快，在t=40 ms左右时，测点R₁~R₄压力值都基本趋于零。从图 7(b)、图 7(c)中可以看出，管廊燃气仓顶部和底部的测点变化规律相同，压力都在短时间内增高达到峰值，随后逐渐降低，与测点正压相比，负压峰值较大，持续时间长，大约等于3倍的超压峰值。

图  7  各节点超压时程曲线

Figure  7.  Time history diagram of node overpressure

下载: 全尺寸图片 幻灯片

由于冲击波在管廊内传播时，在波阵面内部形成负压区，导致结构受到负压影响，燃气仓内墙受到负压影响出现正负压振荡，持续时间较长，结构在受到交变应力作用下容易发生疲劳破坏。虽燃气仓顶部和底部结构测点超压值不高，但是负压值峰值较大，且持续时间长，故负压值对结构的影响不容忽视。在工程结构抗爆设计时，应当考虑此因素对结构的影响。

由图 7可知，离爆源最近的测点超压峰值最大，其余各点的峰值压力基本随偏离爆心距离的增大而滞后，超压峰值随爆心距的增加而降低。测点的超压曲线虽然各有差异，但变化规律相同，都分为3个阶段：(1)燃气被引燃，气体膨胀，压力上升到最大值；(2)压力达到峰值后，经历强烈的正负压振荡；(3)振荡幅度变小，逐渐回到初始值。

图 8为T₁、B₁、L₁和R₁测点的位移时程曲线。从图 8中可以看出, 测点位移随着时间变化不断减小。测点R₁虽然在爆炸初期产生了较大的位移，但随后位移减小。大约在t=40 ms时，测点T₁、B₁和R₁位移趋于稳定，残余变形为1.5 mm左右。在t=10 ms左右时，测点L₁位移达到最大，为10.47 mm，且残余变形为2.5 mm左右。

图  8  节点位移时程曲线变化图

Figure  8.  Time history diagram of node displacement

下载: 全尺寸图片 幻灯片

图 9给出了T₁、B₁、L₁和R₁测点的x方向和y方向的速度时程曲线。由图 9可以看出，测点的速度呈现周期性的变化，且周期内的峰值与上一次周期峰值相比降低许多。L₁、R₁测点的y方向速度较小，基本趋于零。在t=10.48 ms时，L₁测点的x方向的速度最大，达到3.303 m/s；在t=13.19 ms时，R₁测点x方向的速度最大，达到0.908 m/s。B₁和T₁测点在x方向的速度变化周期相同，速度较小。在t=15 ms时，B₁和T₁测点在y方向的速度同时达到最大，且方向相反，速度最大值分别为-0.641和-0.563 m/s。

图  9  节点速度时程曲线变化图

Figure  9.  Time history diagram of node velocity

下载: 全尺寸图片 幻灯片

3.   结论

采用流固耦合及多物质ALE算法，对地下管廊结构进行了爆炸作用下的动力仿真模拟分析。结果表明：

(1) 在爆炸过程中，超压峰值从爆源中心转移到波阵面上，导致爆源附近压力低于波阵面压力，爆源附近形成负压区，这是导致结构产生正负压振荡的主要原因；

(2) 混合气体爆炸后，冲击波遇到管道壁面形成多次发生反射，上下管壁(左右管壁)的反射波向相反方向传播，相互叠加，导致能量消散，致使冲击波强度降低；

(3) 燃气仓内发生爆炸，测点L₁受压最为明显，在t=7.8 ms时，L₁超压峰值达到最大，为18.65 MPa，同燃气仓其他壁面的测点相比，燃气仓内墙上的测点正负压振荡时间持续较长，混凝土结构容易失稳发生破坏；

(4) 在爆炸荷载作用下，管廊燃气仓内墙振动现象最为明显，且速度呈现周期性变化，在t=10 ms左右时，位移和速度达到最大值，分别为10.47 mm、3.303 m/s。