Numerical Investigation of Network-Based Shock Wave Propagation of Designated Methane Explosion Source in Subsurface Mine Ventilation System Using 1D FDM Code

Abstract

In coal mining operations, methane explosions constitute a severe safety risk, endangering miners’ lives and causing substantial economic losses, which, in turn, weaken the production efficiency and economic benefits of the mining industry and hinder the sustainable development of the industry. To address this challenge, this article explores the application of decoupling network-based methods in methane explosion simulation, aiming to optimize underground mine ventilation system design through scientific means and enhance safety protection for miners. We used the one-dimensional finite difference method (FDM) software Flowmaster to simulate the propagation process of shock waves from a gas explosion source in complex underground tunnel networks, covering a wide range of scenarios from laboratory-scale parallel network samples to full-scale experimental mine settings. During the simulation, we traced the pressure loss in the propagation of the shock wave in detail, taking into account the effects of pipeline friction, shock losses caused by bends and obstacles, T-joint branching connections, and cross-sectional changes. The results of these two case studies were presented, leading to the following insights: (1) geometric variations within airway networks exert a relatively minor influence on overpressure; (2) the positioning of the vent positively contributes to attenuation effects; (3) rarefaction waves propagate over greater distances than compression waves; and (4) oscillatory phenomena were detected in the conduits connecting to the surface. This research introduces a computationally efficient method for predicting methane explosions in complex underground ventilation networks, offering reasonable engineering accuracy. These research results provide valuable references for the safe design of underground mine ventilation systems, which can help to create a safer and more efficient mining environment and effectively protect the lives of miners.

1. Introduction

Shock waves emanating from methane explosions directly impact ventilation facilities, leading to significantly reduced airflow and the subsequent accumulation of toxic and flammable gases, posing a considerable risk to miners’ safety and lives [1,2]. In extreme cases, shock waves can inflict bodily harm, resulting in fatalities [3]. Methane explosions and their resulting shock waves pose a serious threat to miners’ safety, impact mining productivity and economic efficiency, and hinder sustainable development [4]. A well-conceived design not only effectively mitigates the impact of methane explosions but also ensures the rapid restoration of ventilation in emergency situations, minimizing the accumulation of toxic gases. This, in turn, maximizes the protection of miners’ lives, fostering healthy, stable, and sustainable development within the mining industry.

Two key factors influencing shock propagation in underground ventilation systems are scaling and geometric change effects [5,6]. The scale effect pertains to the airway dimensions, influencing the pressure evolution characteristics of an explosion. The geometric change effect describes how shock waves can be modified as they traverse bends, obstacles, branches, and cross-sectional variations. These effects jointly determine a shock wave’s destructiveness and propagation distance [7]. Conducting full-scale experiments in underground spaces is hazardous and costly [8], necessitating numerical predictions. However, numerical models are not without uncertainties and require validation. Van Wingerden’s research [9] indicates that the scale effect is significantly correlated with the normalized flame propagation velocity, though this influence diminishes in turbulent scenarios. Catlin and Johnson proposed compensating for scale effects by enriching oxygen concentration [10], a theoretically valid approach for Reynolds numbers below 10,000, which are generally lower than those for practical cases. Zhang et al. [11] utilized the CFD code AutoReaGas to investigate scale effects across three scales (1:1, 1:10, 1:100) [11,12]. Building upon our previous work [13], which delved into scaling effects in detail, we used this understanding as the explosion source (initial and boundary conditions) for our current one-dimensional (1D) numerical study.

Research indicates that geometric similarity in explosion parameters is achieved when the length-to-diameter ratio exceeds 80, with explosion overpressure positively correlated with explosion diameter [14]. Beyond scaling effects, geometric variations within airways significantly impact blast wave migration, necessitating thorough study and quantification for subsurface ventilation system explosion descriptions. Typical geometric variations in longwall and room-and-pillar mining operations include curvatures, bifurcations, obstructions, and modifications in cross-sectional profiles [3]. For low-speed duct water flow, extensive experiments have addressed minor and shock losses due to factors such as duct expansions, contractions, valves, and elbows, though data are primarily confined to low Mach numbers and incompressible flow [15,16]. Limited experimental studies offer qualitative insights into how geometric changes affect compressible, discontinuous shock wave propagation [17,18,19]. Researchers assert that bends can either amplify or mitigate gas explosion overpressure, depending on the presence of fuel-filled spaces [20,21]. Branching and expansion generally dampen methane explosion flame propagation and overpressure [22]. Meanwhile, obstacles exhibit behavior similar to that of bends, enhancing overpressure in fuel-filled areas while hindering shock wave propagation [23].

To illustrate the influence of scaling and geometric variations on methane explosions within subsurface ventilation networks, a comprehensive three-dimensional (3D) explosion database is necessary, serving as the foundation for initializing, with bounding conditions in a one-dimensional (1D) network-centric shock wave propagation model, which is meticulously explored in this article. The handling of pressure losses stemming from friction and geometric alterations requires thorough discussion. Furthermore, this article presents two illustrative case studies that employ a decoupled methane explosion prediction methodology. These case studies demonstrate the applicability of the methodology in predicting methane explosions within both a representative parallel network model and a full-scale experimental mine facility. Both studies leverage numerical simulation technology to uncover insights into how underground mine ventilation system designs can be optimized for safe and sustainable development.

2. Numerical Modeling

In order to identify the unique characteristics of the two families within the methane explosion zone, we divided the domain into two distinct components: a driver section (characterized by field function analysis) and a blast wave section (defined by network function dynamics). To model the explosion within the driver section, we employed a 3D field-based numerical approach, leveraging the Large Eddy Simulation (LES) turbulent model and an adapted c-equation chemical reaction model. Comprehensive details on the utilized boundary and initial conditions, spatial–temporal resolution assessments, combustion model benchmarks, and validation procedures were extensively covered in the author’s prior research [13]. Furthermore, a comprehensive database of explosion sources, incorporating diverse scaling factors and equivalence ratios, was established.

This study focuses on the numerical simulation of a planar shock wave’s propagation within a blast wave segment, abstracting turbulence and chemical reactions for simplicity. Simplifying the geometric model to a one-dimensional (1D) representation allowed for significant reductions in computational time, while still enabling adherence to the requirements of ventilation engineering simulations within intricate networks. The previously developed explosion source database [13] served as a valuable resource, furnishing the 1D model with a broad range of initial and boundary conditions grounded in the foundational 3D field simulations. Additionally, the analysis accounted for geometric variations by adjusting parameters specific to the components within the one-dimensional network geometry.

For the numerical analysis of blast wave propagation, Flowmaster (Version number V7), a widely recognized commercial 1D Computational Fluid Dynamics (CFD) software, was selected as the computational platform for network-based simulations. This software employs the implicit finite difference method (FDM) and is well known for its use in modeling fluid dynamics in gas and water pipeline systems. Nevertheless, its performance in predicting blast wave propagation necessitated thorough examination, and its attenuation models necessitated customization to suit the specific research objectives.

To accommodate the requirements of Flowmaster, the formulation of the governing one-dimensional equations was tailored. These equations, grounded in the assumptions of one-dimensionality, inviscidity, and compressibility, comprehensively encapsulated the conservation principles of momentum, energy, and mass, as tailored for the blast wave propagation domain. This detailed formulation is referenced in [24], providing a comprehensive basis for the numerical simulations undertaken in this study.

$${\displaystyle \frac{\partial p}{\partial t}}+\rho V^2{\displaystyle \frac{\partial V}{\partial x}}+\rho a^2{\displaystyle \frac{\partial V}{\partial x}}={\displaystyle \frac{a^2}{C_{P}T}}\left[1+{\displaystyle \frac{T}{Z}}{({\displaystyle \frac{\partial Z}{\partial T}})}_T\right]{\displaystyle \frac{\mathsf{\Omega }+WV}{A}}$$

(1)

$${\displaystyle \frac{\partial V}{\partial t}}+V{\displaystyle \frac{\partial V}{\partial x}}+{\displaystyle \frac{1}{\rho }}{\displaystyle \frac{\partial p}{\partial x}}={\displaystyle \frac{W}{A\rho }}-gsin\theta$$

(2)

$${\displaystyle \frac{\partial T}{\partial t}}+V{\displaystyle \frac{\partial T}{\partial x}}+{\displaystyle \frac{a^2}{C_P}}\left[1+{\displaystyle \frac{T}{Z}}{({\displaystyle \frac{\partial Z}{\partial T}})}_p\right]{\displaystyle \frac{\partial V}{\partial x}}={\displaystyle \frac{a^2}{C_{P}p}}\left[1-{\displaystyle \frac{p}{Z}}{({\displaystyle \frac{\partial Z}{\partial p}})}_T\right]{\displaystyle \frac{\mathsf{\Omega }+WV}{A}}$$

(3)

$$W=f{\displaystyle \frac{A}{D}}\rho {\displaystyle \frac{V\left|V\right|}{2}}$$

(4)

Here, $V$ represents velocity in meters per second, while $x$ denotes axial distance measured in meters (m). The velocity of sound is given by $a$, also in m/s, and $C_p$ signifies specific heat, measured in joules per gram kelvin (J/(g*K)). $Z$ corresponds to gas compressibility, $A$ represents the flow area in square meters (m²), and $\mathsf{\Omega }$ denotes heat input per unit length. The gravitational acceleration is g and $\theta$ is the pipe’s angle to the horizontal, measured in degrees. Additionally, $f$ represents the Darcy Frictional Factor, which will be elaborated upon in Section 3.2. It is important to highlight that gravitational acceleration and the pipe’s bearing angle were included in the calculations of Equation (2). Nevertheless, these elements were not factored in during the creation of the explosion source database, especially in cases with horizontal or near-horizontal airways, where the effects induced by gravitational acceleration are considered insignificant. Additionally, the blast wave propagation was assumed to follow an adiabatic process, meaning no heat exchange with the environment took place. As a result, in the formulation of Equations (1) and (3), the heat input term, denoted by $\mathsf{\Omega }$, was omitted.

3. Pressure Losses

3.1. Pressure Losses

Based on one-dimensional Computational Fluid Dynamics (CFD) software Flowmaster, we modeled the spread of explosion waves in underground networks. During this process, pressure losses arose, categorized by Flowmaster into two main types: frictional losses and losses attributable to geometric variations. The Darcy Friction Factor ($f$) was incorporated into the analysis to account for frictional losses in the pipe flow. This dimensionless parameter took into consideration the geometric characteristics of each pipe segment, specifically its length and absolute roughness.

Regarding pressure losses caused by geometric changes, Flowmaster employs specific governing equations tailored to different network components. The attenuation factor data, previously validated and published in our earlier work [25], were directly incorporated into the relevant component equations. In subsequent sections, we will detail the equations for T-junctions, cross-sectional alterations, bends, and obstacles, facilitating a comprehensive understanding of their impact on blast wave propagation.

3.2. Frictional Loss

Shear heating along a wall causes a decrease in pressure known as friction loss. This process involves the conversion of a blast wave’s mechanical energy into heat, which is subsequently dissipated. The Darcy Friction Factor ($f$) in Equation (4) was a dynamic parameter that varied in accordance with the Reynolds number (Re), which differed across three distinct flow scenarios. The dynamic viscosity ($\mu$) of air at 25 °C was constant, equal to 1.983 × 10⁻⁵, and played a pivotal role in determining Re. Additionally, Colebrook–White’s approximation [26] offered a practical method for estimating f, further enhancing the accuracy of the analysis.

For laminar flow (Re < 2000):

$$f=f_l={\displaystyle \frac{64}{\mathrm{Re}}}$$

(5)

For transition zone (2000 < Re < 4000):

$$f=x{f}_t+\left(1-x\right)f_l; x={\displaystyle \frac{\mathrm{Re}-2000}{2000}}$$

(6)

For turbulent flow (Re > 4000):

$$f=f_t={\displaystyle \frac{0.25}{{\left[\log \left(\frac{k}{0.37D}+\frac{5.74}{{\mathrm{Re}}^{0.9}}\right)\right]}^2}}$$

(7)

where $f_l$ represents the laminar friction factor, $f_t$ denotes the turbulent friction factor, and $k$ symbolizes the absolute roughness, measured in millimeters (mm).

3.3. Bends and Obstacles

As blast waves propagate, there is also a dissipation of energy caused by a reduction in pressure due to geometric changes. In our prior research [25], we explored the attenuation effects of four common geometric variations in underground coal mines: changes in cross-sections, branches, obstacles, and bends. Notably, while bends are often overlooked in long-distance water pipeline networks modeled with Flowmaster, their influence on methane explosion attenuation cannot be disregarded. Consequently, we incorporated a Discrete Loss component into our analytical framework to accurately capture this effect.

Equation (8) shows that Flowmaster accounted for pressure losses caused by geometric changes using the $C_d$ Discrete Loss component. In order to improve the accuracy of the simulation, we tried to modify the $C_d$ coefficient, and attenuation factors were smoothly integrated into the equation.

$$\Delta p={\displaystyle \frac{8Lf}{{\rho }^2\pi d^5}}{\left({\displaystyle \frac{C_d·A_t·P_{t1}·\psi ·\sqrt{\frac{2}{R_s·T_{t1}·Z_1}}}{m_{t1}}}\right)}^2$$

(8)

where $T_{t1}$, $m_{t1}$, and $P_{t1}$ are the total temperature upstream of the Cd Discrete Loss, mass flow rate in kg/s, and total pressure in Pa, respectively. Equation (7) explains how the Darcy Friction Factor $f$ was calculated. $R_s$ is a gas constant, $Z_1$ is compressibility coefficient for upstream flow, and $\psi$ is the flow function, which could be calculated by $\psi ={\displaystyle \frac{\rho {\int }_V\phi \left(x_{SGS}\right)d{x}_{SGS}}{V\rho }}$. Three $V$ is the mesh size for the Finite Volume Method (FVM) used by ANSYS Fluent, $\rho$ is a filtered or time-averaged value of the product of density, and $\phi \left(x_{SGS}\right)$ is the Sub-Grade value. Such a relationship was maintained between the attenuation factor $\eta$ and the pressure loss $\delta p$ due to bending, and the calculation formula was $\delta p=P_{t1}\left(1-\eta \right)/\eta$. If $\eta$ was determined, the change in the discharge factor $C_d\prime$ due to bending could be reversed using the determined $\delta p$. The flow property in Discrete Loss was expressed by the total $C_d$. It must be pointed out that $C_d$ needed to be calculated by $\eta$ and $P_{t1}$, where $P_{t1}$ represents the peak overpressure input from the source of the explosion. The formula for calculating $\eta$ was as follows: $\eta =P_0/P_1$, where $P_0$ and $P_1$ represent the peak overpressures before and after the geometric change.

Based on experimental observations, the pressure preceding a bend is consistently comparable in magnitude to the initial peak overpressure. Consequently, obstacles were addressed similarly to bends in this study, with the primary distinction lying in their unique attenuation factors. Further elaboration on the inclusion of obstacles is omitted for brevity.

3.4. T-Branches

Blast wave propagation in the T-branch has been analyzed in detail in the literature [25]. The length of the three arms of the T-branches selected in this article was 0.35 m. Analogous to bends, the T-branch element (also termed T-joint) exhibited negligible influence on lengthy water pipelines modeled within the software. Consequently, when simulating the propagation of blast waves, it became imperative to adjust the attenuation coefficient to accurately reflect its impact.

Flowmaster calculated the pressure loss in T-branches through the following equation:

$$\Delta p={\displaystyle \frac{C_{\mathrm{Re}}K\stackrel{·}{m_c}\left|\stackrel{·}{m_c}\right|}{2\rho A_c}}$$

(9)

where $C_{Re}$ represents the correction for the Reynolds number, while $K$ denotes the loss coefficient for one of the arms and $\dot{m_c}$ is the mass flow rate per unit area in kg/(m*s). The total pressure loss $\Delta p^{\prime }$ for an arm of a T-branch, accounting for the additional pressure loss due to geometric change, could be expressed as $\Delta p^{\prime }=\Delta p+P_{t1}\left(1-\eta \right)/\eta V$.

Attenuation factors were predicted for T-branches with flows originating from both the main arm and the branch arm. $\eta$ and $P_{t1}$ were used to calculate the loss coefficient $K$.

3.5. Cross-Sectional Change

According to Bernoulli’s principle, when a duct enlarges, the blast wave overpressure increases while the velocity decreases. This transformation sees kinetic energy converted into volumetric energy, thereby intensifying the internal overpressure. To encapsulate this phenomenon, Flowmaster utilizes a Transition Component, the behavior of which is dictated by the subsequent equation:

$${\displaystyle \frac{P_{t1}}{P_{t2}}}=1-K\left(1-{\displaystyle \frac{P_{s2}}{P_{t2}}}\right)$$

(10)

where $P_{s2}$ denotes the static pressure downstream of the expanded region, prior to its subsequent contraction. Notably, the chosen expand–contract configuration exerted minimal influence on the propagation of a blast wave. The primary consideration revolved around the heightened overpressure within the expanded zone, stemming from reflections off the duct walls. This interaction resulted in a pronounced elevation of overpressure within the expanded area.

4. Numerical Case Studies

4.1. Parallel Network

4.1.1. Problem Statement

A one-dimensional Computational Fluid Dynamics (CFD) code, Flowmaster, was used to model a methane explosion in a sample parallel airway. Figure 1a illustrates the network schematic.

A pressure source was used to simulate the explosion source on the left side of the network. In this experimental setup, pipes were standardized at 1 m in length, with exceptions being C8 at 1.31 m and C13 at 1.84 m. All components within the network model featured a diameter of 0.08 m. Near the left-dead end, we simulated a methane explosion, resulting in a blast wave traversing the duct, encountering a T-branch, and subsequently dividing into two parallel ducts traversing through bends. The bending angles of the C5-C8-C3 series were set at 130° and 50°, respectively. The difference was that the C2-C9-C4 branch had 90° bends throughout. $\eta =P_0/P_1$ was used to calculate the attenuation coefficient for each bend [27]. At the T-branch, the two ducts converged and vented to the atmosphere, with a constant 1-bar pressure source placed on the right side.

Figure 1a illustrates the geometric model of the Flowmaster parallel network, with component numbers denoted alongside each element as Ci.

Figure 1b presents the Flowmaster-specific geometric model, mirroring the network depicted in Figure 1a. Pipe lengths were set according to Figure 1a, with select Cd Discrete Loss components to replace the bends, a total of four, as described in Equation (8). The Cd values, calculated from the angles (one at 50°, one at 130°, and two at 90°), were utilized in the Discrete Loss components to detail the variations in the attenuation effects caused by the varying bend angles, as seen in Figure 1a.

4.1.2. Results

The monitoring of pressure histories (measured in bars) along pipe components with varying individual lengths (spanning from 0 m to 1 m in this context) was employed to provide a clear description of the propagation of the explosion wave. The peak pressures observed in each plot accurately coincided with the moment when blast waves traversed specific pipe segments. The waves experienced attenuation during their journey owing to frictional dissipation within the pipes and geometrical changes encountered along their path.

In the initial configuration of the 1D simulation, the overpressure history stemming from an 8% volumetric methane explosion within an experimental airway of dimensions of both width and height were 0.08 m and 4.25 m in length, retrieved from a pre-existing, comprehensive database. This database encompasses field-validated simulation outcomes across diverse scenarios. Leveraging the pressure evolution profile depicted in Figure 2, an assumed explosion origin was positioned in C16. The recorded overpressure maxima, situated 0.3 m from the ignition point, were noted as 0.664 bar, 0.599 bar, and 0.221 bar, respectively. Figure 3 offers a comprehensive examination of the pressure dynamics, synchronized with the respective local pipe lengths (all pipe segments were measured end-to-end) and temporal correlations across all pipe components. For clarity, the numbering and spatial arrangement of a selected pipe can be traced back to Figure 1b.

The intricate relationships between pressure (in bars) and time (in seconds) of different pipe components (the numbers from small to large are C2, C4, C5, C8, C9, C13, C14, and C17) identified in Figure 1 and the length of these pipe sections (in meters) are clearly recorded in Figure 3a–h. These illustrations vividly exhibit the progression of blast waves through individual pipe segments, characterized by notable pressure fluctuations over time and subtle pressure gradients along the pipe lengths. The dynamic visualizations underscore the wave’s traversal through each component, offering a precise portrayal of pressure dynamics within the system.

Figure 4 showcases a snapshot of the pressure distribution across all pipe components at 0.065 s, with pressure magnitudes denoted in bars (rounded to two decimal places) adjacent to each pipe. Upon detailed examination, a minute disparity in pressure was discernible between the upper and lower branches, exhibiting a comparable pressure decrement of roughly 0.1 bar for both. This marginal attenuation could be attributed to the deflagrative character of the simulated explosion and the relatively low peak overpressure (less than 0.7 bar) assigned to the explosion source, consequently leading to a negligible reduction in pressure during its propagation.

4.1.3. Discussion

Figure 3 underscores a pivotal observation: at approximately 0.065 s, the peak overpressures reached their zenith. Notably, C14, which was exposed to the open atmosphere, exhibited a distinctive rarefaction wave with a magnitude of approximately 0.91 bar, marking a slight depression of 0.09 bar below the ambient pressure level. Over time, the pressure within the pipe gradually harmonized with the surrounding environment, approaching ambient conditions. To further clarify, the maximum pressure levels achieved by each pipe component, quantified precisely in bars, were systematically documented in

Table 1. Example of maximum overpressure in a pipe assembly within a parallel network.

Component #	C17	C2	C5	C8	C9	C13	C4	C14
Peak Overpressure (bar)	0.20862	0.18832	0.21908	0.20935	0.13642	0.05882	0.06313	0
Arrival Time (s)	0.065	0.065	0.065	0.065	0.065	0.065	0.065	Na

# The number of the selected pipe.

, offering a comprehensive and concise summary of the pressure dynamics.

Table 1 highlights a distinct contrast: when compared to measurements in isolated ducts, the effect of geometric modifications on overpressure diminished. This phenomenon of attenuation arose due to the complex interplay between the varying geometries within the underground air passageways. In particular, downstream bends counteracted the attenuation triggered by upstream bends, as demonstrated in [28].

Disregarding temporal pressure fluctuations, Figure 3 portrays a marginal decrease in pressure as the pipe lengthened, with the most significant decline observed in C17, situated closest to the explosion epicenter. This pressure decrement was attributed to the augmentation of frictional forces, which intensified in accordance with pipe elongation.

Furthermore, the positioning of vents played a pivotal role in modulating the attenuation of geometric variations, particularly for pipes in close vicinity to vents. Notably, C13 and C4, positioned nearest to the vent, exhibited a marked reduction in overpressure, emphasizing the importance of vent location in influencing pressure dynamics.

4.2. Experimental Mine

4.2.1. Problem Statement

Using the decoupled simulation approach, a supplementary case study was conducted to model a methane explosion in a full-scale experimental mine. This primary mine, situated at 12350 Private Drive 7002 Rolla, MO, is employed by the Missouri University of Science and Technology for both educational and research endeavors (as illustrated in Figure 5). The network within the experimental mine exhibits significantly greater complexity in comparison to the sample parallel network model. Furthermore, its airway dimensions offer a much more accurate representation than those found in the lab-scale network model.

This study focused on an underground room-and-pillar mine located at the Mining and Nuclear Engineering Department of Missouri University of Science and Technology in Rolla, Missouri, known as the experimental mine. The simulation of the pressure distribution and blast wave arrival time during a methane explosion in this mine was conducted using the one-dimensional CFD tool, Flowmaster (Viasoft Limited, Southampton, UK).

The Flowmaster network model was established using a geometric model generated with VentSim, a ventilation simulation tool developed by Chasm Consulting in Queensland, Australia (the tool, accessible at http://www.ventsim.com/contact/, was accessed on 28 August 2024). The angles between the underground airways are indicated in the network schematic shown in Figure 5. The airway layout and dimensions for the geometric model are provided in the literature [29].

To ensure the convergence of calculations, reasonable simplifications were necessary due to the complexity of the network at the experimental mine. Figure 5 shows that the eastern region was linked to the central and western regions by a narrow passage measuring 2.9 m in width and 2.8 m in height. Consequently, the methane explosion designated in the mine’s central area had minimal impact on the eastern part, justifying its exclusion from the simulation for simplification purposes. Furthermore, when two airways formed an angle of more than 165° and were adjacent, they merged into a single airway, and the length was recorded as the sum of the two original airway lengths. Figure 6 depicts the geometric model that was simplified for Flowmaster, encompassing the western and middle regions.

Figure 6 shows that the airway circuit was used as the basis for the division of the middle and western regions, with a total of eight sub-regions, but region 6 was an exception, which connected the two regions via three series pipes. Airways belonged to two sub-regions if the airways were shared by two sub-regions. One airway was selected in each sub-region to assess the impact of the explosion: C59 in sub-region 1, C9 in sub-region 2, C24 in sub-region 3, C11 in sub-region 4, C31 in sub-region 5, C29 in sub-region 6, C43 in sub-region 7, and C50 in sub-region 8.

Between pipes C59 and C21 was an explosion simulation point. The corresponding bending angles were given next to each Cd Discrete Loss, and the selected airway and their regions were highlighted. Two additional measurements were taken from airways C53 and C2, which were directly connected to the portal and completely independent of the previously selected airways. Data were collected from blast waves exiting through two portals, each denoted by a “P” in red: portal number 1 (constant 1 bar pressure) in the top left corner and portal number 2 (also at 1 bar) at the bottom.

In Figure 6, the pressure source (denoted by the red “P” near region 1) indicates the explosion source’s location. This source had a methane volumetric concentration of 8% and dimensions of 8.5 m in length and 2.56 m by 2.56 m, aligning with the selected part of the experimental mine. The chosen methane explosion source is detailed in the pre-developed database [13] and will not be reiterated here. The inflating area was 2.56 m wide and high and 8.5 m long, with a hydraulic diameter similar to that of the C59 and C21 airways. A time step of 0.0013 s was selected to match the 3D simulations.

Discrete pressure loss components were used to model bends in the experimental mine’s network, similar to the sample parallel network. For joints with more than three arms, the double T-branch structure was employed to represent such configurations, accommodating smooth cross-sectional transitions by altering the dimensions of the pertaining airways. Since abrupt cross-sectional changes were absent, no transition component was required.

4.2.2. Results

The results of the simulation are depicted as three-dimensional surfaces, correlating with the local airway temporal dimensions, pressure, and length. Figure 7 illustrates the simulation outcomes of selected airways from eight distinct regions and their proximities to the portals.

The temporal and spatial pressure fluctuations of blast waves within specific components are depicted in Figure 7, showcasing the relationship between pressure, time, and local pipe length. Notably, pressure oscillations were observed when the simulation time was less than 0.5 s. Additionally, the pressure gradient along the pipe length, as observed in the pressure–length plane, exhibited a less significant decline compared to the parallel models presented in Figure 3. This was attributed to friction playing a lesser role in attenuating blast waves in larger ducts.

The network pressure distribution for 0.039 and 0.159 s is recorded in Figure 8. The detection of three blast wave crests underscored the impact of the explosion source on the entire mine complex.

The color map in Figure 8a on the left displays the upper pressure range in bar. The highlighted colors indicate the peak pressure for corresponding components (Cd discrete losses, pipes, and T-junctions).

The impact of the imaginary explosion on the middle and western regions of the network at 0.039 s is depicted in Figure 8, corresponding to the arrival of the first pressure peak. These regions are displayed separately (top and bottom). Small figures of different colors (red, yellow, and green) are placed next to the components to indicate their peak pressure. In Figure 8a, there was no significant pressure value change from region 1 to region 5, while the blast wave had a limited impact on regions 6 and 7 and nearly no effect on region 8.

Figure 8b shows the pressure distribution at 0.195 s as the third blast wave peaks arrived. At this point, four airways in region 7 reached 1.1 bar, indicating a stronger impact from the methane explosion compared to other times and suggesting that the blast wave reached regions 7 and 8 later than other areas. The pressure of about 1.1 bar in the western region remained significantly lower compared to the pressure of 1.4–1.5 bar in the middle region.

4.2.3. Discussion

This paper focuses on clarifying the propagation process of shock waves in a mine through numerical simulations, which helps optimize mine design and ensure personnel and property safety. Wang et al. [30] used FLACS software to simulate the characteristics of methane explosions in interconnected container systems. They found a 6.97% error between the simulation results and experimental data. Liang et al. [31] used GASFLOW-MPI for the numerical simulation of a full-scale pipeline network methane explosion, based on experimental results from the U.S. Lynn Lake Experimental Mine No. 501. They found that the error between the simulation and experimental results in complex networks was only 4.25%. This indicates that numerical simulation is a promising and reliable method for studying methane explosions in a full-scale mine network. In addressing the complexity of the network, this study considered factors such as the attenuation factor $\eta$ and Cd Discrete Loss to simulate the effects caused by pipeline bends, branches, and obstacles. Figure 7 shows that as the pipeline length increased, overpressure gradually decreased. Additionally, overpressure dropped at pipeline bends, which was attributed to shock wave reflection and interference and was consistent with previous studies. Prior research has shown that factors like pipeline length [30], pipeline bending angles and branches [32], and obstacles [31] can lead to pressure variations within a pipeline, which is related to shock wave reflections, interference, and vortices during propagation [31]. In Figure 8, the overpressure in the pipeline near the explosion point was significantly higher than that in pipelines farther away, which was likely due to the weakening effect of the pipeline branches [30], an important factor for interconnected systems.

A hypothetical methane explosion was simulated for the experimental mine network. The relationship between the length, pressure, and time of the pipe components selected from the eight regions is described in Figure 7, along with two pipes, shaft 1 and portal 2, which were directly connected to the surface. Figure 7 shows that similar pressure oscillation patterns occurred in regions 1 through 6. However, there were about 0.2 bar oscillations found in regions 7 and 8 (see Figure 7g,h). Negative impulses, around 0.1 bar, were larger than positive ones, indicating that rarefaction waves traveled farther than compressive waves. Negative pulses were also observed in airways near shaft 1 and portal 2 (see Figure 7i,j), likely due to rarefaction waves from the blast wave’s exhaust process near the surface.

Peak pressures at 0.039 s and 0.195 s are shown in Figure 8. At 0.195 s, the parallel-connected regions 2 through 5 were still strongly impacted by the explosion, as shown in Figure 8b. However, the explosion wave was greatly reduced in regions 6 through 8, connected in series. Three airways in region 6 linked the explosion source directly to the western region. Friction and geometric changes played key roles in attenuating the explosion waves.

Table 2. Peak overpressures and arrival times for each pipe in the experimental mine.

Region	Component #	Peak Overpressure (bar)	Arrival Time (s)
R1	C59	0.522	0.039
R2	C9	0.492	0.039
R3	C24	0.512	0.039
R4	C11	0.375	0.117
R5	C31	0.443	0.195
R6	C29	0.144	0.156
R7	C43	0.151	0.156
R8	C50	0.04	0.195
Shaft 1	C53	0.001	Na
Portal 2	C2	0	Na

# The number of the selected pipe.

shows the peak overpressures and arrival times for each pipe component.

Table 2 indicates a decrease in peak overpressure with increasing distance from the explosion source. When multiple airways linked cells, as seen in regions 2 through 5, the attenuation effect of geometric changes was negligible. However, when only a single airway connected two regions, geometric changes played a crucial role in blast wave propagation. The drop in overpressure was noticeable between regions 6 and 7, as well as between regions 7 and 8.

5. Conclusions

This article introduced a decoupled network method applied to methane explosion simulation, emphasizing how these insights contribute to designing underground mine ventilation systems that better protect miners’ safety. The Microsoft SQL Express server’s pre-constructed explosion source database provided initial and boundary conditions. Pressure losses during the spread of the explosion wave were recorded through Flowmaster’s various components, such as transition components for changes in cross-sections, T-junctions for branching, Cd Discrete Loss for bends and obstacles, and pipe frictional losses. To predict the overpressure distribution from a gas explosion in a complex subsurface opening, a 1D FDM-based code, Flowmaster, was employed to simulate network-based shock propagation in the explosion wave section. Hypothetical methane explosions were modeled through both a full-scale experimental mine and a lab-scale parallel network.

The main conclusions of this paper are as follows:

(1)

The influence of geometric changes in a pipeline network on overpressure is significantly smaller than that in a single pipeline. The effect of geometric changes on overpressure reduction is more significant in series pipelines compared to parallel ones. It can be seen that changes in the geometric structure of a pipeline network lead to the interaction of shock waves within the pipelines, resulting in complex effects on overpressure. Therefore, when predicting the explosion process in a complex pipeline network, the results from a single pipeline should not be simply extrapolated. Instead, the network effect should be taken into account.

(2)

The attenuation effect of the explosion wave at the vent position is significant. With the attenuation caused by geometric changes in the pipeline network, the peak overpressure and destructiveness of the explosion wave can be effectively reduced. The location and quantity of vents can significantly influence the propagation and attenuation of the explosion wave within the pipeline. Adjusting the location, size, and other parameters of the vents can enhance the ability to prevent and control the explosion wave.

(3)

The negative pressure generated by a gas explosion in a full-size roadway is significant, and the pressure action time is prolonged. Therefore, its destructive effect cannot be ignored. The pressure-time history curve can be divided into three stages: pressure growth, pressure decay, and pressure fluctuation. When the compression wave propagates (overpressure), the gas density increases, whereas, during the propagation of the rarefaction wave (negative pressure), the gas density decreases. These differences in physical characteristics lead to variations in the propagation process of the two waves. Additionally, geometric changes in the mine channel can cause waveform distortion and energy attenuation during compression wave propagation, while they have little influence on the propagation of the rarefaction wave. As a result, the rarefaction wave can propagate over longer distances than the compression wave.

(4)

In a pipe exposed to the air, the propagation of the pressure wave can reflect and interfere within the pipe, leading to a shock phenomenon. This results in the pressure inside the pipe initially increasing, then decreasing, and, eventually, approaching the ambient air pressure. The occurrence of this phenomenon may be related to the velocity gradient within the pipe, and the bending or obstruction of the pipe can affect the vibration. Excessive velocity gradients lead to complex flow structures and increase the pressure drop in the pipeline. Pipeline vibration is a major factor contributing to structural damage and accidents.

(5)

The varying bending angles and T-shaped branches in the pipeline network result in different waveform distortions and pressure attenuations as the explosion wave passes through these areas. Bifurcations in the pipeline cause the shunting of the pressure wave, leading to a weakening of energy and affecting the wave’s propagation. Different pipe bending angles lead to inconsistent interactions within the pipe, impacting both the propagation of the explosion wave and pressure changes.

Overall, this study enhances the efficiency of predicting gaseous explosions, and a decoupling numerical method was employed to simulate the full-scale underground network, offering robust technical support and a scientific foundation for the safe design of underground mine ventilation systems.