Numerical Investigation on the One-Way Coupling of Gas Leakage-Explosion and a New Quantitative Overpressure Attenuation Law in Underground Culverts

1  Introduction

With adjustments to the national energy structure and changes in energy consumption patterns, the use of natural gas has been increasing year by year. As the main energy transportation route, natural gas pipelines are increasingly attracting public attention for safety issues. In recent years, there have been multiple natural gas pipeline leakage and explosion accidents worldwide. Compared with open spaces, the consequences of natural gas leaks and gas cloud explosions in enclosed spaces are more severe.

To study the leakage of natural gas pipelines, Jo and Ahn [1] and Woodward and Mudan [2] used a small-hole leakage model to calculate the leakage rate. De Almeida et al. [3] and Montiel et al. [4] proposed a new model for pipeline leakage and small-hole leakage. Cai et al. [5] developed a 3D full-scale coupled CFD-EnKF model to predict the dynamic spatial-temporal evolution of gas leakage and dispersion within utility tunnels. Liao et al. [6] conducted a comprehensive parametric numerical study to quantify the influence of varied boundary conditions on natural gas diffusion behavior in tunnel environments. Zhang et al. [7] investigated gas explosion overpressure characteristics in underground squares and evaluated the mitigation efficacy of active water spray mechanisms.

To investigate the actual leakage and diffusion of gas, Denisenko et al. [8] found that when the leakage rate was low, the leaked gas first gathered at the top of the confined space and gradually diffused downward; when the leakage rate was high, the leaked gas was evenly distributed near the leakage point. Zhou et al. [9] constructed a 10-m full-size tunnel model in an underground space and analyzed the diffusion characteristics of CO2 as a tracer gas in that space. Li et al. [10] conducted experiments in a room with the same proportions to verify the accuracy of the simulation results for the leakage and diffusion of difluoromethane indoors. Tan et al. [11] studied the dispersion mechanism of hazardous gases in complex urban environments using numerical and experimental methods. Wang et al. [12] evaluated the structural vulnerability of overhead pipeline systems under gas explosion loads in utility tunnels through dynamic structural response modeling.

To study the law of explosion overpressure, Sato et al. [13] conducted explosion tests on hydrogen-air mixtures in polyethylene films and found that the explosion load was almost constant at positions with the same proportion of detonation distance. Gieras et al. [14] studied the variation of peak overpressure of methane-air mixture explosions at different temperatures and concentrations based on multiple sets of methane-air mixture explosion experiments. Mercx et al. [15] conducted experiments on a methane-air mixture cloud, and the results showed that the presence of obstacles significantly increased the explosive shock load. Zheng et al. [16] implemented a fluid-solid coupling simulation framework to analyze the dynamic interaction between methane explosion shock waves and the surrounding rock mass in restricted underground roadways. Ciccarelli et al. [17] conducted natural gas air explosion tests in a confined pipeline and found that the flame experienced multiple accelerations from initial formation to final steady-state propagation after ignition. The coupling effect of overpressure and obstacles significantly accelerated the spread of the flame.

In numerical simulation research, Chen and Wu [18] used ANSYS/LS-DYNA finite element software to explore the propagation process and attenuation law of explosive shock waves of combustible gases in narrow underground confined spaces, and fitted the attenuation calculation formula of explosive shock waves at different depths of underground spaces. Zhang et al. [19] used finite element software to simulate the explosion process of a methane-air mixture at different scales in tunnels, and the results showed that the scale effect gradually decreased with the increase of spatial distance containing methane air mixture. Wang et al. [20] used numerical simulation software to establish a physical model based on actual explosion accidents, and analyzed the most suitable filling scheme and gas reaction equivalence ratio. Que et al. [21] utilized a 1D finite difference method code to numerically investigate the network-based propagation characteristics of shock waves induced by a methane explosion within a subsurface ventilation system.

To study the damage patterns of explosions on buildings and personnel, Schenker et al. [22] conducted full-scale explosion tests on concrete slabs with and without protective measures. Xu and Lu [23], Low and Hao [24], Zhou et al. [25], Tai et al. [26], Yan et al. [27], Li et al. [28], and others used finite element software to analyze the dynamic response of reinforced concrete slabs under explosive impact. By analyzing the factors that affect the failure of reinforced concrete slabs, the degree of damage to reinforced concrete slabs was evaluated. Baker et al. [29] obtained a criterion for determining the damage caused by explosive shock waves to the human body within 3–5 ms based on a large amount of experimental data, and provided a survival rate curve for lung injury according to the overpressure impulse criterion. Bai et al. [30] performed a dynamic resilience assessment of natural gas compartments in utility tunnels using computational fluid dynamics (CFD) simulations to optimize quantitative risk-based emergency strategies.

However, despite extensive previous efforts, existing methodologies still face significant challenges in accurately evaluating gas explosion risks in urban underground infrastructures. Most current studies decouple the leakage-diffusion process from the subsequent explosion mechanics, often assuming a simplified, uniform gas distribution that fails to reflect the complex concentration gradients found in real-world culvert geometries. Furthermore, damage assessments typically treat structural integrity and human vulnerability as isolated metrics, lacking a unified quantitative framework for comprehensive urban safety evaluation. To address these limitations, this study develops a systematic multi-material coupled numerical framework. By implementing a one-way data mapping strategy, the steady-state methane concentration field simulated in FLUENT is integrated into ANSYS/LS-DYNA as the initial material state, ensuring high-fidelity explosion precursors. This research investigates the overpressure propagation characteristics across both enclosed and semi-enclosed culvert conditions. The primary contribution is the derivation of a novel overpressure attenuation law that explicitly accounts for the depth-to-width ratio r of the confined space. This formula, combined with integrated damage criteria for buildings and personnel, provides a robust tool for accurately determining the impact scope of natural gas explosion accidents in underground urban environments.

2  Leakage Diffusion and Explosion Theory

2.1 Theory of Combustible Gas Leakage and Diffusion

In this study, a small-hole leakage model is adopted as the leakage aperture d is less than 20 mm (d/D < 0.2). The subsequent dispersion of natural gas within the underground culvert is governed by fundamental fluid mechanics conservation laws [31]. The transport of mass, momentum, and energy for the methane-air mixture is resolved using the continuity, Navier-Stokes, and energy equations. These governing equations, integrated with the Standard k-ε turbulence model and species transport model, describe the dynamic evolution of gas concentration fields over time. Detailed formulations of these classical conservation laws are well-documented in the literature [32,33].

When the research object is an ideal gas, the state equation of a simple ideal gas can be used as [34]:

PV=nRT(1)

The equation of state for actual gases:

PV=ZnRT(2)

where P represents the absolute pressure of gas (Pa); V represents the gas volume (m3/mol); n represents the amount of substance (mol); Z represents the compressibility factor of gas; R represents the gas constant (J/(mol·K)); T represents the thermodynamic temperature (K).

2.2 Theory of Combustible Gas Explosion

At present, there are five methods for calculating overpressure in gas cloud explosions, including the Henrych model, Mills model, The Netherlands Organization (TNO) multi-energy method, trinitrotoluene (TNT) equivalent method, and the methane-air mixture method. Chi et al. [35] conducted overpressure prediction on full-scale leakage explosion experiments of gas pipelines. It was found that the Henrych model, Mills model, and TNT equivalent numerical simulation method were not suitable for accurately calculating the shock wave overpressure of gas cloud explosions. The TNO multi-energy method is simple to operate and has low computational costs, but some calculated values are seriously disconnected from experimental data. The numerical simulation method for methane-air mixture gas is cumbersome to operate and computationally expensive, but the overall relative error in predicting peak overpressure is small, and the accuracy is the highest. Therefore, the methane-air mixture gas method was chosen as the numerical simulation method for gas cloud explosions in this article.

When numerically simulating the process of mixed gas explosion, it is necessary to determine the parameters of the explosion shock wave front, where the explosion pressure p and the specific volume V of the explosion gas can be used to describe the internal energy e of the combustible gas, expressed as e = e (p, V).

The Rayleigh equation for the explosion wave is:

p−p0=D2V02(V0−V)(3)

The adiabatic Hugoniot equation for explosion waves is:

e−e0=12(p+p0)(V0−V)+QV(4)

where QV represents the heat of reaction (kJ/kg); D represents the propagation speed of C-J detonation wave (m/s); V and V0 represent the specific volume before and after the explosion.

For gas-phase detonation follow the ideal gas law. The internal energy of the gas can be expressed using the specific heat ratio equation containing γ:

e0=P0V0γ−1(or e=PVγ−1)(5)

Substituting the internal energy equation of the gas with specific heat ratio before and after the explosion into (4) yields:

PVγ−1−P0V0γ−1=12(P+P0)(V0−V)+QV(6)

Substitute Eqs. (3) and (4) into Eq. (6) and organize them into a quadratic equation with (p − p0) as the unknown variable:

V02D⋅λ+1λ−1(p−p0)2+[γp0V0(γ−1)D−Dγ−1](p−p0)+DV0⋅QV=0(7)

Solving the equation yields Eq. (8):

p−p0=γ−1γ+1⋅DV0{Dγ−1−γγ−1⋅p0V0D±[γp0V0(γ−1)D−Dγ−1]2−2(γ+1)γ−1⋅QV}(8)

For C-J stable detonation, there are:

[γp0V0(γ−1)DJ−DJγ−1]2=2(γ+1)γ−1⋅QV(9)

Simplifying Eq. (7) yields:

pJ−p0=γ−1γ+1⋅DJV0[DJγ−1−γγ−1⋅p0V0DJ](10)

Further simplifying Eq. (10) and ignoring p0, i.e., p ≫ p0, Eqs. (11) and (12) can be obtained:

DJ=2(γ2−1)QV(11)

pJ=11+γρ0DJ2(12)

By substituting Eq. (12) into Eqs. (1)–(8) and simplifying it, Eq. (14) can be obtained:

ρJ=γ+1γρ0(13)

For the explosion of methane air mixture gas, the state equation of C-J product can be considered as:

pJVJ=RTJMJ(14)

where R represents the universal gas constant, R=8.31J/(mol⋅K); MJ represents the average molar mass of C-J detonation products; TJ represents the gas explosion temperature (K).

Substituting Eqs. (12) and (14) into Eq. (14) yields Eq. (15):

TJ=MJR⋅γDJ2(γ+1)2(15)

The initial density formula for a mixed gas is:

ρ0=∑niMi22.4×10−3∑ni(16)

where Mi represents the molar mass of each component before the explosion of the mixed gas; ni represents the number of moles of each component before the explosion of the mixed gas.

The formula for the average molar mass MJ of gas detonation products is:

MJ=∑iNniMi∑iNni(17)

where Mi represents the molar mass of each component of the detonation product; ni represents the number of moles of each component in the detonation product.

The specific heat ratio formula for detonation products is:

γ=∑niCpi¯∑niCVi¯=∑niCVi¯+∑niR∑niCVi¯(18)

where CVi¯ and Cpi¯ represent the average specific heat capacity at constant volume and pressure of each component of the detonation product (J/(mol·K)); ni represents the number of moles of each component in the detonation product.

The average specific heat capacity of detonation products at constant volume can be obtained using the Caster average heat capacity calculation formula. The Caster average heat capacity calculation formula is as follows:

CV¯=a+bt×10−4t(19)

where t represents the gas explosion temperature (K); a and b represent the constant, the value depends on the gas properties.

3  Model Establishment

3.1 Establishment of Flammable Gas Leakage and Diffusion Model

A finite element model is established for a typical underground space in a certain urban underground culvert. Assuming a small hole leakage occurs in the middle position of the pipeline, with a leakage hole size of 10 mm, assuming the leakage hole is circular and the internal pressure of the pipeline is 0.4 MPa.

There are two types of underground culverts in urban areas: confined and semi confined. The specific working conditions are shown in Fig. 1, where a represents the confined condition. The area ratio of outlet b is 1/2, and the area ratio of outlet c is 1.

Figure 1: Finite element model of underground underdrain (a) Confined working condition (b) Export area ratio of 1/2 (c) Export area ratio of 1.

3.2 Establishment of Natural Gas Explosion Shock Model in Underground Confined Space

Establish a finite element model of an underground culvert with a length, width, and depth of 30 m × 3 m × 2 m, respectively. A concrete cover plate with a thickness of 0.17 m is laid above the culvert in the city. Considering symmetry issues, a 1/4 finite element model is established, and the finite element model of underground culvert pavement air mixed gas explosive is shown in Fig. 2.

Figure 2: Finite element model of underground culvert pavement air mixed gas explosive.

3.3 Data Coupling Strategy

The research adopts a one-way coupling method to integrate the leakage and explosion simulations. First, the natural gas leakage and diffusion process were simulated using FLUENT until a stable methane concentration field was formed. Then, the spatial distribution data of methane concentration (coordinate points and mass fractions) were extracted and discretized. These data were mapped onto the element integration points of the LS-DYNA mesh to define the initial material status of the gas-air mixture. It is assumed that the ignition occurs instantaneously; therefore, the gas cloud geometry and concentration distribution at the moment of ignition are frozen and serve as the initial conditions for the subsequent explosion and shock wave propagation analysis.

3.4 Grid Independence Verification

The indoor natural gas leakage model is a relatively simple numerical model. Therefore, only after using SpaceClaim to establish a suitable geometric model, FLUENT Meshing is used to divide the unstructured grid to meet the calculation requirements. The grid spacing is 6 cm, and the leakage area is encrypted. The total number of grids is 187,904. In order to obtain more reasonable and accurate results, it is necessary to verify the mesh independence of the mesh size of the leakage area. The mesh size of the leakage area is set to 0.1, 0.2 and 0.4 cm, respectively. Taking the underground closed space as an example, the center point is selected as the measuring point to analyze the change of methane concentration.

It can be seen from Fig. 3 that when the grid size of the leakage area is 0.1 and 0.2 cm, the methane concentration at the central measuring point in the room does not change much, and the two lines almost overlap. When the grid size of the leakage area is 0.4 cm, compared with the grid size of the leakage area of 0.2 cm, the methane concentration of the indoor central measuring point is quite different. Therefore, under the premise of satisfying the accuracy and computational efficiency, the grid size of the leakage area is 0.2 cm.

Figure 3: The variation of methane concentration at the central measuring point with different grid sizes.

3.5 Model Parameters and Boundary Conditions

When the finite element software FLUENT is used to simulate the leakage and diffusion of natural gas, the energy equation is turned on. Considering the effects of gravity and buoyancy, the gravity acceleration in the Z direction is set to −9.8 m/s2, and the Full Buoyancy Effects option under the Viscous Model control panel is turned on. The turbulence model adopts the k-ε model, and the species transport model is used to simulate the leakage and diffusion process of natural gas in the air. The ambient pressure is 101,325 Pa and the ambient temperature is 298 K.

In the process of natural gas leakage, the influence of wind speed is not considered. One end of the underground confined space is the leakage inlet, which is set as the boundary condition of the velocity inlet. The leakage of the pipeline is 0.0032 m3/s. Except for the inlet and outlet of natural gas leakage; other surfaces are set as wall boundary conditions.

3.6 Verification of the Numerical Framework Using Experimental Benchmarks

In order to verify the reliability of the numerical simulation method, this section carries out finite element modeling and calculation on the indoor natural gas leakage test done by Guo [36], and compares the finite element simulation results with the test results to verify the model parameters and boundary conditions in the finite element simulation process. The cuboid with a volume of 1.6 m × 0.9 m × 1.6 m was selected as the model space for the indoor natural gas leakage test. The opening on the right side of the space was the exit, and the size was set to 30 cm × 60 cm and 30 cm × 90 cm. In this paper, the working condition of the opening size of 30 cm × 60 cm was selected for verification. The leakage gas is natural gas, the pressure is 2000 Pa, the leakage port is circular, and the nominal diameter is 10 mm. The combustion efficiency tester was used to measure the methane concentration at the center of the model. The measurement frequency was measured every 60 s for 60 consecutive measurements. The specific test equipment is shown in Fig. 4.

Figure 4: Test equipment diagram [36].

According to the test parameters, the geometric model with the same size as the test is established in ANSYS/SpaceClaim, and the parameter setting and model calculation are carried out to obtain the methane concentration at the intermediate measuring point of the numerical model, which is compared with the test data, as shown in Fig. 5.

Figure 5: Comparison of experimental and simulated data results.

200 s before the leakage time, the simulation results are larger than the test results, mainly because there is air in the probe of the measuring point before the start of the test. When the test begins, the methane enters the probe and mixes with the air, diluting the methane, making the methane concentration measured in the early stage of the test smaller. After the leakage time of 1440 s, the simulation results are basically similar to the experimental results and reach a stable trend, and the error between the two is within 10%. It is verified that the numerical simulation method is accurate and reliable to analyze the model parameters and boundary conditions in the process of natural gas leakage and diffusion.

4  Analysis of Results

4.1 Analysis of Diffusion Results of Combustible Gas Leakage

4.1.1 Natural Gas Diffusion Patterns in Underground Culverts

(1)   Confined working condition

A YZ cross-section with X of 1.5 m was selected to analyze the methane concentration profiles at different leak times.

As shown in Fig. 6, after the gas pipeline leak, natural gas quickly gathered at the top of the culvert, and with the increase of the leak time, the natural gas spread to both sides, and the methane concentration decreased from the middle area to both sides. The methane concentration in the culvert was greater than 5% when the leak time was 400 s, and when the leak time was 800 s, the methane concentration in 2/3 of the space in the culvert exceeded 17%.

Figure 6: Distribution of methane concentration at different times.

A ZX section with Y of 15 m was selected to analyze the methane concentration profiles for different leak times.

As shown in Fig. 7, natural gas was injected vertically upward at the leakage hole, and the methane concentration was high at the leakage hole and vertically upward area when the leakage time was 400 s, and the methane concentration at the leakage hole, the top of the culvert, and the two sides of the culvert was stratified obviously. The methane concentration in the culvert was more than 15% when the leakage time was 800 s.

Figure 7: Distribution of methane concentration at different times.

(2)   Export area ratio of 1/2

A YZ cross-section with X of 1.5 m was selected to analyze the methane concentration profiles at different leak times.

As shown in Fig. 8, when the left end of the underground culvert is confined and an outlet exists at the right end, the leakage time is 400 s in addition to the low methane concentration in the outlet area, and the methane concentration in the remaining space is more than 5%, and the leakage time is 800 s when the methane concentration at the outlet end is low and the methane concentration at the confined end is uniformly increasing, the leakage time continues to grow, the methane concentration at the outlet end grows slowly and the methane concentration is stratified obviously.

Figure 8: Distribution of methane concentration at different times.

A ZX section with Y of 15 m was selected to analyze the methane concentration profiles for different leak times.

As shown in Fig. 9, when there is an outlet in the culvert, it has a significant effect on the diffusion of natural gas leakage, the methane concentration at the leakage hole and the vertical upward area is higher at the leakage time of 400 s, and a polygonal area with lower concentration appears on both sides of the leakage hole at the leakage time of 800 s due to the existence of the outlet. When the leak time is 4000 s, the area of the polygonal area with lower concentration decreases gradually.

Figure 9: Distribution of methane concentration at different times.

(3)   Export area ratio of 1

A YZ section with X of 1.5 m was selected to analyze the distribution of methane concentration at different leak times.

As shown in Fig. 10, when the left end of the culvert is confined and the right end is completely open, the methane concentration at the exit end centered on the leakage hole is lower compared to the concentration at the confined end, and the stratification of methane concentration at the exit end is more pronounced.

Figure 10: Distribution of methane concentration at different times.

A ZX section with Y of 15 m was selected to analyze the methane concentration profiles for different leak times.

As shown in Fig. 11, when the right side of the culvert is completely open, the methane concentration in the leakage hole and the vertical upward area is higher when the leakage time is 400 s. When the leakage time is 800 s, the methane concentration on both sides of the leakage hole is not uniformly distributed, and as the leakage time grows, the uneven area of the methane concentration is gradually reduced.

Figure 11: Distribution of methane concentration at different times.

4.1.2 Migration Patterns in Exploded Areas of Underground Culverts

(1)   Confined working condition

Post-processing of the simulation results using ANSYS CFD-Post to analyze the migration pattern of the explosion region at different leak times.

As shown in Fig. 12a,b, when the underground culvert is a confined space, the explosion region is full of the entire underground space in 400 s, and the concentration of methane in the space from the center of the leakage hole to the two sides of the gradual decrease. When the methane concentration in the underground culvert continues to rise, exceeding the upper limit of the explosion of 15%, as shown in Fig. 12c, the explosion region from the center to the ends gradually disappear. As shown in Fig. 12d, when the leakage time of 800 s, the underground culvert is filled with methane concentration of more than 15% of the gas cloud, there is no explosion region in the space.

Figure 12: Map of the distribution of explosion regions at different times.

(2)   Export area ratio of 1/2

ANSYS CFD-Post was used to post-process the simulation results and analyze the migration pattern of the explosion region at different leak times.

As shown in Fig. 13, when the underground culvert is a semi-confined space, the migration pattern of the explosion region with the confined space has obvious differences, the leak time is 200 s, similar to the confined space, the explosion region to the leakage hole as the center to the length of the culvert on both sides of the expansion, the leak time of 600 s, due to the existence of the exit, the part of the space close to the exit of the lower methane concentration, the exit end of the part of the space there is no Explosion region. With the growth of the leak time, the leak time of 2000 s, compared with the confined space near the exit end of the space there is still a certain volume of the explosion region, is due to the exit is connected to the outside world, near the exit end of the space of the methane concentration is maintained in the explosion limit range.

Figure 13: Map of the distribution of explosion regions at different times.

(3)   Export area ratio of 1

Post-processing of the simulation results using ANSYS CFD-Post to analyze the migration pattern of the explosion region at different leak times.

As shown in Fig. 14, the distribution of the explosion region in the first 800 s of the leak time is similar to the working condition b, but with the increase of the proportion of the exit area, the concentration of methane at the same location at the same time is lower than the working condition b. Leakage time of 2000 s, due to the increase of the exit area, near the exit end of the space to maintain the natural gas explosion within the limits of the explosion region is larger in volume.

Figure 14: Map of the distribution of explosion regions at different times.

As shown in Fig. 15, underground culvert for the confined space, leakage time 0~450 s when the volume of the explosion region increased suddenly, leakage time 450~800 s explosion region volume appeared to decrease suddenly, leakage time 800 s after the underground culvert in the gas cloud in the methane concentration of all greater than 15%, so the underground culvert does not exist in the explosion region.

Figure 15: Volume change in explosion region.

Underground culvert belongs to the semi-confined space, the exit area ratio of 1/2, the leakage time 0~440 s explosion region volume increased suddenly, when the leakage time of 750 s when the explosion region volume reaches its maximum value. The leakage time of 750~1000 s explosion region volume decreased suddenly, after the underground culvert explosion region to maintain a certain volume no longer change. With the increase in the exit area, the volume of the explosion zone that continues to exist at the outlet end of the underground culvert becomes larger.

4.2 Analysis of Natural Gas Explosion Impact Results in Underground Confined Spaces

4.2.1 Shock Wave Propagation Law of Natural Gas Explosion in Underground Confined Space

Numerical model calculation time is 4 ms, the entire concrete cover damage failure, the time is 10 ms, the shock wave is almost all washed out of the culvert. Natural gas explosion shock wave initially presented elliptical spherical outward propagation, as shown in Fig. 16a–e, with the increase in time, the explosion impact range increases, as shown in Fig. 16f, time for 100 ms when the gradual development of spherical shock wave.

Figure 16: Schematic diagram of the propagation of the blast wave on the road at different times.

In the numerical model to select a number of measurement points to record the explosion process at different locations in the shock wave overpressure change rule, measurement point schematic shown in Fig. 17.

Figure 17: Schematic diagram of measurement point selection.

Based on the measurement points selected in Fig. 17, the overpressure time-course curves for each position are shown in Fig. 18.

Fisgure 18: Overpressure time histories at different measuring point locations.

Fig. 18 different measurement points of the overpressure time curve can be seen, the culvert gas cloud explosion, 3 ms peak overpressure reached 840 kPa, and in about 10 ms shock wave overpressure value quickly reduced to normal levels, the value of 0. Explosion shock wave in the space of the culvert reflection superposition of the formation of a strong overpressure in the culvert above the position of 2 m above the culvert (Y1) peak overpressure reached 525 kPa, in the height of 16 m from the culvert at the location (Y8) shock wave overpressure peak value fell to 10.8 kPa.

After the explosion shock wave breaks through the culvert cover, it firstly expands outward in the shape of an ellipsoid ball, and then gradually develops into a spherical shock wave. The overpressure value of the shock wave at different positions on the wavefront surface is different, with the largest overpressure value in the Y-axis direction, the second largest in the X-axis direction, and the smallest in the Z-axis direction, and the overall overpressure value of the shock wave shows the trend of ‘gradually becoming smaller with the increase of the range’.

4.2.2 Scope of Impact of Natural Gas Explosions in Underground Confined Spaces

According to the typical explosion overpressure criterion and the damage degree of the building structure [37], four levels are divided, as shown in Table 1.

Four classes are classified according to the typical explosion overpressure criterion and the degree of injury to personnel [38], as shown in Table 2.

In this paper, in the calculation of underground confined space natural gas explosion accident consequences of the scope of selection of overpressure under the action of the building structure, personnel damage level, follow the ‘maximum risk’ principle, the level as the most conservative value, that is, select the minimum of all the damage threshold as the overpressure damage threshold.

Because the underground culvert natural gas explosion shock wave mainly in the horizontal direction (X-axis, Z-axis) on the road to produce damage to the building structure and personnel, through the Fig. 18 measurement point of the overpressure time curve can be seen, the horizontal direction of the explosion shock wave at the same distance in the direction of the X-axis of the overpressure value is greater, so the selection of the X-axis direction of the propagation of the shock wave law as an important basis for judging the scope of the impact of the consequences of the explosion on the data along the X-axis direction of the data of each measurement point for formula fitting. At the same time, in order to study the impact of different underground confined space depth on the consequences of the explosion, the establishment of the depth of the confined space were 1.5, 1.75, 2, 2.25, 2.5, 3 m of the six groups of numerical models, in which the numerical model only change the depth parameter, the other size parameters do not change. Simulation calculations to obtain the explosion shock wave propagation law at different depths and data acquisition, extract the peak shock wave overpressure at different locations, Fig. 18 for the depth of the confined space under the peak shock wave overpressure-distance curve (where r is the ratio of the depth and width of the underground confined space).

The Hopkinson-Cranz proportionality law is applied to obtain the corresponding shock wave attenuation laws for confined space explosions with different widths [39].

As shown in Fig. 19, the larger the size of the depth of the underground confined space, the larger the peak overpressure, with the shock wave propagation process in the road surface energy consumption, the farther along the edge of the width of the culvert shock wave peak overpressure is smaller. With the increase of distance under different r-value of the shock wave peak overpressure basically maintain a certain size ratio, in line with the classical explosion cube root law, but also indirectly verified the validity of the numerical calculation.

Figure 19: Variation of peak blast shock wave overpressure at different distances from the edge of the confined space.

According to the curve trend shown in Fig. 19, a modified exponential function model is used to fit the data, and the functional expression form is shown in Eq. (20), and the fitting results are shown in Fig. 20, and the specific function parameters are shown in Table 3.

Pex=A×eB×R(20)

where Pex represents the peak overpressure (kPa); R represents the distance from the point of location to the edge of the underground confined space (m); A and B represents the parameters to be fitted.

Figure 20: Fitting results of the peak overpressure of the blast shock wave at different distances from the edge of the confined space.

In order to obtain the underground confined space explosion dimensionless distance under the shock wave overpressure attenuation law, the function parameters in Table 3 will be fitted to the underground confined space structure size ratio r, the fitting results are characterised by a linear function of the equation, the fitting curve is shown in Fig. 21, the calculation of the formulae are as follows:

A=79.53+97.02r(21)

B=−0.32+0.11r(22)

Figure 21: Parameter fitting graph.

will fit the results into the Eq. (20), that is, at this time by the underground confined space structure size changes caused by changes in the propagation law of the explosive shock wave uniformly transformed by the depth of the underground confined space and the width of the ratio of r characterisation, and at the same time, according to the law of proportionality of the explosion to take into account the proportion of different widths of the structure of the width of the current model, you can get the different sizes of the rectangular underground confined space explosion shock wave attenuation law equation:

Pex=D×(79.53+97.02r)×e(−0.32+0.11r)×R(23)

where Pex represents the peak overpressure at the calculated position point; D represents the width of the underground confined space; r represents the depth to width ratio of the underground confined space; R represents the distance from the location point to the edge of the underground confined space. This formula is applicable to underground rectangular confined spaces with a depth-to-width ratio r ranging from 0.5 to 1.0.

Based on the explosion shock wave attenuation law formula, can be calculated for the maximum area of the underground confined space explosion, by the Eq. (23) can be transformed into:

R=f(Pex,D,r)=1−0.32+0.11rln⁡(PexD(79.53+97.02r))(24)

From the above formula combined with Tables 1 and 2, the damage level distance between the building structure and the personnel can be found in Table 4.

4.3 Example Analyses

To verify the engineering applicability of the generalized Eq. (23) derived in the previous section, a real-world urban underground culvert case is selected for comprehensive risk zoning and damage assessment, as shown in Fig. 22. The underground culvert has a length of 20 m, a width of 7.1 m and a depth of 3.95 m, and some of its sections are close to residential areas, where the density of people is high. Based on these dimensions, the depth-to-width ratio r is calculated as 0.5566. Assuming that a gas pipeline leaks through the underground culvert section near a residential building, the external factors lead to the formation of a confined space in the section, which is filled with a methane concentration of 9.5% (representing the stoichiometric concentration for worst-case scenario analysis) of the gas cloud.

Figure 22: Underground culvert in a city.

According to Eq. (24), can find the explosion accident on the building damage level impact range, the calculation results are shown in Table 5.

According to the Eq. (24), can be derived from the explosion on the personnel damage level impact range, the results of the calculation shown in Table 6.

Based on the calculation results of the consequences of the underground confined space explosion accident, the damage regions are categorized. The final impact ratio R1–R4 in Fig. 23 is determined by selecting the most conservative values from Tables 5 and 6 according to the maximum risk principle. Following the consequence distance formulas for different levels, the explosion impact zones of the underground confined space are accurately divided, as illustrated in Fig. 23.

Figure 23: Comprehensive risk zoning of the explosion accident based on the maximum risk principle.

5  Conclusion

From the above simulation study, this paper draws the following conclusions:

1.    Developed a multi-material coupled model (leakage-diffusion-explosion) for complex geometries.

2.    Quantified the impact of the depth-to-width ratio r on shock wave attenuation.

3.    Formulated a predictive model for overpressure decay, applicable for safety distance planning in urban underground spaces.

Acknowledgement: The authors would like to thank all individuals and organizations that contributed to this study through administrative and technical assistance, as well as material or equipment support. Their contributions were essential to the successful completion of this research.

Funding Statement: Financial support from the National Natural Science Foundation of China (Grant No. 52374068) is appreciated.

Author Contributions: All authors have made significant contributions to this work. Pengcheng Kang wrote the original draft and performed visualization. Yuanyuan Tian developed the methodology and conducted formal analysis. Ying Liu, Qian Xu, Yuting Chen, Lixin Jia, Shuge Guo and Heng Rong contributed to investigation, resources, data curation, supervision, and software implementation. Taolong Xu supervised the study, acquired funding, and led manuscript revision. All authors reviewed and approved the final version of the manuscript.

Availability of Data and Materials: Data will be made available on request.

Ethics Approval: This study involved exclusively numerical simulations of natural gas leakage and explosion scenarios. No human participants, animal subjects, personal data, or biological samples were involved in the research. Therefore, ethical approval from an institutional review board was not required for this work.

Conflicts of Interest: The authors declare no conflicts of interest.

References

1. Jo YD, Ahn BJ. A simple model for the release rate of hazardous gas from a hole on high-pressure pipelines. J Hazard Mater. 2003;97(1–3):31–46. doi:10.1016/S0304-3894(02)00261-3. [Google Scholar] [PubMed] [CrossRef]

2. Woodward JL, Mudan KS. Liquid and gas discharge rates through holes in process vessels. J Loss Prev Process Ind. 1991;4(3):161–5. doi:10.1016/0950-4230(91)80031-O. [Google Scholar] [CrossRef]

3. De Almeida JC, Antonio Velásquez J, Barbieri R. Development and experimental validation of a computational model for the analysis of transient events in a natural gas distribution network. Can J Chem Eng. 2014;92(10):1776–82. doi:10.1002/cjce.22024. [Google Scholar] [CrossRef]

4. Montiel H, Vílchez JA, Casal J, Arnaldos J. Mathematical modelling of accidental gas releases. J Hazard Mater. 1998;59(2–3):211–33. doi:10.1016/S0304-3894(97)00149-0. [Google Scholar] [CrossRef]

5. Cai J, Wu J, Yuan S, Kong D, Zhang X. Prediction of gas leakage and dispersion in utility tunnels based on CFD-EnKF coupling model: a 3D full-scale application. Sustain Cities Soc. 2022;80(6):103789. doi:10.1016/j.scs.2022.103789. [Google Scholar] [CrossRef]

6. Liao K, Wang Y, Chen D, He G, Huang Y, Zhang S, et al. Parametric study on natural gas leakage and diffusion in tunnels. J Pipeline Syst Eng Pract. 2023;14(2):04023003. doi:10.1061/jpsea2.pseng-1396. [Google Scholar] [CrossRef]

7. Zhang C, Ma J, Shen J, Jiao D, Chen J, Wu X, et al. Gas explosion characteristics and spray control mechanism in underground square. PLoS One. 2024;19(4):e0293421. doi:10.1371/journal.pone.0293421. [Google Scholar] [PubMed] [CrossRef]

8. Denisenko VP, Kirillov IA, Korobtsev S, Nikolaev II. Hydrogen-air explosive envelope behavior in confined space at different leak velocities. In: Proceedings of the 4th International Seminar on Fire and Explosion Hazards; 2003 Sep 8–12; Londonderry, UK. 194 p. [Google Scholar]

9. Zhou K, Li F, Cai H, Yang Y, Peng F, Chen L, et al. Experimental and numerical investigation of gas diffusion under an urban underground construction. Energy Built Environ. 2021;2(4):436–44. doi:10.1016/j.enbenv.2020.09.002. [Google Scholar] [CrossRef]

10. Li YN, Tao JJ, Han YC, Han X, Qin J. Numerical and experimental study on the diffusion property of difluoromethane (HFC-32) in leakage. Procedia Eng. 2014;71:34–43. doi:10.1016/j.proeng.2014.04.006. [Google Scholar] [CrossRef]

11. Tan W, Wang K, Li C, Liu L, Wang Y, Zhu G. Experimental and numerical study on the dispersion of heavier-than-air gas in a complex urban environment. Process Saf Environ Prot. 2018;116:640–6. doi:10.1016/J.PSEP.2018.03.027. [Google Scholar] [CrossRef]

12. Wang D, Dong J, Chen X, Du J, Shu D, Krassowska J. Influence of gas explosions in utility tunnels on the structural safety of overhead pipelines. Buildings. 2025;15(18):3391. doi:10.3390/buildings15183391. [Google Scholar] [CrossRef]

13. Sato Y, Iwabuchi H, Groethe M, Merilo E, Chiba S. Experiments on hydrogen deflagration. J Power Sources. 2006;159(1):144–8. doi:10.1016/j.jpowsour.2006.04.062. [Google Scholar] [CrossRef]

14. Gieras M, Klemens R, Rarata G, Wolański P. Determination of explosion parameters of methane-air mixtures in the chamber of 40dm3 at normal and elevated temperature. J Loss Prev Process Ind. 2006;19(2–3):263–70. doi:10.1016/j.jlp.2005.05.004. [Google Scholar] [CrossRef]

15. Mercx WP, van den Berg AC, Hayhurst CJ, Robertson NJ, Moran KC. Developments in vapour cloud explosion blast modeling. J Hazard Mater. 2000;71(1–3):301–19. doi:10.1016/s0304-3894(99)00085-0. [Google Scholar] [PubMed] [CrossRef]

16. Zheng Q, Ding P, Yan Z, Zhu Y, Zhang J. Dynamic response of methane explosion and roadway surrounding rock in restricted space: a simulation analysis of fluid-solid coupling. Appl Sci. 2025;15(17):9454. doi:10.3390/app15179454. [Google Scholar] [CrossRef]

17. Ciccarelli G, Johansen CT, Parravani M. The role of shock–flame interactions on flame acceleration in an obstacle laden channel. Combust Flame. 2010;157(11):2125–36. doi:10.1016/j.combustflame.2010.05.003. [Google Scholar] [CrossRef]

18. Chen GH, Wu JJ. Propagation laws of gas explosion shock waves in underground confined space. Nat Gas Ind. 2017;37(2):120–5. (In Chinese). [Google Scholar]

19. Zhang Q, Pang L, Liang HM. Effect of scale on the explosion of methane in air and its shockwave. J Loss Prev Process Ind. 2011;24(1):43–8. doi:10.1016/j.jlp.2010.08.011. [Google Scholar] [CrossRef]

20. Wang D, Qian X, Yuan M, Ji T, Xu W, Liu S. Numerical simulation analysis of explosion process and destructive effect by gas explosion accident in buildings. J Loss Prev Process Ind. 2017;49:215–27. doi:10.1016/j.jlp.2017.07.002. [Google Scholar] [CrossRef]

21. Que S, Zeng J, Wang L. Numerical investigation of network-based shock wave propagation of designated methane explosion source in subsurface mine ventilation system using 1D FDM code. Sustainability. 2024;16(22):9935. doi:10.3390/su16229935. [Google Scholar] [CrossRef]

22. Schenker A, Anteby I, Gal E, Kivity Y, Nizri E, Sadot O, et al. Full-scale field tests of concrete slabs subjected to blast loads. Int J Impact Eng. 2008;35(3):184–98. doi:10.1016/j.ijimpeng.2006.12.008. [Google Scholar] [CrossRef]

23. Xu K, Lu Y. Numerical simulation study of spallation in reinforced concrete plates subjected to blast loading. Comput Struct. 2006;84(5–6):431–8. doi:10.1016/j.compstruc.2005.09.029. [Google Scholar] [CrossRef]

24. Low HY, Hao H. Reliability analysis of direct shear and flexural failure modes of RC slabs under explosive loading. Eng Struct. 2002;24(2):189–98. doi:10.1016/S0141-0296(01)00087-6. [Google Scholar] [CrossRef]

25. Zhou X, Kuznetsov V, Waschl J, Hao H. Numerical calculation of concrete slab response to blast loading. Trans Tianjin Univ. 2006;12(10):94–9. doi:10.1016/j.ijimpeng.2008.01.004. [Google Scholar] [CrossRef]

26. Tai YS, Chu TL, Hu HT, Wu JY. Dynamic response of a reinforced concrete slab subjected to air blast load. Theor Appl Fract Mech. 2011;56(3):140–7. doi:10.1016/j.tafmec.2011.11.002. [Google Scholar] [CrossRef]

27. Yan S, Zhang L, Wang D. Failure mode analysis for RC slab under explosive loads. J Shenyang Jianzhu Univ. 2005;3:177–80. [Google Scholar]

28. Li ZX, Shi YC, Shi XS. Damage analysis and assessment of RC slabs under blast load. J Build Struct. 2009;30(6):60–6. (In Chinese). [Google Scholar]

29. Baker WE, Cox PA, Kulesz JJ, Strehlow RA, Westine PS. Explosion hazards and evaluation. Amsterdam, The Netherlands: Elsevier; 2012. [Google Scholar]

30. Bai Y, Wu J, Yuan S, Reniers G, Yang M, Cai J. Dynamic resilience assessment and emergency strategy optimization of natural gas compartments in utility tunnels. Process Saf Environ Prot. 2022;165(3):114–25. doi:10.1016/j.psep.2022.07.008. [Google Scholar] [CrossRef]

31. Niu Y, Jiang L, Li Z, Gong Y, Du B, Mi H. Study on the law of pressure and flame propagation during gas explosion in the gas cabin of the utility tunnel. ACS Omega. 2025;10(16):16236–44. doi:10.1021/acsomega.4c10174. [Google Scholar] [PubMed] [CrossRef]

32. Hou Q. Simulation study of natural gas pipeline leakage and natural gas diffusion in the atmosphere [dissertation]. Harbin, China: Harbin Institute of Technology; 2009. (In Chinese). [Google Scholar]

33. Versteeg HK, Malalasekera W. An introduction to computational fluid dynamics. New York, NY, USA: Pearson Education Ltd; 2007. [Google Scholar]

34. Wang R, Mo L, Yang Y, Shi H, Xiao S, Chen C. A numerical study of surface building damage caused by natural gas explosion in underground spaces. Process Saf Environ Prot. 2025;202:107806. doi:10.1016/j.psep.2025.107806. [Google Scholar] [CrossRef]

35. Chi M, Jiang H, Lan X, Xu T, Jiang Y. Study on overpressure propagation law of vapor cloud explosion under different building layouts. ACS Omega. 2021;6(49):34003–20. doi:10.1021/acsomega.1c05332. [Google Scholar] [PubMed] [CrossRef]

36. Guo Y. Simulation of indoor combustible gas leakage diffusion and explosion effects [dissertation]. Chongqing, China: Chongqing University; 2011. (In Chinese). [Google Scholar]

37. Gu D. Analysis of structural damage to ground level brick building by explosion of buried gas pipeline [dissertation]. Chengdu, China: Southwest Petroleum University; 2019. (In Chinese). [Google Scholar]

38. Wang K, Liu Z, Qian X, Li M, Huang P. Comparative study on blast wave propagation of natural gas vapor cloud explosions in open space based on a full-scale experiment and PHAST. Energy Fuels. 2016;30(7):6143–52. doi:10.1021/acs.energyfuels.6b01293. [Google Scholar] [CrossRef]

39. Eslami M, Keshavarz MirzaMohammadi P, Khalilpour SH, Parsa H, Kodure V. Experimental and numerical investigation of blast wave attenuation by using barriers in different configurations and shapes. J Struct Eng. 2023;149(1):04022224. doi:10.1061/jsendh.steng-11408. [Google Scholar] [CrossRef]