Research on the Dynamic Leaking and Diffusion Law of Hydrogen-Blended Natural Gas under the Soil–Atmosphere Coupled Model

Abstract

With the breakthrough in mixing hydrogen into natural gas pipelines for urban use, the widespread application of hydrogen-blended natural gas (HBNG) in energy delivery is imminent. However, this development also introduces significant safety concerns due to notable disparities in the physical and chemical properties between methane and hydrogen, heightening the risks associated with gas leaks. Current models that simulate the diffusion of leaked HBNG from buried pipelines into the atmosphere often employ fixed average leakage rates, which do not accurately represent the dynamic nature of gas leakage and diffusion. This study uses computational fluid dynamics (CFD) 2024R1 software to build a three-dimensional simulation model under a soil–atmosphere coupling model for HBNG leakage and diffusion. The findings reveal that, in the soil–atmosphere coupling model, the gas diffusion range under a fixed leakage rate is smaller than that under a dynamic leakage rate. Under the same influencing factors in calm wind conditions, the gas primarily diffuses in the vertical direction, whereas under the same influencing factors in windy conditions, the gas mainly diffuses in the horizontal direction.

1. Introduction

At present, countries mainly rely on traditional fossil energy sources to facilitate development, with them playing an important role in promoting economic development but also bringing about problems such as carbon pollution and emissions. This runs counter to the existing concept of sustainable development. Therefore, governments are committed to the development of clean energy, and there is growing demand for sustainable, low-pollution clean energy [1,2,3]. Hydrogen energy is seen as an important carrier to ensure energy security and renewable energy [4,5,6]. In recent years, governments and researchers have been focusing on the development of hydrogen energy, which has the advantages of ‘high calorific value and zero carbon emissions’ [7,8,9].

Currently, many countries including Russia [10], France [11], the United Kingdom [12,13], the Netherlands [14], Australia [15], and others are conducting tests and research on hydrogen delivery via natural gas pipelines. Figure 1 illustrates the blending ratios of hydrogen in the major natural gas pipelines of developed countries, with the current safety limit not exceeding 50% [16]. Additionally, the Chinese government has incorporated ‘hydrogen energy technology’ into its upcoming five-year plan for national key research and development [17,18,19]. Concurrently, it has initiated projects involving hydrogen blending in urban gas pipelines [20].

However, according to relevant researchers and academics [21], mixing hydrogen with natural gas changes the original properties of natural gas and increases the risk of pipeline operation issues, including pipeline leaks and compressor failure. Hydrogen-doped natural gas pipelines are not only affected by pipeline aging and third-party damage but also by significant potential impacts such as hydrogen embrittlement, leading to inevitable leaks in such pipelines [20,22].

In recent years, a number of explosion accidents related to gas leakage from buried pipelines have occurred both at home and abroad, causing significant loss of personnel and property as well as serious damage to the surrounding environment. For example, in 2019, a hydrogen cylinder fire and explosion occurred at a South Korean enterprise; Airgas Corporation in the United States also experienced a hydrogen explosion [23,24]. The compression coefficient, ignition energy, and other parameters of hydrogen and methane are quite different, and hydrogen has a wider range of explosion limits [25]. Comparative analysis of hydrogen and natural gas incidents reveals that hydrogen accidents tend to have a higher fatality rate and more severe consequences than natural gas accidents.

Bu et al. [26] carried out a numerical simulation study on the diffusion law after the leakage of buried natural gas pipelines using the numerical simulation software CFD, focusing on the formation of hazardous area time (FDT) and the farthest diffusion range (GDR) after the gas diffusion. According to the simulation results, the nature of the soil has a significant effect on the diffusion of the gas, and the FDT increases by 5 times when the soil type changes from sandy soil to clay. The increase in pipeline pressure leads to an increase in the hazardous area. Zhang et al. [27] established a three-dimensional simulation model of leakage diffusion for underground hydrogen pipelines, focusing on the effects of leakage holes, soil type, pipeline pressure, and pipeline diameter on hydrogen diffusion. The results of the study show that when the hydrogen pipeline leaks, the hydrogen concentration increases with the increase in the leakage time, showing a symmetrical distribution trend; hydrogen has a large difference in different soil properties, and its diffusion is the fastest in sandy soil, and the diffusion range and concentration are higher in than other soil types. Hu [28] conducted a simulation study using Computational Fluid Dynamics (CFD) software to investigate the diffusion characteristics of buried pipelines containing hydrogen doping or pure hydrogen, influenced by porous media. The study revealed that diffusion of the gas through the soil is notably slow. Moreover, it was found that a higher ratio of pore diameter to pore length and width correlates with an increased risk of hydrogen leakage. Su et al. [29] carried out a numerical simulation study of urban buried hydrogen-doped pipeline leakage, focusing on exploring the extent of the influence of hydrogen doping ratio (5–20%), leakage pressure, leakage aperture, and other factors on gas diffusion. The results showed that the leakage of the hydrogen-doped pipeline had no significant effect on the generation of relevant risk factors. Zhang et al. [30] used the numerical simulation software CFD to carry out a numerical simulation of pinhole leakage in medium-pressure (≤0.4 MPa) buried hydrogen pipelines under different working conditions, focusing on exploring the effects of pipeline pressure, soil properties, pinhole locations, and pinhole diameters on leakage diffusion. According to the simulation results, the correlation between surface disaster radius and time was established. Zhu et al. [31] set up an experimental platform for the leakage of buried hydrogen-doped natural gas pipelines, focusing on different hydrogen doping ratios, the location of leakage holes, and leakage pressure, and measured the leakage of hydrogen-doped natural gas through experiments. The results of the study showed that leakage pressure exhibited a wide diffusion range when a significant proportion of individual components existed in hydrogen-doped natural gas. Peng [16] and others investigated the diffusion problem of hydrogen-doped natural gas in long-distance pipelines. A two-dimensional planar model was established using CFD numerical simulation software, focusing on the effects of hydrogen doping ratio, leakage aperture, temperature, and other factors on the diffusion of hydrogen-doped natural gas leakage. Ye [32] and others used CFD numerical simulation software to investigate gas diffusion in the coupled environment of soil and atmosphere. The results showed that compared to underground leakage, the diffusion distance and height of aboveground leakage under different conditions exhibit errors, with horizontal diffusion distance errors generally being larger.

Currently, most mainstream research simplifies leakage sources to a constant leakage area and leakage rate, which deviates significantly from real-world conditions. As a result, the potential impact zones obtained are conservative and do not provide a more accurate range of hazardous areas. Therefore, this study addresses this issue by fitting the gas leakage velocity and diffusion area and by developing a User-Defined Function (UDF) to control the dynamic changes in leakage velocity and area, thereby reflecting the actual gas leakage conditions. Compared to previously developed fixed leakage source models, this approach better simulates the gas leakage and diffusion patterns under real conditions. Consequently, this research employs the coupled model to investigate the diffusion characteristics of hydrogen-blended natural gas leaking from buried pipelines into the atmosphere, considering factors such as the hydrogen blending ratio (HBR), wind speed, temperature, and soil properties. The findings of this study can inform the development of future emergency decision-making procedures for managing pipeline leaks.

2. Numerical Method

2.1. Physical Model

To investigate the diffusion patterns of natural gas in the atmosphere after soil seepage, a three-dimensional simulation model is employed. Initially, an underground leakage diffusion model was established with dimensions of 4 m × 4 m × 3 m (length × width × height). Upon solving the underground leakage diffusion model, it was observed that the velocity of hydrogen-doped natural gas seeping from the soil to the surface approximates an elliptical distribution. Consequently, the gas diffusion inlet into the atmospheric domain model is configured as a circular region located at the center of the model’s bottom. To comprehensively examine the coupled diffusion of hydrogen-doped natural gas from soil to atmosphere, the dimensions of the atmospheric domain model were set to 450 m × 300 m × 150 m (length × width × height). The geometric configuration is illustrated in Figure 2.

2.2. Mathematical Model

The following assumptions were made for this study:

(a)

The soil is set up as an isotropic homogeneous porous medium and the spatial structure of the soil does not change during mass transfer.

(b)

There is no chemical reaction between the leaking gas and the surrounding soil.

(c)

Heat transfer between the gas and the soil is ignored. Only the mass transfer process occurs.

(d)

Assume that the soil pores are filled with air and ignore moisture in the soil.

When gas flows and diffuses, the fluid dynamics conservation equations, namely mass conservation, momentum conservation, and energy conservation, are usually followed. Therefore, the corresponding continuity, momentum, and energy equations can be derived. In addition, there are gas state equations and component transport equations [33].

The continuity conservation equation is shown in Equation (1):

$${\displaystyle \frac{\partial (\rho \epsilon )}{\partial t}}+\nabla \cdot \left(\rho v_i\right)=0$$

(1)

where $v_i$ is the gas velocity in x, y, and z directions, m/s.

The momentum conservation equation is shown in Equation (2):

$${\displaystyle \frac{\partial \left(\rho v_i\right)}{\partial t}}+{\displaystyle \frac{\partial \left(\rho v_i{v}_j\right)}{\partial x_j}}=-{\displaystyle \frac{\partial P}{\partial x_i}}+{\displaystyle \frac{\partial {\theta }_{ij}}{\partial x_j}}+F_i$$

(2)

where $\rho$ is the fluid density, kg/m³; ${\theta }_{ij}$ is the shear stress, N/m²; $F_i$ is $i$ the directional mass force, m/s².

The substance transport equation is shown in Equation (3):

$${\displaystyle \frac{\partial }{\partial t}}(\epsilon \rho {\omega }_i)+\nabla \cdot \left(\rho {\omega }_i\nu \right)=\nabla \cdot (\rho D\nabla {\omega }_i)+S_i$$

(3)

where v is the velocity vector, m/s; D is the diffusion coefficient; $\nabla$ is the Hamiltonian operator; $\epsilon$ is the porosity; $S_i$ is the user-determined rate of the i substance, m/s; $\rho$ is the fluid density, kg/m³; ${\omega }_i$ is the volume fraction of component i.

The turbulent flow of gases presents a highly intricate issue necessitating the selection of appropriate turbulence models tailored to specific scenarios. Factors influencing model selection typically include whether the fluid behaves as an ideal gas, the requisite computational precision, and the computational capabilities available. In this study, the $k$-$\epsilon$ model was employed, with its specific equations detailed in Equations (4) and (5) below:

$k$ equation:

$${\displaystyle \frac{\partial (\rho k)}{\partial t}}+{\displaystyle \frac{\partial (\rho k{u}_i)}{\partial x_i}}={\displaystyle \frac{\partial }{\partial x_j}}\left[\left(\mu +{\displaystyle \frac{{\mu }_t}{{\sigma }_k}}\right){\displaystyle \frac{\partial k}{\partial x_j}}\right]+G_k-\rho \epsilon$$

(4)

$\epsilon$ equation:

$${\displaystyle \frac{\partial (\rho \epsilon )}{\partial t}}+{\displaystyle \frac{\partial (\rho \epsilon u_i)}{\partial x_i}}={\displaystyle \frac{\partial }{\partial x_j}}\left[\left(\mu +{\displaystyle \frac{{\mu }_t}{{\sigma }_{\epsilon }}}\right){\displaystyle \frac{\partial \epsilon }{\partial x_j}}\right]+\rho C_{1}E\epsilon -\rho C_2{\displaystyle \frac{{\epsilon }^2}{k+\sqrt{\nu \epsilon }}}$$

(5)

2.3. Implementation of the Soil Atmosphere Coupled Diffusion Dynamic Leakage Source

Firstly, according to the calculation of Le Chatelier’s law in Formula (6) for the explosive limit of mixed gases [34,35,36,37], the lower explosion limit (LEL) of the gas mixture at various HBRs is determined, as shown in

Table 1. Lower explosive limit of HBNG.

HBR (%)	0	5	10	15	20	25	30
LEL (vol%)	5	4.94	4.88	4.82	4.76	4.7	4.65

. The potential explosion zone of the gas mixture is defined by its LEL, while for methane and hydrogen, their respective lower explosion limits mark their potential explosion zones.

$${\displaystyle \frac{1}{{\mathrm{LEL}}_{\mathrm{mix}}}}={\displaystyle \sum }_{i=1}^N{\displaystyle \frac{{\mathrm{Y}}_{\mathrm{i}}}{{\mathrm{LEL}}_i}}$$

(6)

where $LE{L}_{mix}$ is the lower explosive limit of the gas mixture; $i$ and $mix$ are gas component i and the gas mixture; $Y_i$ is the concentration of component i, %; $LE{L}_i$ is the lower explosive limit of component i, %; $N$ is the amount of combustible gas components.

Based on the diffusion law of gas under the influence of porous media, we fit the diffusion range and diffusion rate of gas on the ground surface over time. According to GB/T 20936.1-2022 “Gas detectors for explosive environments—Part 1: Performance requirements for combustible gas detectors” and references [33,38], the maximum alarm concentration of gas shall not exceed 20% of the LEL. Therefore, relevant regulations are also adopted in the UDF compilation. Taking HBR = 15% as an example, if the LEL of the mixed gas is 4.82%, the alarm concentration should not exceed 0.96%. The diffusion range and gas diffusion rate when the HBNG concentration on the ground reaches 0.96% will be recorded, and the values will be fitted. Gas leakage is a process that changes dynamically with the leakage time. The dynamic change rate of gas leakage is obtained by integrating the dynamic change curve, which truly reflects the situation of gas leakage. We obtained btain Formulas (2) and (3) as the formulas for compilation in the UDF. Figure 3 shows the gas diffusion range and leakage rate under loam soil conditions at HBR = 15%.

$$R_0=2.1-[3.4/(1+(\exp (t+39.95)/453.89))]$$

(7)

$${\mathrm{V}}_0=9.44\times {10}^{-5}+[0.015/(1+\exp ((t+5871.52)/1046.43))]$$

(8)

where R₀ is the radius of gas diffusion, m; ${\mathrm{V}}_0$ is the gas diffusion velocity, m/s; $t$ is the leakage time, s.

In this study, Equations (7) and (8) were fitted and integrated into a user-defined function (UDF). This UDF was developed to facilitate its importation into the numerical simulation solution via the C editor in Visual Studio 2019. This integration aims to dynamically adjust both the diffusion range and diffusion rate over time within the specified region. Figure 4 illustrates the flowchart of the UDF compilation process. The essential code snippets are detailed below:

begin _f_ loop(f, t)

{F_CENTROID(xc, f, t);

temp = sqrt((xc[0]-pX)*(xc[0]-pX) + (xc[1]-pY)*(xc[1]-pY));// Determine the coordinates of the point of the center of diffusion of the dynamic circle leakage.

if (temp < R0)

{F_PROFILE(f, t, i) = flow_V0;}}

After compilation via UDF, Figure 5 illustrates the variation in leakage velocities across different time points within the dynamic circular domain for HBR = 15%. In this figure, leakage velocities are centered at the circle’s midpoint and progressively expand over time.

2.4. Simulation Parameters and Boundary Conditions

This article chooses a three-dimensional model with a velocity inlet set on the bottom wall and a wind speed inlet set on the left wall. Typically, researchers set the wind speed to a constant value to simplify model calculations, but this differs greatly from reality. The present article employs variable wind speed and utilizes an exponential function to depict the near-ground wind speed. The mathematical expression is shown in Equation (9) [39]. Other walls are set as pressure outlets to simulate the free flow of gas. The ground is set as “Wall”, and the specific boundary conditions are shown in

Table 2. Setting of boundary conditions.

Boundary Name	Boundary Type	Conditional Settings
Left wall	Velocity inlet	Wind Speed Setting
Lower wall	Wall	Preserve the default
Other wall	Pressure outlet	Preserve the default
Dynamic faceted regions	Velocity inlet	Mounting the UDF

.

$$u_z=u_1{({\displaystyle \frac{z}{z_1}})}^{\alpha }$$

(9)

where $u_z$ is the average wind speed at a near-surface height of z; $u_1$ is the standard wind speed (m/s); $\alpha$ is the surface roughness, which generally takes the value of 0.1–0.4; in this paper, it takes 0.4 [40]; $z_1$ is the height from the ground and $z_1$ is taken as 10 m.

The soil can be roughly classified into sandy and loamy soils by the average grain size of the soil [41,42]. China has abundant and diverse soil resources, which can be roughly divided into sandy soil and loam soil based on their average particle size.

Table 3. Physical parameters of different soil types.

Soil Type	Average Particle Size/mm	Porosity	$\mathbf{Viscous} \mathbf{Drag} \mathbf{Coefficient} (1/\mathit{\alpha })$/m−2	$\mathbf{Inertial} \mathbf{Drag} \mathbf{Coefficient} (\mathit{C}$2)/m−1
Sandy soil	0.50	0.25	2.16 × 1010	3.36 × 105
Loam soil	0.05	0.43	2.45 × 1011	5.02 × 105

provides the physical parameters of different soil properties.

According to the common pressure range of medium- and low-pressure pipelines in GB 55009-2021 [43] Gas Engineering Project Specification and referencing relevant literature [28,30,44,45,46,47,48,49], we set the pipeline pressure to 0.4 MPa and selected three hydrogen blending ratios of 5%, 15%, and 30%. Due to the characteristic of pipeline leakage mostly being small hole leakage [50], in this study, we mainly explore the leakage of small holes with a leakage aperture of 10 mm. Based on the research of relevant researchers and relevant references [40,51], the wind speed was selected based on the annual average wind speed in a certain area of southwestern Chongqing, which ranges from 0.6 m/s to 7 m/s, and we selected the temperatures of 278 K, 288 K, and 300 K.

2.5. Model Reliability and Grid Independence Verification

Despite the correlation analysis of the mathematical model, it is still necessary to verify the model to ensure the reliability of the simulation. Liang [52] calculated the leakage and diffusion law of natural gas in the atmospheric domain through a simulation model of 50 m × 50 m × 50 m, in which the mass flow rate was defined as 0.07 kg/s and the leakage time was 60 s. In this paper, the physical model, boundary conditions, and parameter values consistent with the beam are selected, and the numerical simulation results of this study are compared with those of the beam, as shown in Figure 6 below. Through comparison, it can be observed that the gas diffusion law in the two cloud images after the gas leakage is roughly similar, showing an upward diffusion plume, and the development direction of the jet is basically the same. Therefore, it can be considered that the numerical simulation model proposed in this paper is feasible.

The computational domain was meshed using unstructured meshing with Fluent Meshing software, as illustrated in Figure 7. The grid size was locally refined to 0.01 m at the leakage points. After verification, the grid quality exceeded 0.7, ensuring sufficient calculation accuracy. For simulating gas leakage in the atmosphere, the model was tested with four different grid densities: 580,000, 1,210,000, 2,020,000, and 3,100,000 grid cells. A monitoring point positioned 10 m above the leakage source tracked the methane concentration over time, depicted in Figure 8. After evaluating calculation accuracy and computer performance, a grid of 1,210,000 cells was selected for the final simulation.

3. Results and Discussion

To ensure the safety of personnel in the vicinity following a gas leak, this study adheres to relevant national standards [53] and references [54,55,56]. As a precautionary measure, the concentration limit is set at 1/50 of the 10% Lower Explosive Limit (LEL). Areas where gas concentrations exceed this limit are identified as potentially hazardous and require safety investigation. Consequently, subsequent sections will present diffusion distribution cloud maps of hydrogen-blended natural gas (HBNG) under varying operational conditions. These maps will delineate potentially hazardous areas bounded by 1/50 of the 10% LEL of the gas mixture’s concentration.

3.1. Comparative Analysis of the Leakage Fixed Source and Dynamic Leakage Source

The leakage velocity of hydrogen-doped natural gas from the soil to the ground surface typically forms an elliptical distribution. There exists a significant disparity between the diffusion rates of the gas in soil and in the atmospheric environment. Consequently, many researchers use the stabilized surface leakage velocity as the initial condition. To determine this average leakage velocity, one must first integrate the leakage velocity over time, as demonstrated in Equation (10).

$${\displaystyle {\int }_{100}^{1800}{v}_0(t)}dt$$

(10)

When considering the average leakage velocity, Figure 9 illustrates contour plots depicting the diffusion distribution of HBNG in the atmosphere at various time intervals. Within 60 s of the leakage, the gas diffusion range expands rapidly vertically due to minimal wind turbulence. By 600 s, the gas concentration reaching hazardous levels extends vertically to the model boundary, covering a distance of 150 m, with a horizontal diffusion of 41 m. Between 600 and 1800 s, the diffusion range stabilizes under the influence of atmospheric buoyancy forces. Over time, the diffusion range of low-concentration areas distant from the leakage point exhibits significant changes, with vertical diffusion distances of HBNG increasing gradually while the growth rate diminishes.

When the leakage rate varies over time (dynamic leakage rate), Figure 10 presents contour plots illustrating the diffusion distribution of HBNG in the atmosphere at different time intervals. The plots indicate that as leakage time increases, the gas diffusion range gradually expands vertically. Although the diffusion patterns of dynamic and fixed leakage rates are generally similar, the dynamic rate results in a larger diffusion range. For instance, at a leakage time of 10 s, the vertical extent of the hazardous area calculated using a fixed velocity is 18.2 m, whereas with dynamic leakage velocity calculation, it extends to 18.45 m vertically.

Figure 11 shows the range of potential ignition and explosion region formed by gas in a vertical direction calculated by using the fixed leakage source and dynamic leakage source at different times. It shows that the use of UDF dynamic leakage velocity results in a greater vertical spread of the potential explosion area compared to using a fixed leakage velocity. For instance, at 5 s of HBNG leakage, the vertical extent of the potential explosive region measures 5.65 m with dynamic leakage velocity and 5.4 m with fixed leakage velocity. Similarly, at 30 s of leakage, the vertical extent is 25.3 m with dynamic leakage velocity compared to 25.22 m with fixed leakage velocity.

3.2. Analysis of Dynamic Leakage Sources

3.2.1. Effects of Different Wind Speeds

(1)

Z-X planar HBNG leakage cloud map

In this section, conditions 1, 2, 3 and 4 in

Table 4. Operating condition setting table.

NO.	Pipeline Pressure/MPa	HBR/%	Soil Type	Air Velocity/m/s	Temperature/K
1	0.4	15	Loam	0	288
2	0.4	15	Loam	1.5	288
3	0.4	15	Loam	3	288
4	0.4	15	Loam	5	288
5	0.4	5	Loam	0	288
6	0.4	5	Loam	1.5	288
7	0.4	5	Loam	3	288
8	0.4	5	Loam	5	288
9	0.4	30	Loam	0	288
10	0.4	30	Loam	1.5	288
11	0.4	30	Loam	3	288
12	0.4	30	Loam	5	288
13	0.4	15	Sandy	0	288
14	0.4	15	Sandy	1.5	288
15	0.4	15	Sandy	3	288
16	0.4	15	Sandy	5	288
17	0.4	15	Loam	0	278
18	0.4	15	Loam	0	300

are used to study the diffusion law of HBNG in the atmosphere under different wind speeds, and the influence of different wind speeds on gas diffusion after leakage is discussed. Simulations are conducted under four wind speed conditions: 0 m/s, 1.5 m/s, 3 m/s, and 5 m/s. Figure 12 depicts the hazardous area’s cloud diagram resulting from HBNG diffusion at different wind speeds in the Z-X plane. From Figure 12, it is evident that wind speed significantly influences the diffusion path of HBNG. At 0 m/s wind speed, the gas primarily diffuses vertically. As wind speed increases, the vertical diffusion range decreases, showing a negative correlation. Conversely, the horizontal diffusion range increases with wind speed, indicating a positive correlation. Hence, greater safety risks exist for populations residing in downwind directions.

(2)

HBNG diffusion range

Different wind speeds affect the diffusion of HBNG in potentially hazardous and explosive areas, as illustrated in Figure 13. It is evident from the figure that both the vertical diffusion range of potentially hazardous areas and explosive areas of HBNG decreases with increasing wind speed. Conversely, the horizontal diffusion range increases with wind speed. The maximum vertical diffusion range of the potentially hazardous area of HBNG ranges from 24 m to 150 m. At a wind speed of 0 m/s, the vertical diffusion range of the potential explosive area of HBNG is greatest. At a wind speed of 5 m/s, the vertical diffusion range decreases to a minimum of 2.82 m.

3.2.2. Effect of Different HBR

(1)

Z-X planar HBNG leakage cloud map

The study investigated the impact of varying HBR ratios on gas diffusion following leakage. Calculations were conducted for three scenarios: HBR = 5%, HBR = 15%, and HBR = 30%. Figure 14 illustrates the gas diffusion distances for HBNG under these conditions. From the figure, it can be observed that at HBR = 5%, the vertical diffusion extends up to 24 m within the potential explosive zone, whereas at HBR = 30%, it extends to 28 m, indicating a significant increase with higher HBR ratios. In the horizontal direction, the diffusion range at HBR = 5% spans 4.05 m, while at HBR = 30%, it slightly increases to 4.15 m. This minor change suggests that higher HBR ratios do not significantly enhance gas diffusion in the horizontal plane.

(2)

HBNG diffusion range

In this section, conditions 1–12 in Table 4 are selected, the ranges of potentially hazardous and explosive areas formed by gases under different HBRs are shown in Figure 15, from which it can be observed that the vertical diffusion range of the potentially hazardous and explosive areas of HBNG decreases with the increase in wind force, while the horizontal diffusion range increases with the increase in the wind speed. The maximum vertical diffusion range of the potentially hazardous area of HBNG is 150 m, and the minimum diffusion range is 23 m; the maximum vertical diffusion range of the potential explosive area is 28 m, and the minimum diffusion height is only 2 m.

3.2.3. Effect of Different Temperature

(1)

Z-X planar HBNG leakage cloud map

In order to study the influence of different temperatures on the gas diffusion after leakage, a city in southwest China was taken as an example, and the typical temperatures of T = 278 K, T = 288 K and T = 300 K were analyzed by selecting working conditions 1, 17 and 18 in Table 4. The gas diffusion ranges of HBNG at these temperatures are illustrated in Figure 16. Analysis of Figure 16 indicates that temperature minimally affects the gas diffusion range in the atmospheric environment post leakage. After 1800 s of leakage, the vertical diffusion heights of the potentially explosive gas region measure 21.5 m at T = 278 K, 25.41 m at T = 288 K, and 26.1 m at T = 300 K. The vertical diffusion heights of potentially hazardous regions have fully extended to the model’s boundary.

(2)

HBNG diffusion range

The range of potentially hazardous and explosive areas formed by HBNG at various temperatures is depicted in Figure 17. The figure illustrates that the extent of these areas increases slightly with rising temperature, although the overall change is minimal. At 278 K, the vertical diffusion height of the potentially explosive gas region measures 21.5 m, with a horizontal diffusion range of 3.95 m. When the temperature rises to 300 K, these measurements change to a vertical diffusion height of 26 m and a horizontal diffusion range of 4.15 m.

3.2.4. Effect of Different Soil Properties

(1)

Z-X planar HBNG leakage cloud map

In previous studies, we have observed significant effects of soil properties on gas diffusion. Thus, this section examines gas diffusion patterns in the atmosphere by considering gas diffusion velocities under different soil conditions as the initial boundary conditions. Prior research indicates that gas diffusion concentrations from clay soil on the ground surface do not reach hazardous levels. This study focuses on comparative analysis between loamy and sandy soils. Figure 18 illustrates the variation in HBNG gas concentrations over time following leaks under different soil conditions. The figure demonstrates that the gas concentration at ground level significantly influences diffusion. Specifically, after 1800 s of leakage, HBNG reaches 33.8 m at the lower explosive limit contour in sandy soil compared to 25.41 m in loamy soil, marking a 33% difference between the two.

(2)

HBNG diffusion range

Working conditions 1–4 and 13–16 in Table 4 were selected for research, the distribution of potentially hazardous and explosive zones created by HBNG under varying soil properties is illustrated in Figure 19. The figure highlights significant differences in the extent of explosion and hazard zones resulting from gas leakage under sandy and loamy soil conditions. In sandy soil, the vertical diffusion height of the potential explosive zone at wind speeds of 0 m/s reaches 33.8 m, whereas in loamy soil, it only extends to 25.41 m. As the wind speed increases, the vertical difference between the two zones gradually diminishes, and the horizontal extent of the potential explosive zone expands approximately linearly with wind speed. Regarding the hazardous area’s diffusion range, at a vertical wind speed of 0 m/s, gas diffuses up to 150 m from the top boundary of the model after 1800 s of leakage. As the wind speed increases, vertical diffusion diminishes while horizontal spread increases; at wind speeds of 1.5 m/s and above, gas diffuses to the model’s boundary under wind influence.

4. Conclusions

To characterize the similarities and differences in the diffusion patterns of leaking gases from the soil to the atmosphere coupled, a model using a fixed leakage rate and a dynamic leakage rate was used. In this study, based on CFD numerical simulation and the three theoretical equations of fluid dynamics, a functional relationship between leakage velocity and time is proposed and compiled and solved using the UDF. On this basis, a three-dimensional model with a length × width × height of 450 m × 300 m × 150 m was established to study the diffusion characteristics of the gas, and the effects of different wind speeds, HBR, temperature, and soil properties on the gas leakage and diffusion characteristics were mainly analyzed, and the main conclusions were drawn as follows:

(1)

After comparative analysis of the leakage source and dynamic leakage source, the diffusion range of the potential explosion region formed by gas in the dynamic leakage source is larger than that obtained by the leakage source used by previous researchers;

(2)

Wind speed has a significant impact on the hydrogen-doped natural gas leakage; with the increase in ambient wind speed, the vertical diffusion height of hydrogen-doped natural gas is gradually reduced, and the two show a negative correlation. The horizontal diffusion distance of the gas increases with the increase in wind speed, and the two show a positive correlation;

(3)

Different hydrogen doping ratios (HBR) have a certain effect on the diffusion range of the gas in the vertical region and have less effect on the diffusion range of the gas in the horizontal direction;

(4)

At atmospheric temperatures of 278 K, 288 K and 300 K, the diffusion pattern of the leaking gas and the range of potentially hazardous areas and explosive areas formed in the atmosphere are not greatly changed, and the atmospheric temperature has less influence on the diffusion pattern of the gas after the leakage and the range of potential influence areas formed;

(5)

At the same moment, the gas diffusion velocity of different soil properties has a significant effect on the diffusion of gas in the atmosphere; the vertical diffusion height of hydrogen-doped natural gas from sandy soil to the atmospheric domain after diffusion for 1800 s is 33.8 m and that of loamy soil is 25.41 m, which is a 33% difference between the two.