# Analysis of the Characteristics and Influencing Factors of Gas Explosion in Heading Face

Authors: Xue-bo Zhang; Jian-liang Gao; Jing-zhang Ren; Chun-xia Wang
Journal: Shock and Vibration
 (2020)

---

## Abstract

In order to accurately grasp the characteristics and influencing factors of gas explosion in heading face, the mathematical model of gas explosion was determined. According to the actual size of a heading face of a coal mine, a 3D geometric model with a length of 100 m was established, and the effects of ignition energy and gas explosion equivalent on the gas explosion characteristics of the heading face were analyzed. The results show the following. (1) The mathematical models for numerical simulation of gas explosion can accurately simulate the gas explosion and its propagation process. The time-space step size has a great influence on the simulation results. The grid spacing for numerical simulation of mine gas explosion is determined to be 0.1 m and the time step length is determined to be 0.001 s. (2) The ignition energy has a limited effect on gas explosion characteristics. It only has a certain influence on the gas explosion process, but has little influence on the overpressure of shock wave. The larger the ignition energy is, the faster the explosion reaction speed is, and the maximum overpressure increases slightly. When the ignition energy increases to a certain value, the time of peak shock wave and the maximum overpressure both tend to be stable. The ignition energy has little effect on gas explosion characteristics when an explosion accident occurs underground with a large amount of gas accumulation. (3) The gas explosion equivalent has a great influence on the overpressure of gas explosion shock wave. The higher the explosion equivalent is, the greater the pressure is, and the peak value of the shock wave overpressure increases with the explosion equivalent as a power function. The research results have important guiding significance for the research and development of new technology for prevention and control of gas explosion.


---

## Body

## 1. Introduction

Coal mine gas explosion accident is one of the most serious disaster accidents [1], especially the major and large gas explosion accident, resulting in a large number of casualties and serious equipment damage. For example, on September 27, 2000, a gas explosion occurred in the Muchonggou Coal Mine of Guizhou Shuicheng Mining Bureau, resulting in 162 deaths and more than 12 million yuan of direct economic losses. In 2016, huge gas explosion accidents occurred in the Jinshangou Coal Mine in Chongqing and Bauma Coal Mine in Inner Mongolia, which killed 33 and 32 people, respectively [2, 3]. In the same year, 36 people died in Severnaya coal mine in Russia [4]. Mine gas explosion can be regarded as a gas explosion process in which methane–air mixed gas is ignited by external fire source, which can be divided into two stages: ignition and propagation. The destructiveness of explosion is mainly reflected in the propagation stage [5], that is, the combustion of combustible and the mechanical damage in the propagation of explosion shock wave caused by the propagation of explosion flame [6]. The destructiveness of explosion is mainly reflected in the propagation stage [5], that is, the combustion of combustible materials caused by the propagation of explosive flame and the mechanical damage caused by the propagation of explosive shock wave [6]. Through studying the propagation law of gas explosion in underground roadway and judging and predicting the propagation range of shock wave and the losses caused by it, it can provide theoretical guidance for the research and development of new technologies for gas explosion prevention and control [7, 8].Many scholars at home and abroad have carried out a lot of researches on the characteristics of gas explosion and the propagation law of shock wave. The characteristics of gas explosion and the propagation of shock wave are affected by the volume, composition, and concentration of premixed gas, the shape, energy, temperature, and space size of ignition source. It is very difficult to study the theory of gas explosion, and experimental study and numerical simulation are the main means to study gas explosion. In terms of experimental research, in view of the strong destructiveness of gas explosion and the risk of test, most scholars study the propagation and destruction characteristics of gas explosion in small-size pipelines. Zhao et al. [9], Jiang et al. [10, 11], and Zhu et al. [12] experimentally studied the propagation law and influencing factors of gas explosion flame and shock wave in a square straight pipe with centimetre-level diameter. Yu et al. [13] studied the influence of pressure relief port on gas explosion shock wave and flame propagation in a square straight pipe with dimension of 10 cm × 10 cm × 1 m. However, some scholars carried out gas explosion tests in real scale of mine. For example, Chongqing Coal Research Institute of China carried out a gas explosion test with methane–air premixed volume of 100–200 m3 in a square single-head straight lane with a section of 7.2 m2 and a length of 900 m [14]. The National Institute for Occupational Safety and Health (NIOSH) conducted an experiment on the damage characteristics of closed walls caused by methane gas explosion with a maximum methane–air premixed volume of 210 m3 in a straight mine lane with a section of about 12.5 m2 and a length of about 480 m [15]. In terms of numerical simulation, Scott et al. [16] simulated the gas explosion accident in the Upper Big Branch (UBB) coal mine in the United States. Jiang et al. [17, 18] simulated and studied the propagation characteristics of gas explosion shock wave in parallel roadway as well as the propagation law of gas explosion shock wave to adjacent coal mining face in heading face. Zhu et al. [19] used numerical simulation to study the impact and oscillation characteristics of the explosion wave in the closed system and the variation law of its characteristic parameters. Cheng et al. [20] studied the influence of pipe bifurcation and turning on gas explosion propagation characteristics. Qiu [21] used the k−ε turbulence model and laminar velocity/eddy-dissipation combustion model to simulate the shock wave propagation laws under pipe turning, roadway bifurcation, and abrupt change of cross section of pipe. Meng et al. [22] used RNG k−ε turbulence model and EDM combustion model to simulate and analyze the influence of obstacles on the propagation characteristics of gas explosion shock wave in large-size ventilation pipe. Cheng et al. [23–26] used numerical simulation to study the modelling mine gas explosive pattern in underground mine gob and overlying strata and effects of explosion impact load on underground mine seal. Zhu et al. [27] studied the bending angles of the three angles. The result showed that the smaller the bending angle, the stronger the reflection, the stronger the turbulence, and the greater the peak overpressure. Ma et al. [28] studied the propagation characteristics of methane explosion in the pipeline network by numerical simulation method. The result showed that, in the parallel pipeline network, the peak of overpressure increased significantly at the intersection point, while the flame velocity decreased at the intersection point.For the current research results, no matter the experimental research or numerical simulation, most of them study the gas explosion characteristics and shock wave propagation law in laboratory small-size pipelines or local roadway, so they cannot objectively and quantitatively reflect the propagation law of gas explosion shock wave in real scale of mine. Considering the risk and cost of gas explosion test research at the existing mine scale, methane–air premix volume is small, and roadway configuration and sensor layout are relatively simple, which are not enough to fully reflect the propagation law of gas explosion shock wave at the real mine scale. The numerical simulation method provides the possibility to study the propagation law of gas explosion shock wave in mine scale, but the following three problems need to be studied. (1) Many scholars have used different theoretical models to numerically simulate the characteristics and propagation law of gas explosion. Further studies are needed to determine which theoretical model is appropriate and how much influence the time-space step size has on the simulation results. (2) There are many ignition sources that cause gas explosion in underground coal mines, such as electric spark, friction and collision spark, open flame and thermal spontaneous combustion, etc. Their ignition size and ignition energy are different, and their influence characteristics on gas explosion need to be studied. (3) Under the influence of coal seam gas occurrence and system ventilation, the gas accumulation amount in roadway is different, and the influence of gas accumulation amount (gas explosion equivalent) on the gas explosion characteristics and shock wave propagation law in heading face has not been studied. It has important guiding significance for accurately predicting the propagation law of gas explosion shock wave and effectively controlling the destructive effect of gas explosion in mine.


## 2. Mathematical Model of Mine Gas Explosion and Key Parameter Determination

### 2.1. Determination of Mathematical Model of Mine Gas Explosion

#### 2.1.1. Establishment of Mathematical Model of Gas Explosion

Gas explosion can be assumed to be an ideal gas heated accelerated expansion process, which can be described by a mathematical model composed of kinetic process control equation (mass conservation equation, momentum conservation equation, energy conservation equation, and component balance equation), turbulence model, and gas chemical reaction model (combustion model and turbulence flame velocity model).The turbulence model includesk−ε model, LES model, and DES model. Gas chemical reaction models include laminar flow finite rate model, eddy-dissipation model, laminar flow finite rate/eddy-dissipation model, and EDC model. The laminar finite rate model uses the Arrhenius formula to calculate the chemical source terms but ignores the influence of turbulence fluctuation, so it is suitable for simulating laminar flame propagation. Eddy-dissipation model (EDM) can simulate the rapid combustion of most fuels, and the overall reaction rate is controlled by turbulent mixing. When LES turbulence model is used (LES + eddy-dissipation model), turbulent mixing rate is replaced by subgrid-scale mixing rate, which can better simulate the propagation process of gas explosion shock wave. Eddy-dissipation concept (EDC) model is an extension of eddy-dissipation model, which can simulate the turbulent reaction flows with chemical reaction mechanisms. However, typical mechanisms have different rigidity, their numerical integration calculation is very expensive, and the model can only be used when fast chemical reaction is assumed to be invalid. Therefore, LES turbulence model and eddy-dissipation model (EDM) are finally adopted to simulate the propagation law of gas explosion shock wave.


### 2.2. Combustion Model and Turbulent Flame Velocity Model

(1)Rfu=−minRi,fu,Rj,fu,(2)Ri,fu=vi,fu′Mω,iAρτsgs−1minYi,fuvi,fuMω,fu,(3)Rj,fu=vj,fu′Mω,jABρτsgs−1∑jYj,fu∑jNvj,fu″Mω,j,(4)τsgs−1=2SijSij,where Yi,fu is the mass fraction of any production species; Yj,fu is the mass fraction of a particular reactant; A, B is empirical constant equal to 4.0 and 0.5, respectively; τsgs−1 is subgrid-scale mixing rate; and Sij is strain rate tensor.The wall function method is adopted to deal with the area near the wall, and the standard wall function given by Lauder and Splading is adopted to deal with the fluid flow near the wall of a confined space [29].


#### 2.2.1. Determination of the Mathematical Model of Gas Explosion

As shown in Figure1, the geometric model was established according to the pipe size of the gas explosion experiment in Dr. Jia Zhiwei’s thesis. The pipe sections were all 0.08 m × 0.08 m, the straight pipe length was 19.2 m, and the vertical pipe length was 5 m. Two measuring points were set up in the experiment. The measuring point 1, being 19 m away from the leftmost end, was located on the central line of the horizontal pipeline; the measuring point 2, being 1 m away from vertical pipeline, was located on the central line of the vertical pipeline.Figure 1
Geometric model diagram.In order to verify the correctness of the selected mathematical models, four commonly used models, namely, LES turbulence model + eddy-dissipation model (EDM),k−ε turbulence model + laminar flow/eddy-dissipation model, RNG k−ε + eddy-dissipation model (EDM), and k−ε turbulence model + eddy-dissipation conceptual (EDC) model, were used for numerical simulation. In the simulation of the four models, the grid spacing of the geometric model was set to 5 mm, the time step length was set to 0.001 s, and the maximum number of iterations was 30 steps. When the length of gas filling area was 4.0 m, 5.5 m, and 7.0 m, respectively, the overpressure results measured at measuring points 1 and 2 by numerical simulation of four mathematical models were obtained, which were then compared with the test results of Jia Zhiwei’s gas explosion experiment, as shown in Figure 2.Figure 2
Comparison between simulation results and experimental results. (a) Comparison of results at measuring point 1. (b) Comparison of results at measuring point 2.
(a)(b)It can be seen from Figure2 that the numerical simulation results of LES turbulence model and v eddy-dissipation model (EDM) were closer to the experimental results, and the relative errors were less than 10%. Therefore, the established mathematical model of gas explosion can be used to simulate the propagation law of gas explosion shock wave.


### 2.3. Determination of Key Parameters of Numerical Simulation of Mine Gas Explosion

According to the actual size of fully mechanized excavation face of a coal mine return air channeling, a 3D geometric model with a length of 100 m was established, as shown in Figure3. Assume that the area 10 m from the heading face is filled with gas with a concentration of 9.5%, and the ignition sources are arranged in the center of the heading face. The open boundary is set on the right side of the model, the closed boundary is set on the left side, and the roadway wall is on all sides.Figure 3
Geometric model diagram of driving surface. (a) Schematic diagram of geometric model of driving surface. (b) 3D model diagram of driving surface.


#### 2.3.1. Determination of Grid Size

The grid density has great influence on the numerical results, so it is very important to determine the appropriate grid size. The error sources of numerical calculation mainly include approximate error of physical model, truncation error of differential equation, discrete error of solving region, iterative error, and rounding error. Generally speaking, the truncation error and the discrete error decrease as the grid becomes thinner, but as the grid becomes thinner, the number of discrete points increases and the rounding error also increases. As grids become denser, computing scale and storage space will increase and computing efficiency will decrease. Under the conditions of balancing computation efficiency, storage space, and accuracy, coarse grid can be used to solve the problem first, and then the grid can be continuously refined until the increase of the grid number has little influence on the calculation result. Typical grid-refinement study results are shown in Figure4.Figure 4
Relationship between solution error and grid size.The curve graph of solution error and grid size can be divided into three areas in the condition of ensuring the convergence of the solution results. (1) When the grid size is large, the calculation error is mainly truncation error. It is the error caused by the Taylor series expansion of the differential operator, which is mainly affected by the order of the element. (2) As the mesh size gradually decreases, it gradually enters the asymptotic range, and the error in this area is mainly discretization error. Generally speaking, the more dense the mesh is, the smaller the calculation error is, and the closer the calculation result is to the real solution; that is, the mesh size is located in the above two regions. (3) When the grid size is very small, the calculation error is mainly dominated by floating-point error. At this time, the smaller the grid size, the larger the error.In order to determine the grid size, the gas explosion process in heading face was numerically simulated when the grid spacing was set to 0.05 m, 0.1 m, 0.2 m, and 0.4 m (the number of grids was 14800000, 1850000, 237500, and 26000, resp.), the step length of iteration time was 0.001 s, and the maximum number of iterations was 30 steps. The simulated maximum overpressure obtained at the measuring point (80 m away from the heading face) is shown in Table1.Table 1
Simulation results of different grid sizes.
NumberGrid size (m)Grid numberMaximum overpressure time (s)Maximum overpressure (kPa)Note10.05148000000.138355.920.118500000.145351.630.22375000.171337.640.4260000.703315.2As can be seen from Table1, the grid size has a great influence on the simulation results. With the continuous refinement of the grid, the maximum overpressure at the measuring point increased gradually (the increase rate slowed down gradually). When the grid spacing was refined from 0.4 m, 0.2 m, and 0.1 m to 0.05 m, the maximum overpressure increased by 6.6%, 2.7%, and 0.6%, respectively. The duration time of maximum overpressure of shock wave was continuously shortened (the shortening rate was slowing down gradually). When the grid spacing was refined from 0.4 m, 0.2 m, and 0.1 m to 0.05 m, the duration time of maximum overpressure was shortened by 80.4%, 19.3%, and 4.8%, respectively. According to the simulation results, the results were very close when the grid spacing was 0.1 m and 0.05 m. Considering the calculation efficiency, storage space, and accuracy, the grid spacing was finally determined to be 0.1 m.


#### 2.3.2. Time Step Determination

The proper selection of time step can not only ensure the convergence of simulation results, but also ensure the accuracy of numerical simulation results and improve the calculation efficiency. In order to determine the appropriate time step, the gas explosion characteristics in heading surface were simulated when the time step was set to 0.01 s, 0.001 s, and 0.0005 s, respectively, the grid spacing was 0.1 m, and the maximum number of iterations was 30 steps. The overpressure attenuation law obtained at the measured points arranged at 40 m, 50 m, 60 m, 70 m, 80 m, and 90 m from the heading face is shown in Figure5.Figure 5
The maximum overpressure and its occurrence time at different positions.As can be seen from Figure5, the time step has a great influence on the simulation results. The smaller the time step is, the greater the simulated overpressure is and the faster the peak overpressure appears. The simulation result with time step size of 0.01 s differed greatly from that with time step size of 0.001 s. The simulation result with time step size of 0.001 s was close to that with time step size of 0.0005 s, and the difference between them was less than 6%. Taking into account the calculation efficiency and accuracy, the time step was finally determined to be 0.001 s.


## 2.1. Determination of Mathematical Model of Mine Gas Explosion

### 2.1.1. Establishment of Mathematical Model of Gas Explosion

Gas explosion can be assumed to be an ideal gas heated accelerated expansion process, which can be described by a mathematical model composed of kinetic process control equation (mass conservation equation, momentum conservation equation, energy conservation equation, and component balance equation), turbulence model, and gas chemical reaction model (combustion model and turbulence flame velocity model).The turbulence model includesk−ε model, LES model, and DES model. Gas chemical reaction models include laminar flow finite rate model, eddy-dissipation model, laminar flow finite rate/eddy-dissipation model, and EDC model. The laminar finite rate model uses the Arrhenius formula to calculate the chemical source terms but ignores the influence of turbulence fluctuation, so it is suitable for simulating laminar flame propagation. Eddy-dissipation model (EDM) can simulate the rapid combustion of most fuels, and the overall reaction rate is controlled by turbulent mixing. When LES turbulence model is used (LES + eddy-dissipation model), turbulent mixing rate is replaced by subgrid-scale mixing rate, which can better simulate the propagation process of gas explosion shock wave. Eddy-dissipation concept (EDC) model is an extension of eddy-dissipation model, which can simulate the turbulent reaction flows with chemical reaction mechanisms. However, typical mechanisms have different rigidity, their numerical integration calculation is very expensive, and the model can only be used when fast chemical reaction is assumed to be invalid. Therefore, LES turbulence model and eddy-dissipation model (EDM) are finally adopted to simulate the propagation law of gas explosion shock wave.


## 2.1.1. Establishment of Mathematical Model of Gas Explosion

Gas explosion can be assumed to be an ideal gas heated accelerated expansion process, which can be described by a mathematical model composed of kinetic process control equation (mass conservation equation, momentum conservation equation, energy conservation equation, and component balance equation), turbulence model, and gas chemical reaction model (combustion model and turbulence flame velocity model).The turbulence model includesk−ε model, LES model, and DES model. Gas chemical reaction models include laminar flow finite rate model, eddy-dissipation model, laminar flow finite rate/eddy-dissipation model, and EDC model. The laminar finite rate model uses the Arrhenius formula to calculate the chemical source terms but ignores the influence of turbulence fluctuation, so it is suitable for simulating laminar flame propagation. Eddy-dissipation model (EDM) can simulate the rapid combustion of most fuels, and the overall reaction rate is controlled by turbulent mixing. When LES turbulence model is used (LES + eddy-dissipation model), turbulent mixing rate is replaced by subgrid-scale mixing rate, which can better simulate the propagation process of gas explosion shock wave. Eddy-dissipation concept (EDC) model is an extension of eddy-dissipation model, which can simulate the turbulent reaction flows with chemical reaction mechanisms. However, typical mechanisms have different rigidity, their numerical integration calculation is very expensive, and the model can only be used when fast chemical reaction is assumed to be invalid. Therefore, LES turbulence model and eddy-dissipation model (EDM) are finally adopted to simulate the propagation law of gas explosion shock wave.


## 2.2. Combustion Model and Turbulent Flame Velocity Model

(1)Rfu=−minRi,fu,Rj,fu,(2)Ri,fu=vi,fu′Mω,iAρτsgs−1minYi,fuvi,fuMω,fu,(3)Rj,fu=vj,fu′Mω,jABρτsgs−1∑jYj,fu∑jNvj,fu″Mω,j,(4)τsgs−1=2SijSij,where Yi,fu is the mass fraction of any production species; Yj,fu is the mass fraction of a particular reactant; A, B is empirical constant equal to 4.0 and 0.5, respectively; τsgs−1 is subgrid-scale mixing rate; and Sij is strain rate tensor.The wall function method is adopted to deal with the area near the wall, and the standard wall function given by Lauder and Splading is adopted to deal with the fluid flow near the wall of a confined space [29].


### 2.2.1. Determination of the Mathematical Model of Gas Explosion

As shown in Figure1, the geometric model was established according to the pipe size of the gas explosion experiment in Dr. Jia Zhiwei’s thesis. The pipe sections were all 0.08 m × 0.08 m, the straight pipe length was 19.2 m, and the vertical pipe length was 5 m. Two measuring points were set up in the experiment. The measuring point 1, being 19 m away from the leftmost end, was located on the central line of the horizontal pipeline; the measuring point 2, being 1 m away from vertical pipeline, was located on the central line of the vertical pipeline.Figure 1
Geometric model diagram.In order to verify the correctness of the selected mathematical models, four commonly used models, namely, LES turbulence model + eddy-dissipation model (EDM),k−ε turbulence model + laminar flow/eddy-dissipation model, RNG k−ε + eddy-dissipation model (EDM), and k−ε turbulence model + eddy-dissipation conceptual (EDC) model, were used for numerical simulation. In the simulation of the four models, the grid spacing of the geometric model was set to 5 mm, the time step length was set to 0.001 s, and the maximum number of iterations was 30 steps. When the length of gas filling area was 4.0 m, 5.5 m, and 7.0 m, respectively, the overpressure results measured at measuring points 1 and 2 by numerical simulation of four mathematical models were obtained, which were then compared with the test results of Jia Zhiwei’s gas explosion experiment, as shown in Figure 2.Figure 2
Comparison between simulation results and experimental results. (a) Comparison of results at measuring point 1. (b) Comparison of results at measuring point 2.
(a)(b)It can be seen from Figure2 that the numerical simulation results of LES turbulence model and v eddy-dissipation model (EDM) were closer to the experimental results, and the relative errors were less than 10%. Therefore, the established mathematical model of gas explosion can be used to simulate the propagation law of gas explosion shock wave.


## 2.2.1. Determination of the Mathematical Model of Gas Explosion

As shown in Figure1, the geometric model was established according to the pipe size of the gas explosion experiment in Dr. Jia Zhiwei’s thesis. The pipe sections were all 0.08 m × 0.08 m, the straight pipe length was 19.2 m, and the vertical pipe length was 5 m. Two measuring points were set up in the experiment. The measuring point 1, being 19 m away from the leftmost end, was located on the central line of the horizontal pipeline; the measuring point 2, being 1 m away from vertical pipeline, was located on the central line of the vertical pipeline.Figure 1
Geometric model diagram.In order to verify the correctness of the selected mathematical models, four commonly used models, namely, LES turbulence model + eddy-dissipation model (EDM),k−ε turbulence model + laminar flow/eddy-dissipation model, RNG k−ε + eddy-dissipation model (EDM), and k−ε turbulence model + eddy-dissipation conceptual (EDC) model, were used for numerical simulation. In the simulation of the four models, the grid spacing of the geometric model was set to 5 mm, the time step length was set to 0.001 s, and the maximum number of iterations was 30 steps. When the length of gas filling area was 4.0 m, 5.5 m, and 7.0 m, respectively, the overpressure results measured at measuring points 1 and 2 by numerical simulation of four mathematical models were obtained, which were then compared with the test results of Jia Zhiwei’s gas explosion experiment, as shown in Figure 2.Figure 2
Comparison between simulation results and experimental results. (a) Comparison of results at measuring point 1. (b) Comparison of results at measuring point 2.
(a)(b)It can be seen from Figure2 that the numerical simulation results of LES turbulence model and v eddy-dissipation model (EDM) were closer to the experimental results, and the relative errors were less than 10%. Therefore, the established mathematical model of gas explosion can be used to simulate the propagation law of gas explosion shock wave.


## 2.3. Determination of Key Parameters of Numerical Simulation of Mine Gas Explosion

According to the actual size of fully mechanized excavation face of a coal mine return air channeling, a 3D geometric model with a length of 100 m was established, as shown in Figure3. Assume that the area 10 m from the heading face is filled with gas with a concentration of 9.5%, and the ignition sources are arranged in the center of the heading face. The open boundary is set on the right side of the model, the closed boundary is set on the left side, and the roadway wall is on all sides.Figure 3
Geometric model diagram of driving surface. (a) Schematic diagram of geometric model of driving surface. (b) 3D model diagram of driving surface.


### 2.3.1. Determination of Grid Size

The grid density has great influence on the numerical results, so it is very important to determine the appropriate grid size. The error sources of numerical calculation mainly include approximate error of physical model, truncation error of differential equation, discrete error of solving region, iterative error, and rounding error. Generally speaking, the truncation error and the discrete error decrease as the grid becomes thinner, but as the grid becomes thinner, the number of discrete points increases and the rounding error also increases. As grids become denser, computing scale and storage space will increase and computing efficiency will decrease. Under the conditions of balancing computation efficiency, storage space, and accuracy, coarse grid can be used to solve the problem first, and then the grid can be continuously refined until the increase of the grid number has little influence on the calculation result. Typical grid-refinement study results are shown in Figure4.Figure 4
Relationship between solution error and grid size.The curve graph of solution error and grid size can be divided into three areas in the condition of ensuring the convergence of the solution results. (1) When the grid size is large, the calculation error is mainly truncation error. It is the error caused by the Taylor series expansion of the differential operator, which is mainly affected by the order of the element. (2) As the mesh size gradually decreases, it gradually enters the asymptotic range, and the error in this area is mainly discretization error. Generally speaking, the more dense the mesh is, the smaller the calculation error is, and the closer the calculation result is to the real solution; that is, the mesh size is located in the above two regions. (3) When the grid size is very small, the calculation error is mainly dominated by floating-point error. At this time, the smaller the grid size, the larger the error.In order to determine the grid size, the gas explosion process in heading face was numerically simulated when the grid spacing was set to 0.05 m, 0.1 m, 0.2 m, and 0.4 m (the number of grids was 14800000, 1850000, 237500, and 26000, resp.), the step length of iteration time was 0.001 s, and the maximum number of iterations was 30 steps. The simulated maximum overpressure obtained at the measuring point (80 m away from the heading face) is shown in Table1.Table 1
Simulation results of different grid sizes.
NumberGrid size (m)Grid numberMaximum overpressure time (s)Maximum overpressure (kPa)Note10.05148000000.138355.920.118500000.145351.630.22375000.171337.640.4260000.703315.2As can be seen from Table1, the grid size has a great influence on the simulation results. With the continuous refinement of the grid, the maximum overpressure at the measuring point increased gradually (the increase rate slowed down gradually). When the grid spacing was refined from 0.4 m, 0.2 m, and 0.1 m to 0.05 m, the maximum overpressure increased by 6.6%, 2.7%, and 0.6%, respectively. The duration time of maximum overpressure of shock wave was continuously shortened (the shortening rate was slowing down gradually). When the grid spacing was refined from 0.4 m, 0.2 m, and 0.1 m to 0.05 m, the duration time of maximum overpressure was shortened by 80.4%, 19.3%, and 4.8%, respectively. According to the simulation results, the results were very close when the grid spacing was 0.1 m and 0.05 m. Considering the calculation efficiency, storage space, and accuracy, the grid spacing was finally determined to be 0.1 m.


### 2.3.2. Time Step Determination

The proper selection of time step can not only ensure the convergence of simulation results, but also ensure the accuracy of numerical simulation results and improve the calculation efficiency. In order to determine the appropriate time step, the gas explosion characteristics in heading surface were simulated when the time step was set to 0.01 s, 0.001 s, and 0.0005 s, respectively, the grid spacing was 0.1 m, and the maximum number of iterations was 30 steps. The overpressure attenuation law obtained at the measured points arranged at 40 m, 50 m, 60 m, 70 m, 80 m, and 90 m from the heading face is shown in Figure5.Figure 5
The maximum overpressure and its occurrence time at different positions.As can be seen from Figure5, the time step has a great influence on the simulation results. The smaller the time step is, the greater the simulated overpressure is and the faster the peak overpressure appears. The simulation result with time step size of 0.01 s differed greatly from that with time step size of 0.001 s. The simulation result with time step size of 0.001 s was close to that with time step size of 0.0005 s, and the difference between them was less than 6%. Taking into account the calculation efficiency and accuracy, the time step was finally determined to be 0.001 s.


## 2.3.1. Determination of Grid Size

The grid density has great influence on the numerical results, so it is very important to determine the appropriate grid size. The error sources of numerical calculation mainly include approximate error of physical model, truncation error of differential equation, discrete error of solving region, iterative error, and rounding error. Generally speaking, the truncation error and the discrete error decrease as the grid becomes thinner, but as the grid becomes thinner, the number of discrete points increases and the rounding error also increases. As grids become denser, computing scale and storage space will increase and computing efficiency will decrease. Under the conditions of balancing computation efficiency, storage space, and accuracy, coarse grid can be used to solve the problem first, and then the grid can be continuously refined until the increase of the grid number has little influence on the calculation result. Typical grid-refinement study results are shown in Figure4.Figure 4
Relationship between solution error and grid size.The curve graph of solution error and grid size can be divided into three areas in the condition of ensuring the convergence of the solution results. (1) When the grid size is large, the calculation error is mainly truncation error. It is the error caused by the Taylor series expansion of the differential operator, which is mainly affected by the order of the element. (2) As the mesh size gradually decreases, it gradually enters the asymptotic range, and the error in this area is mainly discretization error. Generally speaking, the more dense the mesh is, the smaller the calculation error is, and the closer the calculation result is to the real solution; that is, the mesh size is located in the above two regions. (3) When the grid size is very small, the calculation error is mainly dominated by floating-point error. At this time, the smaller the grid size, the larger the error.In order to determine the grid size, the gas explosion process in heading face was numerically simulated when the grid spacing was set to 0.05 m, 0.1 m, 0.2 m, and 0.4 m (the number of grids was 14800000, 1850000, 237500, and 26000, resp.), the step length of iteration time was 0.001 s, and the maximum number of iterations was 30 steps. The simulated maximum overpressure obtained at the measuring point (80 m away from the heading face) is shown in Table1.Table 1
Simulation results of different grid sizes.
NumberGrid size (m)Grid numberMaximum overpressure time (s)Maximum overpressure (kPa)Note10.05148000000.138355.920.118500000.145351.630.22375000.171337.640.4260000.703315.2As can be seen from Table1, the grid size has a great influence on the simulation results. With the continuous refinement of the grid, the maximum overpressure at the measuring point increased gradually (the increase rate slowed down gradually). When the grid spacing was refined from 0.4 m, 0.2 m, and 0.1 m to 0.05 m, the maximum overpressure increased by 6.6%, 2.7%, and 0.6%, respectively. The duration time of maximum overpressure of shock wave was continuously shortened (the shortening rate was slowing down gradually). When the grid spacing was refined from 0.4 m, 0.2 m, and 0.1 m to 0.05 m, the duration time of maximum overpressure was shortened by 80.4%, 19.3%, and 4.8%, respectively. According to the simulation results, the results were very close when the grid spacing was 0.1 m and 0.05 m. Considering the calculation efficiency, storage space, and accuracy, the grid spacing was finally determined to be 0.1 m.


## 2.3.2. Time Step Determination

The proper selection of time step can not only ensure the convergence of simulation results, but also ensure the accuracy of numerical simulation results and improve the calculation efficiency. In order to determine the appropriate time step, the gas explosion characteristics in heading surface were simulated when the time step was set to 0.01 s, 0.001 s, and 0.0005 s, respectively, the grid spacing was 0.1 m, and the maximum number of iterations was 30 steps. The overpressure attenuation law obtained at the measured points arranged at 40 m, 50 m, 60 m, 70 m, 80 m, and 90 m from the heading face is shown in Figure5.Figure 5
The maximum overpressure and its occurrence time at different positions.As can be seen from Figure5, the time step has a great influence on the simulation results. The smaller the time step is, the greater the simulated overpressure is and the faster the peak overpressure appears. The simulation result with time step size of 0.01 s differed greatly from that with time step size of 0.001 s. The simulation result with time step size of 0.001 s was close to that with time step size of 0.0005 s, and the difference between them was less than 6%. Taking into account the calculation efficiency and accuracy, the time step was finally determined to be 0.001 s.


## 3. Influence of Ignition Energy on the Gas Explosive Characteristics in Heading Face

The ignition source size is a very important parameter in numerical simulation of gas explosion, which directly affects the accuracy of simulation results. The size of ignition sources that may cause gas explosions in coal mines is also different. In order to analyze the influence of point ignition source size on gas explosion, the gas explosion process in the heading surface with gas accumulation volume of 145.52 m3 was simulated under ignition source volume of 0.52 m3, 2.14 m3, 4.19 m3, 7.23 m3, and 10.30 m3.The ignition source generates energy to heat the local gas near the ignition source, making its temperature rise rapidly to the ignition temperature. Then the gas is spread by the flame to set the whole mixture on fire. According to the hot spot fire theory, the high temperature burned gas was set as the ignition source in the simulation calculation and arranged in the center of the heading face with the ignition temperature of 1600 K. Considering that different size of ignition sources will lead to change in gas accumulation area, different lengths of gas accumulation zones were set to ensure the same volume of gas explosion when the ignition source volumes were different. Figures6 and 7 show the relationship between the peak value, duration time, and the rise rate of overpressure peak value at the measuring point (arranged 80 m away from the heading face) and the size of the ignition source.Figure 6
Relationship of the overpressure and peak time with the ignition source.Figure 7
Relationship of rise rate of the overpressure with the ignition source.From Figures6 and 7, it can be concluded that (1) the larger the size of the ignition source, the greater the ignition energy, and the more gas involved in the initial reaction, and the more the free radicals generated, the faster the gas combustion rate, the shorter the initiation time, the shorter the duration time of shock wave peak, the higher the rise rate of the maximum overpressure, and the larger the value of maximum overpressure. When the size of the ignition source increased from 0.52 m3 to 2.14 m3, 4.19 m3, and 7.23 m3 (the volume of the ignition source increased to 4.12, 8.06 and 13.90 times, resp.), the duration time of shock wave peak gradually shortened to 65.2%, 65.2%, and 67.8% of the original level, and the rise rate of shock wave peak increased by 37.2%, 50.8%, and 57.5%, respectively, and the maximum overpressure increased by 1.4%, 2.3%, and 2.7%, respectively. (2) When the size of the ignition source increased to a certain value, the duration time of shock wave peak, the rise rate of the maximum overpressure, and the maximum overpressure itself all tended to be stable. When the size of the ignition source increased from 7.23 m3 to 10.3 m3 (the volume of the ignition source increased to 1.42 times of the original), the duration time of shock wave peak remained unchanged, while the rise rate of shock wave peak and the maximum overpressure only increased by 0.06% and 0.06%, respectively.When the volume of ignition source was set too small, unsuccessful ignition occurred and shock waves failed to form, resulting in the failure of simulation. The reasons may be as follows: (1) There was a minimum energy source for gas ignition. According to the spark ignition theory, the spark instantly established a small volume of high-temperature gas, the heat generated by the spark flowed rapidly to the surrounding unburned gas, the temperature in the spark volume dropped rapidly, and the temperature of the surrounding gas quickly increased to trigger chemical reactions, and then the approximately spherical combustion wave formed and spread outward, as shown in Figure8. Whether the combustion wave can develop steadily and finally form shock wave depends on the volume of ignition gas when the initial temperature drops to about the normal flame temperature. If the expansion volume of the ignition gas is too small, the rate of heat release within the approximately spherical chemical reaction zone is insufficient to compensate for the rate of heat loss to the preheated unburned gas zone. Thus, the loss of heat to the unburned gas continuously exceeds the heat generated from the chemical reaction, the temperature in the entire volume of the reaction decreases, the reaction gradually stops, and eventually the combustion wave is extinguished after only a small amount of gas is burned around the original spark. (2) From the perspective of chain reaction, the ignition process of gas explosion is basically a chain reaction process formed by the combination of multiple elementary reactions. Arrhenius pointed out that only the molecules with energy exceeding a certain activation energy can produce chemical reactions. The reaction chain of gas explosion requires a certain amount of energy to break the chemical bond and produce free radicals. To propagate the flame, the reaction rate must be fairly fast, which requires a higher ignition energy. When the ignition energy is low, it is difficult to make more gases in the internal ignition core participate in the reaction and release less heat, so the gas in the external unburned area cannot be preheated and the gas in the external unburned area cannot participate in the reaction, thus leading to the failure of the chemical reaction and the failure of forming explosion shock wave.Figure 8
Minimum flame model.In a word, the ignition energy has a certain influence on the gas explosion process, but has little influence on the overpressure value of shock wave. The ignition energy has little effect on gas explosion characteristics when the explosion accident occurs with a large amount of gas accumulation underground.


## 4. Influence of Gas Explosive Equivalent on Gas Explosive Characteristics in Heading Surface

The blast wave of gas explosion in underground roadway propagates only along the direction of roadway, which is different from free surface explosion in air. To simplify the problem, the roadway wall is assumed to be a rigid body, and the shock wave roadway wall is 100% reflected. According to the law of energy similarity (equal overpressure of explosion shock wave), it can be obtained that the mass of 1 m3 gas converted into TNT is 1.0582 kg [30].In order to study the influence of gas explosion equivalent on the explosion propagation law in heading face, the propagation law of gas explosion shock wave was simulated with gas accumulation of 31.16 m3, 77.91 m3, 155.82 m3, 311.64 m3, 467.46 m3, and 779.10 m3. When gas accumulation amount in the heading face is 467.46 m3, the propagation law of gas explosion shock wave and the distribution law of overpressure along the length direction of the driving roadway at different moments are shown in Figures 9 and 10, respectively.Figure 9
Propagation law of gas explosion shock wave at heading face (gas accumulation 467.46 m3). (a) Law of gas explosion shock wave propagation in heading face. (b) Law of gas explosion shock wave propagation in heading face (X = 2.45 m).
(a)(b)Figure 10
Distribution law of overpressure on driving surface (X = 2.45 m, Y = 1.6 m).It can be seen from Figure9 that, in the early stage, the gas combustion and explosion in the gas accumulation area in the heading face formed shock wave, and the shock wave overpressure in the gas accumulation area gradually increased. With the extension of time, the shock wave in the air zone propagated continuously along the tunneling roadway, the overpressure attenuated continuously, and the shock wave gradually changed from a curved surface wave to a plane wave. The larger the volume of gas accumulation, the greater the explosion shock wave energy, the faster the propagation speed of shock wave, the greater the maximum overpressure at various points, and the farther the distance of shock wave changing from surface wave to plane wave.Figure10 shows the distribution law of maximum overpressure of gas explosion shock wave in the heading face along the length of the driving roadway obtained by simulation with different gas accumulation volume. Figure 11 shows the relationship between the maximum overpressure at measuring points with different distances away from the heading face and the explosive equivalent.Figure 11
Relationship between overpressure peak and explosive equivalent at different positions in the heading surface.The following can be seen from Figure10. (1) In the gas accumulation area, the gas kept burning violently, and the shock wave gradually formed and spread forward. The overpressure of the shock wave increased roughly as a quadratic curve function along the axial direction of the tunneling roadway. With the increase of the propagation distance, the increasing speed of the shock wave overpressure gradually slowed down, especially near the boundary of the gas accumulation area where the overpressure dropped sharply. When the volume of gas involved in the explosion was different, the overpressure in other areas remained basically unchanged, except that the overpressure near the boundary of the gas accumulation area dropped sharply. (2) In the general air zone, the overpressure peak value of explosion shock wave tended to decline as a power function along the axial direction of driving roadway. With the continuous forward propagation of shock wave, the attenuation velocity of overpressure decreased with the increase of propagation distance. The larger the volume of gas involved in the explosion, the faster the attenuation speed of shock wave overpressure along the tunneling roadway. When 467.46 m3 and 155.82 m3 of gas accumulated in the heading face, the shock wave propagated from 40 m to 50 m away from the heading face, and the overpressure decreased by 98.2 kPa and 39.6 kPa, respectively. When the shock wave propagated from 80 m to 90 m away from the heading face, the overpressure decreased by 83.5 kPa and 14.8 kPa.It can be seen from Figure11 that the volume of gas involved in the explosion has a great influence on the overpressure of the explosion shock wave. The larger the volume of gas involved in the explosion, that is, the greater the explosive equivalent, the greater the overpressure of the explosion shock wave. The maximum overpressure of shock wave increased with the explosion equivalent as a power function, and the relation between them can be fitted as Δp=A⋅ωB. Therefore, it is necessary to timely monitor the gas concentration in places where gas accumulation is easy, such as the heading face, increase air volume, and develop gas extraction approaches, so as to a large amount accumulation of high-concentration gas, reduce the risk of human casualties and property losses and ensure the safe and stable mine production.


## 5. Conclusions

(1)
The mathematical model and space-time step for gas explosion numerical simulation are determined. The determined mathematical models (LES and EDM) for numerical simulation of gas explosion can accurately simulate the gas explosion and its propagation process. The time-space step size has a great influence on the simulation results, and the time-space step size is determined to be 0.1 m and 0.001 s.(2)
The ignition energy has a certain influence on the gas explosion process, but has little influence on the overpressure of shock wave. Ignition energy can speed up gas explosion to some extent, but it cannot change shock wave overpressure. The ignition energy has little effect on gas explosion characteristics when gas explosion occurs with a large volume of gas accumulation.(3)
The gas explosion equivalent has great influence on the overpressure of gas explosion shock wave. The maximum overpressure of the shock wave increases with the explosive equivalent as a power function. In the gas accumulation area, the overpressure of the shock wave rises along the axial direction of the tunneling roadway as a conic function. When the volume of gas involved in the explosion is different, the overpressure in other areas is basically the same except that the overpressure near the boundary of the gas accumulation area drops sharply. In general air zone, the overpressure peak value of explosion shock wave attenuates along the axial direction of driving roadway as a power function. The research results provide theoretical guidance for accurately grasping the propagation law of gas explosion in underground roadway and for the research and development of new technologies for prevention and control of gas explosion.



---

*Source: 8871865-2020-10-29.xml*