# Visual Simulation and Case Inversion of Gas Explosion in Underground Mine

Authors: Qiu Jinwei; Jiang Bingyou; Tang Mingyun; Zhou Liang; Ren Bo; Yang Yingdi
Journal: Advances in Civil Engineering
 (2021)

---

## Abstract

To truly represent the developing changes and spatiotemporal effects of gas explosion in underground mine, this paper discussed the field modeling and visual simulation technology of mine gas explosion. Through analyzing the characteristics of various field parameters, space geometric field models, physical attribute field models, and meteorological parameter field models were established to describe mine engineering, ventilation network, and explosion disaster process, respectively. Based on the generated model class and simulation timing flow of gas explosion disaster field, OpenGL technology was adopted to visualize the model of gas explosion disaster field, which simulated the overpressure propagation process of gas explosion shockwave and dynamic spread process of disaster gas, thereby to reveal the developing changes of gas explosion parameters in time and space. Taking the gas explosion accident that occurred in a certain mine in China as an example, the disaster process was inferred to examine the availability of the mine gas explosion disaster field model. The research achievements in the paper not only embody the great engineering value of visual simulation technology in the field of safety engineering but also provide references for accident impact prediction, assessment, and emergency plan.


---

## Body

## 1. Introduction

The simulation of mine gas explosion is involved in multiple fields and multiple disciplines such as computer technology, ventilation network, and fluid dynamics. With the improvement of computer technology, the emergence of network solution algorithm [1] and visualization technology provides an effective means and method for the disaster simulation of mine gas explosion. Owing to the complexity of natural conditions of underground space and the abruptness of the gas explosion process, conducting the large-scale mine gas explosion simulation tests turns out to be impractical. Consequently, it appears to be particularly significant to figure out how the simulation computing technology can be used effectively to visualize the developing process of mine gas explosion, to truly represent the evolution process of mine gas explosion, and to provide support for mine disaster relief as well as avoiding the disasters.The shockwave, high-temperature airflow, and the interaction of disaster gas and ventilation system, which are generated by gas explosion disaster, expand the scope of the disaster. The dynamic developing process of the gas explosion is an issue of complicated system engineering; thus professional software must be adopted to achieve the simulation of its disaster process. Regarding the methods, it is relatively mature to use the numerical simulation software, such as AutoReaGas, FLUENT, and FLCAS, to achieve the research on mine gas explosion. Using AutoReaGas, Jiang et al. [2] simulated the multiparameter attenuation characteristics of gas explosions and revealed the effects of changes in fuel volume on the explosion-proof distance, which provided a theoretical basis of research on explosion-proof and suppression devices for the underground gas pipeline. Sun et al. [3] established a pipe model with a length of 100 m and simulated the propagation characteristics of gas deflagrations for pipe cross-sectional areas. Based on FLUENT simulation software, Sun et al. [4] explored the influences of different ignition conditions, such as ignition locations, ignition temperatures, and ignition area, on the deflagration characteristics of the premixed H2/air mixture in a closed pipe with length of 1000 mm. Luo et al. [5] simulated the gas explosion process of different concentrations of methane using the software FLACS. Nevertheless, a large number of numerical simulation software programs are widely used to simulate gas explosion, but as a result of the limitation to the capacity of calculation, the numerical simulation software can only simulate the gas explosion in a partial roadway in the underground mine. Few studies on numerical simulation of the whole mine have been reported.In terms of software engineering, most of the software programs are designed for mine fire simulation, of which the representative ones include MFIRE, POZAR, uQEM, CFIRE, and MFRDSS, developed by the scientific research institutes from the United States, Poland, Australia, and China [6–9]. Currently, the software programs of fire simulation are mainly reflected on the differential equation method and the time interval method [10], where the former obtains the arithmetic solution by establishing the differential equations and the finite difference method, while the latter carries out the solution by adopting the “Hattie-Claus Method” or “Newton Method” [11, 12] and using the virtual simulation technology [13, 14] and 3D technology [15, 16] to describe the fire propagation process. Nevertheless, few studies have been conducted on the macroscopic ventilation dynamic system, the wind turbulence situation, the gas explosion disaster field, and the spatiotemporal evolution relationship of the gas explosion process caused by a gas explosion; in particular, there is a lack in the analysis of the 3D mine gas explosion disaster simulation and the development of visualization simulation software.Therefore, based on the predecessor’s research foundation of the numerical simulation in the field of mine gas explosion, the analysis and modeling on the characteristics of each field variable factor in the process of mine gas explosion disaster were conducted, and a simulation system of the gas explosion disaster process that was integrated of simulation calculation, data management, and a spatiotemporal process by adopting OpenGL three-dimensional visualization method was developed to achieve the visual expression of the mine gas explosion disaster process and represent its disaster process.


## 2. Establishment of the Field Variable Model for Gas Explosion Disaster Process

Based on the theory of mine compound field [17], the digital mine not only includes geometric information such as the spatial shape of the ore body and the spatial relationship of the mine roadway, the physical attribute information such as the material composition, and mechanical parameters of the ore body and the wall rock but also covers meteorological attribute information such as the temperature, humidity, and noxious gas concentration in underground mine.Consequently, regarding the spatiotemporal influence of the occurrence period of mine gas explosion and the postdisaster developing changes on the underground environment, the disaster field of underground mine gas explosion is proposed, which is a complex model composed of multi-information and multiple attributes [10, 18], as is shown in Figure 1. The field variable parameters of the gas explosion in an underground mine include spatial geometric field variable, physical attribute field variable, and meteorological attribute field variable.Figure 1
Disaster field variable model of gas explosion.


### 2.1. The Analysis of Characteristics of Each Field Variable Factor

From the analysis of the characteristics of each field variable parameter, it can be concluded that the gas explosion disaster field is reflected in different characteristics, mainly including the spatial geometric characteristics, attribute distribution characteristics, running status characteristics, vector direction characteristics, and temporal characteristics.The spatial geometric characteristics are used to describe the position, shape, and spatial relationship of each mine roadway, node, and related entities (structures, fans, and explosion sources) in the field variable. Additionally, they also include the relatively complex topological relationship between entities.The attribute distribution characteristics are used to describe attribute parameters that are closely related to the field variable attribute and correspond to a certain time and space position, such as wind volume, wind speed, cross-sectional area, length, wind resistance, and wind pressure, of each mine roadway.The running status characteristics are used to describe the operation condition of fan equipment, ventilation door, and firewall in the gas explosion disaster field, mainly referring to the opening and closing condition of the fan and the intact and ineffective condition of the ventilation structure.The vector direction characteristics are used to describe the direction information about airflow, shock wave, disaster gas, disaster avoidance, and disaster relief personnel movement in the field variables. While carrying out the analysis of vector direction characteristics, it is necessary to describe the spatiotemporal characteristics of the gas explosion process, such as the propagation direction of gas explosion shockwave overpressure and disaster gas in a certain space position.The temporal characteristics mean that all the properties in the gas explosion disaster field are time-sequenced, and the dynamic change process appears with the occurrence of the explosion and the development of the shockwave and the disaster gas. Currently, the motion and propagation of disaster gas in the roadway are usually simulated by the “method of smoke plume frontal surface” based on the temporal characteristics [19], but the period of propagation of the shockwave overpressure in the mine roadway is as short as several seconds to the longest or several milliseconds to the shortest. Therefore, it is difficult to subdivide the period again.


### 2.2. Modeling of the Spatial Geometric Field Variable

Modeling of spatial geometric field variable uses the 3D structural model theory and computerized 3D visualization technology to organize and draw the volume element of gas explosion field, which is used to express the spatial position, spatial relationship, and geometric shape of the field variable. It usually uses points, lines, faces, and other graphics to represent the spatial structure of the ventilation network. During the modeling process of the spatial geometric field variable, it is necessary to consider the topological adjacency between nodes and nodes, nodes and roadways, roadways and structures, and explosion sources and roadways.


### 2.3. Modeling of the Physical Attribute Field Variable

The model of a physical attribute during the gas explosion disaster is mainly used to express the spatial distribution of the physical attributes of the wind flow in the ventilation network. The modeling of physical attribute field variable is based on the topology of the ventilation network, and the attribute parameters (e.g., wind volume, wind pressure, wind resistance, resistance coefficient, etc.) of the roadways are added to the model of geometric field variable to form a ventilation shaft model and vectorize each roadway. A topological relationship model of the ventilation network is established based on the topological relationship between roadways and roadways. Meanwhile, according to the ventilation network model, such as the ventilation network solution model, the overpressure propagation model, and the ventilation dynamic model, the physical attribute field variable of the whole network structure is modeled by adding the field variable information of each directed edge to complete the ventilation network solution.


### 2.4. Modeling of the Meteorological Attribute Field Variable

Gas explosion shockwaves and high-speed airflow propagate in the ventilation network for only a few tens of seconds, while the second round of damage is the propagation of high-temperature airflow in the network, which has a long propagation period. As a result of the high temperature, fire pressure is generated during air propagation. Modeling of meteorological attribute field variable is based on the unsteady network solution calculation, combined with the models of wind flow temperature propagation, fire wind pressure, and disaster gas propagation. Through the calculation of airflow temperature, disaster gas concentration, fire wind pressure, shockwave overpressure, and so forth, it completes the representation of the spatial-temporal variation in various condition parameters following the disaster. As a result, the modeling of the meteorological field variable during the mine gas explosion disaster is achieved.


#### 2.4.1. Modeling of Wind Temperature Propagation and Mixing Field Variable

The main models of airflow temperature after the gas explosion are the wind temperature calculation model and the wind temperature mixing model [20].Wind temperature calculation of horizontal roadways is as follows:(1)T=T0+T1−T0e−abUx/mCp,where T0 is wind temperature before disaster, K; T1 is wind temperature of disaster gas at roadway entry end, K; αb is unstable heat transfer coefficient of disaster gas and roadway wall; m is mass flow of disaster gas, kg/s; x is distance from the air inlet of the explosion source, m; Cp is constant pressure specific heat, J/kg⋅K; and U is perimeter of roadway, m.Wind temperature calculation of inclined roadways is as follows:(2)T=T0−mLgαbU+T1−T0+mLgαbUexp−αbUxmCp,where Lis slope of roadway, m and g is gravity acceleration, m/s2.When calculating the wind-temperature mixing of disaster gas at the junction of roadways, it is necessary to treat the high-temperature disaster gas after the gas explosion as a constant specific heat gas. The wind-temperature mixing calculation at the junction mainly adopts the conservation law of energy, and its mixing formula is as follows:(3)Th=∑i=1nmiT2i∑i=1nmi,where T2 is wind temperature at return air end of roadway, K; n is number of branches flowing into the node; and mi is mass flow of disaster gas in roadway i, kg/s.


#### 2.4.2. Modeling of Fire Wind Pressure Field Variable of Explosive High-Temperature Gas

At present, there is no consensus on the calculation method of fire wind pressure in academia [14, 21, 22]. In the paper, the temperature change before and after the explosion is considered and derived from the gas equation:(4)ΔH=zγ0Ta−TbTa,where ΔH is fire wind pressure value of explosive high-temperature gas, Pa; z is the height difference between the beginning and the end of the roadway, m; γ0 is average air weight rate in the roadway before and after the explosion, kg/m2⋅s2; Tb and Ta are average absolute air temperatures in the roadway before and after a gas explosion, respectively, K.


#### 2.4.3. Modeling of Disaster Gas Concentration Field Variable

Gas explosion disaster gas enters the network nodes and mixes with the airflow of other roadways, the disaster gas is diluted, and the pollution area is enlarged. Currently, the disaster gas transmission theory mainly includes gradient transport theory and turbulence statistical theory [23, 24]. Regarding the free space, the models that are widely used are Gaussian model, BM model, Sutton model, FEM finite element model, and so forth [25–27]. The mine gas explosion ends in an instant; thus the generation process of the disaster gas can be assumed to be an instantaneous gas leakage process. The main disaster gas that is generated in the gas explosion is CO, and the Gaussian model is mainly suitable for the light gas or the gas whose density is not much different from air density. Consequently, this research uses the Gaussian model to study the diffusion model of CO disaster gas.Model of disaster gas propagation in ventilation network. Suppose that the ventilation network is GV,E(V=mE=n), and the disastrous gas is generated in roadways ej at the initial moment; then the distribution of the disastrous gas concentration in roadways is as follows:(5)cjx,t=cj0+cj′x,t,where cj0 is a steady gas concentration in the roadway; x represents location coordinates with the origin of concentrated disaster gas; cj′x,t is concentrated disaster gas concentration, which can be expressed as(6)cj′x,t=Mj4πExjte−x−ujt2/4Exjt,where Mj, uj, and Exj are mass flow rate, wind speed, and dispersion coefficient of disastrous gas produced by a gas explosion in the roadway ej, respectively.When the disaster gas spreads to the node, it mixes with the airflow in other roadways and then distributes to the airflow in other roadways. As shown in Figure2, if the gas from a gas explosion is evenly mixed at the nodes, the concentration of the disaster gas at the nodes is as follows:(7)Dit=∑j=1nbijQjcjxj,t∑j=1nbijQj.where bij is the coefficient, when the endpoint of the roadway ej is node i: bij = 1; when the endpoint of the roadway ej is node j, bij = 0;Dit is disaster gas concentration at nodeiat timet; Qjj=1,2,3,…,n is air volume of the roadway j, m3/s; xj is the distance between node i and origin, m.Figure 2
Toxic gas diffusion of the explosion in network nodes.After convergence, the hazard gas spreads downwind with the concentration value of nodei, and the gas concentration in the roadway is as follows:(8)cikx,t=cik0x,t+kicikx,t,where ki is the node dilution factor of disaster gas and cik is disaster gas concentration in the roadway ek with node i as the start point, as shown as follows:(9)ki=Ditcjxj,t,cik=Mik4πExkte−x−uktk2/4Ekjt.where uk is the wind speed of the roadway ek; Exk is the dispersion coefficient of the roadway ek; tk=t−xj/uj is propagation time in the roadway ek.The amount of disastrous gas entering the roadwayek is calculated as(10)Mik=MjQk∑bijQj.Through formulas (5)–(10) and combined with unsteady ventilation network theory, the dissemination range of disastrous gas generated by a gas explosion in mine ventilation network can be calculated.


#### 2.4.4. Modeling of Shockwave Overpressure Field Variable

In the ventilation network, there are not only many straight roadways but also a large number of turning roadways, bifurcate roadways, and abrupt roadways. These roadways form a complex ventilation network. When the shockwave overpressure propagates in the complex network, the approximate formula method is used to solve the problem.Suppose that a selected point is a pointA, and the shockwave overpressure value between point B to be sought and point A isΔp (regardless of the partial attenuation coefficient). If there is a local turning, bifurcation, or abrupt structure between B and A, the local attenuation coefficient needs to be considered. The local overpressure attenuation coefficient is solved by the formula [28], and the product of the local overpressure attenuation coefficient between the two points is fk,δ,θ,Ω,r:(11)Δp′=Δpfk,δ,θ,Ω,r,where Δp′represents the shockwave overpressure of point B with a distance of r from the selected point A, Pa; Δp denotes the shockwave overpressure between point B and point A without considering the local attenuation coefficient, Pa； and fk,δ,θ,Ω,r represents the product of local overpressure attenuation coefficient between B and A.No matter how many times the local resistance is encountered, the shockwave attenuation coefficient can be substituted for equation (11) to obtain the shockwave overpressure value of point B.


## 2.1. The Analysis of Characteristics of Each Field Variable Factor

From the analysis of the characteristics of each field variable parameter, it can be concluded that the gas explosion disaster field is reflected in different characteristics, mainly including the spatial geometric characteristics, attribute distribution characteristics, running status characteristics, vector direction characteristics, and temporal characteristics.The spatial geometric characteristics are used to describe the position, shape, and spatial relationship of each mine roadway, node, and related entities (structures, fans, and explosion sources) in the field variable. Additionally, they also include the relatively complex topological relationship between entities.The attribute distribution characteristics are used to describe attribute parameters that are closely related to the field variable attribute and correspond to a certain time and space position, such as wind volume, wind speed, cross-sectional area, length, wind resistance, and wind pressure, of each mine roadway.The running status characteristics are used to describe the operation condition of fan equipment, ventilation door, and firewall in the gas explosion disaster field, mainly referring to the opening and closing condition of the fan and the intact and ineffective condition of the ventilation structure.The vector direction characteristics are used to describe the direction information about airflow, shock wave, disaster gas, disaster avoidance, and disaster relief personnel movement in the field variables. While carrying out the analysis of vector direction characteristics, it is necessary to describe the spatiotemporal characteristics of the gas explosion process, such as the propagation direction of gas explosion shockwave overpressure and disaster gas in a certain space position.The temporal characteristics mean that all the properties in the gas explosion disaster field are time-sequenced, and the dynamic change process appears with the occurrence of the explosion and the development of the shockwave and the disaster gas. Currently, the motion and propagation of disaster gas in the roadway are usually simulated by the “method of smoke plume frontal surface” based on the temporal characteristics [19], but the period of propagation of the shockwave overpressure in the mine roadway is as short as several seconds to the longest or several milliseconds to the shortest. Therefore, it is difficult to subdivide the period again.


## 2.2. Modeling of the Spatial Geometric Field Variable

Modeling of spatial geometric field variable uses the 3D structural model theory and computerized 3D visualization technology to organize and draw the volume element of gas explosion field, which is used to express the spatial position, spatial relationship, and geometric shape of the field variable. It usually uses points, lines, faces, and other graphics to represent the spatial structure of the ventilation network. During the modeling process of the spatial geometric field variable, it is necessary to consider the topological adjacency between nodes and nodes, nodes and roadways, roadways and structures, and explosion sources and roadways.


## 2.3. Modeling of the Physical Attribute Field Variable

The model of a physical attribute during the gas explosion disaster is mainly used to express the spatial distribution of the physical attributes of the wind flow in the ventilation network. The modeling of physical attribute field variable is based on the topology of the ventilation network, and the attribute parameters (e.g., wind volume, wind pressure, wind resistance, resistance coefficient, etc.) of the roadways are added to the model of geometric field variable to form a ventilation shaft model and vectorize each roadway. A topological relationship model of the ventilation network is established based on the topological relationship between roadways and roadways. Meanwhile, according to the ventilation network model, such as the ventilation network solution model, the overpressure propagation model, and the ventilation dynamic model, the physical attribute field variable of the whole network structure is modeled by adding the field variable information of each directed edge to complete the ventilation network solution.


## 2.4. Modeling of the Meteorological Attribute Field Variable

Gas explosion shockwaves and high-speed airflow propagate in the ventilation network for only a few tens of seconds, while the second round of damage is the propagation of high-temperature airflow in the network, which has a long propagation period. As a result of the high temperature, fire pressure is generated during air propagation. Modeling of meteorological attribute field variable is based on the unsteady network solution calculation, combined with the models of wind flow temperature propagation, fire wind pressure, and disaster gas propagation. Through the calculation of airflow temperature, disaster gas concentration, fire wind pressure, shockwave overpressure, and so forth, it completes the representation of the spatial-temporal variation in various condition parameters following the disaster. As a result, the modeling of the meteorological field variable during the mine gas explosion disaster is achieved.


### 2.4.1. Modeling of Wind Temperature Propagation and Mixing Field Variable

The main models of airflow temperature after the gas explosion are the wind temperature calculation model and the wind temperature mixing model [20].Wind temperature calculation of horizontal roadways is as follows:(1)T=T0+T1−T0e−abUx/mCp,where T0 is wind temperature before disaster, K; T1 is wind temperature of disaster gas at roadway entry end, K; αb is unstable heat transfer coefficient of disaster gas and roadway wall; m is mass flow of disaster gas, kg/s; x is distance from the air inlet of the explosion source, m; Cp is constant pressure specific heat, J/kg⋅K; and U is perimeter of roadway, m.Wind temperature calculation of inclined roadways is as follows:(2)T=T0−mLgαbU+T1−T0+mLgαbUexp−αbUxmCp,where Lis slope of roadway, m and g is gravity acceleration, m/s2.When calculating the wind-temperature mixing of disaster gas at the junction of roadways, it is necessary to treat the high-temperature disaster gas after the gas explosion as a constant specific heat gas. The wind-temperature mixing calculation at the junction mainly adopts the conservation law of energy, and its mixing formula is as follows:(3)Th=∑i=1nmiT2i∑i=1nmi,where T2 is wind temperature at return air end of roadway, K; n is number of branches flowing into the node; and mi is mass flow of disaster gas in roadway i, kg/s.


### 2.4.2. Modeling of Fire Wind Pressure Field Variable of Explosive High-Temperature Gas

At present, there is no consensus on the calculation method of fire wind pressure in academia [14, 21, 22]. In the paper, the temperature change before and after the explosion is considered and derived from the gas equation:(4)ΔH=zγ0Ta−TbTa,where ΔH is fire wind pressure value of explosive high-temperature gas, Pa; z is the height difference between the beginning and the end of the roadway, m; γ0 is average air weight rate in the roadway before and after the explosion, kg/m2⋅s2; Tb and Ta are average absolute air temperatures in the roadway before and after a gas explosion, respectively, K.


### 2.4.3. Modeling of Disaster Gas Concentration Field Variable

Gas explosion disaster gas enters the network nodes and mixes with the airflow of other roadways, the disaster gas is diluted, and the pollution area is enlarged. Currently, the disaster gas transmission theory mainly includes gradient transport theory and turbulence statistical theory [23, 24]. Regarding the free space, the models that are widely used are Gaussian model, BM model, Sutton model, FEM finite element model, and so forth [25–27]. The mine gas explosion ends in an instant; thus the generation process of the disaster gas can be assumed to be an instantaneous gas leakage process. The main disaster gas that is generated in the gas explosion is CO, and the Gaussian model is mainly suitable for the light gas or the gas whose density is not much different from air density. Consequently, this research uses the Gaussian model to study the diffusion model of CO disaster gas.Model of disaster gas propagation in ventilation network. Suppose that the ventilation network is GV,E(V=mE=n), and the disastrous gas is generated in roadways ej at the initial moment; then the distribution of the disastrous gas concentration in roadways is as follows:(5)cjx,t=cj0+cj′x,t,where cj0 is a steady gas concentration in the roadway; x represents location coordinates with the origin of concentrated disaster gas; cj′x,t is concentrated disaster gas concentration, which can be expressed as(6)cj′x,t=Mj4πExjte−x−ujt2/4Exjt,where Mj, uj, and Exj are mass flow rate, wind speed, and dispersion coefficient of disastrous gas produced by a gas explosion in the roadway ej, respectively.When the disaster gas spreads to the node, it mixes with the airflow in other roadways and then distributes to the airflow in other roadways. As shown in Figure2, if the gas from a gas explosion is evenly mixed at the nodes, the concentration of the disaster gas at the nodes is as follows:(7)Dit=∑j=1nbijQjcjxj,t∑j=1nbijQj.where bij is the coefficient, when the endpoint of the roadway ej is node i: bij = 1; when the endpoint of the roadway ej is node j, bij = 0;Dit is disaster gas concentration at nodeiat timet; Qjj=1,2,3,…,n is air volume of the roadway j, m3/s; xj is the distance between node i and origin, m.Figure 2
Toxic gas diffusion of the explosion in network nodes.After convergence, the hazard gas spreads downwind with the concentration value of nodei, and the gas concentration in the roadway is as follows:(8)cikx,t=cik0x,t+kicikx,t,where ki is the node dilution factor of disaster gas and cik is disaster gas concentration in the roadway ek with node i as the start point, as shown as follows:(9)ki=Ditcjxj,t,cik=Mik4πExkte−x−uktk2/4Ekjt.where uk is the wind speed of the roadway ek; Exk is the dispersion coefficient of the roadway ek; tk=t−xj/uj is propagation time in the roadway ek.The amount of disastrous gas entering the roadwayek is calculated as(10)Mik=MjQk∑bijQj.Through formulas (5)–(10) and combined with unsteady ventilation network theory, the dissemination range of disastrous gas generated by a gas explosion in mine ventilation network can be calculated.


### 2.4.4. Modeling of Shockwave Overpressure Field Variable

In the ventilation network, there are not only many straight roadways but also a large number of turning roadways, bifurcate roadways, and abrupt roadways. These roadways form a complex ventilation network. When the shockwave overpressure propagates in the complex network, the approximate formula method is used to solve the problem.Suppose that a selected point is a pointA, and the shockwave overpressure value between point B to be sought and point A isΔp (regardless of the partial attenuation coefficient). If there is a local turning, bifurcation, or abrupt structure between B and A, the local attenuation coefficient needs to be considered. The local overpressure attenuation coefficient is solved by the formula [28], and the product of the local overpressure attenuation coefficient between the two points is fk,δ,θ,Ω,r:(11)Δp′=Δpfk,δ,θ,Ω,r,where Δp′represents the shockwave overpressure of point B with a distance of r from the selected point A, Pa; Δp denotes the shockwave overpressure between point B and point A without considering the local attenuation coefficient, Pa； and fk,δ,θ,Ω,r represents the product of local overpressure attenuation coefficient between B and A.No matter how many times the local resistance is encountered, the shockwave attenuation coefficient can be substituted for equation (11) to obtain the shockwave overpressure value of point B.


## 2.4.1. Modeling of Wind Temperature Propagation and Mixing Field Variable

The main models of airflow temperature after the gas explosion are the wind temperature calculation model and the wind temperature mixing model [20].Wind temperature calculation of horizontal roadways is as follows:(1)T=T0+T1−T0e−abUx/mCp,where T0 is wind temperature before disaster, K; T1 is wind temperature of disaster gas at roadway entry end, K; αb is unstable heat transfer coefficient of disaster gas and roadway wall; m is mass flow of disaster gas, kg/s; x is distance from the air inlet of the explosion source, m; Cp is constant pressure specific heat, J/kg⋅K; and U is perimeter of roadway, m.Wind temperature calculation of inclined roadways is as follows:(2)T=T0−mLgαbU+T1−T0+mLgαbUexp−αbUxmCp,where Lis slope of roadway, m and g is gravity acceleration, m/s2.When calculating the wind-temperature mixing of disaster gas at the junction of roadways, it is necessary to treat the high-temperature disaster gas after the gas explosion as a constant specific heat gas. The wind-temperature mixing calculation at the junction mainly adopts the conservation law of energy, and its mixing formula is as follows:(3)Th=∑i=1nmiT2i∑i=1nmi,where T2 is wind temperature at return air end of roadway, K; n is number of branches flowing into the node; and mi is mass flow of disaster gas in roadway i, kg/s.


## 2.4.2. Modeling of Fire Wind Pressure Field Variable of Explosive High-Temperature Gas

At present, there is no consensus on the calculation method of fire wind pressure in academia [14, 21, 22]. In the paper, the temperature change before and after the explosion is considered and derived from the gas equation:(4)ΔH=zγ0Ta−TbTa,where ΔH is fire wind pressure value of explosive high-temperature gas, Pa; z is the height difference between the beginning and the end of the roadway, m; γ0 is average air weight rate in the roadway before and after the explosion, kg/m2⋅s2; Tb and Ta are average absolute air temperatures in the roadway before and after a gas explosion, respectively, K.


## 2.4.3. Modeling of Disaster Gas Concentration Field Variable

Gas explosion disaster gas enters the network nodes and mixes with the airflow of other roadways, the disaster gas is diluted, and the pollution area is enlarged. Currently, the disaster gas transmission theory mainly includes gradient transport theory and turbulence statistical theory [23, 24]. Regarding the free space, the models that are widely used are Gaussian model, BM model, Sutton model, FEM finite element model, and so forth [25–27]. The mine gas explosion ends in an instant; thus the generation process of the disaster gas can be assumed to be an instantaneous gas leakage process. The main disaster gas that is generated in the gas explosion is CO, and the Gaussian model is mainly suitable for the light gas or the gas whose density is not much different from air density. Consequently, this research uses the Gaussian model to study the diffusion model of CO disaster gas.Model of disaster gas propagation in ventilation network. Suppose that the ventilation network is GV,E(V=mE=n), and the disastrous gas is generated in roadways ej at the initial moment; then the distribution of the disastrous gas concentration in roadways is as follows:(5)cjx,t=cj0+cj′x,t,where cj0 is a steady gas concentration in the roadway; x represents location coordinates with the origin of concentrated disaster gas; cj′x,t is concentrated disaster gas concentration, which can be expressed as(6)cj′x,t=Mj4πExjte−x−ujt2/4Exjt,where Mj, uj, and Exj are mass flow rate, wind speed, and dispersion coefficient of disastrous gas produced by a gas explosion in the roadway ej, respectively.When the disaster gas spreads to the node, it mixes with the airflow in other roadways and then distributes to the airflow in other roadways. As shown in Figure2, if the gas from a gas explosion is evenly mixed at the nodes, the concentration of the disaster gas at the nodes is as follows:(7)Dit=∑j=1nbijQjcjxj,t∑j=1nbijQj.where bij is the coefficient, when the endpoint of the roadway ej is node i: bij = 1; when the endpoint of the roadway ej is node j, bij = 0;Dit is disaster gas concentration at nodeiat timet; Qjj=1,2,3,…,n is air volume of the roadway j, m3/s; xj is the distance between node i and origin, m.Figure 2
Toxic gas diffusion of the explosion in network nodes.After convergence, the hazard gas spreads downwind with the concentration value of nodei, and the gas concentration in the roadway is as follows:(8)cikx,t=cik0x,t+kicikx,t,where ki is the node dilution factor of disaster gas and cik is disaster gas concentration in the roadway ek with node i as the start point, as shown as follows:(9)ki=Ditcjxj,t,cik=Mik4πExkte−x−uktk2/4Ekjt.where uk is the wind speed of the roadway ek; Exk is the dispersion coefficient of the roadway ek; tk=t−xj/uj is propagation time in the roadway ek.The amount of disastrous gas entering the roadwayek is calculated as(10)Mik=MjQk∑bijQj.Through formulas (5)–(10) and combined with unsteady ventilation network theory, the dissemination range of disastrous gas generated by a gas explosion in mine ventilation network can be calculated.


## 2.4.4. Modeling of Shockwave Overpressure Field Variable

In the ventilation network, there are not only many straight roadways but also a large number of turning roadways, bifurcate roadways, and abrupt roadways. These roadways form a complex ventilation network. When the shockwave overpressure propagates in the complex network, the approximate formula method is used to solve the problem.Suppose that a selected point is a pointA, and the shockwave overpressure value between point B to be sought and point A isΔp (regardless of the partial attenuation coefficient). If there is a local turning, bifurcation, or abrupt structure between B and A, the local attenuation coefficient needs to be considered. The local overpressure attenuation coefficient is solved by the formula [28], and the product of the local overpressure attenuation coefficient between the two points is fk,δ,θ,Ω,r:(11)Δp′=Δpfk,δ,θ,Ω,r,where Δp′represents the shockwave overpressure of point B with a distance of r from the selected point A, Pa; Δp denotes the shockwave overpressure between point B and point A without considering the local attenuation coefficient, Pa； and fk,δ,θ,Ω,r represents the product of local overpressure attenuation coefficient between B and A.No matter how many times the local resistance is encountered, the shockwave attenuation coefficient can be substituted for equation (11) to obtain the shockwave overpressure value of point B.


## 3. Simulation Technology Framework and Function Module of the Gas Explosion Disaster Process

The simulation system uses Windows as the platform, selects the VB.net object-oriented language, uses ODBC and ADO technology to access the Access database, and adopts the internationally recognized 3D graphics industry-standard OpenGL as the technical support [29] to realize the 3D visualization simulation of mine gas explosion disaster process.


### 3.1. Establishment of Simulation Model Class for the Disaster Process

According to the aforementioned gas explosion disaster field model, the simulation model of the disaster process is generated. The most important classes are wind roadway class, wind network class, roadway condition class, explosion source class, structure class, gas explosion information class, shockwave overpressure propagation class, and spreading of gas explosion disaster gas class. The relationships between classes and classes are shown in Figure3. The figure shows the basic data structure of the model classes involved in the gas explosion simulation process. The basic data members such as the space location, network topology relationship, and attribute information are listed in the given data structure.Figure 3
Visual model class diagram of gas explosion disaster process.


### 3.2. Time Sequence Flow of Disaster Process Simulation

Based on the establishment of the visualization model class of gas explosion disaster process and the analysis structure of time sequence characteristics of field variable, the timing flow diagrams of gas explosion disaster process simulation are given. Figures4 and 5 are the simulation timing flow diagram of shockwave overpressure propagation and disaster gas simulation spreading.Figure 4
Simulation time sequence diagram of shockwave propagation of gas explosion.Figure 5
Simulation time sequence diagram of gas explosion in gas explosion.The specialized time sequence flow of gas explosion shockwave propagation simulation is as follows: first invoke ExpPos class and ExpPre class to determine the location in the explosion source, the highest pressure of the explosion source, and other related parameters and then transfer these parameters through the GetExpPos and GetExpPre methods. Then use the received parameters to call WaveProSimu class to simulate the shockwave propagation, e.x. the propagation attenuation value of the shockwave overpressure in ventilation network; meanwhile, the StrFailureCal class is invoked to determine the failure mode of the ventilation structure; finally, the WaveProVisu class is invoked to realize the visual display of the shockwave overpressure propagation process.In the same way, the time sequence flow of the disaster gas simulation spreading has the following processes: first invoke the ExpPos, ExpTem, and ExpCon classes to determine the related parameters such as the location of explosion source, the maximum temperature of explosion source, and the concentration of explosion disaster gas; then transfer these parameters through the GetExpPos, GetExpTem, and GetExpCon methods; use the received parameters and the failure condition of the previously determined ventilation structure and invoke the UpdateNet method to update the mine ventilation network structure; then, based on it, the GasSprSimu class is invoked to simulate the spread of disaster gas in the ventilation network; finally, invoke the GasSprVisu class to visualize the disaster gas spread process.


### 3.3. Visualization System Function Module

The simulation system of the gas explosion disaster process is based on OpenGL three-dimensional technology and database technology. With the simulation of the gas explosion disaster process as the core, it realizes a comprehensive platform integrating calculation analysis, data management, and visual display and serves mine disaster prediction and prevention and formulates the corresponding decision-making measures. Its main system functions and technical architecture are shown in Figure6.Figure 6
System function and technical architecture.The main functional modules of the system include the data management module, 3D visualization module, and analysis decision module. The data management module mainly uses the relation-object-oriented method to store spatial data and attribute data in an integrated manner. It is mainly responsible for data organization, storage, and management and can realize data addition, deletion, modification, query, and export. The 3D visualization module is the basic module of the whole system. The module is used to dynamically display the simulation data and its simulation results on the 3D visualization platform by means of point, line, surface, and so forth and visually interact, which can realize the scaling, rotation, translation, interactive editing, picking, and positioning of graphics. The analysis decision module is the core module of the software. According to the network solution model, the simulation of the gas explosion disaster process is achieved on the basis of the collected network basic data and gas explosion disaster data. The module can mainly achieve the analysis of the influential range of shockwave overpressure, the determination of structure damage, gas diffusion simulation of explosion disasters, wind flow simulation of disaster process, the selection of optimal escape and disaster relief routes, and so forth.


## 3.1. Establishment of Simulation Model Class for the Disaster Process

According to the aforementioned gas explosion disaster field model, the simulation model of the disaster process is generated. The most important classes are wind roadway class, wind network class, roadway condition class, explosion source class, structure class, gas explosion information class, shockwave overpressure propagation class, and spreading of gas explosion disaster gas class. The relationships between classes and classes are shown in Figure3. The figure shows the basic data structure of the model classes involved in the gas explosion simulation process. The basic data members such as the space location, network topology relationship, and attribute information are listed in the given data structure.Figure 3
Visual model class diagram of gas explosion disaster process.


## 3.2. Time Sequence Flow of Disaster Process Simulation

Based on the establishment of the visualization model class of gas explosion disaster process and the analysis structure of time sequence characteristics of field variable, the timing flow diagrams of gas explosion disaster process simulation are given. Figures4 and 5 are the simulation timing flow diagram of shockwave overpressure propagation and disaster gas simulation spreading.Figure 4
Simulation time sequence diagram of shockwave propagation of gas explosion.Figure 5
Simulation time sequence diagram of gas explosion in gas explosion.The specialized time sequence flow of gas explosion shockwave propagation simulation is as follows: first invoke ExpPos class and ExpPre class to determine the location in the explosion source, the highest pressure of the explosion source, and other related parameters and then transfer these parameters through the GetExpPos and GetExpPre methods. Then use the received parameters to call WaveProSimu class to simulate the shockwave propagation, e.x. the propagation attenuation value of the shockwave overpressure in ventilation network; meanwhile, the StrFailureCal class is invoked to determine the failure mode of the ventilation structure; finally, the WaveProVisu class is invoked to realize the visual display of the shockwave overpressure propagation process.In the same way, the time sequence flow of the disaster gas simulation spreading has the following processes: first invoke the ExpPos, ExpTem, and ExpCon classes to determine the related parameters such as the location of explosion source, the maximum temperature of explosion source, and the concentration of explosion disaster gas; then transfer these parameters through the GetExpPos, GetExpTem, and GetExpCon methods; use the received parameters and the failure condition of the previously determined ventilation structure and invoke the UpdateNet method to update the mine ventilation network structure; then, based on it, the GasSprSimu class is invoked to simulate the spread of disaster gas in the ventilation network; finally, invoke the GasSprVisu class to visualize the disaster gas spread process.


## 3.3. Visualization System Function Module

The simulation system of the gas explosion disaster process is based on OpenGL three-dimensional technology and database technology. With the simulation of the gas explosion disaster process as the core, it realizes a comprehensive platform integrating calculation analysis, data management, and visual display and serves mine disaster prediction and prevention and formulates the corresponding decision-making measures. Its main system functions and technical architecture are shown in Figure6.Figure 6
System function and technical architecture.The main functional modules of the system include the data management module, 3D visualization module, and analysis decision module. The data management module mainly uses the relation-object-oriented method to store spatial data and attribute data in an integrated manner. It is mainly responsible for data organization, storage, and management and can realize data addition, deletion, modification, query, and export. The 3D visualization module is the basic module of the whole system. The module is used to dynamically display the simulation data and its simulation results on the 3D visualization platform by means of point, line, surface, and so forth and visually interact, which can realize the scaling, rotation, translation, interactive editing, picking, and positioning of graphics. The analysis decision module is the core module of the software. According to the network solution model, the simulation of the gas explosion disaster process is achieved on the basis of the collected network basic data and gas explosion disaster data. The module can mainly achieve the analysis of the influential range of shockwave overpressure, the determination of structure damage, gas diffusion simulation of explosion disasters, wind flow simulation of disaster process, the selection of optimal escape and disaster relief routes, and so forth.


## 4. A Case Study Analysis

In the study, taking the gas explosion accident that occurred in a coal mine in China as an example, the 3D simulation of shockwave propagation and disaster gas dynamic propagation in the gas explosion was conducted to identify the simulation influence.


### 4.1. Accident Profiles

The minefield is a monoclinic structure, bounded by the coal seam outcrop in the east, the lower SHIJIE minefield in the north, and the coal seam boundary in the south and west. It is natural sedimentation, heading northeast, with the length of 5.5 km, northwest inclination, the width of 3.7 km, and the area of 20.4 km2. The shallow section of the minefield has a dip angle of 20 to 30° and a depth of 6 to 9°. The minefield has five coal seams of 1, 2, 3, 4, and 5, where coal seams 4 and 5 are the main mining layer with a thickness of 8.12 m. There are four mining areas in the whole mine and an intensive mining mode on one side of the mine. The entire mines adopt the method of 8 forward 1 backward and the extraction ventilation mode in the mining area.Before the “11.28” gas explosion accident, the belts in the no. 4 mining area were drilled down to 180 m, and excavation 415 was prepared to transport the working surface. Since the 1# contact roadway was closed to form a blind alley, a large amount of gas that accumulated in the 1# contact roadway was poured into the rear area of the tail beam of lower corner bracket in 415 working face, where the roof explosion occurred when the top shot was placed. The gas explosion disaster wave involved all areas of the four mining areas below the main platform, including 2 working face systems (415 working face and 416 working face), 1 driving working face (417 return roadway), and 3 extension roadways (4 downward railways, 4 downward belt-ways, and 4 main return airways), the whole chambers in the no. 4 mining area, drilling areas, and contact roadways, as well as ANZIGOU ventilation shaft in the no. 4 mining area, and so forth. The entire process threatened 293 people, of which there were 166 deaths and 45 injuries, as shown in Figures7 and 8 showing the layout of the ventilation facilities before the explosion.Figure 7
Distribution of gas explosion casualties in mine.Figure 8
Position map of mine ventilation facilities.The mine was equipped with KJ66 monitoring and with 50 sensors. Based on the mine monitoring system and the sensor system and combined with the scene investigation after the gas explosion, the main damage evidence was obtained, as shown in Figure9, and the gas explosion arrival position, the shockwave propagation range, and the underground personnel injury situation were judged.Figure 9
Main damage evidence in field investigation.
(a)(b)Based on the accident investigation result, the explosion source was located at the rear side of the corner stern beam at the lower surface of the 415 working faces, as is shown in Figure10. Because the no. 1 connecting roadway has been sealed before the gas explosion accident, the accumulated high concentration gas flows into the lower corner and goaf through the intersection and surrounding fissure [30]. It is estimated that a total of 3007 m3 of gas was emitted in 87 hours before the explosion in the no. 1 connecting roadway. According to the gas concentration extracted from the no. 1–no. 6 gas drainage holes in goaf 415, which was 8.4%–50%, it can be seen that the gas concentration within the scope of the goaf is within the explosion limit. It is estimated that the total amount of gas involved in a gas explosion in the goaf is about 4813 m3. The total amount of gas explosion at the explosion source is 7820 m3.Figure 10
Explosion source point and 415 working face.In general, the equivalent TNT is used as the basic parameter of gas explosion intensity, so the explosion intensity of gas explosion accident in the mine ism=0.945VCH4 = 7.39 t, which is equivalent to the explosion effect of 7.39 t TNT explosive.


### 4.2. Damage Analysis of Shock Wave to Structures

The damage judgment method of shock wave overpressure on ventilation facilities is shown in the literature [28]. Figure 8 is the location map of mine ventilation facilities, taking ventilation facilities damper F1 as an example of calculation.The overpressure value of shock waves propagating to damperF1 is 0.787 MPa. Assuming that the damper F1 is a wooden damper, the thickness of the rectangular plate in the middle of the damper is h = 0.1 M, width is b = 0.3 m, the spacing between damper columns is L = 1 m, the cross section of tunnel where the damper is located is s = 13.2 m2, the height of roadway is 2.4 m, and the coefficient of wall friction is β = 0.045.According to formulas (5–11)–(5–17) in [28], the calculated damage value of damper F1 caused by shock wave overpressure is obtained, as shown in Table 1. As the yield strength of wood is equal to 50 MPa, it can be seen from the judgment criterion that the damper is damaged by the shock wave. Using the same method, the damage to other air doors can be judged.Table 1
The calculation value of the failure of the shock wave overpressure onF1.
OverpressureΔp (MPa)Reflection pressureΔp1 (MPa)Maximum bending momentMmax (kN.m)Standard stressσmax (MPa)0.7874.0611.0152030


### 4.3. Analysis of the Scope of the Gas Explosion

Before the simulation, the 3D spatial geometric data of the mine and the basic data of ventilation networks need to be input, and the 3D models of the mine are generated by the simulation system. The 3D model scale of the mine is as follows: 116 coordinate points, 145 roadways, and 20 dampers and adjustment facilities.


#### 4.3.1. Simulation of Gas Explosion Shockwave Propagation Process

According to the simulation time sequence flow of gas explosion shockwave propagation, the visual simulation of gas explosion shockwave overpressure is conducted on the basis of a certain time step, and the dynamic changes of explosion shockwave overpressure in the ventilation network at different moments are obtained. Figure11 shows the propagation of the shockwave overpressure in 415 working faces at four different times at t = 0.15 s, 0.48 s, 1.25 s, and 3.52 s. In the figure, yellow indicates the roadways not affected by the explosion shockwave, red dots represent the location of the explosion source, and blue represents the roadways affected by the explosion shockwave overpressure.Figure 11
Visual simulation of overpressure propagation process of shockwave. (a) 0.15 s. (b) 0.48 s. (c) 1.25 s. (d) 3.52 s.
(a)(b)(c)(d)It can be seen from Figure11 that when t = 0.15 s, the explosion shockwave overpressure propagates to the 415 working face, 415 conveying roadway, and no. 1 contact roadway; when it reaches the no. 1 contact roadway air door F1, the overpressure is 0.787 MPa. The air door is damaged at this time; then, the shockwave is transmitted to the closed wall with the overpressure of 0.612 MPa, and the closed wall is also damaged. When t = 0.48 s, the explosion shockwave overpressure propagates to 417 air return way, four downward railway extensions, four downward belt-way extensions, 415 conveying roadway, and 415 high position roadway. Meanwhile, a total of four air doors are damaged. When t = 1.25 s, the explosion shockwave overpressure enters the 416 working face by the 416 conveying roadway, 416 air return way, and 416 grouting roadway. Meanwhile, it enters the total return airway extension, and a total of 13 air doors suffer from the damage. When t = 3.52 s, the explosion shockwave overpressure completely covers the 415 working face and the 416 working face, and the explosion shockwave propagates backward along the air inlet side to the four downward railway extensions and the four downward belt-way extensions, which also causes the damage of the partial air doors at the air inlet side. At this time, a total of 17 air doors suffer from the damage. Table 2 shows the time sequence of the damage to the air doors at four different times.Table 2
Damage to the air door at different times.
Time0.15 s0.48 s1.25 s3.52 sDamage of air doorF1F2, F4, F5, F6, F14F3, F7, F8, F11, F12, F15F9, F10, F13, F18, F19According to the simulation results, there are 98 nodes and 89 roadways in the entire ventilation network, and 85% of the structures are damaged by different degrees. The impact time of explosion shockwave overpressure is about 5 s, and the affected wind location mainly has 415 working face, 417 air return way working face, 416 preparatory working face, and 17 air doors.


#### 4.3.2. Simulation of Disaster Gas Spread Process

Before simulating the disastrous gas spread process, the structure of the mine ventilation network is analyzed. Figure12 is a ventilation network diagram before the gas explosion occurs, and Figure 12 is a ventilation network diagram after the gas explosion occurs. Contrasting Figures 12(a) and 12(b), the structure of the mine ventilation network has undergone major changes. The wind flows at multiple locations are short-circuited, in which air door F16 is destroyed, 407 working face is short-circuited, air door F7 is destroyed, and 416 working face area is short-circuited. After 416 working face and 407 working face are short-circuited, the wind flow is in a windless or breezy condition, which causes the gas accumulation on the working face. Before the explosion, the wind volume of 415 working face was 18.7 m3/s, the wind volume of 416 conveying roadway and grouting roadway was 25 m3/s and 10.7 m3/s, respectively, and the wind volume of 407 working face was 6 m3/s. After the explosion, the wind volume of 415 working face was about 10 m3/s, and the wind volume of 416 conveying roadway and grouting roadway is 0 m3/s and 0 m3/s, respectively, and the wind volume of 407 working face is 0 m3/s.Figure 12
Evolution of ventilation network before and after a disaster of the gas explosion. (a) Before the explosion. (b) After the explosion.Based on the ventilation network structure after the explosion and the time sequence process of gas explosion simulation, the disaster gas diffusion process is simulated according to a certain time step. Figure13 shows the propagation process of 415 working face explosion disaster gas at four different times at t = 20 s, 135 s, 312 s, and 500 s. In the figure, yellow indicates the roadways not affected by the gas spread of the explosion disaster, red dots represent the location of the explosion source, and blue represents the roadways affected by the gas spread of the explosion disaster. It can be seen from Figure 13 that when t = 20 s, the explosion shockwave propagation ends, and the disaster gas CO spreads with the wind flow from the no. 1 contact roadway and 415 working face to 415 grouting roadway and 415 air return way. When t = 135 s, the CO gas covers 415 air return way, 415 high position roadway, and 415 grouting roadway and begins to spread to the total return airway. When t = 312 s, the CO gas spreads to the downward belt-way extension and the downward railway extension, and some CO gas spreads to 416 along the 416 grouting roadway. When t = 500 s, the CO gas spreads to the main return airways and is transferred from the return air well to the ground.Figure 13
Visual simulation of gas spreading process of explosion disaster. (a) 20 s. (b) 135 s. (c) 312 s. (d) 500 s.
(a)(b)(c)(d)Figure14 shows the variation curve of flue gas concentration and temperature on the simulation node. The closer the distance to the explosion source, the more direct the fire smoke spread and the more close the shape of the monitoring value curve to the combustion change process model of the explosion source point. It can be seen from Figures 14(a) and 14(b) that, after the explosion (t = 98 s), the maximum concentration of flue gas at the explosion source point reaches 11.8%, and the maximum temperature of flue gas reaches 1150°C. In the process of smoke concentration and temperature spreading, the temperature distribution gradually attenuates along the direct smoke invasion area downstream of the airflow. From the analysis of the death toll of disasters, the main harm of mine after the explosion is the invasion of smoke. Table 3 shows the time when disaster gas spreads to each node.Figure 14
Time series change of flue gas concentration and temperature at each simulated monitoring point. (a) Simulation of flue gas concentration change. (b) Simulation of flue gas temperature change.
(a)(b)Table 3
Disaster gas spread to each node time.
Serial numberTime (s)Location nodeSerial numberTime (s)Location node1069520473, 7422063, 65626780345677312101, 103413568, 1148500111It is known from the simulation analysis of the disaster gas diffusion process that 30 nodes and 39 roadways in the entire ventilation network are affected. Gas explosion disaster gas spread time is about 500 s. Influenced wind locations mainly have 415 working face and 417 air return way working face. Additionally, airflow stagnation occurs in the entire 416 working face area and 406 working face area, and the wind volume of the 415 working face is reduced.It can be seen from the above analysis that the gas explosion shockwave overpressure covers the entire area of 415 working face and 416 working face, and the disaster gas covers 415 working face area and the air return way areas such as main return airway extension, downward belt-way extension, and downward railway extension. Additionally, as a result of the destruction of the ventilation facilities and the change of the structure of the ventilation network, the entire areas of 416 working face and 407 working face are in a breeze or no-wind condition for a long time. The shockwave overpressure propagation and disaster gas coverage obtained by the simulation are basically consistent with the distribution of the number of casualties in the explosion accident in Figure7; a total of 17 air doors were damaged in the simulation calculation, where there are some errors in the damage of the 20 air doors obtained after the accident. The major reason perhaps is that the air doors data that are evaluated for the calculation are different from the actual values when these are simulated.


## 4.1. Accident Profiles

The minefield is a monoclinic structure, bounded by the coal seam outcrop in the east, the lower SHIJIE minefield in the north, and the coal seam boundary in the south and west. It is natural sedimentation, heading northeast, with the length of 5.5 km, northwest inclination, the width of 3.7 km, and the area of 20.4 km2. The shallow section of the minefield has a dip angle of 20 to 30° and a depth of 6 to 9°. The minefield has five coal seams of 1, 2, 3, 4, and 5, where coal seams 4 and 5 are the main mining layer with a thickness of 8.12 m. There are four mining areas in the whole mine and an intensive mining mode on one side of the mine. The entire mines adopt the method of 8 forward 1 backward and the extraction ventilation mode in the mining area.Before the “11.28” gas explosion accident, the belts in the no. 4 mining area were drilled down to 180 m, and excavation 415 was prepared to transport the working surface. Since the 1# contact roadway was closed to form a blind alley, a large amount of gas that accumulated in the 1# contact roadway was poured into the rear area of the tail beam of lower corner bracket in 415 working face, where the roof explosion occurred when the top shot was placed. The gas explosion disaster wave involved all areas of the four mining areas below the main platform, including 2 working face systems (415 working face and 416 working face), 1 driving working face (417 return roadway), and 3 extension roadways (4 downward railways, 4 downward belt-ways, and 4 main return airways), the whole chambers in the no. 4 mining area, drilling areas, and contact roadways, as well as ANZIGOU ventilation shaft in the no. 4 mining area, and so forth. The entire process threatened 293 people, of which there were 166 deaths and 45 injuries, as shown in Figures7 and 8 showing the layout of the ventilation facilities before the explosion.Figure 7
Distribution of gas explosion casualties in mine.Figure 8
Position map of mine ventilation facilities.The mine was equipped with KJ66 monitoring and with 50 sensors. Based on the mine monitoring system and the sensor system and combined with the scene investigation after the gas explosion, the main damage evidence was obtained, as shown in Figure9, and the gas explosion arrival position, the shockwave propagation range, and the underground personnel injury situation were judged.Figure 9
Main damage evidence in field investigation.
(a)(b)Based on the accident investigation result, the explosion source was located at the rear side of the corner stern beam at the lower surface of the 415 working faces, as is shown in Figure10. Because the no. 1 connecting roadway has been sealed before the gas explosion accident, the accumulated high concentration gas flows into the lower corner and goaf through the intersection and surrounding fissure [30]. It is estimated that a total of 3007 m3 of gas was emitted in 87 hours before the explosion in the no. 1 connecting roadway. According to the gas concentration extracted from the no. 1–no. 6 gas drainage holes in goaf 415, which was 8.4%–50%, it can be seen that the gas concentration within the scope of the goaf is within the explosion limit. It is estimated that the total amount of gas involved in a gas explosion in the goaf is about 4813 m3. The total amount of gas explosion at the explosion source is 7820 m3.Figure 10
Explosion source point and 415 working face.In general, the equivalent TNT is used as the basic parameter of gas explosion intensity, so the explosion intensity of gas explosion accident in the mine ism=0.945VCH4 = 7.39 t, which is equivalent to the explosion effect of 7.39 t TNT explosive.


## 4.2. Damage Analysis of Shock Wave to Structures

The damage judgment method of shock wave overpressure on ventilation facilities is shown in the literature [28]. Figure 8 is the location map of mine ventilation facilities, taking ventilation facilities damper F1 as an example of calculation.The overpressure value of shock waves propagating to damperF1 is 0.787 MPa. Assuming that the damper F1 is a wooden damper, the thickness of the rectangular plate in the middle of the damper is h = 0.1 M, width is b = 0.3 m, the spacing between damper columns is L = 1 m, the cross section of tunnel where the damper is located is s = 13.2 m2, the height of roadway is 2.4 m, and the coefficient of wall friction is β = 0.045.According to formulas (5–11)–(5–17) in [28], the calculated damage value of damper F1 caused by shock wave overpressure is obtained, as shown in Table 1. As the yield strength of wood is equal to 50 MPa, it can be seen from the judgment criterion that the damper is damaged by the shock wave. Using the same method, the damage to other air doors can be judged.Table 1
The calculation value of the failure of the shock wave overpressure onF1.
OverpressureΔp (MPa)Reflection pressureΔp1 (MPa)Maximum bending momentMmax (kN.m)Standard stressσmax (MPa)0.7874.0611.0152030


## 4.3. Analysis of the Scope of the Gas Explosion

Before the simulation, the 3D spatial geometric data of the mine and the basic data of ventilation networks need to be input, and the 3D models of the mine are generated by the simulation system. The 3D model scale of the mine is as follows: 116 coordinate points, 145 roadways, and 20 dampers and adjustment facilities.


### 4.3.1. Simulation of Gas Explosion Shockwave Propagation Process

According to the simulation time sequence flow of gas explosion shockwave propagation, the visual simulation of gas explosion shockwave overpressure is conducted on the basis of a certain time step, and the dynamic changes of explosion shockwave overpressure in the ventilation network at different moments are obtained. Figure11 shows the propagation of the shockwave overpressure in 415 working faces at four different times at t = 0.15 s, 0.48 s, 1.25 s, and 3.52 s. In the figure, yellow indicates the roadways not affected by the explosion shockwave, red dots represent the location of the explosion source, and blue represents the roadways affected by the explosion shockwave overpressure.Figure 11
Visual simulation of overpressure propagation process of shockwave. (a) 0.15 s. (b) 0.48 s. (c) 1.25 s. (d) 3.52 s.
(a)(b)(c)(d)It can be seen from Figure11 that when t = 0.15 s, the explosion shockwave overpressure propagates to the 415 working face, 415 conveying roadway, and no. 1 contact roadway; when it reaches the no. 1 contact roadway air door F1, the overpressure is 0.787 MPa. The air door is damaged at this time; then, the shockwave is transmitted to the closed wall with the overpressure of 0.612 MPa, and the closed wall is also damaged. When t = 0.48 s, the explosion shockwave overpressure propagates to 417 air return way, four downward railway extensions, four downward belt-way extensions, 415 conveying roadway, and 415 high position roadway. Meanwhile, a total of four air doors are damaged. When t = 1.25 s, the explosion shockwave overpressure enters the 416 working face by the 416 conveying roadway, 416 air return way, and 416 grouting roadway. Meanwhile, it enters the total return airway extension, and a total of 13 air doors suffer from the damage. When t = 3.52 s, the explosion shockwave overpressure completely covers the 415 working face and the 416 working face, and the explosion shockwave propagates backward along the air inlet side to the four downward railway extensions and the four downward belt-way extensions, which also causes the damage of the partial air doors at the air inlet side. At this time, a total of 17 air doors suffer from the damage. Table 2 shows the time sequence of the damage to the air doors at four different times.Table 2
Damage to the air door at different times.
Time0.15 s0.48 s1.25 s3.52 sDamage of air doorF1F2, F4, F5, F6, F14F3, F7, F8, F11, F12, F15F9, F10, F13, F18, F19According to the simulation results, there are 98 nodes and 89 roadways in the entire ventilation network, and 85% of the structures are damaged by different degrees. The impact time of explosion shockwave overpressure is about 5 s, and the affected wind location mainly has 415 working face, 417 air return way working face, 416 preparatory working face, and 17 air doors.


### 4.3.2. Simulation of Disaster Gas Spread Process

Before simulating the disastrous gas spread process, the structure of the mine ventilation network is analyzed. Figure12 is a ventilation network diagram before the gas explosion occurs, and Figure 12 is a ventilation network diagram after the gas explosion occurs. Contrasting Figures 12(a) and 12(b), the structure of the mine ventilation network has undergone major changes. The wind flows at multiple locations are short-circuited, in which air door F16 is destroyed, 407 working face is short-circuited, air door F7 is destroyed, and 416 working face area is short-circuited. After 416 working face and 407 working face are short-circuited, the wind flow is in a windless or breezy condition, which causes the gas accumulation on the working face. Before the explosion, the wind volume of 415 working face was 18.7 m3/s, the wind volume of 416 conveying roadway and grouting roadway was 25 m3/s and 10.7 m3/s, respectively, and the wind volume of 407 working face was 6 m3/s. After the explosion, the wind volume of 415 working face was about 10 m3/s, and the wind volume of 416 conveying roadway and grouting roadway is 0 m3/s and 0 m3/s, respectively, and the wind volume of 407 working face is 0 m3/s.Figure 12
Evolution of ventilation network before and after a disaster of the gas explosion. (a) Before the explosion. (b) After the explosion.Based on the ventilation network structure after the explosion and the time sequence process of gas explosion simulation, the disaster gas diffusion process is simulated according to a certain time step. Figure13 shows the propagation process of 415 working face explosion disaster gas at four different times at t = 20 s, 135 s, 312 s, and 500 s. In the figure, yellow indicates the roadways not affected by the gas spread of the explosion disaster, red dots represent the location of the explosion source, and blue represents the roadways affected by the gas spread of the explosion disaster. It can be seen from Figure 13 that when t = 20 s, the explosion shockwave propagation ends, and the disaster gas CO spreads with the wind flow from the no. 1 contact roadway and 415 working face to 415 grouting roadway and 415 air return way. When t = 135 s, the CO gas covers 415 air return way, 415 high position roadway, and 415 grouting roadway and begins to spread to the total return airway. When t = 312 s, the CO gas spreads to the downward belt-way extension and the downward railway extension, and some CO gas spreads to 416 along the 416 grouting roadway. When t = 500 s, the CO gas spreads to the main return airways and is transferred from the return air well to the ground.Figure 13
Visual simulation of gas spreading process of explosion disaster. (a) 20 s. (b) 135 s. (c) 312 s. (d) 500 s.
(a)(b)(c)(d)Figure14 shows the variation curve of flue gas concentration and temperature on the simulation node. The closer the distance to the explosion source, the more direct the fire smoke spread and the more close the shape of the monitoring value curve to the combustion change process model of the explosion source point. It can be seen from Figures 14(a) and 14(b) that, after the explosion (t = 98 s), the maximum concentration of flue gas at the explosion source point reaches 11.8%, and the maximum temperature of flue gas reaches 1150°C. In the process of smoke concentration and temperature spreading, the temperature distribution gradually attenuates along the direct smoke invasion area downstream of the airflow. From the analysis of the death toll of disasters, the main harm of mine after the explosion is the invasion of smoke. Table 3 shows the time when disaster gas spreads to each node.Figure 14
Time series change of flue gas concentration and temperature at each simulated monitoring point. (a) Simulation of flue gas concentration change. (b) Simulation of flue gas temperature change.
(a)(b)Table 3
Disaster gas spread to each node time.
Serial numberTime (s)Location nodeSerial numberTime (s)Location node1069520473, 7422063, 65626780345677312101, 103413568, 1148500111It is known from the simulation analysis of the disaster gas diffusion process that 30 nodes and 39 roadways in the entire ventilation network are affected. Gas explosion disaster gas spread time is about 500 s. Influenced wind locations mainly have 415 working face and 417 air return way working face. Additionally, airflow stagnation occurs in the entire 416 working face area and 406 working face area, and the wind volume of the 415 working face is reduced.It can be seen from the above analysis that the gas explosion shockwave overpressure covers the entire area of 415 working face and 416 working face, and the disaster gas covers 415 working face area and the air return way areas such as main return airway extension, downward belt-way extension, and downward railway extension. Additionally, as a result of the destruction of the ventilation facilities and the change of the structure of the ventilation network, the entire areas of 416 working face and 407 working face are in a breeze or no-wind condition for a long time. The shockwave overpressure propagation and disaster gas coverage obtained by the simulation are basically consistent with the distribution of the number of casualties in the explosion accident in Figure7; a total of 17 air doors were damaged in the simulation calculation, where there are some errors in the damage of the 20 air doors obtained after the accident. The major reason perhaps is that the air doors data that are evaluated for the calculation are different from the actual values when these are simulated.


## 4.3.1. Simulation of Gas Explosion Shockwave Propagation Process

According to the simulation time sequence flow of gas explosion shockwave propagation, the visual simulation of gas explosion shockwave overpressure is conducted on the basis of a certain time step, and the dynamic changes of explosion shockwave overpressure in the ventilation network at different moments are obtained. Figure11 shows the propagation of the shockwave overpressure in 415 working faces at four different times at t = 0.15 s, 0.48 s, 1.25 s, and 3.52 s. In the figure, yellow indicates the roadways not affected by the explosion shockwave, red dots represent the location of the explosion source, and blue represents the roadways affected by the explosion shockwave overpressure.Figure 11
Visual simulation of overpressure propagation process of shockwave. (a) 0.15 s. (b) 0.48 s. (c) 1.25 s. (d) 3.52 s.
(a)(b)(c)(d)It can be seen from Figure11 that when t = 0.15 s, the explosion shockwave overpressure propagates to the 415 working face, 415 conveying roadway, and no. 1 contact roadway; when it reaches the no. 1 contact roadway air door F1, the overpressure is 0.787 MPa. The air door is damaged at this time; then, the shockwave is transmitted to the closed wall with the overpressure of 0.612 MPa, and the closed wall is also damaged. When t = 0.48 s, the explosion shockwave overpressure propagates to 417 air return way, four downward railway extensions, four downward belt-way extensions, 415 conveying roadway, and 415 high position roadway. Meanwhile, a total of four air doors are damaged. When t = 1.25 s, the explosion shockwave overpressure enters the 416 working face by the 416 conveying roadway, 416 air return way, and 416 grouting roadway. Meanwhile, it enters the total return airway extension, and a total of 13 air doors suffer from the damage. When t = 3.52 s, the explosion shockwave overpressure completely covers the 415 working face and the 416 working face, and the explosion shockwave propagates backward along the air inlet side to the four downward railway extensions and the four downward belt-way extensions, which also causes the damage of the partial air doors at the air inlet side. At this time, a total of 17 air doors suffer from the damage. Table 2 shows the time sequence of the damage to the air doors at four different times.Table 2
Damage to the air door at different times.
Time0.15 s0.48 s1.25 s3.52 sDamage of air doorF1F2, F4, F5, F6, F14F3, F7, F8, F11, F12, F15F9, F10, F13, F18, F19According to the simulation results, there are 98 nodes and 89 roadways in the entire ventilation network, and 85% of the structures are damaged by different degrees. The impact time of explosion shockwave overpressure is about 5 s, and the affected wind location mainly has 415 working face, 417 air return way working face, 416 preparatory working face, and 17 air doors.


## 4.3.2. Simulation of Disaster Gas Spread Process

Before simulating the disastrous gas spread process, the structure of the mine ventilation network is analyzed. Figure12 is a ventilation network diagram before the gas explosion occurs, and Figure 12 is a ventilation network diagram after the gas explosion occurs. Contrasting Figures 12(a) and 12(b), the structure of the mine ventilation network has undergone major changes. The wind flows at multiple locations are short-circuited, in which air door F16 is destroyed, 407 working face is short-circuited, air door F7 is destroyed, and 416 working face area is short-circuited. After 416 working face and 407 working face are short-circuited, the wind flow is in a windless or breezy condition, which causes the gas accumulation on the working face. Before the explosion, the wind volume of 415 working face was 18.7 m3/s, the wind volume of 416 conveying roadway and grouting roadway was 25 m3/s and 10.7 m3/s, respectively, and the wind volume of 407 working face was 6 m3/s. After the explosion, the wind volume of 415 working face was about 10 m3/s, and the wind volume of 416 conveying roadway and grouting roadway is 0 m3/s and 0 m3/s, respectively, and the wind volume of 407 working face is 0 m3/s.Figure 12
Evolution of ventilation network before and after a disaster of the gas explosion. (a) Before the explosion. (b) After the explosion.Based on the ventilation network structure after the explosion and the time sequence process of gas explosion simulation, the disaster gas diffusion process is simulated according to a certain time step. Figure13 shows the propagation process of 415 working face explosion disaster gas at four different times at t = 20 s, 135 s, 312 s, and 500 s. In the figure, yellow indicates the roadways not affected by the gas spread of the explosion disaster, red dots represent the location of the explosion source, and blue represents the roadways affected by the gas spread of the explosion disaster. It can be seen from Figure 13 that when t = 20 s, the explosion shockwave propagation ends, and the disaster gas CO spreads with the wind flow from the no. 1 contact roadway and 415 working face to 415 grouting roadway and 415 air return way. When t = 135 s, the CO gas covers 415 air return way, 415 high position roadway, and 415 grouting roadway and begins to spread to the total return airway. When t = 312 s, the CO gas spreads to the downward belt-way extension and the downward railway extension, and some CO gas spreads to 416 along the 416 grouting roadway. When t = 500 s, the CO gas spreads to the main return airways and is transferred from the return air well to the ground.Figure 13
Visual simulation of gas spreading process of explosion disaster. (a) 20 s. (b) 135 s. (c) 312 s. (d) 500 s.
(a)(b)(c)(d)Figure14 shows the variation curve of flue gas concentration and temperature on the simulation node. The closer the distance to the explosion source, the more direct the fire smoke spread and the more close the shape of the monitoring value curve to the combustion change process model of the explosion source point. It can be seen from Figures 14(a) and 14(b) that, after the explosion (t = 98 s), the maximum concentration of flue gas at the explosion source point reaches 11.8%, and the maximum temperature of flue gas reaches 1150°C. In the process of smoke concentration and temperature spreading, the temperature distribution gradually attenuates along the direct smoke invasion area downstream of the airflow. From the analysis of the death toll of disasters, the main harm of mine after the explosion is the invasion of smoke. Table 3 shows the time when disaster gas spreads to each node.Figure 14
Time series change of flue gas concentration and temperature at each simulated monitoring point. (a) Simulation of flue gas concentration change. (b) Simulation of flue gas temperature change.
(a)(b)Table 3
Disaster gas spread to each node time.
Serial numberTime (s)Location nodeSerial numberTime (s)Location node1069520473, 7422063, 65626780345677312101, 103413568, 1148500111It is known from the simulation analysis of the disaster gas diffusion process that 30 nodes and 39 roadways in the entire ventilation network are affected. Gas explosion disaster gas spread time is about 500 s. Influenced wind locations mainly have 415 working face and 417 air return way working face. Additionally, airflow stagnation occurs in the entire 416 working face area and 406 working face area, and the wind volume of the 415 working face is reduced.It can be seen from the above analysis that the gas explosion shockwave overpressure covers the entire area of 415 working face and 416 working face, and the disaster gas covers 415 working face area and the air return way areas such as main return airway extension, downward belt-way extension, and downward railway extension. Additionally, as a result of the destruction of the ventilation facilities and the change of the structure of the ventilation network, the entire areas of 416 working face and 407 working face are in a breeze or no-wind condition for a long time. The shockwave overpressure propagation and disaster gas coverage obtained by the simulation are basically consistent with the distribution of the number of casualties in the explosion accident in Figure7; a total of 17 air doors were damaged in the simulation calculation, where there are some errors in the damage of the 20 air doors obtained after the accident. The major reason perhaps is that the air doors data that are evaluated for the calculation are different from the actual values when these are simulated.


## 5. Conclusions

(1)
Based upon the theory of mine compound field, the disaster field of mine gas explosion facing the development and change of mine gas explosion and its influence on the underground environment in time and space is put forward. The parameters of the field variables include the roadway passage for underground ventilation and disaster spread, spatial geometric field variables of physical facilities such as fan and air ventilation, the physical attribute field variables of underground ventilation such as wind volume, wind resistance, wind speed, and resistance coefficient, and the field variables of meteorological parameters such as wind flow temperature, explosion shockwave pressure, disaster gas concentration, and fire wind pressure.(2)
By analyzing the mine gas explosion disaster field’s characteristics of spatial geometry, attribute distribution, running status, vector direction, and time sequence, an integrated model structure that is suitable to the mine gas explosion disaster field is established, and the models of the geometric attribute, physical attribute, and meteorological parameters of the explosion disaster process are constructed, and the expression of the logical structure and topological structure of the field information including the mine roadway engineering, ventilation network, and explosive disaster are realized.(3)
By establishing the model class of mine gas explosion disaster field and analyzing the time sequence of field, the simulation system of mine gas explosion disaster process is developed by using OpenGL visualization method. The visual expression of process parameters (including air volume, overpressure, and disaster gas concentration) of gas explosion shockwave propagation and disaster gas dynamic propagation is realized. Through inversion of gas explosion accident cases, the disaster process of underground gas explosion can be truly reproduced, which can provide support for coal mine disaster relief and avoidance work.(4)
Taking a coal mine gas explosion accident in China as the research object, this paper analyzes the explosion source intensity, the destruction of ventilation structures, and the dynamic evolution process of ventilation network structure. The propagation of shock wave overpressure in the network and the gas spreading process of gas explosion disaster are simulated by using the developed gas explosion disaster process simulation system. The correctness and practicability of the simulation system are verified.



---

*Source: 6644263-2021-01-15.xml*