Simulation Study of Gas Seepage in Goaf Based on Fracture–Seepage Coupling Field

Abstract

In order to solve the problem of gas overrun in the fully mechanized caving face and the upper corner of high gas and extra-thick coal seam, the fracture and caving process of the roof in the goaf is analyzed and studied by using the relevant theories of fracture mechanics and seepage mechanics. The mathematical model of fracture and caving of the immediate roof and main roof in the goaf is established. Combined with ANSYS Fluent 6.3.26, the seepage process of gas in coal and rock accumulation in the goaf under different ventilation modes is simulated. The distribution law of gas concentration in the goaf is obtained, and the application scope of different ventilation modes is determined. In addition, the influence of the tail roadway application and the wind speed size on the gas concentration in the goaf and the upper corner of the fully mechanized caving face is also explored. The results show that, affected by wind speed and rock porosity, along the strike of the goaf, about 30 m near the working face, the gas concentration is low and growth is slow. In the range of 30~160 m, the gas concentration increases rapidly and reaches a higher value. After 160 m, the gas concentration tends to be stable. Along with the tendency of the working face, the gas concentration in the goaf increases gradually from the inlet side to the return side, and the gas concentration increases noticeably near the return air roadway. Along the vertical direction of the goaf, the gas concentration gradually increases, and the concentration of the fracture zone basically reaches 100%. Different ventilation modes have different application scopes. The U-type ventilation mode is suitable for the scenario of less desorption gas in the coal seam, while U + I and U + L-type ventilation modes are suitable for the scenario of more desorption gas in coal seam or higher mining intensity. The application of the tail roadway can reduce the gas concentration in the upper corner to a certain extent, but it has limited influence on the overall gas concentration distribution in the goaf. In addition, when the wind speed of the working face should be controlled at 2.0~3.5 m/s, it is more conducive to the discharge of gas, the method of reducing the gas concentration in the upper corner by increasing the wind speed of the working face is more suitable for the case where the absolute gas emission of the fully mechanized caving face is low, and the effect is limited when the absolute gas emission is high. The above conclusions provide a reference for solving the problem of gas overrun in the goaf and the upper corner of a fully mechanized caving face.

1. Introduction

As an important basic energy source for economic and social development, coal has experienced great progress in mining scope and intensity. It plays a vital role in national economic development and modernization.

With the development of the economy and the progress of science and technology, large and super-large coal mines have been put into operation one after another. Coal mining has also changed from traditional blasting and general mining to fully mechanized mining with a higher efficiency. The record of mining coal seam depth has been consecutively broken. The high gas pressure and high content of coal seam have become important characteristics of modern coal mining. With the increase in coal mining depth, the intensity of gas explosions increases yearly [1,2]. Accompanied by gas explosions, especially serious gas explosion accidents still occur frequently. For example, from 2003 to 2013, there were 16 gas explosion accidents with more than 10 victims in Anhui Shengluling Coal Mine, Henan Daping Coal Mine, Shanxi Chenjiashan Coal Mine, Heilongjiang Dongfeng Coal Mine, Jilin Babao Coal Industry Company, Jilin Qingxing Coal Mine, Guizhou Dashan Coal Mine, Sichuan Taozigou Coal Mine, Hunan Simachong Coal Mine, and so on. There were seven gas explosion accidents with more than 100 deaths in each accident. And with the increase in coal mining intensity, gas disaster accidents are more prominent. Therefore, the prevention and control of gas disasters in the fully mechanized mining face of the high gas and extra-thick coal seam have become the focus of coal mine safety management. In order to explore the correlation between gas pressure and disaster prevention and control, researchers have carried out many related experiments, indicating that gas seepage and desorption in coal seams are the key causes of gas outbursts and explosions [3,4,5,6,7,8,9,10,11,12,13,14,15,16]. However, due to the particularity of its occurrence conditions, the seepage process of gas is affected by various external forces such as the in situ stress field, geoelectric field, and geothermal field, and various factors will have different degrees of influence on its flow in porous media. In this regard, experts and scholars have made active attempts to summarize the migration law of the coal seam and goaf gas under the action of a multi-physical field.

As early as 1943, Terzaghi K [17] studied the coupling of rock–soil–fluid. He first treated the flow of saturated fluid in one-dimensional elastic porous media as a flow–deformation coupling problem and proposed the famous effective stress formula. A one-dimensional consolidation model was established. Since then, many experts and scholars at home and abroad have made more valuable attempts to summarize the migration law of gas in the coal seam and goaf under the action of a multi-physical field [18,19,20,21,22,23,24,25]. In recent decades, Zhao et al. [26] proposed a nonlinear coupling mathematical model of solid deformation and gas seepage and established a finite element model for the numerical analysis to explain the migration and exchange of methane in the rock matrix and fractures. In 2006, Yang et al. [27] discussed the influence of fracture pore pressure on methane and carbon dioxide gas seepage in coal seams. It was found that the gas permeability coefficient in rock fractures changed parabolically with fracture pore pressure, and the gas permeability was negatively exponential with normal deformation and tangential deformation. Based on the influence mechanism of permeability, Connell (2009) [28] established a coupling model of gas flow and geomechanics and explored the applicability of the rock mechanics hypothesis to coal seam gas extraction. In 2010, Zhao et al. [29] studied the seepage law of gas in the goaf and the overlying strata. By establishing the mathematical model of the seepage field in the goaf, the seepage field equation was solved by the LBM digital method. Ronny et al. (2011) [30] established a one-dimensional mathematical coupling model that includes gas flow, adsorption, and geomechanical properties. In 2012, Li et al. [31] studied the influence of deep coal seam mining on gas seepage. In 2015, Wang et al. [32] studied the gas seepage characteristics of the coal seam by using the self-developed gas seepage triaxial experimental system, analyzed the influence of the Klinkenberg effect on gas seepage, and proposed a method for calculating the permeability of the coal seam considering the influence of gas dynamic viscosity, the compression coefficient, and the Klinkenberg effect. Danesh et al. (2016, 2017) [33,34] established a multi-field coupling model considering creep, simulated the creep-seepage evolution of the coal seam during gas recovery, and discussed the influence of creep on coal seam permeability under different conditions. Bertrand et al. (2017) [35] developed a hydraulic model for CBM production and analyzed a range of parameters that can highlight the impact of production scenarios or key reservoir-related parameters. In 2017, Zhang et al. [36] studied the influence of loading rate on gas seepage and temperature in coal and found that the seepage velocity of gas had critical deceleration characteristics under load. In 2019, Qiu et al. [37] established a microscopic pore structure model of coal and simulated the seepage behavior of gas in coal under 3 MPa pressure. It was found that the heterogeneity of the coal structure would lead to irregular changes in the cross-sectional area of pore channels along the seepage direction, which made the average seepage velocity of gas fluctuate greatly. With the increase in seepage length, the average seepage pressure and gas mass flow gradually decreased. In the same year, Xue et al. [38] established a hydraulic–mechanicaldamage coupling model considering the coupling relationship between coal damage, gas seepage, and coal deformation. They obtained the migration law of gas with coal seam failure. In 2021, Yan et al. [39] used the lattice Boltzmann method (LBM) at the representative elementary volume (REV) scale to simulate the gas seepage of the reconstructed pore structure and obtain the gas velocity distribution. The effects of porosity and pore structure on gas seepage were studied. It was found that the gas seepage velocity and apparent permeability increased with the increase in porosity. In 2024, based on the smooth particle hydrodynamics (SPH) method, Mu et al. [40] proposed a meshless numerical model of gas seepage in fractured coal bodies, studied the seepage law of gas in homogeneous and heterogeneous coal seams and coal bodies with single prefabricated defects, and revealed the seepage and stress coupling mechanism of fractured coal bodies during gas extraction.

Reviewing the existing research shows that the influence of various factors, including in situ stress, pore pressure, heterogeneity of the coal seam, and so on, on gas seepage has been discussed, and the law of gas seepage under different factors has been obtained. However, in the existing research field, few studies consider the influence of the ventilation mode and wind speed on gas seepage. In addition, due to the low permeability of the coal seam and the short time of gas extraction, the gas concentration in the working face, especially in the upper corner, is still high in the mining process of high gas and extra-thick coal seam. At present, the commonly used treatment methods in mines are buried pipe gas drainage in the goaf [41], using high drainage roadways to extract gas in the fracture zone, excavating directional long drilling groups [42], and excavating tail roadways to control gas seepage in the goaf. Among them, tail roadway extraction technology is one of the most important methods of controlling gas in the upper corner. It is important to facilitate the extraction of gas sources such as the goaf, adjacent layer, and surrounding rock in the working face by excavating the special gas discharge roadway so as to reduce the gas concentration in the upper corner of the working face, thus eliminating the gas overrun in the upper corner. This technology has been widely used, but there are few studies on the mechanisms of gas overrun in the upper corner of the tail roadway, especially for analyzing the gas seepage law. Starting from the theoretical basis, this paper fully analyzes the formation and evolution law of the goaf in the process of fully mechanized top coal caving mining, and it establishes the mathematical model of the initial fracture of the coal seam roof and the separation and caving model of the immediate roof and main roof, which lays a foundation for the study of the gas migration law in the goaf by a numerical simulation method. The influence of the working face wind speed and working face ventilation mode on gas concentration in the working face, goaf, and upper corner is emphatically analyzed. It also makes up for the shortcomings of other experts and scholars in this field.

In summary, this paper takes the caving time effect of the roof surrounding rock as the access point and uses the actual situation of gas control in the goaf of the N2306 working face of a mine in Shanxi Province. By studying the fracture and caving process of the roof, the ANSYS Fluent 6.3.26 is used to simulate the fluid–solid coupling process of gas in the collapsed coal and rock mass under different ventilation modes to analyze the distribution law of gas concentration in the goaf and to explore the influence of the ventilation mode and wind speed on gas seepage in the goaf. In addition to them, by focusing on the analysis of the influence of different ventilation methods (i.e., whether to set up tail roadways and different types of tail roadways) on the law of gas seepage, the application of outburst tail roadways has a positive effect on solving the problem of gas accumulation in the upper corner.

2. Materials and Methods

When the fully mechanized caving mining is adopted, the goaf is formed when the upper coal body is released. There are the pressure arch hypothesis, cantilever beam hypothesis, hinged rock block hypothesis, pre-formed fracture hypothesis, and masonry beam hypothesis in the failure process of the goaf roof. Although the contents and forms expressed by each hypothesis differ, they all show a delay in roof fracture and caving in the goaf. That is, with the advancement of the working face, the upper rock will form a certain length of probe in the goaf under the influence of its strength. The formation and existence of this probe create favorable conditions for the tail roadway to extract gas from the goaf and reduce gas concentration in the upper corner.

2.1. Initial Fracture Step Distance of Coal Seam Roof

Under the condition that the coal seam roof is simply supported on one side of the goaf and fixed on three sides, when it is in the limit suspension state, the three fixed edges form a negative bending moment area and the maximum main bending moment value, $M_a,$ is in the middle of the long fixed edge. The positive bending moment area is formed in the center of the goaf, and the maximum main bending moment, $M_c=\sqrt{M_x^2+M_y^2}$. According to the modified Marcus solution, the following can be obtained:

$$\left|M_a\right|={\displaystyle \frac{\left(1-{\mu }^2\right)\left(1+\mu {\lambda }_1^2\right)}{12\left(1+{\lambda }_1^4\right)}}·q{a}_1^2$$

(1)

$$\left|M_a\right|={\displaystyle \frac{\left(1-{\mu }^2\right){\lambda }_1^2\left(\mu +{\mu }_1^2\right)}{12\left(1+{\lambda }_1^4\right)}}·q{b}_1^2$$

(2)

According to the following relationship between bending moment and stress:

$$M_a={\displaystyle \frac{h^2{\sigma }_s}{6}}$$

(3)

By substituting Equation (3) into Equation (2), the relationship between the first fracture step distance of the coal seam roof under the condition of simple support and three-sided fixed support on one side of the goaf can be obtained as follows:

$$a={\displaystyle \frac{2h}{1-{\mu }^2}}\sqrt{{\displaystyle \frac{{\sigma }_s}{q}}·{\displaystyle \frac{2+{\lambda }^4}{4+3\mu {\lambda }^2}}}$$

(4)

Let $l_m={\displaystyle \frac{h}{1-{\mu }^2}}\sqrt{{\displaystyle \frac{{\sigma }_s}{q}}}$,

Then, Equation (4) is modified as follows:

$$a=\left\{\begin{array}{c}{b}^4\sqrt{{\displaystyle \frac{l_m^2}{2\left(b^2-l_m^2\right)}}}\\ {\displaystyle \frac{b}{l_m}}\sqrt{b^2-\sqrt{b^4-2l_m^4}}\end{array}\right.$$

(5)

where $\mu$ is the Poisson’s ratio of the rock stratum; $q$ is the weight of the rock stratum and its load; and ${\lambda }_1$ is the geometric shape coefficient of the goaf, ${\lambda }_1=a_1/b$, $a_1$ it is the advancing distance of the working face, and $b$ is the length of the working face. Moreover, $h$ is the thickness of the main roof strata, and ${\sigma }_s$ is the tensile strength of the main roof strata.

2.2. The Separation Layer Collapse of the Immediate and Main Roofs

Whether the immediate roof and the main roof fall off from the layer can be judged by the maximum disturbance.

The maximum disturbance of the main roof is given as follows:

$$y_{max}={\displaystyle \frac{\left(\gamma h_1+q_1\right)L_1^4}{384E_1{J}_1}}$$

(6)

where $q_1$ is the load applied to the main roof, $\gamma$ is the load per unit length of the main roof itself, $L_1$ is the initial caving step distance, $h_1$ is the thickness of the main roof, $E_1$ is the elastic modulus of the main roof, and $J_1$ is the section inertia moment of the main roof.

The maximum disturbance of the immediate roof is given as follows:

$$y_{max}={\displaystyle \frac{\sum h\gamma L_1^4}{384E_2{J}_2}}$$

(7)

where $\sum h$ is the thickness of the immediate roof, $E_2$ is the elastic modulus of the immediate roof, and $J_2$ is the section moment of inertia of the immediate roof.

According to the conditions of the separation phenomenon between the immediate roof and the main roof, that the following can be observed:

When ${\displaystyle \frac{\left(\gamma h_1+q_1\right)L_1^4}{384E_1{J}_1}}<{\displaystyle \frac{\sum h\gamma L_1^4}{384E_2{J}_2}}$, the immediate roof and the main roof are prone to separation and collapse.

By simplifying Equation (7), it can be approximated that when the thickness of the immediate roof is less than or equal to the thickness of the main roof, the immediate roof separation layer is prone to collapse.

3. Numerical Simulation

3.1. Overview of Fully Mechanized Caving Face

The N2306 working face of a mine in Shanxi has a length of 275 m, an average height of 6.31 m, an inclination angle of 2~6°, a mining height of 3.5 m, a caving height of about 2.8 m, and a mining/caving ratio of l:1.04. The cycle occurs six times a day, and the cycle progress is 0.8 m. The gas content of the coal seam is 1.81~9.68 m³/t, with an average of 5.49 m³/t. On the top of the coal seam is a false roof, the lithology is generally black carbonaceous mudstone and mudstone; the thickness is not stable, generally 0.25~0.44 m, only in local development; the immediate roof is composed of gray-black mudstone and sandy mudstone, and it is sometimes interbedded with sand and mud. The thickness is generally 1.0~7.0 m, with an average thickness of 5.35 m. The lithology of the main roof is medium-fine gray-white quartz feldspar sandstone, locally coarse sandstone, with an average thickness of 5.8 m, and the height of the caving zone in the goaf is about 30 m. The relevant parameters of the working face are shown in

Table 1. Parameters of N2306 working face coal seam roof.

Strata Name	Strata Thickness (m)	Bulking Coefficient	Poisson Ratio	Elastic Modulus (MPa)	Tensile Strength (MPa)	Rock Load (MPa)	Coal Seam Thickness (m)	Fracture Interval (m)	Falling Body Height (m)
Immediate roof	5.35	1.2	0.25	1.5 × 104	2.2	13.5	6.31	0.91	6.4
Main roof	5.8	1.2	0.25	3.5 × 104	7	13.7	6.31	6.35	13.4

. Combined with Equation (6) and Equation (7), it can be seen that in advancing the working face, the immediate and main roofs of the coal seam will separate and collapse. According to Equation (5), the length of the immediate roof suspension beam is 0.91 m, and the accumulation height of the immediate roof caving rock layer is 6.4 m. The length of the main roof suspension beam is 6.35 m, the accumulation height of the direct roof and the main roof caving rock mass is 13.4 m, and the height of the goaf fracture zone is about 30 m.

3.2. The Establishment of the Geometric Model and the Determination of Simulation Parameters

3.2.1. Establishment of the Geometric Model of Goaf

According to the theory of spatial zone division of rock mass deformation around the stope, the goaf can be divided into ‘three zones’ in both vertical and horizontal directions, shown in Figure 1.

Due to the limitations of the site, it is not suitable to carry out the actual measurement, so the height of the caving zone and the fracture zone are estimated by the empirical formula. The height of the caving zone in the goaf is affected by the thickness of the mining coal seam and the occurrence conditions and bulkiness of the roof rock, and it is usually three to five times the thickness of the coal seam. In the actual mining process, the height of the air-conducting zone and the caving zone can be expressed by Equation (8) and Equation (9), as follows:

$$h={\displaystyle \frac{100m}{1.6m+3.6}}\pm 5.6$$

(8)

where $m$ is the thickness of the mining coal seam.

$$h={\displaystyle \frac{m}{\left(k_p-1\right)cos\alpha }}$$

(9)

where $\alpha$ is the dip angle of coal seam and $k_p$ represents the bulkiness of the rock.

In the N2306 fully mechanized caving face, the thickness of the coal seam is 6.31 m, and the dip angle is 2~6°. According to Equation (9), the height of the caving zone is given as follows:

$$h={\displaystyle \frac{6.11}{\left(1.4-1\right)cos\left(4/180\right)}}=15.28 \mathrm{m}$$

(10)

In the N2306 working face, the roof and floor are based on mudstone and contain a small amount of fine-grained sandstone. The characteristics of the rock strata are treated according to the medium-hard rock strata, and the height of the air-conducting zone (the sum of the caving zone and the fractured zone) can be obtained as follows:

$$h={\displaystyle \frac{100\times 6.11}{1.6\times 6.11+3.6}}\pm 5.6=45.68\pm 5.6 \mathrm{m}$$

(11)

Therefore, the height of the caving zone in the N2306 working face is 15.28 m, the fracture zone is 30.4 m, and the air conduction zone is 45.68 m.

According to the actual situation of the site, combined with the need for numerical simulation, the geometric model of the goaf in the N2306 working face will be simplified. The details are as follows:

• The working face, goaf, intake and return air roadways, and tail roadway are all regarded as cuboids without considering the equipment of the working face;
• The goaf is divided into caving zone 1, caving zone 2, caving zone 3, fracture zone 1, and fracture zone 2 to distinguish different porosities;
• The length of the working face is 275 m, the width is 5.5 m, and the height is 3.5 m. The strike length of the goaf is 200 m, the dip length is 285 m (including the intake and return air roadways), and the vertical height is 40 m;
• The length, width, and height of the intake and return air roadways are 15 m, 5 m, and 3.5 m, respectively, and the length, width, and height of the tail roadway are 20 m, 4.4 m, and 3.5 m, respectively.

The geometric model of the working face is shown in Figure 2.

3.2.2. Determination of Simulation Parameters

Through the detailed investigation of the layout of the fully mechanized caving face, ventilation, gas occurrence, and the operation of the main equipment on the spot, based on reality, the relevant simulation parameters, including the rock porosity, the gas emission of each part of the goaf, the viscous resistance coefficient, and the boundary conditions are set as follows:

(1) Determination of porosity

As the fully mechanized caving face advances, the roof of the goaf begins to collapse, the rock falling behind the goaf is gradually compacted, and the expansion coefficient of the goaf from shallow to deep gradually decreases. According to the empirical formula, the bulking coefficient of the three parts of the caving zone can be taken as 1.5, 1.3, and 1.2. The porosity of the caving zone can be obtained from the bulking coefficient of 0.333, 0.231, and 0.167. The porosity of fracture zones 1 and 2 can be 0.05 and 0.02.

(2) Determination of gas emission (source term) in each part of the goaf

Using the source prediction method, it is predicted that the gas emission in the goaf of the N2306 working face is about 44.89 m³/min, the gas emission in the fracture zone is about 18 m³/min, and the gas emission in the caving zone is about 27 m³/min. Among them, gas emissions of fracture zones 1 and 2 are 13.54 m³/min and 4.45 m³/min, respectively. The gas source term of each part of the caving zone is calculated according to the broken expansion coefficient, in which the gas density is 0.716 kg/m³.

The total gas source term of the caving zone is given as follows:

$${\displaystyle \frac{2.78÷60\times 0.716}{200\times 285\times 16}}=3.51\times {10}^{-7} \mathrm{kg}/\left({\mathrm{m}}^3·\mathrm{s}\right)$$

(12)

According to the length and bulking coefficient of each part of the caving zone, the source terms of the three parts of the caving zone can be calculated according to the following weighted average method:

The gas source term of caving zone 1 is as follows:

$${\displaystyle \frac{1.5\times 20}{1.27\times 200}}\times 3.51\times {10}^{-7}=0.419\times {10}^{-7} \mathrm{kg}/\left({\mathrm{m}}^3·\mathrm{s}\right)$$

(13)

The gas source term of caving zone 2 is as follows:

$${\displaystyle \frac{1.3\times 80}{1.27\times 200}}\times 3.51\times {10}^{-7}=1.433\times {10}^{-7} \mathrm{kg}/\left({\mathrm{m}}^3·\mathrm{s}\right)$$

(14)

The gas source term of caving zone 3 is as follows:

$${\displaystyle \frac{1.2\times 100}{1.27\times 200}}\times 3.51\times {10}^{-7}=1.658\times {10}^{-7} \mathrm{kg}/\left({\mathrm{m}}^3·\mathrm{s}\right)$$

(15)

The gas source term of fracture zone 1 is as follows:

$${\displaystyle \frac{13.54÷60\times 0.716}{200\times 275\times 12}}=1.75\times {10}^{-7} \mathrm{kg}/\left({\mathrm{m}}^3·\mathrm{s}\right)$$

(16)

The gas source term of fracture zone 2 is as follows:

$${\displaystyle \frac{4.45÷60\times 0.716}{200\times 275\times 12}}=0.8\times {10}^{-7} \mathrm{kg}/\left({\mathrm{m}}^3·\mathrm{s}\right)$$

(17)

(3) Determination of viscous resistance coefficient

Due to the different degrees of rock collapse in the goaf, the permeability is also quite different. However, since measuring the expansion coefficient everywhere in the goaf is difficult, it is regarded as heterogeneous isotropic in the numerical simulation process and calculated according to the empirical formula.

The empirical formula for the permeability of each partition in the goaf is as follows:

$$k={\displaystyle \frac{D^2{n}^3}{180{\left(1-n\right)}^2}}$$

(18)

where $D$ represents the average particle size of the rock and $n$ is the porosity.

According to the empirical formula, the permeability of the natural accumulation zone, the load affected zone, and the compaction stable zone can be set to 1.5 × 10⁻⁷ m², 0.7 × 10⁻⁷ m², and 0.2 × 10⁻⁷ m². The permeability of fracture zone 1 and fracture zone 2 can be considered as 1.0 × 10⁻⁵ m² and 0.5 × 10⁻⁵ m², respectively. The viscous resistance coefficient is the reciprocal of the permeability of the goaf.

(4) Setting boundary conditions

The setting of boundary conditions mainly considers the intake wind speed and outlet pressure. The intake wind speeds were 1.0, 1.5, 2.0, 2.5, and 3.0 m/s, respectively. The turbulent k-epsilon and species transport models are selected to open the energy equation. The numerical simulation parameters of the goaf gas seepage law are shown in

Table 2. Simulation parameters of gas seepage law in goaf.

Name	Caving Zone 1	Caving Zone 2	Caving Zone 3	Fracture Zone 1	Fracture Zone 2
Porosity	0.333	0.231	0.167	0.05	0.02
Coefficient of viscous resistance (107 m2)	0.7	1.4	5	100	200
Methane source term (10−7 Kg/m3·s)	0.419	1.433	1.658	1.75	0.8

.

4. Numerical Results and Analysis

4.1. Numerical Simulation of Gas Seepage Law in Goaf Under U-Type Ventilation Mode

4.1.1. Before the Initial Weighting

From the beginning of mining to the formation of the initial weighting before the collapse of the main roof of the fully mechanized caving face, the main roof of the coal seam does not collapsed under the influence of its strength. The cracks in the roof and floor of the goaf are not developed, and effective communication with the adjacent coal seam is not formed. The mining face is less affected by the gas of the adjacent coal seam. At this time, the main source of gas in the working face is the gas emitted from the coal wall of the fully mechanized caving face, the gas desorbed from the coal falling in the working face and the goaf, and a small amount of gas from the seepage of the surrounding rock. Before the initial weighting, only the gas emitted from the coal wall and desorbed from the falling coal is considered the gas source.

According to the technical parameters of the N2306 working face coal cutter, the roof control distance of the working face is 5.16~5.96 m, the model is 5.5 m, the mining height of the working face is 3.5 m, and the tendency length is 285 m. In the model, the intake and return air roadway length is 15 m, and the section is 5 × 3.5 m². The first caving step distance of the main roof is 25~40 m, and the goaf is 30 m along the strike. The height of the goaf is the sum of the mining coal seam’s thickness and the immediate roof’s thickness, which is 8.5 m. The three-dimensional geometric model of the goaf in the fully mechanized caving face is shown in Figure 3.

The numerical simulation is carried out according to the parameters and boundary conditions of the goaf before the first weighting of the fully mechanized caving face. The velocity cloud map of the goaf is shown in Figure 4, the velocity vector map in Figure 5, and the gas concentration distribution map in Figure 6.

From Figure 4 and Figure 5, it can be seen that under the conditions of U-type ventilation, the ventilation airflow mainly flows through the working face from the return air roadway, and a very small part of the airflow bypasses the support between the working face and the goaf into the goaf. Affected by the support and the caving coal and rock mass, the airflow velocity entering the goaf decreases rapidly and tends to stop as the wind speeds away from the working face. When the airflow moves in the working face, the wind speed near the roof and floor is low, which is also one of the reasons for the high gas concentration near the roof of the working face.

According to Figure 6, before the initial weighting, the maximum gas concentration in the working face appears near the coal wall near the return airway. The maximum value is about 0.85% (blue circle area in Figure 6), the gas concentration in the upper corner of the working face is about 0.70% (orange circle area in Figure 6), and the gas concentration in the return airway is 0.80%. The gas concentration is higher in the range of about 130 m from the return air roadway along the tendency of the goaf, and the farther away from the return air roadway, the lower the gas concentration; the gas concentration near the roof and floor is high, and the gas concentration in the middle of the working face is low.

4.1.2. After the Initial Weighting

In the N2306 fully mechanized caving face, after the initial weighting, the coal seam roof falls under gravity, and the bending subsidence zone produces many cracks. The caving coal and rock mass is under the influence of fragmentation, self-weight, and external force, so the goaf of the working face presents a strip distribution. According to the distribution law of each strip, a three-dimensional geometric model of the goaf of the working face with a tendency of 285 m, a strike of 200 m, and a height of 40 m is established, as shown in Figure 7.

(1) Simulation of air leakage in goaf after initial weighting under U-type ventilation mode

The intake velocity is set to 2.85 m/s, the outlet is the free outlet, and the goaf is the porous medium. The multi-component transport model is used to simulate the goaf after the first roof pressure, and the velocity cloud map of the goaf is shown in Figure 8.

It can be seen from the figure that under the ventilation airflow in the working face, there is also a weak airflow in the goaf. Its influence range is concentrated in the lower and upper corners, and the influence range is greater in the upper corner. When the wind speed of the working face is 2.85 m/s, the wind speed of the goaf at the upper corner is 0.1 m/s in the area of 5~10 m from the working face, and the wind speed is 0.05 m/s in the area of 15~18 m from the working face. The farther away from the working face, the lower the wind speed. The wind speed also shows a decreasing trend in the vertical height. At 20 m from the floor, the wind speed has been below 0.01 m/s.

(2) Numerical simulation of gas concentration distribution law in goaf under different wind speed conditions of working face

By adjusting the boundary conditions, the distribution law of gas concentration in the goaf is simulated under wind speed conditions of 1.0 m/s, 1.5 m/s, 2.0 m/s, 2.5 m/s, and 3.0 m/s in the fully mechanized caving face. Figure 9, Figure 10 and Figure 11 depict 1.5 m/s, 2.0 m/s, 2.5 m/s of wind speed, respectively, from the floor z = 3 m, and the 10 m goaf gas concentration distribution cloud.

It can be seen from the figures that the gas concentration can be divided into three stages in the whole goaf trend as follows: the gas concentration increases slowly, and the gas concentration is low in the first 30 m near the working face (location of the blue line in the picture); in the range of 30 m~160 m (blue line and orange line area in the picture), the gas concentration increases rapidly, and there is obvious stratification; after 160 m (orange line position in the picture), the growth of gas concentration gradually slows down and eventually stabilizes. The main reason for this phenomenon is the effect of the wind speed of the working face. The air leakage in the goaf near the working face is high, and the air leakage gradually decreases with the depth increase. On the other hand, because the porosity in the goaf gradually decreases from shallow to deep, the bulking coefficient in the re-compacted area is close to 1, which leads to its accumulation in the deep part of the goaf, and the concentration can reach more than 90%.

Along the dip direction of the working face, the gas concentration distribution in the goaf shows a trend of increasing gas concentration from the intake side to the return side, especially in the gas concentration on the return side of the working face. The reason for this phenomenon is that in the early stage of coal release, the immediate roof and the main roof are prone to the cantilever beam phenomenon under influence of their strength, and the caving rock in the goaf accumulates freely, the porosity is great, and the air leakage wind speed is high. Part of the gas accumulates on the return side driven by the airflow, and with the increase in the depth of the goaf, the air leakage airflow gradually disappears, and its influence on the gas concentration gradually disappears.

In the vertical direction of the goaf, the gas concentration gradually increases from the bottom plate, and the gas concentration in the fracture zone can even reach 100%. The main reason for this phenomenon is that the density of gas is low relative to the air, and it migrates to the upper space of the goaf under gravity and buoyancy. At the same time, due to the large difference in air leakage flow between the caving zone and the fracture zone in the goaf, the gas in the caving zone is easily taken away by the airflow, while the accumulated gas in the fracture zone cannot be taken away so that the gas forms a gradient distribution in the vertical direction of the goaf.

Figure 12 is the change curve of gas concentration in the upper corner of the working face when the wind speed is 1.0 m/s, 1.5 m/s, 2.0 m/s, 2.5 m/s, and 3.0 m/s, respectively, under the U-type ventilation mode of the working face.

It can be seen from Figure 12 that the change in wind speed in the working face greatly influences the gas concentration in the upper corner. However, in practice, this should be determined by the amount of gas emission from the working face. For the fully mechanized caving face with a high absolute gas emission, reducing the gas concentration in the upper corner and the return air side is limited to only increasing the wind speed (pressure difference). As the pressure difference increases, the air leakage volume also increases, removing more gas from the goaf so that the gas concentration in the return air side and the upper corner does not decrease significantly or even increase.

4.2. Numerical Simulation of Gas Seepage Law in Goaf Under U + I-Type Ventilation Mode

Based on the process parameters of the N2306 working face, the working face is designed in accordance with a U + I-type ventilation mode. A three-dimensional geometric model of the goaf in a fully mechanized caving face is established, which is 15 m in length, 5 m in width, and 3.5 m in the height for the intake and return air roadway; 15 m in length, 4.4 m in width, and 3 m in height for the tail roadway; 10 m from the coal seam floor to the bottom of the tail roadway and 20 m from the horizontal distance for the return air roadway; 275 m in length, 5 m in width, and 3.5 m in height for the working face; and 285 m in length, 200 m in width, and 40 m in height for the goaf, as shown in Figure 13.

Under the condition that other parameters remain unchanged, the intake roadway is selected as the speed intake, the wind speed is 2.5 m/s, and the return air roadway and the tail roadway are free exits. The absolute gas emission in the goaf is 44.89 m³/min, respectively, to simulate the distribution and migration law of gas concentration in the goaf when the U + I-type ventilation mode is simulated.

Figure 14 and Figure 15 show the gas concentration distribution cloud maps in the horizontal and vertical directions of the goaf under the U + I-type ventilation mode.

It can be seen from the figure that compared with the U-type ventilation mode, some gas is discharged with the ‘I’ roadway in the U + I-type, and the gas concentration in the upper corner and the support influence area can be effectively reduced in the horizontal direction. In the vertical direction, the stratification of gas concentration distribution in the deep goaf is more obvious. The reason for this phenomenon is that the gas in the shallow caving zone of the goaf is eliminated with the airflow of the tail roadway. The deep wind speed in the goaf is very low, and the gas mainly depends on gravity and free diffusion. In addition, the viscous resistance coefficient of the fracture zone is high, the porosity is low, and the gas accumulation gradually increases the concentration.

4.3. Numerical Simulation of Gas Seepage Law in Goaf Under U + L-Type Ventilation Mode

Based on the process parameters of the N2306 working face, the working face is designed as a U + L-type ventilation mode. A three-dimensional geometric model of the goaf in a fully mechanized caving face with the same size of intake and return air roadways and tail roadway is established, with a length of 15 m, a width of 5 m, and a height of 3.5 m. The tail roadway is parallel to the return air roadway, and the horizontal spacing is 15 m. The working face has a length of 275 m, a width of 5 m, and a height of 3.5 m. The goaf has a length of 285 m, a width of 200 m, and a height of 40 m, as shown in Figure 16.

Under the condition that other parameters remain unchanged, by changing the ratio of the air volume of the return air roadway and the ‘L’ roadway, the distribution and migration law of the gas concentration in the goaf are simulated when the air volume ratio is 2:1, 3:1, and 4:1 in the U + L-type ventilation mode.

When the air volume ratio of the return air roadway and the ‘L’ roadway is 2:1 (the wind speed of the return air roadway is 1.33 m/s, and the wind speed of the ‘L’ roadway is 0.67 m/s), 3:1 (the wind speed of the return air roadway is 1.5 m/s, and the wind speed of the ‘L’ roadway is 0.5 m/s), and 4:1 (the wind speed of the return air roadway is 1.6 m/s, and the wind speed of the ‘L’ roadway is 0.4 m/s), the gas concentration distribution cloud map of the goaf is shown in Figure 17, Figure 18 and Figure 19.

It can be seen from the figure that due to the existence of the ‘L’ roadway, the gas in the goaf flows to the ‘L’ roadway with the airflow, so the concentration of the coal mining face and the return air roadway is low, and the gas concentration of the outer tail roadway is high. The influence of the tail roadway on the gas flow field in the goaf of the working face is mainly reflected near the upper corner. The farther away from the tail roadway, the smaller the gas concentration change. Therefore, the outer tail roadway can effectively reduce the gas concentration in the upper corner and alleviate the gas overrun problem. By comparing the gas concentration distribution in the goaf with the air volume ratios of 2:1, 3:1, and 4:1 in the return air roadway and the outer staggered tail roadway, it can be seen that the concentration in the upper corner is significantly lower than that in the U-type ventilation. However, the change in air volume ratio has a limited impact on the gas concentration in the upper corner. If conditions permit, increasing the air volume of the outer staggered tail roadway will help reduce the gas concentration in the gas tail roadway and ensure safety.

5. Discussion

Using the three ventilation modes of U, U + I, and U + L, and the three wind speed conditions of 1.5 m/s, 2.0 m/s, and 2.5 m/s in the intake air roadway, the gas seepage law in the goaf was simulated. Figure 20, Figure 21 and Figure 22 depict the curves of gas concentration change along the goaf strike at 0.5 m from the working face under the three ventilation modes.

It can be seen from the figure that compared with the U-type mode, the gas concentration in the upper corner of the goaf is significantly reduced under the U + I and U + L ventilation modes. The existence of the tail roadway is beneficial in reducing the gas concentration in the upper corner of the goaf in the fully mechanized caving face. However, the influence of the tail roadway on the overall gas concentration distribution in the goaf is limited. Therefore, for the hidden danger control of gas disasters in high gas fully mechanized caving faces, the ventilation mode should be determined according to the specific situation of the mine to reduce the occurrence of gas accumulation accidents in the upper corner of the working face.

Figure 23 represents the U + L ventilation mode. Under the three wind speed conditions, the gas concentration change curve is 80 m from the return air roadway along the goaf.

It can be seen from Figure 23 that the gas concentration gradually increases along the goaf, while the change rate gradually decreases, and the change tends to be gentle. When it is close to the ‘asphyxiation zone’, the gas concentration tends to be stable. This trend is mainly because when the airflow in the working face enters the working face from the air intake, part of the airflow penetrates the goaf, and the gas in the goaf continuously oozes from the goaf under the influence of the airflow, thus reducing the gas concentration. As the goaf is away from the working face, the influence of the working face airflow on the goaf gradually weakens until it disappears. Therefore, the influence of the airflow on the goaf is also weakened or may even disappear. As the goaf is away from the working face, the gas concentration gradually increases and fills the entire ‘asphyxiation zone’.

6. Conclusions

According to the technological characteristics of the fully mechanized caving face in a mine in Shanxi Province, the numerical simulation of gas seepage law in the goaf of a fully mechanized caving mining is carried out. The mathematical model of the gas seepage law in a fully mechanized top coal caving mining is verified, and the main factors affecting gas accumulation in the working face are determined. The simulation results are significant in guiding future studies of gas control technology. The main conclusions are as follows:

• The gas concentration distribution in the goaf shows a similar law for different ventilation modes (including U-type, U + I-type, and U + L-type), affected by wind speed and rock porosity. The gas concentration in the whole goaf can be divided into the following three stages: in the first 30 m near the working face, the gas concentration increases slowly, and the gas concentration is low; in the range of 30 m~160 m, the gas concentration increases rapidly, and there is obvious stratification; and after 160 m, the growth of gas concentration gradually slowed down and eventually stabilized. Along the dip direction of the working face, the gas concentration distribution in the goaf shows a trend of an increasing gas concentration from the intake side to the return side, especially in the gas concentration on the return side of the working face. In the vertical direction of the goaf, the gas is distributed in a gradient. The gas concentration gradually increases from the bottom plate, and the gas concentration in the fracture zone can even reach 100%;
• Different ventilation methods have different scopes of application. When there us less desorption gas in the coal seam, it is easy to adopt the U-type ventilation mode and the U + I or U + L ventilation modes when there is more desorbable gas, or the mining intensity is high. Compared with the U-type ventilation mode, the U + I and U + L-type ventilation modes are more conducive to reducing gas seepage from the goaf to the upper corner of the working face. That is, applying the tail roadway can solve the problem of gas accumulation in the upper corner caused by eddy currents. To a certain extent, this can reduce the gas concentration in the upper corner of the goaf of the fully mechanized caving face. However, in general, the influence of the tail roadway on the overall gas concentration distribution in the goaf is limited;
• The wind speed at the working face is one of the most important factors affecting the gas concentration in the upper corner of the working face’s mining area. In the case of U-type ventilation, i.e., when the absolute gas outflow from the workface is low, the change in the wind speed at the workface has a greater impact on the gas concentration in the upper corner, and by increasing the wind speed appropriately, the gas concentration in the upper corner can be significantly reduced. In the case of U + I- and U + L-type ventilation, i.e., when the absolute gas outflow from the workface is high, there are great limitations in reducing the gas concentration in the upper corner and the return side by only increasing the wind speed (pressure difference), because when the pressure difference is increased, the amount of air leakage increases, which will take away more gas from the mining area, thus making the reduction in the gas concentration in the return side and the upper corner insignificant or even increasing.