瓦斯爆炸起爆强度对巷道壁面冲击破坏的影响

Influence of Initiation Intensity of Gas Explosion on Impact Failure of Roadway Wall

摘要: 矿井瓦斯爆炸会对巷道产生严重破坏，然而在分析巷道瓦斯爆炸热冲击破坏时，不同爆炸强度时的巷道壁面破坏特征分析较少。针对现有研究的不足，本研究利用LS-DYNA建立了巷道瓦斯爆炸物理数学模型，设置的爆炸能量分别为基本模型瓦斯爆炸能量的2，4，6，8，10，15，30倍，模拟了不同爆炸强度的瓦斯爆炸热冲击对瓦斯爆炸巷道壁面破坏特性的影响，并从超压变化、速度变化、位移变化和应力变化四个方面进行分析。研究结果表明：随着爆炸强度的增加，壁面超压增大，封闭端受到的超压载荷越强，封闭端壁面损伤破坏也越严重。随着瓦斯爆炸强度的增加，壁面测点最大位移不断增大。研究成果将为煤矿瓦斯爆炸灾害防治和瓦斯爆炸事故调查提供理论依据和技术指导，对丰富和完善瓦斯爆炸理论也将起到积极的促进作用。

Abstract: Gas explosion in underground coal mine will cause serious damage to the roadway. However, when the thermal shock damage of gas explosion in roadway was analyzed, the results on damage characteristics of roadway wall under different explosion intensity are less. In view of the shortcomings of the existing research, LS-DYNA is used to establish the physical and mathematical model of gas explosion in roadway. The explosion energies are set as 2, 4, 6, 8, 10, 15 and 30 times of the gas explosion energy of the basic model. The influence of different explosion intensity on the wall failure characteristics of gas explosion in roadway is simulated. It is analyzed from four aspects: over-pressure change, velocity change, displacement change and stress change. The results show that with the increase of explosion intensity, the wall overpressure increases, and the stronger the overpressure load on the closed end, the more serious the wall damage of the closed end. With the increase of gas explosion intensity, the maximum displacement of wall measurement points increases continuously. The research results can provide theoretical basis and technical guidance for the prevention and control of gas explosion disasters and the investigation of gas explosion accidents in coal mine, and will also play a positive role in enriching and perfecting the theory of gas explosion.

1. 引言

煤矿瓦斯爆炸事故导致众多人员伤亡、巷道设施设备严重破坏和巨大财产损失，造成不良的社会影响 [1]。为了掌握瓦斯爆炸破坏特性以及有效预防瓦斯爆炸事故，国内外许多的学者利用管道对瓦斯爆炸进行了大量研究 [2] - [7]，并取得大量的研究成果。高娜、张延松等从理论、实验和数值模拟三方面介绍了近年来国内外在瓦斯爆炸特性方面的研究现状和研究成果 [8]。聂百胜等 [9] [10] 从内径尺寸、截面形状、开口情况等方面设计实验管道，研究火焰传播过程中火焰形状、传播速度与超压，以及瓦斯爆炸化学动力特性，获得了煤矿瓦斯爆炸过程中冲击波的衰减特性。由于大型爆炸实验的成本较高、难度较大，进行的实验数量有限 [11] [12]，实验结果很难完整地体现出瓦斯爆炸强度对爆炸破坏特性的影响。虽然人们非常关注巷道瓦斯爆炸冲击破坏问题，并进行了大量的研究工作 [13] - [18]，但是不同爆炸强度的巷道壁面破坏特征分析较少，然而LS-DYNA软件能很好解决热冲击对巷道壁面破坏问题。针对现有研究的不足，本研究利用LS-DYNA数值模拟软件，建立巷道瓦斯爆炸的物理数学模型，设置了爆炸能量分别为基本模型瓦斯爆炸能量的2，4，6，8，10，15，30倍，通过测定瓦斯爆炸的超压变化、速度变化、位移变化和应力变化等分析瓦斯爆炸热冲击时的不同爆炸强度对瓦斯爆炸巷道壁面破坏特性的影响。

2. 研究方法

LS-DYNA是通用显式计算软件，兼有隐式求解功能，具有分析功能强大、单元类型众多、材料模型丰富、接触类型齐全、荷载及约束条件全面等特点。在冲击和爆炸、物体碰撞、流体分析、裂纹扩展分析、流固耦合、热固耦合等领域得到了广泛的应用。LS-DYNA兼有拉格朗日(Lagrange)算法、欧拉(Euler)算法和ALE (A Arbitrary_Lagrange_Euler)算法。运用ANSYS/LS-DYNA模拟时，可先利用ANSYS建立模型，然后用ANSYS/LS-DYNA做显式求解，也可在ANSYS和ANSYS/LS-DYNA之间传递几何信息和结果数据，进行连续的隐式–显式分析或显式–隐式分析。气体爆炸数值模拟方法通常有TNT当量法和气体填充法 [19] [20]，本次选择TNT当量法。

2.1. 物理模型

2.1.1. 有限元模型的建立

根据巷道几何模型的中心对称性，将建立四分之一的巷道内瓦斯爆炸模型进行相应的数值模拟研究，既便于观察巷道内部壁面情况又减少了计算量。本次数值模拟所采用的巷道几何参数为：巷道外半径R = 1 m，巷道内半径R = 0.8 m，巷道共长10.2 m，巷道壁面厚0.2 m。物理模型如图1所示。巷道一端封闭，另一端为开口，掘进巷道内密闭端存放体积为10.048 m³，浓度为9.5%瓦斯预混气体，巷道内空气与瓦斯预混气体之间用薄膜隔开，空气填充长度为5 m。本文中所用材料的模型和参数均采用统一单位制(kg/m/s)。巷道模型如图2所示。

Figure 1. Physical model of roadway

图1. 巷道物理模型

Figure 2. Model of roadway

图2. 巷道模型

2.1.2. 网格划分

本次对巷道瓦斯空气采用六面体单元网格划分；瓦斯气体和空气采用Euler网格模拟，巷道壁面采用Lagrange网格。单元长度控制为0.05 m，划分网格后的有限元模型如图3所示，瓦斯爆炸巷道物理模型共划分为171,700个单元。

Figure 3. The finite element model after meshing

图3. 划分网格后的有限元模型

2.1.3. 定义单元类型与材料模型

1) 空气本构模型 [18] [19]

在数值模拟中，通常将空气设置为理想气体，标准状态下空气的初始状态参数：P = 0.1 MPa，ρ = 1.29 kg/m³，T = 298 K，由于冲击波波阵面很薄，所以传播过程中热交换以及与巷道壁面之间的摩擦可忽略不计；瓦斯爆炸可按理想气体膨胀来研究爆炸后冲击波；瓦斯爆炸冲击波的膨胀假设是绝热过程，根据Gama准则，一般采用线性多项式的状态方程来描述其性能。

2) 空气材料模型及状态方程 [18] [19]

空气采用NULL材料模型以及*EOS_LINEAR_POLYNOMIAL 状态方程加以描述。线性多项式状态方程为：

$\begin{array}{l}P=C_0+C_1\mu +C_2{\mu }^2+C_3{\mu }^3+\left(C_4+C_5\mu +C_6{\mu }^2\right)E\\ \mu =\frac{1}{V}-1\end{array}$ (1)

式中，P为爆轰压力；E为单位体积内能；V为相对体积； $C_0~C_6$ 为状态方程参数，为常数。

3) 瓦斯气体材料模型及状态方程 [18] [19]

瓦斯气体采用HIGH_EXPLOSIVE_BURN材料模型以及式(1)所述*EOS_LINEAR_ POLYNOMIAL状态方程加以描述，瓦斯和空气参数如表1所示。

4) 巷道壁面材料参数 [18] [19]

本文巷道模型的壁面采用刚体材料模型，此模型在LS-DYNA中的关键字为*MAT_RIGID。选用的材料如表2所示。

Table 1. Air and gas parameters

表1. 瓦斯和空气参数

Table 2. Material parameters of roadway wall

表2. 巷道壁面材料参数

5) 沙漏控制

如果使用全积分算法会耗费大量CPU时间，通过在模型中采用单点积分可以有效的节约CPU时间，但是很易造成沙漏的零能模式，此时就需对沙漏进行必要的控制。对于实体单元SOLID164沙漏控制，K文件设置如表3所示 [19] [20]。

Table 3. Hourglass control

表3. 沙漏控制

6) 点火位置

在K文件中，模型点火位置设置如表4所示 [19] [20]。

Table 4. Ignition position setting

表4. 点火位置设置

7) 时间控制

求解终止时间设置为0.05 s，时间步长控制为0.67，如表5所示。

Table 5. Time control

表5. 时间控制

2.2. 数学模型

2.2.1. 基本假设

为了方便数值模拟和建立数学模型，本次模拟进行如下假设：① 忽略其化学反应中间产物，但需要考虑气体粘性。② 在标准大气压条件下。③ 瓦斯爆炸过程为单向不可逆过程，爆炸过程为绝热过程。④ 壁面为绝热刚性壁面，不考虑巷道内部空间与外界的热量传递辐射。⑤ 边界设置为无反射边界条件。⑥ 瓦斯混合气体均匀混合，处于常温且比热容遵循混合规则。⑦ 巷道内部无障碍。

2.2.2. 基本控制方程

ANSYS/LS-DYNA软件主要采用了Lagrangian描述增量法。本次取初始时间的质点位置 $X_i\left(i=1,2,3\right)$，在任意t时刻，该质点位置为 $x_i\left(i=1,2,3\right)$。质点的运动方程为 [19] [20]：

$x_i=x_i\left(X_j,t\right)\left(i,j=1,2,3\right)$ (2)

在t = 0时，初始条件为：

$x_i\left(X_j,0\right)=x_i$ (3)

${\stackrel{\dot{}}{x}}_i\left(X_j,0\right)=V_i\left(X_j,0\right)$ (4)

式中： $V_i$ 为初始速度。

1) 动量守恒方程

${\sigma }_{ij,j}+\rho f_i=\rho {\ddot{x}}_i$ (5)

式中： ${\sigma }_{ij}$ 为Cauchy应力； $f_i$ 为单位质量体积力； ${\ddot{x}}_i$ 加速度。

2) 质量守恒方程

$\rho V={\rho }_0$ (6)

式中： $\rho$ 为当前质量密度； ${\rho }_0$ 为初始质量密度； $V=\left|F_{ij}\right|$ 为相对体积； $F_{ij}=\frac{\partial x_i}{\partial x_j}$ 为变形梯度。

3) 能量守恒方程

$\stackrel{\dot{}}{E}=V{S}_{ij}{\stackrel{\dot{}}{\epsilon }}_{ij}-\left(p+q\right)\stackrel{\dot{}}{V}$ (7)

$S_{ij}={\sigma }_{ij}+\left(p+q\right){\delta }_{ij}$ (8)

$p=-\frac{1}{3}{\sigma }_{ij}{\delta }_{ij}-q=-\frac{1}{3}{\sigma }_{kk}-q$ (9)

式中， $\epsilon$ 为应变率；q表示体积粘性阻力；S为偏应力；p为压力； $\delta$ 克罗内克符号。

4) 边界条件

① 面力边界条件：

${\sigma }_{ij}{n}_j=t_i\left(t\right)$ (10)

式中： $n_j\left(j=1,2,3\right)$ 表示边界的外法线方向余弦； $t_i\left(i=1,2,3\right)$ 为面力荷载。

② 位移边界条件：

$x_i\left(X_j,t\right)=K_i\left(t\right)$ (11)

式中： $K_i\left(t\right)$ 为给定位移函数 $\left(i=1,2,3\right)$。

③ 滑动接触面几何间断处的位移条件：

$\left({\sigma }_{ij}^{+}-{\sigma }_{ij}^{-}\right)n_i=0$ (12)

在求解循环的过程中，新的时步长为全部单元时步长的最小值，即

$\Delta t^{n+1}=\alpha \min \left(\Delta t_{e1},\Delta t_{e2},\cdots ,\Delta t_{eN}\right)$ (13)

式中：N为单元数目。

单元的极限时步长可以由式(14)计算可得：

$\Delta t_e=\frac{L_2}{Q+\sqrt{Q^2+c^2}}$ (14)

3. 数值模拟及结果分析

壁面测点选取如图4所示，为便于区分，在将单元测点编号命名为A、B、C号测点，其中A测点为41605号单元，坐标位置为a (0, 0, 0.2)；B测点为25729号单元，坐标位置为b (0, 1, 2.5)，C测点为36097号单元，坐标位置为c (0, 1, 7.5)。

Figure 4. Measurement point settings

图4. 测点设置

3.1. 超压变化情况

图5展示了巷道壁面A、B、C三个测点在不同瓦斯爆炸强度下的超压变化规律。由图5中A测点在不同爆炸强度下超压时程曲线可以发现，在不同爆炸强度下，A测点随时间的变化规律基本一致，在达到最大峰值后均迅速衰减。随着爆炸强度的增加，A测点的最大值在不断增加。当瓦斯爆炸强度为2倍时，A测点超压最大值为1.38E07 Pa，当瓦斯爆炸强度为30倍时，A测点超压最大值为2.01E08 Pa，由于A测点为轴向封闭端的中心，瓦斯爆炸主要沿轴向传播，因此A测点超压测值较大。由此可得巷道内瓦斯爆炸强度越大，封闭端受到的超压载荷越强，封闭端壁面损伤破坏也越严重。随着爆炸强度的增加，B测点的最大值也在不断增加。由于巷道壁面约束，冲击超压在壁面来回反射，瓦斯爆炸强度越大，B测点壁面反射超压也越大。由图5中C测点在不同爆炸强度下超压时程曲线可以得到，其各规律和B测点一致，由于C测点处于空气区，B测点处于瓦斯区，所以C测点超压测值在相同瓦斯爆炸强度上总体上要偏小。例如，在10倍瓦斯爆炸强度时，B测点最大值为4.95E07 Pa，C测点最大值为2.02E07 Pa，B测点是C测点的2倍多。因此B测点周围壁面损伤程度要比C测点严重些。

Figure 5. Time-history curve of wall point overpressure under different explosion intensities

图5. 不同爆炸强度下壁面测点超压时程曲线

3.2. 速度变化情况

图6展示了巷道壁面A、B、C三个测点在不同瓦斯爆炸强度下的速度变化规律。由6图中A、B、C测点在不同爆炸强度下的速度时程曲线可以发现，在不同爆炸强度下A、B、C测点速度的变化规律基本一致，即随着爆炸强度的增加，各测点的速度峰值在不断增加。在A测点，封闭端壁面受到爆炸冲击波后，速度开始增加，并迅速达到峰值，当强度为2倍时，A测点在0.0004 s速度达到0.547 m/s，当强度为30倍时，同时刻A测点速度达到8.184 m/s，随着爆炸反应速度测值迅速衰减，而在B、C测点，速度并没有迅速衰减，而是产生来回振动，其中B测点振动频率高于C测点。总体上A、B、C各测点在相同的爆炸强度下速度峰值依次减小，对应了各区域的壁面损伤破坏程度。

Figure 6. Time-history curve of wall point speed under different explosion intensities

图6. 不同爆炸强度下壁面测点速度时程曲线

3.3. 位移变化情况

图7展示了巷道壁面A、B、C三个测点在不同瓦斯爆炸强度下的位移变化规律。由图7中A测点在Z方向不同爆炸强度下的位移时程曲线可以发现，在不同瓦斯爆炸强度下，封闭端壁面对应着不同的位移量，在0.05 s时位移最终稳定。瓦斯爆炸强度从2倍到30倍，位移量分别为0.00173 m，0.00243 m，0.00297 m，0.00343 m，0.00383 m，0.00468 m，0.00661 m。由图7中B测点和C测点在Y方向不同爆炸强度下的位移时程曲线可以发现，随着瓦斯爆炸强度的增加，壁面测点最大位移不断增大。当瓦斯爆炸强度为30倍时，A、B、C各测点最大位移分别为0.00669 m，0.000624 m，0.000548 m。位移量致使巷道壁面扩张变形直接关系到巷道壁面的损伤破坏情况，总的来说，有效降低瓦斯爆炸强度对巷道壁面损伤破坏至关重要。

Figure 7. Time-history curve of wall point displacement under different explosion intensities

图7. 不同爆炸强度下壁面测点位移时程曲线

3.4. 应力变化情况

图8展示了巷道壁面A、B、C三个测点在不同瓦斯爆炸强度下的应力变化规律。由图8中A、B、C测点在不同爆炸强度下应力时程曲线可以发现，随着爆炸强度的增加，A、B、C测点的最大值在不断增加。以A测点为例，随着爆炸强度的增加，应力最大值测值分别为1.83E07 Pa，3.67E07 Pa，5.50E07 Pa，7.32E07 Pa，9.14E07 Pa，1.36E08 Pa，2.72E08 Pa。瓦斯爆炸冲击载荷加载到壁面，壁面形变产生应力，当应力较大壁面形变不可恢复时达到壁面承受应力极限，若应力不断增加对壁面就会造成不可修复破坏。在B、C测点壁面约束导致应力测值往复循环，应力的来回加载容易导致壁面疲劳损伤，若加载强度不衰减最终壁面会被破坏。在相同爆炸强度下通过比较各测点壁面应力同样满足A > B > C。壁面在爆炸载荷作用下的超压、速度、位移和应力是衡量壁面损伤破坏程度的重要因素，通过模拟得到随着爆炸强度的增加，壁面超压、速度、位移和应力也不断上升。

Figure 8. Time-history curve of wall point Stress under different explosion intensities

图8. 不同爆炸强度下壁面测点应力时程曲线

4. 结论

本文运用ANSYS/LSDYNA对巷道内不同瓦斯爆炸强度时的瓦斯爆炸对巷道壁面破坏特性进行了数值模拟，得出如下结论。

1) 巷道壁面A、B、C三个测点在不同瓦斯爆炸强度下的超压变化规律显示，随着爆炸强度的增加，壁面超压增大，封闭端受到的超压载荷越强，封闭端壁面损伤破坏也越严重。

2) 随着爆炸强度的增加，各测点的速度峰值在不断增加。总体上A、B、C各测点在相同的爆炸强度下的速度峰值依次减小，对应了各区域的壁面损伤破坏程度。

3) 随着瓦斯爆炸强度的增加，壁面测点最大位移不断增大。位移量致使巷道壁面扩张变形直接关系到巷道壁面的损伤破坏情况，总的来说有效降低瓦斯爆炸强度对巷道壁面损伤破坏至关重要。

4) 在不同瓦斯爆炸强度下，封闭端壁面响应随着爆炸强度增加愈加强烈。在轴向壁面上，瓦斯区测点的超压、速度、位移和应力测值比空气区大，瓦斯区壁面云图响应明显，瓦斯区壁面损伤比空气区壁面损伤严重。

基金项目

国家自然科学基金项目“矿井瓦斯爆炸能量释放转化特性及热冲击能量损耗研究”(编号：52174177)、“矿井多爆源瓦斯爆炸传播特性及热冲击动力学机制研究”(编号：52174178)。湖南省教育厅青年项目“煤矿瓦斯爆炸巷道壁面破坏过程热冲击能量耗损特征”(编号：20B204)。贵州省教育厅青年成长项目“黔西南州山地旅游安全评价与综合治理模式建构——以贵州醇风景区为例”(编号：黔教合KY字[2019]22)。