# Recapitulation and Prospect of Research on Flow Field in Coal Mine Gob

Authors: Rui Zhang; Jianwei Cheng; Zui Wang; Zhenlu Shao
Journal: Shock and Vibration
 (2021)

---

## Abstract

Coal mine gob, mined-out areas in underground coal mines, often accumulates explosive methane-air mixtures that pose a deadly hazard to miners. A good understanding of the flow field in a sealed coal mine area is crucial in preventing and minimizing accidents associated with mine combustible gases and also for planning and implementing a mine rescue strategy. In recent years, the research on the flow field in the gob has changed from qualitative research in the past to quantitative research. This paper synthesizes the research results of flow field in gob in recent 40 years, covering the permeability of quarried areas, the airflow simulation in quarried areas, and the influence of ventilation parameters and geohydrological conditions on the flow field. Firstly, the overburden failure mechanism and fracture development characteristics of the mine gob, the distribution of porosity and permeability in the gob, and the relationship between them are introduced. Secondly, the development of research methods and numerical models used to study the flow field in mine gob is discussed. The distribution of the flow field in the gob under different conditions is expounded. Thirdly, the research on the prevention and control of fire and explosion risks in the gob is discussed. Finally, the problems to be solved in such research direction are addressed and suggestions are put forward.


---

## Body

## 1. Introduction

For a long time, the problem of gob disaster has been the focus of most researchers. The gob is the abandoned space after coal mining, and its internal medium is the falling gangue and residual coal, which is rich in pores and cracks. Gas and air are distributed in various areas of gob in the form of mixed gas with different concentrations. Due to the wind pressure of air leakage flow in working face, the impact of roof failure, the power of gas emission, the fire pressure of spontaneous combustion of residual coal and other external forces, the mixed gas flows in the gob, and the concentration of components also change, so the flow field distribution is very complex [1]. In addition, the gob is prone to spontaneous combustion of residual coal, which provides a long-term fire source for gas explosion and causes secondary gas explosion, resulting in serious consequences of accident disasters. In-depth study of the flow field distribution in the gob is not only conducive to mastering the flow of mixed gas but also reveals the mutual transformation law of gas combustion and explosion in the gob, which is of great significance for the prevention and control of coal spontaneous combustion and gas explosion disasters.The viewpoint of flow field in gob was first proposed in 1978 [2], which was formed with the development of computational fluid dynamics and computer technology. Based on rock mechanics, porous media seepage dynamics, and numerical simulation, the research object is the space outside the underground roadway or the coal and rock area, which is aimed at directly revealing the nature and mechanism of disasters. Practice has proved that it is an effective analysis method to explore the gob gas distribution law and prevent disaster accidents of the gob.This paper will start from three aspects: the permeability of the gob, the theory of quarry flow field simulation, and the influence of environmental factors on the flow field distribution. It will explain the methods of setting the porosity and permeability of the gob in numerical simulation, compare different analysis methods of the quarry flow field problem, combine the quarry simulation experiments with relevant materials and the simulation results, and sort out the influence of different environmental factors on the gas flow situation in the gob. The research on the prevention and control of fire and explosion risks in the mined area is also analyzed, and finally, suggestions are made for future research directions of the mined area flow field.


## 2. Permeability of the Gob

### 2.1. Theory of Stope Fracture Formation

From the process of mineral pressure during mining when the rocks on the roof are risen and fallen, to the compaction of the back of the gob, the medium inside the gob is a porous medium with the nonuniform medium. The pore size distribution is related to the working face height, the size and arrangement of the rock, the lithology of the coal bed and adjacent layers, the original stress and the mining stress, etc. There are two different characteristics of voids in the gob, i.e., mined voids and original voids, which differ significantly in their characteristics, and generally speaking, the size and permeability of mined voids are much larger than those of original voids [3]. In general, the size and permeability of the active pore space are much larger than those of the primary pore space. The active pore space is interconnected throughout the entire area of the extraction space and is not only spatially unevenly distributed but is also influenced by the extraction pressure during the extraction process.


#### 2.1.1. The Mine Pressure Hypothesis

The formation of mining overburden fractures is associated with various mine pressures and deformation damage to the surrounding rock within the quarry. From the beginning, mine pressure phenomena have been observed and different explanations have been proposed based on the complex nature of the rock mass. These explanations, which reveal the intrinsic connection between mine pressure phenomena, are called mine pressure hypotheses. The pressure arch hypothesis was proposed by the Germans Haacke and Gulitzel in 1928. In the same vein as the pressure arch hypothesis, the cantilever beam hypothesis was proposed by the German Sledke and later supported by the British Ferriday and the former Soviet Germain, among others [4]. By the 1950s, the articulated block hypothesis and the prefracture hypothesis emerged. The preformed fissure hypothesis was first proposed by a Belgian geologist. In the late 1970s and early 1980s, scholars represented by academician Qian based on the articulated block hypothesis and the rock performed fissure hypothesis, and through on-site observation of the internal movement of rock formations and extensive production practice, the masonry beam theory of rock structure was proposed and developed into a hypothesis. In recent years, as the scientific study of rock formations continued to deepen and to understand the broader issues in rock formation activities, Academician Qian et al. [5] further proposed the key layer theory of rock control based on the masonry beam.


#### 2.1.2. The “O” Circle Theory and the Critical Layer Theory [6–10]

Multiple layers of rock of varying thickness and strength exist above the direct top. It has been proven that one to several of these thick hard rock layers play a major controlling role in the overlying rock activity of the quarry. The rock layers that control the activity of the overlying rock layers in the quarry, either locally or up to the surface, are referred to as key layers. The former is referred to as the subcritical layer and the latter as the main critical layer. The key layer is generally a relatively thick and hard rock layer with the following characteristics: geometric features, with thicker single layers relative to other similar rock formations.Lithological characteristics, relatively hard compared to other rock formations, i.e., higher modulus of elasticity and higher strength; deformation characteristics, where the subsidence deformation of the key layer is synchronized with the amount of subsidence of all or part of the overlying rock layers; fracture features, where the fracture of a key layer will result in the simultaneous fracture of all or 1ocal overlying rock layers, causing movement of the rock over a 1arger area; the bearing characteristics, in the form of a slab structure before the key layer breaks as the bearing body for all or part of the rock, become a masonry beam structure after the break if the S-R stability conditions for rock block structures are met, succeeding as the bearing body.The practice has shown that seam activity leads to the creation of a large number of mining fissures, and for coal seams with low permeability, the formation of mining fissures can lead to an increase in permeability of tens to hundreds of times, which creates conditions for the extraction and transport of gas from the coal seam. Two types of fissures are formed in the overlying rock after coal mining: one is the off-bed fissures, which are along-bed fissures that appear between the bed and the seam as the seam sinks, and which cause the coal seam to expand and deform and unload the gas and cause the unloaded gas to gush out along the off-bed fissures; the other is the vertical fracture fissures, which are through-bed fissures that form as the seam sinks and breaks, and these fissures enable the gas to rise and float continuously, forming a rising float.Along the direction of workface advancement, the dynamic distribution of deviations under the key layer shows a two-stage development pattern. Before the initial breakage of the key layer, the amount of deviations increases continuously as the workface advances, and the maximum amount of deviations is located in the middle of the extraction zone. After the initial breakage of the key layer, the key layer tends to be compacted in the middle of the extraction zone, while a delamination zone is still maintained on each side of the extraction zone. The interlayer zones in the upper and lower chutes are connected through the cuttings and the interlayer zone on the working face side, forming an interconnected ring of interlocking interlayer fracture development around the mining void. Generally speaking, the geometry of the open-cut eye, the upper and lower chute, and the working face side are rectangular; that is, the mining void area forms the shape of the moving rock layer as shown in Figure1, which is called the “O” shaped circle of the mining fissure, the porosity of this area is large, easy for the wind flow, and is the main wind leakage channel at the early stage of the mining of the back mining face. With the back mining workface gradually far away from the position of the cutting eye, cutting eye position of the loose triangle by the back mining workface at both ends of the role of wind pressure gradually reduced, when the workface advance to a certain distance, if there are no other sources of wind leakage in the gob and the role of wind leakage sink, the location is no longer the main wind leakage channel in the gob.Figure 1
Diagram of the “O” ring. (a) Tendency profile. (b) Plan view.
(a)
(b)


#### 2.1.3. The Zoning of Overburden Movement Failure in Gob

When the coal is fully extracted, the overlying rock layer loses its supporting role and gradually breaks up, collapses, delaminates, fractures, and deforms, showing a process of sinking of the overlying rock layer from static to dynamic to static.The theory of “three horizontal zones and three vertical belts” was proposed by our academician Qian [11]. The overburden fractures formed by mine pressure and complex stresses are divided in different directions, forming “three zones” in the vertical direction, i.e., the caving band, the fracture band, and the bending subsidence belt, which are similar to the “three zones” in the longwall working face as suggested by Karmis et al. [12]. The basic structure and zoning of the overburden in the gob is shown in Figure 2.Figure 2
Zoning of the overlying rock layers of the quarry and delineation of the mining void. (A) Zone of influence of coal wall support. (B) Zone of rock departure. (C) Zone of recompaction. (a) Zone of natural accumulation. (b) Zone of load influence. (c) Zone of compaction and stability. (I) Caving band. (II) Fracture band. (III) Bending subsidence belt.(i) Caving Band. The overburdened rock in the extraction area has been affected by complex mining, and the rock body at the direct top is completely fractured, resulting in a random rock block shape size and arrangement. The medium in the free-fall zone is loose, isotropic, and porous, and it is the main medium in the flow field of the extraction zone. It is generally considered that the coefficient of fragmentation and expansion in this area is around 1.5. As the workings continue to advance, the first blocks to fall are compacted by the subsequent collapse and the coefficient of fragmentation is reduced by approximately 31%.(ii) Fracture Band. When the underground workings are mined under certain conditions in time and space, the overlying rock layer emerges in a regular arrangement after fracturing, which is the main characteristic of the fissure zone. The fissures are mainly divided into two types according to the different forces: one is subject to bending and tensile stress, which is approximated as a fissure in the direction of the vertical rock layer; the other is a fissure in which the rock layer and the rock layer are subject to staggered displacement by mineral pressure, which is approximated as a fissure parallel to the rock layer. As mining continues, the gob gradually increases and the fissure zone is extended in the upper part of the overburden. After a certain depth has been reached in the mining space, the extension stops at the maximum height.(iii) Bending Subsidence Belt. The extent of the area is defined by the top of the rift zone as the lower boundary and the ground surface as its upper boundary. The pressure on the rock within the curved subsidence zone is less than that described by the riser zone, and the rock movement appears not as fragmentation or randomness, but as displacement of the entire area. Perpendicular to the surface, there is a small amount of subsidence displacement. However, in harder rock formations, dislocation fractures are very likely to occur due to misalignment movements between rock layers.Gao [13] established a “four-zone” model of rock movement and divides the overlying rock into four zones according to the mechanical structure characteristics after its destruction, namely, fracture zone, separation zone, bending zone, and loose alluvial zone. Singh and Kendorski [14] found that the area of the collapse zone and the fracture zone depends on the thickness of the coal seam, the strength of the roof, and the characteristics of the rock combination.


### 2.2. Study of Porosity Distribution Characteristics in the Gob

In the process of coal mining, due to the influence of mining, the overlying rock layer is damaged, which will not only produce new fissures but also lead to the enlargement of the primary fissures, which will in turn change the porosity of the rock body, and the result leads to changes in the permeability of the rock body. Mastering the characteristics of the porosity distribution in the gob is the key to studying the flow field in the gob. It is impossible to accurately simulate the flow field in the gob by simple porosity analysis alone; therefore, it is important to improve the accuracy of the measurement of the pore distribution in the gob for the study of the flow field in the gob.


#### 2.2.1. Fractional Expansion Coefficients and Distribution Functions in the Gob

As time passes, the overlying rock layer gradually sinks under its own weight, the force on the free-fall zone gradually increases, the rock layer is gradually compacted, the porosity of the free-fall zone gradually decreases, and the initial fragmentation and expansion coefficient is larger compared with that after the fall, resulting in a lower height of the flow field in the gob. The height of the zone can be calculated by the following formula [15]:(1)M=H−SK−1,where M is the height of the fall zone (m); S is the actual settlement value of the lower rock beam (old top) in the fissure zone (m); H is the mining height (m); and K is the coefficient of fragmentation and swelling of the fallen rock body.As can be seen from the above equation, if the fragmentation effect of the rock fall is not taken into account, i.e.,K1, the value of M will tend to infinity, indicating that the overburden of the quarry will fall all the way to the surface, resulting in a wide range of estimates of the fall zone. Therefore, the fragmentation and swelling effect of the rock fall has an important influence on the development pattern of deformation and damage of the overlying rock in the quarry, and the fragmentation and swelling characteristics of the rock fall should be fully considered in numerical calculations to make the numerical results more realistic [16].The volume of the rock after crushing should be increased compared with the whole state, and this property is called the crushing and swelling of the rock. The crushing and swelling coefficient is the ratio between the volume of the rock in the loose state after crushing and the volume of the rock in the whole state before crushing [17]. Practice and research have shown that after the collapse of the roof, broken rocks form a haphazard pile of loose bodies in the gob, and depending on the mechanical properties of the rocks, the size and arrangement of the broken rocks vary, and so does their coefficient of fragmentation and expansion, which is generally 1.10–1.40. Fragmentation coefficient and residual expansion coefficient of common coal measures rock are shown in Table 1. In general, hard rocks are broken in large pieces and arranged neatly, and their coefficient of fragmentation and expansion is small, generally 1.10–1.20. The coefficient of fragmentation and swelling of the soft rock layer is large, generally 1.30–1.40, and the coefficient of fragmentation and swelling of the medium-hard rock layer is between the hard and soft rock layers, generally 1.20–1.30 [18].Table 1
Coefficient of fragmentation and residual fragmentation of common coal rocks.
Lithology
Crushing and swelling factorKP
Residual swelling factorKP′
Crushed coal
<1.2
1.05
Muddy shale
1.4
1.1
Sandy mudstone
1.6∼1.8
1.1∼1.15
Sandstone
1.5∼1.8
1.03∼1.1
Medium-hard sandstone
1.3∼1.5
1.03∼1.08Palchik [19] used vertical boreholes to probe the height of the fallout zone of the medium-weathered, strongly weathered, and deeply buried hard overburden rocks in the shallow buried gob, and the theoretical calculation and observation results showed that the fragmentation and swelling coefficients of the strongly weathered rocks ranged from 1.06 to 1.165, and the fallout ratios ranged from 6.07 to 15.6, while the fragmentation and swelling coefficients of the medium-weathered rocks ranged from 1.09 to 1.24, and the fallout ratios ranged from 4.1 to 11.25. The coefficient of fragmentation and swelling of hard overburden rocks can be 1.38 times of that of strongly weathered rocks, while the riser ratio is only 0.26 times of that of strongly weathered rocks, which leads to the conclusion that the higher the strength of rocks, the higher the coefficient of fragmentation and swelling, and the smaller the height of riser zone.In order to study the compaction characteristics of fractured rocks in the coal bed roof, Su et al. [20] sieved sandstone, sandy mudstone, and mudstone in the roof of the Xin’an coal bed in Yima into 5 block sizes (0∼5, 5∼10, 10∼15, 15∼20, and 20∼25 mm) and then took 20% of each of the above 5 block sizes and mixed them as mixed block sizes (0∼25 mm), and in order the compaction tests were carried out on a hybrid electrohydraulic servo rock mechanic test system using a home-made device. The stress-strain relationships for the three types of crushed rock compaction tests were obtained, and the effects of rock strength, block size, and compaction stress on the compaction characteristics of the crushed rock were analyzed.It can be seen from Figure3 that the breaking expansion coefficient of each rock increases as the block diameter increases. The relationship between the crushing coefficient and the block diameter of the sandy mudstone crushed rock in Figure 5 can be described by a logarithmic relationship, which is different from the linear function between the crushing coefficient and the block diameter of the coal rock in the study of Miao et al. [21], and may be related to the range of the block diameter of the crushed rock; analyzing Figure 4, the residual crushing coefficient of the crushed rock after experiencing compaction is not much related to the rock strength and the block diameter. The reduction of the crushing and swelling coefficient is related to the block diameter, and the larger the block diameter, the greater the reduction.Figure 3
Coefficient of fragmentation and swelling of rock as a function of block diameter.Figure 4
Residual coefficient of fragmentation and swelling of crushed rock as a function of block diameter.Figure 5
Crushing and swelling coefficient of sandy mudstone rubble as a function of block size.Using similar material simulation tests, Zhang [22] developed two different models to study the dynamic 1aw of fracture and swelling of the mined rock mass during longwall mining and partial mining of small strips, respectively. The vertical fracture expansion of the modeled rock mass was approximated to replace the volumetric fracture expansion to analyze the dynamic fracture expansion and compressibility of the mined rock mass. The vertical fragmentation coefficient is defined as the ratio of the distance between two adjacent model points in the vertical direction before and after deformation, i.e.,(2)K=hn−n+1′hn−n+1,where hn−n+1′ and hn−n+1 are the vertical distances between two adjacent measurement points n and n+1 before and after mining, respectively, m.The rock fragmentation coefficient of different areas under full and partial mining conditions was obtained in relation to the distance from the coal seam and the observation time; that is, the rock fragmentation coefficient of the observation line located near the boundary of the gob was the largest during full mining of the longwall; the rock fragmentation coefficient of each observation line varied greatly during the mining process and then leveled off soon after the end of mining; the rock fragmentation coefficient of partial mining of the strip was much smaller than that of full mining of the longwall. The coefficient of rock fragmentation and swelling of partial mining is much smaller than that of longwall mining, the coefficient of rock fragmentation and swelling of partial mining varies with the height of waves, and the compressible amount of coal seam area is larger during loading. These understandings are important to explain the residual settlement of the surface in the gob.The porosity of a fractured rock mass is a function of its coefficient of fracture and expansion, and the study of the coefficient of fracture and expansion is the basis for calculating porosity. The porosity of a fractured rock mass can be expressed by the ratio of the pore volume to the total volume of the rock mass in the fractured state, and according to the definition of the porosity of a fractured rock mass and the coefficient of fracture and swelling, the following relationship exists between the two:(3)n=1−1Kρ,where n is the porosity of the crushed rock and Kρ is the crushing and swelling coefficient of the crushed rock.Considering the original porosity of the coal seam roof rockS0, then the total porosity of the rock in the free stacking state is(4)nm0=1−1Kρ+S0,where nm0 is the total porosity of the emergent rock in the free-stacked state and S0 is the original porosity of the top slab.When the workings are finished, as time passes, the rock in the fall zone tends to compact under its own weight and the overlying load, the coefficient of fragmentation becomes smaller, and the final remaining coefficient of fragmentation is called the residual coefficient of fragmentation.The coefficient of fragmentation and swelling is spatially distributed and follows a certain pattern from the start of the fall to the final compaction of the mined area. Li [23], based on the negative exponential decay law, ignored the influence of the shape of the extraction zone, specifically with respect to the distance of a boundary L, and considered that along the direction perpendicular to L, there is(5)Kp,l=Kρ′+KP0−Kρ′e−ald,where KP0 is the initial crushing and swelling factor and Kρ′ is the central compaction fragmentation and swelling coefficient; al is the decay rate; according to the law of mineral pressure, point d is the distance from L.The geometry of the extraction zone does not affect the distribution of the crushing and swelling coefficient, and the distribution pattern of the rock fall in the extraction zone is complex, so the extraction zone is a nonuniform medium. For the whole three-dimensional space, the rock fall fragmentation coefficient is a function of the location of the extraction zone, i.e., the distribution of the fall fragmentation coefficient is(6)Kpx,y=maxKp,l.In Li’s study [24], based on equations (5) and (6), the distance of the points within the gob from the working face and the boundary is reflected in the spatial distribution function of the fragmentation and swelling coefficient of the gob as(7)Kpx,y=Kp,min+Kp,max—Kp,mine−a1d11−e−δa0d0δ<1,where δ is the adjustment function that controls the distribution shape of the model. Kp,max is the initial crushing and swelling coefficient, i.e., the maximum value. Kp,min is the compaction coefficient of fragmentation and swelling, i.e., the minimum value. a0.a1 is the decay rate, and d0.d1 is the distance between any point and the boundary and the working surface, m.


#### 2.2.2. Porosity Distribution Model

The voids within the fractured rock scatter of a fall zone are usually irregular and disordered [25, 26]. These voids cannot be described and studied using conventional Euclidean geometry. Numerous studies have shown that the shapes of the various sizes of rock masses produced by rock fragmentation have a fractal structure: i.e., individuals of different sizes have statistically significant self-similar characteristics. Further studies have shown that the voids between these broken rock masses and the piles of broken rock masses are also fractal in structure [27].Xia and Huang [28] used the Menger sponge fractal model to study the void distribution characteristics of the rock scatter in the fall zone. The void distribution model of the adventitious scatter within the adventitious zone based on the Menger sponge fractal model is shown in Figure 6, with the similarity ratio of(8)t=yx=1k,where y is the volume of the first level transformed small cubic.Figure 6
Menger sponge model.y can be considered as the generating element of a Koch curve of dimension:(9)D=lnNln1/t.Bringing equation (9) into equation (8), we can obtain the void volume fraction dimension of the rock scatter in the bubble fall zone Dk:(10)Dk=lnk3−nlnk.Based on this model, the total fractal void ratioP of the rock bulk in the fall zone was obtained as(11)P=V−M/ρV=1−ρrmax3−Dρ0rmax3−D−rmin3−D,where ρ0 is the density of the rock fall, kg/m3; V is the total fractal volume of rock bulk within the fall zone, m3; M is the mass of the fallen rock bulk, kg; ρ is the density of the crushed rock bulk, kg/m3; and R is the diameter of individual rock mass, m.Wang et al. [29] divided the extraction zone into free accumulation zone and pressure accumulation zone along the working face advance direction, as shown in Figure 7.Figure 7
Schematic diagram of the state of rock accumulation in the extraction area.Derived from the theory of surrounding rock movement in the stope,(12)hx=m′−Sm1−e−aLf/Lb+1−Kf∑h′,sm=m′+1−KR∑h′,3b−1+5b−1b−12ab=1,where x is the strike distance of a point in the gob from the working face, m; h(x) is the gap between the rock fall and the roof at x, m. Lf is the length of the free accumulation zone, m; ∑h′ is the equivalent direct roof thickness, m; m′ is the equivalent mining height, m. sm is the displacement of the rock seam after it has moved and stabilized, m; L is the distance of the basic stabilization point of the rock seam from the working face, usually about 50∼60m; a and b are coefficients that vary with the distance of the rock seam from the coal seam and the lithological characteristics. Kf is the coefficient of fragmentation of rock in the free stacking zone, and KR is the coefficient of residual fragmentation of rock after compaction.The above equation forms a closed system of equations witha, b, and sm as unknowns. For a given seam, the equivalent direct roof thickness h′ and the equivalent mining height m′ can be considered known, so the length of the free stacking zone Lf, the coefficient of fragmentation of the rock in the free stacking state x, and the coefficient of residual fragmentation of the rock in the compacted state KR can be measured to find a, b, and sm, which in turn allows the sinking displacement of the roof at any x to be determined s(x) at any x.The coefficient of rock fragmentation and swelling in the pressure-bearing accumulation zone can be derived from the following equation based on the top slab subsidence displacement:(13)Kpx=m′−sx+∑h′∑h′,where Kp(x) is the coefficient of fragmentation of the rock fall at x in the pressure zone.Whenx = Lf, we have(14)KpLf=Kf.Whenx ⟶ +∞, there is(15)Kp+∞=KR.This means that the coefficient of fragmentation is continuous from the free zone to the compressive zone. Knowing the coefficient of fragmentation, the void ratio of the fallen rock can be calculated using the following formula:(16)n=1−1K,where n is the void ratio of the rock fall and K is the coefficient of fragmentation and swelling of the rock fall.The derived continuous model of the porosity distribution in the extraction zone is not only continuously derivable within each partition but also continuous at the partition boundary. This property helps to ensure the stability of the solution for seepage in the extraction zone.Other literature sets the porosity of the gob as a constant, or by analyzing the lithology of the roof slab of the comprehensive gob and the crushing and swelling characteristics of the fallen gangue, the overburden of the gob is divided into three zones, namely, the natural accumulation zone, the mining influence zone, and the compaction stability zone, and the porosity is set as a constant in each zone, and then the physical model and boundary conditions are solved according to the set permeability.Gao and Wang [30] used FLUENT, based on the porous media seepage theory, to set the permeability of the gob to be uniformly distributed. The results showed that the flow field of the air leakage in the gob can only be simulated by using FLUENT. There results show that the wind flow pattern in the gob is more realistic only if the permeability distribution that more faithfully reflects the rock emergence and compaction pattern in the gob is used. However, these studies are basically flat studies of the porosity of the extraction zone, in the vertical direction are assumed to have no change in the porosity of the extraction zone, and do not comprehensively consider the influence of pressure, time, and media particles on it; to address the problem, Zhang et al. [31] conducted a spatial three-dimensional analysis of the porosity of the extraction zone, the porosity in the extraction zone as a second-order tensor was studied, and the calculation method of the porosity in each direction was derived. The porosity in the vertical direction of the gob was also analyzed, and the relationship between porosity and rock grain size was obtained through simulation tests, and the whole relationship curve was drawn. At the same time, the correctness of this relationship was proved by using field observations and theoretical analysis. Deng et al. [32] determined that the distribution of porosity in the gob in the strike direction varies nearly negatively exponentially by studying the distribution of mineral pressure and roof lithology in the gob.


### 2.3. Study on the Distribution Characteristics of Permeability in the Gob

In order to further study the specific situation of permeability in each location of the extraction zone after the end of surface recovery, Tu et al. [33] approximated the contour plot of permeability distribution in the extraction zone by using the difference algorithm that comes with Suffer software in order to combine the magnitude of permeability in each drilling stabilization stage. As can be seen from the graph, except for the part just after recovery, the permeability rate is obviously larger than the permeability rate on the cut hole side of the fully compacted workings because the recovery zone is not completely compacted. Within 1000 m of the workface strike length, the permeability is lowest in the middle of the extraction zone and gradually increases toward the edge, with a transverse “O” distribution from the inside to the outside, which is symmetrical from top to bottom and from left to right.Wang et al. [34] considered the characteristics of different porosity in each zone and belt of the extraction zone, combined with the Kozeny–Carman equation, and obtained a three-dimensional inhomogeneous permeability model for the extraction zone, which assumes that the porosity variation obeys the sigmoid function from the working face, both ends to the compaction zone and within the fracture zone. The paper firstly studied the porosity models n(x) and n(y) in the strike direction of the working face, as well as in the tendency direction in the one-dimensional case, and then extended to the two-dimensional model on this basis by replacing the one-dimensional model n(x) as the porosity of the compaction zone on the tendency and the porosity of the fractured rock in n(y) after full compaction with n(x) to obtain the two-dimensional porosity distribution model n(x, y), for the porosity three-dimensional model, the variable nz is introduced to represent the difference between the maximum and minimum porosity values parallel to the xoy coordinate plane, and n(x, y, z) = n(x, y) within the collapse zone and the bottom slab fracture zone, and within the fracture zone, the maximum porosity at the working face in the two-dimensional model n(x, y) is replaced by nmin+nmax−nmin/1+ea3z−zc−b3, indicating that within the fracture zone, the maximum porosity above the working face gradually decreases with increasing height, thus obtaining the three-dimensional inhomogeneous distribution model of porosity in the gob as(17)nx,y,z=nmin+nz1+ea1x−b1+nzea1x−b11+ea2y−b21+ea1x−b1,y<ly2,whennz=nmax−nmin,zd≪z<zc,nmax−nmin1+ea3z−zc−b3,zc≪z<zf.In the formula,a3 related to the height of the vertical three bands, a3=2b3/Zf−Zc, b3 takes the value of 5.Substituting the above equation into the Kozeny–Carman equation, the expression for the inhomogeneous permeability of the extraction zone is(18)k=d02n31801−n2,where k is the permeability, m3; d0 is the effective particle diameter, taken as 0.014 m; and n is the porosity.Ren and Edwards [35] established a three-dimensional geotechnical model of the mining void area of the longwall working face based on finite difference or finite element software, determined the three-dimensional stress or strain distribution in the mining void area, and derived the three-dimensional spatial continuous distribution of permeability in the mining void area through the relationship between stress or strain and permeability. Wang [36] constructed a three-dimensional spatial dynamic distribution model of porosity and permeability in the gob based on Usher’s mathematical function.


### 2.4. Relationship between Permeability and Porosity

Both porosity and permeability are basic parameters of porous media. A large number of experiments have proved that there is a certain relationship between porosity and permeability for a particular porous medium.Kong [37] introduced the permeability of a porous medium that is a single bead filling case made of a single spherical particle filling. According to the Carman–Kozeny empirical formula, the following equation is available:(19)KP=Cφ3τΩ2,where C is the Kozeny constant, a dimensionless constant related only to the geometry of the capillary interface; τ is the tortuosity; φ is the porosity; and Ω is the specific surface.By the above, it can be seen that the permeability is inversely proportional to the tortuosity and specific surface and proportional to the third power of porosity.For narrowly screened granular media, many scholars regard Darcy’s law as a laminar flow law for fluids within porous media, while macroscopic hydrodynamic theories and methods are used to study the relationship between porosity, permeability, and interparticle relationships in granular media stacked beds. The relational equations of porosity and permeability are usually written as a function of the squared geometric scale of the particles. Taking spherical granular media as an example, they have the same form of the equation as follows:(20)KP=d2φ3C1−φ2,where d is the particle diameter.For the coefficientC, Ergun gave C = 150, Bear gave C = 180, and Kozeny–Carmen gave C = 172.8.Lei et al. [38] conducted tests for five different narrow sieved sands with particle size range from 100 to 450 μm and particle size ratio of 1.25 : 1 to analyze the law of permeability and porosity of granular media with particle size, and the results showed that under the same porosity condition, the smaller the particle size, the more curved and narrow the flow channel is, and the smaller the permeability is. The permeability of this sand can be expressed by the following equation:(21)Kμm=1.718×104dmm1.465φ4.69.For the porous media in the extraction area, the porosity and permeability of the extracted pores in the extraction area are much larger than the porosity and permeability of the original pores. Therefore, the main consideration for the porous media in the gob is the mining pore. Influenced by the mining pressure, the porosity and permeability of the porous media in the extraction area are functions of the location of the extraction area. For specific conditions, the relationship between permeability and porosity is not fixed, but a general trend is that permeability increases with increasing porosity.


## 2.1. Theory of Stope Fracture Formation

From the process of mineral pressure during mining when the rocks on the roof are risen and fallen, to the compaction of the back of the gob, the medium inside the gob is a porous medium with the nonuniform medium. The pore size distribution is related to the working face height, the size and arrangement of the rock, the lithology of the coal bed and adjacent layers, the original stress and the mining stress, etc. There are two different characteristics of voids in the gob, i.e., mined voids and original voids, which differ significantly in their characteristics, and generally speaking, the size and permeability of mined voids are much larger than those of original voids [3]. In general, the size and permeability of the active pore space are much larger than those of the primary pore space. The active pore space is interconnected throughout the entire area of the extraction space and is not only spatially unevenly distributed but is also influenced by the extraction pressure during the extraction process.


### 2.1.1. The Mine Pressure Hypothesis

The formation of mining overburden fractures is associated with various mine pressures and deformation damage to the surrounding rock within the quarry. From the beginning, mine pressure phenomena have been observed and different explanations have been proposed based on the complex nature of the rock mass. These explanations, which reveal the intrinsic connection between mine pressure phenomena, are called mine pressure hypotheses. The pressure arch hypothesis was proposed by the Germans Haacke and Gulitzel in 1928. In the same vein as the pressure arch hypothesis, the cantilever beam hypothesis was proposed by the German Sledke and later supported by the British Ferriday and the former Soviet Germain, among others [4]. By the 1950s, the articulated block hypothesis and the prefracture hypothesis emerged. The preformed fissure hypothesis was first proposed by a Belgian geologist. In the late 1970s and early 1980s, scholars represented by academician Qian based on the articulated block hypothesis and the rock performed fissure hypothesis, and through on-site observation of the internal movement of rock formations and extensive production practice, the masonry beam theory of rock structure was proposed and developed into a hypothesis. In recent years, as the scientific study of rock formations continued to deepen and to understand the broader issues in rock formation activities, Academician Qian et al. [5] further proposed the key layer theory of rock control based on the masonry beam.


### 2.1.2. The “O” Circle Theory and the Critical Layer Theory [6–10]

Multiple layers of rock of varying thickness and strength exist above the direct top. It has been proven that one to several of these thick hard rock layers play a major controlling role in the overlying rock activity of the quarry. The rock layers that control the activity of the overlying rock layers in the quarry, either locally or up to the surface, are referred to as key layers. The former is referred to as the subcritical layer and the latter as the main critical layer. The key layer is generally a relatively thick and hard rock layer with the following characteristics: geometric features, with thicker single layers relative to other similar rock formations.Lithological characteristics, relatively hard compared to other rock formations, i.e., higher modulus of elasticity and higher strength; deformation characteristics, where the subsidence deformation of the key layer is synchronized with the amount of subsidence of all or part of the overlying rock layers; fracture features, where the fracture of a key layer will result in the simultaneous fracture of all or 1ocal overlying rock layers, causing movement of the rock over a 1arger area; the bearing characteristics, in the form of a slab structure before the key layer breaks as the bearing body for all or part of the rock, become a masonry beam structure after the break if the S-R stability conditions for rock block structures are met, succeeding as the bearing body.The practice has shown that seam activity leads to the creation of a large number of mining fissures, and for coal seams with low permeability, the formation of mining fissures can lead to an increase in permeability of tens to hundreds of times, which creates conditions for the extraction and transport of gas from the coal seam. Two types of fissures are formed in the overlying rock after coal mining: one is the off-bed fissures, which are along-bed fissures that appear between the bed and the seam as the seam sinks, and which cause the coal seam to expand and deform and unload the gas and cause the unloaded gas to gush out along the off-bed fissures; the other is the vertical fracture fissures, which are through-bed fissures that form as the seam sinks and breaks, and these fissures enable the gas to rise and float continuously, forming a rising float.Along the direction of workface advancement, the dynamic distribution of deviations under the key layer shows a two-stage development pattern. Before the initial breakage of the key layer, the amount of deviations increases continuously as the workface advances, and the maximum amount of deviations is located in the middle of the extraction zone. After the initial breakage of the key layer, the key layer tends to be compacted in the middle of the extraction zone, while a delamination zone is still maintained on each side of the extraction zone. The interlayer zones in the upper and lower chutes are connected through the cuttings and the interlayer zone on the working face side, forming an interconnected ring of interlocking interlayer fracture development around the mining void. Generally speaking, the geometry of the open-cut eye, the upper and lower chute, and the working face side are rectangular; that is, the mining void area forms the shape of the moving rock layer as shown in Figure1, which is called the “O” shaped circle of the mining fissure, the porosity of this area is large, easy for the wind flow, and is the main wind leakage channel at the early stage of the mining of the back mining face. With the back mining workface gradually far away from the position of the cutting eye, cutting eye position of the loose triangle by the back mining workface at both ends of the role of wind pressure gradually reduced, when the workface advance to a certain distance, if there are no other sources of wind leakage in the gob and the role of wind leakage sink, the location is no longer the main wind leakage channel in the gob.Figure 1
Diagram of the “O” ring. (a) Tendency profile. (b) Plan view.
(a)
(b)


### 2.1.3. The Zoning of Overburden Movement Failure in Gob

When the coal is fully extracted, the overlying rock layer loses its supporting role and gradually breaks up, collapses, delaminates, fractures, and deforms, showing a process of sinking of the overlying rock layer from static to dynamic to static.The theory of “three horizontal zones and three vertical belts” was proposed by our academician Qian [11]. The overburden fractures formed by mine pressure and complex stresses are divided in different directions, forming “three zones” in the vertical direction, i.e., the caving band, the fracture band, and the bending subsidence belt, which are similar to the “three zones” in the longwall working face as suggested by Karmis et al. [12]. The basic structure and zoning of the overburden in the gob is shown in Figure 2.Figure 2
Zoning of the overlying rock layers of the quarry and delineation of the mining void. (A) Zone of influence of coal wall support. (B) Zone of rock departure. (C) Zone of recompaction. (a) Zone of natural accumulation. (b) Zone of load influence. (c) Zone of compaction and stability. (I) Caving band. (II) Fracture band. (III) Bending subsidence belt.(i) Caving Band. The overburdened rock in the extraction area has been affected by complex mining, and the rock body at the direct top is completely fractured, resulting in a random rock block shape size and arrangement. The medium in the free-fall zone is loose, isotropic, and porous, and it is the main medium in the flow field of the extraction zone. It is generally considered that the coefficient of fragmentation and expansion in this area is around 1.5. As the workings continue to advance, the first blocks to fall are compacted by the subsequent collapse and the coefficient of fragmentation is reduced by approximately 31%.(ii) Fracture Band. When the underground workings are mined under certain conditions in time and space, the overlying rock layer emerges in a regular arrangement after fracturing, which is the main characteristic of the fissure zone. The fissures are mainly divided into two types according to the different forces: one is subject to bending and tensile stress, which is approximated as a fissure in the direction of the vertical rock layer; the other is a fissure in which the rock layer and the rock layer are subject to staggered displacement by mineral pressure, which is approximated as a fissure parallel to the rock layer. As mining continues, the gob gradually increases and the fissure zone is extended in the upper part of the overburden. After a certain depth has been reached in the mining space, the extension stops at the maximum height.(iii) Bending Subsidence Belt. The extent of the area is defined by the top of the rift zone as the lower boundary and the ground surface as its upper boundary. The pressure on the rock within the curved subsidence zone is less than that described by the riser zone, and the rock movement appears not as fragmentation or randomness, but as displacement of the entire area. Perpendicular to the surface, there is a small amount of subsidence displacement. However, in harder rock formations, dislocation fractures are very likely to occur due to misalignment movements between rock layers.Gao [13] established a “four-zone” model of rock movement and divides the overlying rock into four zones according to the mechanical structure characteristics after its destruction, namely, fracture zone, separation zone, bending zone, and loose alluvial zone. Singh and Kendorski [14] found that the area of the collapse zone and the fracture zone depends on the thickness of the coal seam, the strength of the roof, and the characteristics of the rock combination.


## 2.1.1. The Mine Pressure Hypothesis

The formation of mining overburden fractures is associated with various mine pressures and deformation damage to the surrounding rock within the quarry. From the beginning, mine pressure phenomena have been observed and different explanations have been proposed based on the complex nature of the rock mass. These explanations, which reveal the intrinsic connection between mine pressure phenomena, are called mine pressure hypotheses. The pressure arch hypothesis was proposed by the Germans Haacke and Gulitzel in 1928. In the same vein as the pressure arch hypothesis, the cantilever beam hypothesis was proposed by the German Sledke and later supported by the British Ferriday and the former Soviet Germain, among others [4]. By the 1950s, the articulated block hypothesis and the prefracture hypothesis emerged. The preformed fissure hypothesis was first proposed by a Belgian geologist. In the late 1970s and early 1980s, scholars represented by academician Qian based on the articulated block hypothesis and the rock performed fissure hypothesis, and through on-site observation of the internal movement of rock formations and extensive production practice, the masonry beam theory of rock structure was proposed and developed into a hypothesis. In recent years, as the scientific study of rock formations continued to deepen and to understand the broader issues in rock formation activities, Academician Qian et al. [5] further proposed the key layer theory of rock control based on the masonry beam.


## 2.1.2. The “O” Circle Theory and the Critical Layer Theory [6–10]

Multiple layers of rock of varying thickness and strength exist above the direct top. It has been proven that one to several of these thick hard rock layers play a major controlling role in the overlying rock activity of the quarry. The rock layers that control the activity of the overlying rock layers in the quarry, either locally or up to the surface, are referred to as key layers. The former is referred to as the subcritical layer and the latter as the main critical layer. The key layer is generally a relatively thick and hard rock layer with the following characteristics: geometric features, with thicker single layers relative to other similar rock formations.Lithological characteristics, relatively hard compared to other rock formations, i.e., higher modulus of elasticity and higher strength; deformation characteristics, where the subsidence deformation of the key layer is synchronized with the amount of subsidence of all or part of the overlying rock layers; fracture features, where the fracture of a key layer will result in the simultaneous fracture of all or 1ocal overlying rock layers, causing movement of the rock over a 1arger area; the bearing characteristics, in the form of a slab structure before the key layer breaks as the bearing body for all or part of the rock, become a masonry beam structure after the break if the S-R stability conditions for rock block structures are met, succeeding as the bearing body.The practice has shown that seam activity leads to the creation of a large number of mining fissures, and for coal seams with low permeability, the formation of mining fissures can lead to an increase in permeability of tens to hundreds of times, which creates conditions for the extraction and transport of gas from the coal seam. Two types of fissures are formed in the overlying rock after coal mining: one is the off-bed fissures, which are along-bed fissures that appear between the bed and the seam as the seam sinks, and which cause the coal seam to expand and deform and unload the gas and cause the unloaded gas to gush out along the off-bed fissures; the other is the vertical fracture fissures, which are through-bed fissures that form as the seam sinks and breaks, and these fissures enable the gas to rise and float continuously, forming a rising float.Along the direction of workface advancement, the dynamic distribution of deviations under the key layer shows a two-stage development pattern. Before the initial breakage of the key layer, the amount of deviations increases continuously as the workface advances, and the maximum amount of deviations is located in the middle of the extraction zone. After the initial breakage of the key layer, the key layer tends to be compacted in the middle of the extraction zone, while a delamination zone is still maintained on each side of the extraction zone. The interlayer zones in the upper and lower chutes are connected through the cuttings and the interlayer zone on the working face side, forming an interconnected ring of interlocking interlayer fracture development around the mining void. Generally speaking, the geometry of the open-cut eye, the upper and lower chute, and the working face side are rectangular; that is, the mining void area forms the shape of the moving rock layer as shown in Figure1, which is called the “O” shaped circle of the mining fissure, the porosity of this area is large, easy for the wind flow, and is the main wind leakage channel at the early stage of the mining of the back mining face. With the back mining workface gradually far away from the position of the cutting eye, cutting eye position of the loose triangle by the back mining workface at both ends of the role of wind pressure gradually reduced, when the workface advance to a certain distance, if there are no other sources of wind leakage in the gob and the role of wind leakage sink, the location is no longer the main wind leakage channel in the gob.Figure 1
Diagram of the “O” ring. (a) Tendency profile. (b) Plan view.
(a)
(b)


## 2.1.3. The Zoning of Overburden Movement Failure in Gob

When the coal is fully extracted, the overlying rock layer loses its supporting role and gradually breaks up, collapses, delaminates, fractures, and deforms, showing a process of sinking of the overlying rock layer from static to dynamic to static.The theory of “three horizontal zones and three vertical belts” was proposed by our academician Qian [11]. The overburden fractures formed by mine pressure and complex stresses are divided in different directions, forming “three zones” in the vertical direction, i.e., the caving band, the fracture band, and the bending subsidence belt, which are similar to the “three zones” in the longwall working face as suggested by Karmis et al. [12]. The basic structure and zoning of the overburden in the gob is shown in Figure 2.Figure 2
Zoning of the overlying rock layers of the quarry and delineation of the mining void. (A) Zone of influence of coal wall support. (B) Zone of rock departure. (C) Zone of recompaction. (a) Zone of natural accumulation. (b) Zone of load influence. (c) Zone of compaction and stability. (I) Caving band. (II) Fracture band. (III) Bending subsidence belt.(i) Caving Band. The overburdened rock in the extraction area has been affected by complex mining, and the rock body at the direct top is completely fractured, resulting in a random rock block shape size and arrangement. The medium in the free-fall zone is loose, isotropic, and porous, and it is the main medium in the flow field of the extraction zone. It is generally considered that the coefficient of fragmentation and expansion in this area is around 1.5. As the workings continue to advance, the first blocks to fall are compacted by the subsequent collapse and the coefficient of fragmentation is reduced by approximately 31%.(ii) Fracture Band. When the underground workings are mined under certain conditions in time and space, the overlying rock layer emerges in a regular arrangement after fracturing, which is the main characteristic of the fissure zone. The fissures are mainly divided into two types according to the different forces: one is subject to bending and tensile stress, which is approximated as a fissure in the direction of the vertical rock layer; the other is a fissure in which the rock layer and the rock layer are subject to staggered displacement by mineral pressure, which is approximated as a fissure parallel to the rock layer. As mining continues, the gob gradually increases and the fissure zone is extended in the upper part of the overburden. After a certain depth has been reached in the mining space, the extension stops at the maximum height.(iii) Bending Subsidence Belt. The extent of the area is defined by the top of the rift zone as the lower boundary and the ground surface as its upper boundary. The pressure on the rock within the curved subsidence zone is less than that described by the riser zone, and the rock movement appears not as fragmentation or randomness, but as displacement of the entire area. Perpendicular to the surface, there is a small amount of subsidence displacement. However, in harder rock formations, dislocation fractures are very likely to occur due to misalignment movements between rock layers.Gao [13] established a “four-zone” model of rock movement and divides the overlying rock into four zones according to the mechanical structure characteristics after its destruction, namely, fracture zone, separation zone, bending zone, and loose alluvial zone. Singh and Kendorski [14] found that the area of the collapse zone and the fracture zone depends on the thickness of the coal seam, the strength of the roof, and the characteristics of the rock combination.


## 2.2. Study of Porosity Distribution Characteristics in the Gob

In the process of coal mining, due to the influence of mining, the overlying rock layer is damaged, which will not only produce new fissures but also lead to the enlargement of the primary fissures, which will in turn change the porosity of the rock body, and the result leads to changes in the permeability of the rock body. Mastering the characteristics of the porosity distribution in the gob is the key to studying the flow field in the gob. It is impossible to accurately simulate the flow field in the gob by simple porosity analysis alone; therefore, it is important to improve the accuracy of the measurement of the pore distribution in the gob for the study of the flow field in the gob.


### 2.2.1. Fractional Expansion Coefficients and Distribution Functions in the Gob

As time passes, the overlying rock layer gradually sinks under its own weight, the force on the free-fall zone gradually increases, the rock layer is gradually compacted, the porosity of the free-fall zone gradually decreases, and the initial fragmentation and expansion coefficient is larger compared with that after the fall, resulting in a lower height of the flow field in the gob. The height of the zone can be calculated by the following formula [15]:(1)M=H−SK−1,where M is the height of the fall zone (m); S is the actual settlement value of the lower rock beam (old top) in the fissure zone (m); H is the mining height (m); and K is the coefficient of fragmentation and swelling of the fallen rock body.As can be seen from the above equation, if the fragmentation effect of the rock fall is not taken into account, i.e.,K1, the value of M will tend to infinity, indicating that the overburden of the quarry will fall all the way to the surface, resulting in a wide range of estimates of the fall zone. Therefore, the fragmentation and swelling effect of the rock fall has an important influence on the development pattern of deformation and damage of the overlying rock in the quarry, and the fragmentation and swelling characteristics of the rock fall should be fully considered in numerical calculations to make the numerical results more realistic [16].The volume of the rock after crushing should be increased compared with the whole state, and this property is called the crushing and swelling of the rock. The crushing and swelling coefficient is the ratio between the volume of the rock in the loose state after crushing and the volume of the rock in the whole state before crushing [17]. Practice and research have shown that after the collapse of the roof, broken rocks form a haphazard pile of loose bodies in the gob, and depending on the mechanical properties of the rocks, the size and arrangement of the broken rocks vary, and so does their coefficient of fragmentation and expansion, which is generally 1.10–1.40. Fragmentation coefficient and residual expansion coefficient of common coal measures rock are shown in Table 1. In general, hard rocks are broken in large pieces and arranged neatly, and their coefficient of fragmentation and expansion is small, generally 1.10–1.20. The coefficient of fragmentation and swelling of the soft rock layer is large, generally 1.30–1.40, and the coefficient of fragmentation and swelling of the medium-hard rock layer is between the hard and soft rock layers, generally 1.20–1.30 [18].Table 1
Coefficient of fragmentation and residual fragmentation of common coal rocks.
Lithology
Crushing and swelling factorKP
Residual swelling factorKP′
Crushed coal
<1.2
1.05
Muddy shale
1.4
1.1
Sandy mudstone
1.6∼1.8
1.1∼1.15
Sandstone
1.5∼1.8
1.03∼1.1
Medium-hard sandstone
1.3∼1.5
1.03∼1.08Palchik [19] used vertical boreholes to probe the height of the fallout zone of the medium-weathered, strongly weathered, and deeply buried hard overburden rocks in the shallow buried gob, and the theoretical calculation and observation results showed that the fragmentation and swelling coefficients of the strongly weathered rocks ranged from 1.06 to 1.165, and the fallout ratios ranged from 6.07 to 15.6, while the fragmentation and swelling coefficients of the medium-weathered rocks ranged from 1.09 to 1.24, and the fallout ratios ranged from 4.1 to 11.25. The coefficient of fragmentation and swelling of hard overburden rocks can be 1.38 times of that of strongly weathered rocks, while the riser ratio is only 0.26 times of that of strongly weathered rocks, which leads to the conclusion that the higher the strength of rocks, the higher the coefficient of fragmentation and swelling, and the smaller the height of riser zone.In order to study the compaction characteristics of fractured rocks in the coal bed roof, Su et al. [20] sieved sandstone, sandy mudstone, and mudstone in the roof of the Xin’an coal bed in Yima into 5 block sizes (0∼5, 5∼10, 10∼15, 15∼20, and 20∼25 mm) and then took 20% of each of the above 5 block sizes and mixed them as mixed block sizes (0∼25 mm), and in order the compaction tests were carried out on a hybrid electrohydraulic servo rock mechanic test system using a home-made device. The stress-strain relationships for the three types of crushed rock compaction tests were obtained, and the effects of rock strength, block size, and compaction stress on the compaction characteristics of the crushed rock were analyzed.It can be seen from Figure3 that the breaking expansion coefficient of each rock increases as the block diameter increases. The relationship between the crushing coefficient and the block diameter of the sandy mudstone crushed rock in Figure 5 can be described by a logarithmic relationship, which is different from the linear function between the crushing coefficient and the block diameter of the coal rock in the study of Miao et al. [21], and may be related to the range of the block diameter of the crushed rock; analyzing Figure 4, the residual crushing coefficient of the crushed rock after experiencing compaction is not much related to the rock strength and the block diameter. The reduction of the crushing and swelling coefficient is related to the block diameter, and the larger the block diameter, the greater the reduction.Figure 3
Coefficient of fragmentation and swelling of rock as a function of block diameter.Figure 4
Residual coefficient of fragmentation and swelling of crushed rock as a function of block diameter.Figure 5
Crushing and swelling coefficient of sandy mudstone rubble as a function of block size.Using similar material simulation tests, Zhang [22] developed two different models to study the dynamic 1aw of fracture and swelling of the mined rock mass during longwall mining and partial mining of small strips, respectively. The vertical fracture expansion of the modeled rock mass was approximated to replace the volumetric fracture expansion to analyze the dynamic fracture expansion and compressibility of the mined rock mass. The vertical fragmentation coefficient is defined as the ratio of the distance between two adjacent model points in the vertical direction before and after deformation, i.e.,(2)K=hn−n+1′hn−n+1,where hn−n+1′ and hn−n+1 are the vertical distances between two adjacent measurement points n and n+1 before and after mining, respectively, m.The rock fragmentation coefficient of different areas under full and partial mining conditions was obtained in relation to the distance from the coal seam and the observation time; that is, the rock fragmentation coefficient of the observation line located near the boundary of the gob was the largest during full mining of the longwall; the rock fragmentation coefficient of each observation line varied greatly during the mining process and then leveled off soon after the end of mining; the rock fragmentation coefficient of partial mining of the strip was much smaller than that of full mining of the longwall. The coefficient of rock fragmentation and swelling of partial mining is much smaller than that of longwall mining, the coefficient of rock fragmentation and swelling of partial mining varies with the height of waves, and the compressible amount of coal seam area is larger during loading. These understandings are important to explain the residual settlement of the surface in the gob.The porosity of a fractured rock mass is a function of its coefficient of fracture and expansion, and the study of the coefficient of fracture and expansion is the basis for calculating porosity. The porosity of a fractured rock mass can be expressed by the ratio of the pore volume to the total volume of the rock mass in the fractured state, and according to the definition of the porosity of a fractured rock mass and the coefficient of fracture and swelling, the following relationship exists between the two:(3)n=1−1Kρ,where n is the porosity of the crushed rock and Kρ is the crushing and swelling coefficient of the crushed rock.Considering the original porosity of the coal seam roof rockS0, then the total porosity of the rock in the free stacking state is(4)nm0=1−1Kρ+S0,where nm0 is the total porosity of the emergent rock in the free-stacked state and S0 is the original porosity of the top slab.When the workings are finished, as time passes, the rock in the fall zone tends to compact under its own weight and the overlying load, the coefficient of fragmentation becomes smaller, and the final remaining coefficient of fragmentation is called the residual coefficient of fragmentation.The coefficient of fragmentation and swelling is spatially distributed and follows a certain pattern from the start of the fall to the final compaction of the mined area. Li [23], based on the negative exponential decay law, ignored the influence of the shape of the extraction zone, specifically with respect to the distance of a boundary L, and considered that along the direction perpendicular to L, there is(5)Kp,l=Kρ′+KP0−Kρ′e−ald,where KP0 is the initial crushing and swelling factor and Kρ′ is the central compaction fragmentation and swelling coefficient; al is the decay rate; according to the law of mineral pressure, point d is the distance from L.The geometry of the extraction zone does not affect the distribution of the crushing and swelling coefficient, and the distribution pattern of the rock fall in the extraction zone is complex, so the extraction zone is a nonuniform medium. For the whole three-dimensional space, the rock fall fragmentation coefficient is a function of the location of the extraction zone, i.e., the distribution of the fall fragmentation coefficient is(6)Kpx,y=maxKp,l.In Li’s study [24], based on equations (5) and (6), the distance of the points within the gob from the working face and the boundary is reflected in the spatial distribution function of the fragmentation and swelling coefficient of the gob as(7)Kpx,y=Kp,min+Kp,max—Kp,mine−a1d11−e−δa0d0δ<1,where δ is the adjustment function that controls the distribution shape of the model. Kp,max is the initial crushing and swelling coefficient, i.e., the maximum value. Kp,min is the compaction coefficient of fragmentation and swelling, i.e., the minimum value. a0.a1 is the decay rate, and d0.d1 is the distance between any point and the boundary and the working surface, m.


### 2.2.2. Porosity Distribution Model

The voids within the fractured rock scatter of a fall zone are usually irregular and disordered [25, 26]. These voids cannot be described and studied using conventional Euclidean geometry. Numerous studies have shown that the shapes of the various sizes of rock masses produced by rock fragmentation have a fractal structure: i.e., individuals of different sizes have statistically significant self-similar characteristics. Further studies have shown that the voids between these broken rock masses and the piles of broken rock masses are also fractal in structure [27].Xia and Huang [28] used the Menger sponge fractal model to study the void distribution characteristics of the rock scatter in the fall zone. The void distribution model of the adventitious scatter within the adventitious zone based on the Menger sponge fractal model is shown in Figure 6, with the similarity ratio of(8)t=yx=1k,where y is the volume of the first level transformed small cubic.Figure 6
Menger sponge model.y can be considered as the generating element of a Koch curve of dimension:(9)D=lnNln1/t.Bringing equation (9) into equation (8), we can obtain the void volume fraction dimension of the rock scatter in the bubble fall zone Dk:(10)Dk=lnk3−nlnk.Based on this model, the total fractal void ratioP of the rock bulk in the fall zone was obtained as(11)P=V−M/ρV=1−ρrmax3−Dρ0rmax3−D−rmin3−D,where ρ0 is the density of the rock fall, kg/m3; V is the total fractal volume of rock bulk within the fall zone, m3; M is the mass of the fallen rock bulk, kg; ρ is the density of the crushed rock bulk, kg/m3; and R is the diameter of individual rock mass, m.Wang et al. [29] divided the extraction zone into free accumulation zone and pressure accumulation zone along the working face advance direction, as shown in Figure 7.Figure 7
Schematic diagram of the state of rock accumulation in the extraction area.Derived from the theory of surrounding rock movement in the stope,(12)hx=m′−Sm1−e−aLf/Lb+1−Kf∑h′,sm=m′+1−KR∑h′,3b−1+5b−1b−12ab=1,where x is the strike distance of a point in the gob from the working face, m; h(x) is the gap between the rock fall and the roof at x, m. Lf is the length of the free accumulation zone, m; ∑h′ is the equivalent direct roof thickness, m; m′ is the equivalent mining height, m. sm is the displacement of the rock seam after it has moved and stabilized, m; L is the distance of the basic stabilization point of the rock seam from the working face, usually about 50∼60m; a and b are coefficients that vary with the distance of the rock seam from the coal seam and the lithological characteristics. Kf is the coefficient of fragmentation of rock in the free stacking zone, and KR is the coefficient of residual fragmentation of rock after compaction.The above equation forms a closed system of equations witha, b, and sm as unknowns. For a given seam, the equivalent direct roof thickness h′ and the equivalent mining height m′ can be considered known, so the length of the free stacking zone Lf, the coefficient of fragmentation of the rock in the free stacking state x, and the coefficient of residual fragmentation of the rock in the compacted state KR can be measured to find a, b, and sm, which in turn allows the sinking displacement of the roof at any x to be determined s(x) at any x.The coefficient of rock fragmentation and swelling in the pressure-bearing accumulation zone can be derived from the following equation based on the top slab subsidence displacement:(13)Kpx=m′−sx+∑h′∑h′,where Kp(x) is the coefficient of fragmentation of the rock fall at x in the pressure zone.Whenx = Lf, we have(14)KpLf=Kf.Whenx ⟶ +∞, there is(15)Kp+∞=KR.This means that the coefficient of fragmentation is continuous from the free zone to the compressive zone. Knowing the coefficient of fragmentation, the void ratio of the fallen rock can be calculated using the following formula:(16)n=1−1K,where n is the void ratio of the rock fall and K is the coefficient of fragmentation and swelling of the rock fall.The derived continuous model of the porosity distribution in the extraction zone is not only continuously derivable within each partition but also continuous at the partition boundary. This property helps to ensure the stability of the solution for seepage in the extraction zone.Other literature sets the porosity of the gob as a constant, or by analyzing the lithology of the roof slab of the comprehensive gob and the crushing and swelling characteristics of the fallen gangue, the overburden of the gob is divided into three zones, namely, the natural accumulation zone, the mining influence zone, and the compaction stability zone, and the porosity is set as a constant in each zone, and then the physical model and boundary conditions are solved according to the set permeability.Gao and Wang [30] used FLUENT, based on the porous media seepage theory, to set the permeability of the gob to be uniformly distributed. The results showed that the flow field of the air leakage in the gob can only be simulated by using FLUENT. There results show that the wind flow pattern in the gob is more realistic only if the permeability distribution that more faithfully reflects the rock emergence and compaction pattern in the gob is used. However, these studies are basically flat studies of the porosity of the extraction zone, in the vertical direction are assumed to have no change in the porosity of the extraction zone, and do not comprehensively consider the influence of pressure, time, and media particles on it; to address the problem, Zhang et al. [31] conducted a spatial three-dimensional analysis of the porosity of the extraction zone, the porosity in the extraction zone as a second-order tensor was studied, and the calculation method of the porosity in each direction was derived. The porosity in the vertical direction of the gob was also analyzed, and the relationship between porosity and rock grain size was obtained through simulation tests, and the whole relationship curve was drawn. At the same time, the correctness of this relationship was proved by using field observations and theoretical analysis. Deng et al. [32] determined that the distribution of porosity in the gob in the strike direction varies nearly negatively exponentially by studying the distribution of mineral pressure and roof lithology in the gob.


## 2.2.1. Fractional Expansion Coefficients and Distribution Functions in the Gob

As time passes, the overlying rock layer gradually sinks under its own weight, the force on the free-fall zone gradually increases, the rock layer is gradually compacted, the porosity of the free-fall zone gradually decreases, and the initial fragmentation and expansion coefficient is larger compared with that after the fall, resulting in a lower height of the flow field in the gob. The height of the zone can be calculated by the following formula [15]:(1)M=H−SK−1,where M is the height of the fall zone (m); S is the actual settlement value of the lower rock beam (old top) in the fissure zone (m); H is the mining height (m); and K is the coefficient of fragmentation and swelling of the fallen rock body.As can be seen from the above equation, if the fragmentation effect of the rock fall is not taken into account, i.e.,K1, the value of M will tend to infinity, indicating that the overburden of the quarry will fall all the way to the surface, resulting in a wide range of estimates of the fall zone. Therefore, the fragmentation and swelling effect of the rock fall has an important influence on the development pattern of deformation and damage of the overlying rock in the quarry, and the fragmentation and swelling characteristics of the rock fall should be fully considered in numerical calculations to make the numerical results more realistic [16].The volume of the rock after crushing should be increased compared with the whole state, and this property is called the crushing and swelling of the rock. The crushing and swelling coefficient is the ratio between the volume of the rock in the loose state after crushing and the volume of the rock in the whole state before crushing [17]. Practice and research have shown that after the collapse of the roof, broken rocks form a haphazard pile of loose bodies in the gob, and depending on the mechanical properties of the rocks, the size and arrangement of the broken rocks vary, and so does their coefficient of fragmentation and expansion, which is generally 1.10–1.40. Fragmentation coefficient and residual expansion coefficient of common coal measures rock are shown in Table 1. In general, hard rocks are broken in large pieces and arranged neatly, and their coefficient of fragmentation and expansion is small, generally 1.10–1.20. The coefficient of fragmentation and swelling of the soft rock layer is large, generally 1.30–1.40, and the coefficient of fragmentation and swelling of the medium-hard rock layer is between the hard and soft rock layers, generally 1.20–1.30 [18].Table 1
Coefficient of fragmentation and residual fragmentation of common coal rocks.
Lithology
Crushing and swelling factorKP
Residual swelling factorKP′
Crushed coal
<1.2
1.05
Muddy shale
1.4
1.1
Sandy mudstone
1.6∼1.8
1.1∼1.15
Sandstone
1.5∼1.8
1.03∼1.1
Medium-hard sandstone
1.3∼1.5
1.03∼1.08Palchik [19] used vertical boreholes to probe the height of the fallout zone of the medium-weathered, strongly weathered, and deeply buried hard overburden rocks in the shallow buried gob, and the theoretical calculation and observation results showed that the fragmentation and swelling coefficients of the strongly weathered rocks ranged from 1.06 to 1.165, and the fallout ratios ranged from 6.07 to 15.6, while the fragmentation and swelling coefficients of the medium-weathered rocks ranged from 1.09 to 1.24, and the fallout ratios ranged from 4.1 to 11.25. The coefficient of fragmentation and swelling of hard overburden rocks can be 1.38 times of that of strongly weathered rocks, while the riser ratio is only 0.26 times of that of strongly weathered rocks, which leads to the conclusion that the higher the strength of rocks, the higher the coefficient of fragmentation and swelling, and the smaller the height of riser zone.In order to study the compaction characteristics of fractured rocks in the coal bed roof, Su et al. [20] sieved sandstone, sandy mudstone, and mudstone in the roof of the Xin’an coal bed in Yima into 5 block sizes (0∼5, 5∼10, 10∼15, 15∼20, and 20∼25 mm) and then took 20% of each of the above 5 block sizes and mixed them as mixed block sizes (0∼25 mm), and in order the compaction tests were carried out on a hybrid electrohydraulic servo rock mechanic test system using a home-made device. The stress-strain relationships for the three types of crushed rock compaction tests were obtained, and the effects of rock strength, block size, and compaction stress on the compaction characteristics of the crushed rock were analyzed.It can be seen from Figure3 that the breaking expansion coefficient of each rock increases as the block diameter increases. The relationship between the crushing coefficient and the block diameter of the sandy mudstone crushed rock in Figure 5 can be described by a logarithmic relationship, which is different from the linear function between the crushing coefficient and the block diameter of the coal rock in the study of Miao et al. [21], and may be related to the range of the block diameter of the crushed rock; analyzing Figure 4, the residual crushing coefficient of the crushed rock after experiencing compaction is not much related to the rock strength and the block diameter. The reduction of the crushing and swelling coefficient is related to the block diameter, and the larger the block diameter, the greater the reduction.Figure 3
Coefficient of fragmentation and swelling of rock as a function of block diameter.Figure 4
Residual coefficient of fragmentation and swelling of crushed rock as a function of block diameter.Figure 5
Crushing and swelling coefficient of sandy mudstone rubble as a function of block size.Using similar material simulation tests, Zhang [22] developed two different models to study the dynamic 1aw of fracture and swelling of the mined rock mass during longwall mining and partial mining of small strips, respectively. The vertical fracture expansion of the modeled rock mass was approximated to replace the volumetric fracture expansion to analyze the dynamic fracture expansion and compressibility of the mined rock mass. The vertical fragmentation coefficient is defined as the ratio of the distance between two adjacent model points in the vertical direction before and after deformation, i.e.,(2)K=hn−n+1′hn−n+1,where hn−n+1′ and hn−n+1 are the vertical distances between two adjacent measurement points n and n+1 before and after mining, respectively, m.The rock fragmentation coefficient of different areas under full and partial mining conditions was obtained in relation to the distance from the coal seam and the observation time; that is, the rock fragmentation coefficient of the observation line located near the boundary of the gob was the largest during full mining of the longwall; the rock fragmentation coefficient of each observation line varied greatly during the mining process and then leveled off soon after the end of mining; the rock fragmentation coefficient of partial mining of the strip was much smaller than that of full mining of the longwall. The coefficient of rock fragmentation and swelling of partial mining is much smaller than that of longwall mining, the coefficient of rock fragmentation and swelling of partial mining varies with the height of waves, and the compressible amount of coal seam area is larger during loading. These understandings are important to explain the residual settlement of the surface in the gob.The porosity of a fractured rock mass is a function of its coefficient of fracture and expansion, and the study of the coefficient of fracture and expansion is the basis for calculating porosity. The porosity of a fractured rock mass can be expressed by the ratio of the pore volume to the total volume of the rock mass in the fractured state, and according to the definition of the porosity of a fractured rock mass and the coefficient of fracture and swelling, the following relationship exists between the two:(3)n=1−1Kρ,where n is the porosity of the crushed rock and Kρ is the crushing and swelling coefficient of the crushed rock.Considering the original porosity of the coal seam roof rockS0, then the total porosity of the rock in the free stacking state is(4)nm0=1−1Kρ+S0,where nm0 is the total porosity of the emergent rock in the free-stacked state and S0 is the original porosity of the top slab.When the workings are finished, as time passes, the rock in the fall zone tends to compact under its own weight and the overlying load, the coefficient of fragmentation becomes smaller, and the final remaining coefficient of fragmentation is called the residual coefficient of fragmentation.The coefficient of fragmentation and swelling is spatially distributed and follows a certain pattern from the start of the fall to the final compaction of the mined area. Li [23], based on the negative exponential decay law, ignored the influence of the shape of the extraction zone, specifically with respect to the distance of a boundary L, and considered that along the direction perpendicular to L, there is(5)Kp,l=Kρ′+KP0−Kρ′e−ald,where KP0 is the initial crushing and swelling factor and Kρ′ is the central compaction fragmentation and swelling coefficient; al is the decay rate; according to the law of mineral pressure, point d is the distance from L.The geometry of the extraction zone does not affect the distribution of the crushing and swelling coefficient, and the distribution pattern of the rock fall in the extraction zone is complex, so the extraction zone is a nonuniform medium. For the whole three-dimensional space, the rock fall fragmentation coefficient is a function of the location of the extraction zone, i.e., the distribution of the fall fragmentation coefficient is(6)Kpx,y=maxKp,l.In Li’s study [24], based on equations (5) and (6), the distance of the points within the gob from the working face and the boundary is reflected in the spatial distribution function of the fragmentation and swelling coefficient of the gob as(7)Kpx,y=Kp,min+Kp,max—Kp,mine−a1d11−e−δa0d0δ<1,where δ is the adjustment function that controls the distribution shape of the model. Kp,max is the initial crushing and swelling coefficient, i.e., the maximum value. Kp,min is the compaction coefficient of fragmentation and swelling, i.e., the minimum value. a0.a1 is the decay rate, and d0.d1 is the distance between any point and the boundary and the working surface, m.


## 2.2.2. Porosity Distribution Model

The voids within the fractured rock scatter of a fall zone are usually irregular and disordered [25, 26]. These voids cannot be described and studied using conventional Euclidean geometry. Numerous studies have shown that the shapes of the various sizes of rock masses produced by rock fragmentation have a fractal structure: i.e., individuals of different sizes have statistically significant self-similar characteristics. Further studies have shown that the voids between these broken rock masses and the piles of broken rock masses are also fractal in structure [27].Xia and Huang [28] used the Menger sponge fractal model to study the void distribution characteristics of the rock scatter in the fall zone. The void distribution model of the adventitious scatter within the adventitious zone based on the Menger sponge fractal model is shown in Figure 6, with the similarity ratio of(8)t=yx=1k,where y is the volume of the first level transformed small cubic.Figure 6
Menger sponge model.y can be considered as the generating element of a Koch curve of dimension:(9)D=lnNln1/t.Bringing equation (9) into equation (8), we can obtain the void volume fraction dimension of the rock scatter in the bubble fall zone Dk:(10)Dk=lnk3−nlnk.Based on this model, the total fractal void ratioP of the rock bulk in the fall zone was obtained as(11)P=V−M/ρV=1−ρrmax3−Dρ0rmax3−D−rmin3−D,where ρ0 is the density of the rock fall, kg/m3; V is the total fractal volume of rock bulk within the fall zone, m3; M is the mass of the fallen rock bulk, kg; ρ is the density of the crushed rock bulk, kg/m3; and R is the diameter of individual rock mass, m.Wang et al. [29] divided the extraction zone into free accumulation zone and pressure accumulation zone along the working face advance direction, as shown in Figure 7.Figure 7
Schematic diagram of the state of rock accumulation in the extraction area.Derived from the theory of surrounding rock movement in the stope,(12)hx=m′−Sm1−e−aLf/Lb+1−Kf∑h′,sm=m′+1−KR∑h′,3b−1+5b−1b−12ab=1,where x is the strike distance of a point in the gob from the working face, m; h(x) is the gap between the rock fall and the roof at x, m. Lf is the length of the free accumulation zone, m; ∑h′ is the equivalent direct roof thickness, m; m′ is the equivalent mining height, m. sm is the displacement of the rock seam after it has moved and stabilized, m; L is the distance of the basic stabilization point of the rock seam from the working face, usually about 50∼60m; a and b are coefficients that vary with the distance of the rock seam from the coal seam and the lithological characteristics. Kf is the coefficient of fragmentation of rock in the free stacking zone, and KR is the coefficient of residual fragmentation of rock after compaction.The above equation forms a closed system of equations witha, b, and sm as unknowns. For a given seam, the equivalent direct roof thickness h′ and the equivalent mining height m′ can be considered known, so the length of the free stacking zone Lf, the coefficient of fragmentation of the rock in the free stacking state x, and the coefficient of residual fragmentation of the rock in the compacted state KR can be measured to find a, b, and sm, which in turn allows the sinking displacement of the roof at any x to be determined s(x) at any x.The coefficient of rock fragmentation and swelling in the pressure-bearing accumulation zone can be derived from the following equation based on the top slab subsidence displacement:(13)Kpx=m′−sx+∑h′∑h′,where Kp(x) is the coefficient of fragmentation of the rock fall at x in the pressure zone.Whenx = Lf, we have(14)KpLf=Kf.Whenx ⟶ +∞, there is(15)Kp+∞=KR.This means that the coefficient of fragmentation is continuous from the free zone to the compressive zone. Knowing the coefficient of fragmentation, the void ratio of the fallen rock can be calculated using the following formula:(16)n=1−1K,where n is the void ratio of the rock fall and K is the coefficient of fragmentation and swelling of the rock fall.The derived continuous model of the porosity distribution in the extraction zone is not only continuously derivable within each partition but also continuous at the partition boundary. This property helps to ensure the stability of the solution for seepage in the extraction zone.Other literature sets the porosity of the gob as a constant, or by analyzing the lithology of the roof slab of the comprehensive gob and the crushing and swelling characteristics of the fallen gangue, the overburden of the gob is divided into three zones, namely, the natural accumulation zone, the mining influence zone, and the compaction stability zone, and the porosity is set as a constant in each zone, and then the physical model and boundary conditions are solved according to the set permeability.Gao and Wang [30] used FLUENT, based on the porous media seepage theory, to set the permeability of the gob to be uniformly distributed. The results showed that the flow field of the air leakage in the gob can only be simulated by using FLUENT. There results show that the wind flow pattern in the gob is more realistic only if the permeability distribution that more faithfully reflects the rock emergence and compaction pattern in the gob is used. However, these studies are basically flat studies of the porosity of the extraction zone, in the vertical direction are assumed to have no change in the porosity of the extraction zone, and do not comprehensively consider the influence of pressure, time, and media particles on it; to address the problem, Zhang et al. [31] conducted a spatial three-dimensional analysis of the porosity of the extraction zone, the porosity in the extraction zone as a second-order tensor was studied, and the calculation method of the porosity in each direction was derived. The porosity in the vertical direction of the gob was also analyzed, and the relationship between porosity and rock grain size was obtained through simulation tests, and the whole relationship curve was drawn. At the same time, the correctness of this relationship was proved by using field observations and theoretical analysis. Deng et al. [32] determined that the distribution of porosity in the gob in the strike direction varies nearly negatively exponentially by studying the distribution of mineral pressure and roof lithology in the gob.


## 2.3. Study on the Distribution Characteristics of Permeability in the Gob

In order to further study the specific situation of permeability in each location of the extraction zone after the end of surface recovery, Tu et al. [33] approximated the contour plot of permeability distribution in the extraction zone by using the difference algorithm that comes with Suffer software in order to combine the magnitude of permeability in each drilling stabilization stage. As can be seen from the graph, except for the part just after recovery, the permeability rate is obviously larger than the permeability rate on the cut hole side of the fully compacted workings because the recovery zone is not completely compacted. Within 1000 m of the workface strike length, the permeability is lowest in the middle of the extraction zone and gradually increases toward the edge, with a transverse “O” distribution from the inside to the outside, which is symmetrical from top to bottom and from left to right.Wang et al. [34] considered the characteristics of different porosity in each zone and belt of the extraction zone, combined with the Kozeny–Carman equation, and obtained a three-dimensional inhomogeneous permeability model for the extraction zone, which assumes that the porosity variation obeys the sigmoid function from the working face, both ends to the compaction zone and within the fracture zone. The paper firstly studied the porosity models n(x) and n(y) in the strike direction of the working face, as well as in the tendency direction in the one-dimensional case, and then extended to the two-dimensional model on this basis by replacing the one-dimensional model n(x) as the porosity of the compaction zone on the tendency and the porosity of the fractured rock in n(y) after full compaction with n(x) to obtain the two-dimensional porosity distribution model n(x, y), for the porosity three-dimensional model, the variable nz is introduced to represent the difference between the maximum and minimum porosity values parallel to the xoy coordinate plane, and n(x, y, z) = n(x, y) within the collapse zone and the bottom slab fracture zone, and within the fracture zone, the maximum porosity at the working face in the two-dimensional model n(x, y) is replaced by nmin+nmax−nmin/1+ea3z−zc−b3, indicating that within the fracture zone, the maximum porosity above the working face gradually decreases with increasing height, thus obtaining the three-dimensional inhomogeneous distribution model of porosity in the gob as(17)nx,y,z=nmin+nz1+ea1x−b1+nzea1x−b11+ea2y−b21+ea1x−b1,y<ly2,whennz=nmax−nmin,zd≪z<zc,nmax−nmin1+ea3z−zc−b3,zc≪z<zf.In the formula,a3 related to the height of the vertical three bands, a3=2b3/Zf−Zc, b3 takes the value of 5.Substituting the above equation into the Kozeny–Carman equation, the expression for the inhomogeneous permeability of the extraction zone is(18)k=d02n31801−n2,where k is the permeability, m3; d0 is the effective particle diameter, taken as 0.014 m; and n is the porosity.Ren and Edwards [35] established a three-dimensional geotechnical model of the mining void area of the longwall working face based on finite difference or finite element software, determined the three-dimensional stress or strain distribution in the mining void area, and derived the three-dimensional spatial continuous distribution of permeability in the mining void area through the relationship between stress or strain and permeability. Wang [36] constructed a three-dimensional spatial dynamic distribution model of porosity and permeability in the gob based on Usher’s mathematical function.


## 2.4. Relationship between Permeability and Porosity

Both porosity and permeability are basic parameters of porous media. A large number of experiments have proved that there is a certain relationship between porosity and permeability for a particular porous medium.Kong [37] introduced the permeability of a porous medium that is a single bead filling case made of a single spherical particle filling. According to the Carman–Kozeny empirical formula, the following equation is available:(19)KP=Cφ3τΩ2,where C is the Kozeny constant, a dimensionless constant related only to the geometry of the capillary interface; τ is the tortuosity; φ is the porosity; and Ω is the specific surface.By the above, it can be seen that the permeability is inversely proportional to the tortuosity and specific surface and proportional to the third power of porosity.For narrowly screened granular media, many scholars regard Darcy’s law as a laminar flow law for fluids within porous media, while macroscopic hydrodynamic theories and methods are used to study the relationship between porosity, permeability, and interparticle relationships in granular media stacked beds. The relational equations of porosity and permeability are usually written as a function of the squared geometric scale of the particles. Taking spherical granular media as an example, they have the same form of the equation as follows:(20)KP=d2φ3C1−φ2,where d is the particle diameter.For the coefficientC, Ergun gave C = 150, Bear gave C = 180, and Kozeny–Carmen gave C = 172.8.Lei et al. [38] conducted tests for five different narrow sieved sands with particle size range from 100 to 450 μm and particle size ratio of 1.25 : 1 to analyze the law of permeability and porosity of granular media with particle size, and the results showed that under the same porosity condition, the smaller the particle size, the more curved and narrow the flow channel is, and the smaller the permeability is. The permeability of this sand can be expressed by the following equation:(21)Kμm=1.718×104dmm1.465φ4.69.For the porous media in the extraction area, the porosity and permeability of the extracted pores in the extraction area are much larger than the porosity and permeability of the original pores. Therefore, the main consideration for the porous media in the gob is the mining pore. Influenced by the mining pressure, the porosity and permeability of the porous media in the extraction area are functions of the location of the extraction area. For specific conditions, the relationship between permeability and porosity is not fixed, but a general trend is that permeability increases with increasing porosity.


## 3. Current Status of Theoretical Research on Flow Field in Gob

### 3.1. Methodology for Studying Flow Field in Gob

#### 3.1.1. Graphical Method and Filtered Flow Field Model

Early studies on the gas flow law in the gob usually analyzed the air leakage in the gob by considering the gob and the retreating working face as two parallel wind paths, and the model was applicable to the study of the U-type ventilation method without leakage sources and sinks after the retreating face, but the theory was not applicable to the situation where the number of leakage sources was too many or the geometry of the gob was complicated, which had obvious limitations. In 1978, Chinese scholars Huang and Wang [39], etc., set the quarry with the long-walled backward coal mining method as a standard rectangular body, set the air-gob as a homogeneous porous medium, considered that the air flow conditions on each plane parallel to the quarry floor are similar, and cited the theory of fluid mechanics about the vortex-free steady motion in the plane, in the case of one source and one sink, one-source and one-sink quarry means that the air is fed from the inlet down-slot. In the case of one-source and one-sink quarry, the air is fed from the inlet down-slot and discharged from the return down-slot, and there is no air leakage. The graphical method is used to describe the gas flow law in the rectangular gob, and the mine ventilation system is adjusted according to the graphical results. However, once the number of air leakage sources is too large or the geometry of the quarry is complicated, the graphical method is difficult to apply and even technically infeasible.In the same period, Polish and Soviet scholars put forward the filter flow field theory, which regarded the medium in the gob as a longitudinal and horizontal pipe network and believed that the flow of gas in the gob was similar to the flow in the pipe network. Shi et al. [40] used a computer to divide the fallout zone in the gob into a number of square grids, and each square grid was replaced by a filter flow branch network. The filter flow velocity and air leakage volume of the fallen strip in the gob were determined, and the computer simulation of the technical measures to reduce the air leakage volume in the gob was carried out according to the theory of filter flow field, and the results showed that the air leakage reduction ranged from 17.08% to 24.26%. Ding et al. [41] gridded the bubble fall zone in the gob, as shown in Figure 8, and established the mathematical model of the filter flow field in the gob, calculated the wind resistance of each branch of the filter flow field in the gob, and programmed a C language program to solve the filter flow field according to the law of wind volume balance and the law of wind pressure balance, as well as the iterative method of Scott-Hensley, to obtain the wind speed of each branch of the filter flow in the gob:(22)h=R′Q+R″Q2,where h is the resistance of the filtered flow field branches, Pa; R′ is the filtration flow field laminar wind resistance, kg/(m4•s); R″ is the filtration flow turbulent wind resistance, kg/m7; and Q is the air volume of filter flow branch, m3/s.Figure 8
Air leakage branches and node numbers in the extraction area.


#### 3.1.2. Similarity Simulation Method

Similarity simulation (including physical similarity simulation and electrical simulation) is based on similarity theory and plays an important role in the early research stage of flow field in the gob. The so-called similarity principle means that when the mathematical models of two physical fields, the same or different, are isomorphic, the distribution laws of the two field quantities are the same, but with different signs such as functions or coefficients, and when the solution of one field has been solved or is easy to solve, the other unknown and relatively difficult field can be solved based on it. The similarity simulation method is relatively simple to use, and researchers can study different physical fields under various conditions by creating similar simulations in the laboratory.Electrical simulation is one of the simulation methods with more advantages. Zhu and Hui [42] introduced the electric simulation method into the study of mine ventilation, comparing the wind field with the electric field, and proved by practice that the introduction of the electric simulation method into coal mines is useful for the study of problems such as ventilation, but there are some problems with the electric simulation, such as the unsatisfactory treatment of the junction conditions in places where the field changes drastically, the poor mastery of the model dissection technique, the poor orthogonality in some parts of the flow network, and the existence of some errors.Xu et al. [43] used the physical similarity simulation method to realize the qualitative study of wind flow state by analyzing the smoke flow state. Yang et al. [44] established a model of coal seam permeability coefficient based on similar material experiments simulating coal mining with fractal dimension of mining fractures.


#### 3.1.3. CFD Numerical Simulation Method

Numerical simulation, also known as computer simulation, is the study of engineering and physical problems by means of numerical computation and image display using an electronic computer as a carrier. CFD is a branch of computational fluid dynamics, where the mathematical description of the overall flow field of interest is obtained by numerically solving a system of equations to reflect the fluid flow, heat and mass transfer laws and related phenomena.With the continuous progress of computer science and technology, a lot of CFD commercial software has been created, such as PHOENICS and FLUENT. PHOENICS is the world’s first set of commercial software for computational fluids and computational heat transfer, which was developed by the main founder of international computational fluids and computational heat transfer, Professor D. B. Spalding, a member of the Royal Academy of Engineering, and more than 40 Ph. Although PHOENICS was introduced into the field of coal mine flow field earlier, the software is completely applied to the calculation of the flow field of rock fall media in the gob, there are still many problems, such as the adaptation and debugging of specific problems, with large differences and limitations, and there are still some problems to overcome in the integration of analysis and calculation with this professional problem. FLUENT is currently international more popular commercial CFD software package, with a high market share. It has rich physical models, advanced numerical methods, and powerful pre- and postprocessing functions and has a wide range of applications. It uses the finite volume method based on a completely unstructured grid and has a gradient algorithm based on grid nodes and grid cells.In recent years, CFD numerical simulation of wind flow motion and gas transport and its distribution pattern in mining sites has been carried out in China. Zhou et al. [45] used PHOENICS software to simulate the distribution law of the wind flow and temperature field at the back mining workings and concluded that the wind flow temperature at the back mining workings decreases with the increase of wind speed in a negative power function, increases linearly with the increase of incoming wind flow temperature, and varies with the length of the workings and the location of the coal miner. The numerical simulation program (G3) of the field flow in the gob prepared by Li [23] has the following characteristics: (1) it can adapt to various complex boundary conditions and treat the gob as a bubbling nonhomogeneous medium; (2) it takes into account the natural coupling effect when gas gushes out and can reflect the relationship between the influence of factors such as working face advancement and ventilation volume; (3) it is operable and can make arbitrary simulation test faces for various situations, which is convenient. Theoretically, it can depict the flow state of air leakage in the gob, dynamically depict the distribution state of gas, oxygen, carbon monoxide concentration and temperature, and their change process, and give the graphical distribution solution of each phase of gas and temperature change in the gob under the mining conditions of the working face; (4) it realizes the integration of the joint solution of “a multiphase problem” in the gob, which can meet the needs of the study of the problem. G3 can assist in analyzing many safety problems such as gas emission from the gob, spontaneous combustion and fire, nitrogen injection and fire extinguishing, and large area pressure equalization and pressure regulation. It is powerful and comprehensive and represents the development direction of domestic application software in the field of numerical simulation of field flow in the gob.


### 3.2. Numerical Model Equation Study

#### 3.2.1. Linear Gas Seepage Theory

The numerical model equations were initially simplified according to Darcy’s law, which is a typical representative of linear gas seepage theory, and it is believed that the gas flow law in coal seams is basically in line with the linear percolation law. In the international aspect, the former Soviet Union scholars used Darcy’s law for the first time to analyze and study the problem of linear gas seepage in coal seams by using the nature of gas adsorption, which opened a new chapter for the study of linear gas seepage theory. In China, experts such as Zhou and Sun [46] agreed with the linear gas seepage theory, and on this basis, they analyzed and studied the flow of gas in coal seams for the first time, and regarded coal seams as continuous and uniformly distributed porous media, which laid the theoretical foundation for the analysis and study of gas flow theory by domestic experts and scholars. The linear seepage theory believes that the flow rate of gas in coal and the pressure difference obey a linear mathematical relationship; that is, the flow rate is proportional to the pressure difference of gas. Since the 1980s, scholars in China have carried out research on the numerical equations of linear flow theory and have modified and improved them on this basis.In 1986, Tan and Yuan [47] combined the theory of seepage mechanics and thermodynamics, treated coal seam gas as real gas, proposed the seepage equation of real gas in the coal seam of the mine, and simplified the equation. In 1989, Yu et al. [48] proposed the idea that the amount of gas involved in seepage in coal seam is the partial amount of gas content in coal body, and under the assumption that the process of gas adsorption and desorption in coal body is completely reversible, the controlling equation of gas seepage in the coal seam was established; in 2010, Yang and Li [49] established a mathematical model of gas flow in the coal body around the borehole based on Darcy’s law and the assumptions of coal seam gas flow theory and applied the law of mass conservation to obtain the equation of gas flow around the borehole and its fixed solution conditions.


#### 3.2.2. Nonlinear Gas Seepage Theory

Darcy’s law has a good fit for laminar flow motion that conforms to the linear drag relationship, but in practice the flow of wind flow in the gob is often very complex, with a large portion of the area being transitional seepage or turbulent flow. In these regions, fluid motion does not conform to Darcy’s law. Therefore, many scholars have proposed nonlinear seepage laws based on this problem and conducted research based on nonlinear seepage theory.In 1984, Higuchi Sumitoshi, a professor at Hokkaido University, Japan, determined the methane permeability of coal samples by varying the pressure difference and found that Darcy’s law was not quite consistent with the coal seam gas flow law. Based on a large number of experimental studies, he proposed the basic law of coal seam gas flow that is more consistent with the power law:(23)VN=−Adpdxm,where VN is the uncaused gas flow rate in the corresponding standard state; A is the uncaused gas permeability coefficient; m is the state constant; and dp/dx is the uncaused gas pressure gradient along the X-axis.Sun [50] considered that the mathematical model of gas flow proposed by Professor Higuchi was not rigorous and concluded that the gas flow in coal seams is actually a mixed unsteady flow of compressible fluid by permeation-diffusion in an anisotropic and nonhomogeneous pore-fissure twofold medium. He established the first partial differential equation for compressible gas flow in coal seams under the conditions of homogeneous and nonhomogeneous coal seams according to the generalized form of power law and made practical calculations to verify the flow of gas flow field in homogeneous coal seams based on the actual measured gas flow parameters in 23051 working face of the Jiaozuo Zhongmacun Mine, and the results showed that it was closer to the actual situation than Darcy’s law. In 1994, Liu [51] concluded that the coal seam gas flow model derived in the literature [41] on the basis of the power law and the application of the invariant had errors in the process, and the numerical results obtained by it were unreliable, and derived the correct flow model when the gas flow followed the power law in the paper.Deng et al. [52] analyzed a large amount of experimental data, on this basis generalized the universal characteristic equation of motion for seepage curves, and proposed a nonlinear equation of motion for gas seepage in low-permeability pore-fracture media. However, this type of model mainly adopts the method of phenomenology, which is only a generalization and distillation of the experimental phenomena and does not have the function of in-depth explanation, which is manifested in this study mainly by scholars fitting mathematical equations according to the physical phenomena of nonlinear seepage, and it is difficult to reflect the intrinsic causes and influencing factors that cause the nonlinear flow of gas in coal.To fill this research gap, Zhang and Cheng [53] took nonlinear seepage mechanism as the theoretical basis, combined with actual experimental research results, and established a nonlinear seepage equation describing the characteristics of gas seepage based on mechanical equilibrium equation; the physical meaning of which is as follows: the seepage resistance consists of two parts, the first part is viscous resistance, which is proportional to the primary side of the seepage velocity, and the second part is the additional seepage resistance caused by adsorption, which is proportional to the adsorbed gas content in the coal. It is precisely because of the existence of the intermolecular adsorption between coal and gas, which reduces the flow rate of gas in the coal body. The equation is the Darcy flow equation when the adsorbed gas content in the coal is zero. While most scholars currently use heat shrink tubing to seal the coal sample, using oil pressure to offset the air pressure to implement an annular seal on the sample so that the specimen is partially subjected to the oil pressure, this literature utilizes a sealing cylinder sleeve combined with a sealing reagent to implement a coal sample seal under constrained annular strain, fundamentally eliminating the influence of this factor on the experimental results.


#### 3.2.3. Ground-Field Effect Flow-Solid Coupling Gas Flow Theory

With the in-depth research on the flow of gas in porous media, domestic and foreign experts and scholars also gradually realize that the effect of ground-field effect on gas flow cannot be ignored, and based on Darcy’s law, a large number of analytical studies are conducted in gas percolation, and a mathematical model of gas percolation under the conditions of gas-solid coupling action is established.By varying the pressure difference, scholars determined the permeability of gas-bearing coal under different conditions, integrated a large amount of experimental data, and found that the permeability of the coal body gradually increased with the decrease of gas pressure under low-pressure conditions, and they defined this physical phenomenon as the Klinkenberg effect. It is also suggested that the mechanism underlying this effect is that the gas slips on the solid surface under low-pressure conditions, leading to an accelerated flow of gas in the coal body, which is macroscopically manifested as an increase in the permeability of the coal body. However, the literature [54] argues that this effect does not actually exist, and it is not the slip of gas molecules on the solid surface, but the adsorption between coal and gas molecules that leads to the change of permeability. Liu et al. [55] considered the Klinkenberg effect, the effective stress, and desorption shrinkage on gas seepage and coal deformation, established a fluid-solid coupling model describing gas seepage and coal skeleton deformability, and applied it to practical engineering. The results show that the Klinkenberg effect is more rapid in the vicinity of the extraction borehole than in the absence of the Klinkenberg effect, and the farther away from the extraction borehole, the less effect of the Klinkenberg effect.Experts and scholars at home and abroad have also conducted a lot of systematic analyses and studies on the deformation laws, mechanical properties, and their rheological characteristics of coal body samples containing gas. Wang et al. [56] studied the effect of electric field on the seepage properties of gas in coal using a triaxial seepage experiment device and an electric field implementation device, refined Darcy’s seepage law under the original conditions, proposed the basic law of gas seepage under the condition of ground electric field action, and established the corresponding mathematical model of gas seepage. Based on the research results of other scholars, Sun [57] and others created a gas-solid coupling theory of gas transmissive flow in the multicoal seam system.


## 3.1. Methodology for Studying Flow Field in Gob

### 3.1.1. Graphical Method and Filtered Flow Field Model

Early studies on the gas flow law in the gob usually analyzed the air leakage in the gob by considering the gob and the retreating working face as two parallel wind paths, and the model was applicable to the study of the U-type ventilation method without leakage sources and sinks after the retreating face, but the theory was not applicable to the situation where the number of leakage sources was too many or the geometry of the gob was complicated, which had obvious limitations. In 1978, Chinese scholars Huang and Wang [39], etc., set the quarry with the long-walled backward coal mining method as a standard rectangular body, set the air-gob as a homogeneous porous medium, considered that the air flow conditions on each plane parallel to the quarry floor are similar, and cited the theory of fluid mechanics about the vortex-free steady motion in the plane, in the case of one source and one sink, one-source and one-sink quarry means that the air is fed from the inlet down-slot. In the case of one-source and one-sink quarry, the air is fed from the inlet down-slot and discharged from the return down-slot, and there is no air leakage. The graphical method is used to describe the gas flow law in the rectangular gob, and the mine ventilation system is adjusted according to the graphical results. However, once the number of air leakage sources is too large or the geometry of the quarry is complicated, the graphical method is difficult to apply and even technically infeasible.In the same period, Polish and Soviet scholars put forward the filter flow field theory, which regarded the medium in the gob as a longitudinal and horizontal pipe network and believed that the flow of gas in the gob was similar to the flow in the pipe network. Shi et al. [40] used a computer to divide the fallout zone in the gob into a number of square grids, and each square grid was replaced by a filter flow branch network. The filter flow velocity and air leakage volume of the fallen strip in the gob were determined, and the computer simulation of the technical measures to reduce the air leakage volume in the gob was carried out according to the theory of filter flow field, and the results showed that the air leakage reduction ranged from 17.08% to 24.26%. Ding et al. [41] gridded the bubble fall zone in the gob, as shown in Figure 8, and established the mathematical model of the filter flow field in the gob, calculated the wind resistance of each branch of the filter flow field in the gob, and programmed a C language program to solve the filter flow field according to the law of wind volume balance and the law of wind pressure balance, as well as the iterative method of Scott-Hensley, to obtain the wind speed of each branch of the filter flow in the gob:(22)h=R′Q+R″Q2,where h is the resistance of the filtered flow field branches, Pa; R′ is the filtration flow field laminar wind resistance, kg/(m4•s); R″ is the filtration flow turbulent wind resistance, kg/m7; and Q is the air volume of filter flow branch, m3/s.Figure 8
Air leakage branches and node numbers in the extraction area.


### 3.1.2. Similarity Simulation Method

Similarity simulation (including physical similarity simulation and electrical simulation) is based on similarity theory and plays an important role in the early research stage of flow field in the gob. The so-called similarity principle means that when the mathematical models of two physical fields, the same or different, are isomorphic, the distribution laws of the two field quantities are the same, but with different signs such as functions or coefficients, and when the solution of one field has been solved or is easy to solve, the other unknown and relatively difficult field can be solved based on it. The similarity simulation method is relatively simple to use, and researchers can study different physical fields under various conditions by creating similar simulations in the laboratory.Electrical simulation is one of the simulation methods with more advantages. Zhu and Hui [42] introduced the electric simulation method into the study of mine ventilation, comparing the wind field with the electric field, and proved by practice that the introduction of the electric simulation method into coal mines is useful for the study of problems such as ventilation, but there are some problems with the electric simulation, such as the unsatisfactory treatment of the junction conditions in places where the field changes drastically, the poor mastery of the model dissection technique, the poor orthogonality in some parts of the flow network, and the existence of some errors.Xu et al. [43] used the physical similarity simulation method to realize the qualitative study of wind flow state by analyzing the smoke flow state. Yang et al. [44] established a model of coal seam permeability coefficient based on similar material experiments simulating coal mining with fractal dimension of mining fractures.


### 3.1.3. CFD Numerical Simulation Method

Numerical simulation, also known as computer simulation, is the study of engineering and physical problems by means of numerical computation and image display using an electronic computer as a carrier. CFD is a branch of computational fluid dynamics, where the mathematical description of the overall flow field of interest is obtained by numerically solving a system of equations to reflect the fluid flow, heat and mass transfer laws and related phenomena.With the continuous progress of computer science and technology, a lot of CFD commercial software has been created, such as PHOENICS and FLUENT. PHOENICS is the world’s first set of commercial software for computational fluids and computational heat transfer, which was developed by the main founder of international computational fluids and computational heat transfer, Professor D. B. Spalding, a member of the Royal Academy of Engineering, and more than 40 Ph. Although PHOENICS was introduced into the field of coal mine flow field earlier, the software is completely applied to the calculation of the flow field of rock fall media in the gob, there are still many problems, such as the adaptation and debugging of specific problems, with large differences and limitations, and there are still some problems to overcome in the integration of analysis and calculation with this professional problem. FLUENT is currently international more popular commercial CFD software package, with a high market share. It has rich physical models, advanced numerical methods, and powerful pre- and postprocessing functions and has a wide range of applications. It uses the finite volume method based on a completely unstructured grid and has a gradient algorithm based on grid nodes and grid cells.In recent years, CFD numerical simulation of wind flow motion and gas transport and its distribution pattern in mining sites has been carried out in China. Zhou et al. [45] used PHOENICS software to simulate the distribution law of the wind flow and temperature field at the back mining workings and concluded that the wind flow temperature at the back mining workings decreases with the increase of wind speed in a negative power function, increases linearly with the increase of incoming wind flow temperature, and varies with the length of the workings and the location of the coal miner. The numerical simulation program (G3) of the field flow in the gob prepared by Li [23] has the following characteristics: (1) it can adapt to various complex boundary conditions and treat the gob as a bubbling nonhomogeneous medium; (2) it takes into account the natural coupling effect when gas gushes out and can reflect the relationship between the influence of factors such as working face advancement and ventilation volume; (3) it is operable and can make arbitrary simulation test faces for various situations, which is convenient. Theoretically, it can depict the flow state of air leakage in the gob, dynamically depict the distribution state of gas, oxygen, carbon monoxide concentration and temperature, and their change process, and give the graphical distribution solution of each phase of gas and temperature change in the gob under the mining conditions of the working face; (4) it realizes the integration of the joint solution of “a multiphase problem” in the gob, which can meet the needs of the study of the problem. G3 can assist in analyzing many safety problems such as gas emission from the gob, spontaneous combustion and fire, nitrogen injection and fire extinguishing, and large area pressure equalization and pressure regulation. It is powerful and comprehensive and represents the development direction of domestic application software in the field of numerical simulation of field flow in the gob.


## 3.1.1. Graphical Method and Filtered Flow Field Model

Early studies on the gas flow law in the gob usually analyzed the air leakage in the gob by considering the gob and the retreating working face as two parallel wind paths, and the model was applicable to the study of the U-type ventilation method without leakage sources and sinks after the retreating face, but the theory was not applicable to the situation where the number of leakage sources was too many or the geometry of the gob was complicated, which had obvious limitations. In 1978, Chinese scholars Huang and Wang [39], etc., set the quarry with the long-walled backward coal mining method as a standard rectangular body, set the air-gob as a homogeneous porous medium, considered that the air flow conditions on each plane parallel to the quarry floor are similar, and cited the theory of fluid mechanics about the vortex-free steady motion in the plane, in the case of one source and one sink, one-source and one-sink quarry means that the air is fed from the inlet down-slot. In the case of one-source and one-sink quarry, the air is fed from the inlet down-slot and discharged from the return down-slot, and there is no air leakage. The graphical method is used to describe the gas flow law in the rectangular gob, and the mine ventilation system is adjusted according to the graphical results. However, once the number of air leakage sources is too large or the geometry of the quarry is complicated, the graphical method is difficult to apply and even technically infeasible.In the same period, Polish and Soviet scholars put forward the filter flow field theory, which regarded the medium in the gob as a longitudinal and horizontal pipe network and believed that the flow of gas in the gob was similar to the flow in the pipe network. Shi et al. [40] used a computer to divide the fallout zone in the gob into a number of square grids, and each square grid was replaced by a filter flow branch network. The filter flow velocity and air leakage volume of the fallen strip in the gob were determined, and the computer simulation of the technical measures to reduce the air leakage volume in the gob was carried out according to the theory of filter flow field, and the results showed that the air leakage reduction ranged from 17.08% to 24.26%. Ding et al. [41] gridded the bubble fall zone in the gob, as shown in Figure 8, and established the mathematical model of the filter flow field in the gob, calculated the wind resistance of each branch of the filter flow field in the gob, and programmed a C language program to solve the filter flow field according to the law of wind volume balance and the law of wind pressure balance, as well as the iterative method of Scott-Hensley, to obtain the wind speed of each branch of the filter flow in the gob:(22)h=R′Q+R″Q2,where h is the resistance of the filtered flow field branches, Pa; R′ is the filtration flow field laminar wind resistance, kg/(m4•s); R″ is the filtration flow turbulent wind resistance, kg/m7; and Q is the air volume of filter flow branch, m3/s.Figure 8
Air leakage branches and node numbers in the extraction area.


## 3.1.2. Similarity Simulation Method

Similarity simulation (including physical similarity simulation and electrical simulation) is based on similarity theory and plays an important role in the early research stage of flow field in the gob. The so-called similarity principle means that when the mathematical models of two physical fields, the same or different, are isomorphic, the distribution laws of the two field quantities are the same, but with different signs such as functions or coefficients, and when the solution of one field has been solved or is easy to solve, the other unknown and relatively difficult field can be solved based on it. The similarity simulation method is relatively simple to use, and researchers can study different physical fields under various conditions by creating similar simulations in the laboratory.Electrical simulation is one of the simulation methods with more advantages. Zhu and Hui [42] introduced the electric simulation method into the study of mine ventilation, comparing the wind field with the electric field, and proved by practice that the introduction of the electric simulation method into coal mines is useful for the study of problems such as ventilation, but there are some problems with the electric simulation, such as the unsatisfactory treatment of the junction conditions in places where the field changes drastically, the poor mastery of the model dissection technique, the poor orthogonality in some parts of the flow network, and the existence of some errors.Xu et al. [43] used the physical similarity simulation method to realize the qualitative study of wind flow state by analyzing the smoke flow state. Yang et al. [44] established a model of coal seam permeability coefficient based on similar material experiments simulating coal mining with fractal dimension of mining fractures.


## 3.1.3. CFD Numerical Simulation Method

Numerical simulation, also known as computer simulation, is the study of engineering and physical problems by means of numerical computation and image display using an electronic computer as a carrier. CFD is a branch of computational fluid dynamics, where the mathematical description of the overall flow field of interest is obtained by numerically solving a system of equations to reflect the fluid flow, heat and mass transfer laws and related phenomena.With the continuous progress of computer science and technology, a lot of CFD commercial software has been created, such as PHOENICS and FLUENT. PHOENICS is the world’s first set of commercial software for computational fluids and computational heat transfer, which was developed by the main founder of international computational fluids and computational heat transfer, Professor D. B. Spalding, a member of the Royal Academy of Engineering, and more than 40 Ph. Although PHOENICS was introduced into the field of coal mine flow field earlier, the software is completely applied to the calculation of the flow field of rock fall media in the gob, there are still many problems, such as the adaptation and debugging of specific problems, with large differences and limitations, and there are still some problems to overcome in the integration of analysis and calculation with this professional problem. FLUENT is currently international more popular commercial CFD software package, with a high market share. It has rich physical models, advanced numerical methods, and powerful pre- and postprocessing functions and has a wide range of applications. It uses the finite volume method based on a completely unstructured grid and has a gradient algorithm based on grid nodes and grid cells.In recent years, CFD numerical simulation of wind flow motion and gas transport and its distribution pattern in mining sites has been carried out in China. Zhou et al. [45] used PHOENICS software to simulate the distribution law of the wind flow and temperature field at the back mining workings and concluded that the wind flow temperature at the back mining workings decreases with the increase of wind speed in a negative power function, increases linearly with the increase of incoming wind flow temperature, and varies with the length of the workings and the location of the coal miner. The numerical simulation program (G3) of the field flow in the gob prepared by Li [23] has the following characteristics: (1) it can adapt to various complex boundary conditions and treat the gob as a bubbling nonhomogeneous medium; (2) it takes into account the natural coupling effect when gas gushes out and can reflect the relationship between the influence of factors such as working face advancement and ventilation volume; (3) it is operable and can make arbitrary simulation test faces for various situations, which is convenient. Theoretically, it can depict the flow state of air leakage in the gob, dynamically depict the distribution state of gas, oxygen, carbon monoxide concentration and temperature, and their change process, and give the graphical distribution solution of each phase of gas and temperature change in the gob under the mining conditions of the working face; (4) it realizes the integration of the joint solution of “a multiphase problem” in the gob, which can meet the needs of the study of the problem. G3 can assist in analyzing many safety problems such as gas emission from the gob, spontaneous combustion and fire, nitrogen injection and fire extinguishing, and large area pressure equalization and pressure regulation. It is powerful and comprehensive and represents the development direction of domestic application software in the field of numerical simulation of field flow in the gob.


## 3.2. Numerical Model Equation Study

### 3.2.1. Linear Gas Seepage Theory

The numerical model equations were initially simplified according to Darcy’s law, which is a typical representative of linear gas seepage theory, and it is believed that the gas flow law in coal seams is basically in line with the linear percolation law. In the international aspect, the former Soviet Union scholars used Darcy’s law for the first time to analyze and study the problem of linear gas seepage in coal seams by using the nature of gas adsorption, which opened a new chapter for the study of linear gas seepage theory. In China, experts such as Zhou and Sun [46] agreed with the linear gas seepage theory, and on this basis, they analyzed and studied the flow of gas in coal seams for the first time, and regarded coal seams as continuous and uniformly distributed porous media, which laid the theoretical foundation for the analysis and study of gas flow theory by domestic experts and scholars. The linear seepage theory believes that the flow rate of gas in coal and the pressure difference obey a linear mathematical relationship; that is, the flow rate is proportional to the pressure difference of gas. Since the 1980s, scholars in China have carried out research on the numerical equations of linear flow theory and have modified and improved them on this basis.In 1986, Tan and Yuan [47] combined the theory of seepage mechanics and thermodynamics, treated coal seam gas as real gas, proposed the seepage equation of real gas in the coal seam of the mine, and simplified the equation. In 1989, Yu et al. [48] proposed the idea that the amount of gas involved in seepage in coal seam is the partial amount of gas content in coal body, and under the assumption that the process of gas adsorption and desorption in coal body is completely reversible, the controlling equation of gas seepage in the coal seam was established; in 2010, Yang and Li [49] established a mathematical model of gas flow in the coal body around the borehole based on Darcy’s law and the assumptions of coal seam gas flow theory and applied the law of mass conservation to obtain the equation of gas flow around the borehole and its fixed solution conditions.


### 3.2.2. Nonlinear Gas Seepage Theory

Darcy’s law has a good fit for laminar flow motion that conforms to the linear drag relationship, but in practice the flow of wind flow in the gob is often very complex, with a large portion of the area being transitional seepage or turbulent flow. In these regions, fluid motion does not conform to Darcy’s law. Therefore, many scholars have proposed nonlinear seepage laws based on this problem and conducted research based on nonlinear seepage theory.In 1984, Higuchi Sumitoshi, a professor at Hokkaido University, Japan, determined the methane permeability of coal samples by varying the pressure difference and found that Darcy’s law was not quite consistent with the coal seam gas flow law. Based on a large number of experimental studies, he proposed the basic law of coal seam gas flow that is more consistent with the power law:(23)VN=−Adpdxm,where VN is the uncaused gas flow rate in the corresponding standard state; A is the uncaused gas permeability coefficient; m is the state constant; and dp/dx is the uncaused gas pressure gradient along the X-axis.Sun [50] considered that the mathematical model of gas flow proposed by Professor Higuchi was not rigorous and concluded that the gas flow in coal seams is actually a mixed unsteady flow of compressible fluid by permeation-diffusion in an anisotropic and nonhomogeneous pore-fissure twofold medium. He established the first partial differential equation for compressible gas flow in coal seams under the conditions of homogeneous and nonhomogeneous coal seams according to the generalized form of power law and made practical calculations to verify the flow of gas flow field in homogeneous coal seams based on the actual measured gas flow parameters in 23051 working face of the Jiaozuo Zhongmacun Mine, and the results showed that it was closer to the actual situation than Darcy’s law. In 1994, Liu [51] concluded that the coal seam gas flow model derived in the literature [41] on the basis of the power law and the application of the invariant had errors in the process, and the numerical results obtained by it were unreliable, and derived the correct flow model when the gas flow followed the power law in the paper.Deng et al. [52] analyzed a large amount of experimental data, on this basis generalized the universal characteristic equation of motion for seepage curves, and proposed a nonlinear equation of motion for gas seepage in low-permeability pore-fracture media. However, this type of model mainly adopts the method of phenomenology, which is only a generalization and distillation of the experimental phenomena and does not have the function of in-depth explanation, which is manifested in this study mainly by scholars fitting mathematical equations according to the physical phenomena of nonlinear seepage, and it is difficult to reflect the intrinsic causes and influencing factors that cause the nonlinear flow of gas in coal.To fill this research gap, Zhang and Cheng [53] took nonlinear seepage mechanism as the theoretical basis, combined with actual experimental research results, and established a nonlinear seepage equation describing the characteristics of gas seepage based on mechanical equilibrium equation; the physical meaning of which is as follows: the seepage resistance consists of two parts, the first part is viscous resistance, which is proportional to the primary side of the seepage velocity, and the second part is the additional seepage resistance caused by adsorption, which is proportional to the adsorbed gas content in the coal. It is precisely because of the existence of the intermolecular adsorption between coal and gas, which reduces the flow rate of gas in the coal body. The equation is the Darcy flow equation when the adsorbed gas content in the coal is zero. While most scholars currently use heat shrink tubing to seal the coal sample, using oil pressure to offset the air pressure to implement an annular seal on the sample so that the specimen is partially subjected to the oil pressure, this literature utilizes a sealing cylinder sleeve combined with a sealing reagent to implement a coal sample seal under constrained annular strain, fundamentally eliminating the influence of this factor on the experimental results.


### 3.2.3. Ground-Field Effect Flow-Solid Coupling Gas Flow Theory

With the in-depth research on the flow of gas in porous media, domestic and foreign experts and scholars also gradually realize that the effect of ground-field effect on gas flow cannot be ignored, and based on Darcy’s law, a large number of analytical studies are conducted in gas percolation, and a mathematical model of gas percolation under the conditions of gas-solid coupling action is established.By varying the pressure difference, scholars determined the permeability of gas-bearing coal under different conditions, integrated a large amount of experimental data, and found that the permeability of the coal body gradually increased with the decrease of gas pressure under low-pressure conditions, and they defined this physical phenomenon as the Klinkenberg effect. It is also suggested that the mechanism underlying this effect is that the gas slips on the solid surface under low-pressure conditions, leading to an accelerated flow of gas in the coal body, which is macroscopically manifested as an increase in the permeability of the coal body. However, the literature [54] argues that this effect does not actually exist, and it is not the slip of gas molecules on the solid surface, but the adsorption between coal and gas molecules that leads to the change of permeability. Liu et al. [55] considered the Klinkenberg effect, the effective stress, and desorption shrinkage on gas seepage and coal deformation, established a fluid-solid coupling model describing gas seepage and coal skeleton deformability, and applied it to practical engineering. The results show that the Klinkenberg effect is more rapid in the vicinity of the extraction borehole than in the absence of the Klinkenberg effect, and the farther away from the extraction borehole, the less effect of the Klinkenberg effect.Experts and scholars at home and abroad have also conducted a lot of systematic analyses and studies on the deformation laws, mechanical properties, and their rheological characteristics of coal body samples containing gas. Wang et al. [56] studied the effect of electric field on the seepage properties of gas in coal using a triaxial seepage experiment device and an electric field implementation device, refined Darcy’s seepage law under the original conditions, proposed the basic law of gas seepage under the condition of ground electric field action, and established the corresponding mathematical model of gas seepage. Based on the research results of other scholars, Sun [57] and others created a gas-solid coupling theory of gas transmissive flow in the multicoal seam system.


## 3.2.1. Linear Gas Seepage Theory

The numerical model equations were initially simplified according to Darcy’s law, which is a typical representative of linear gas seepage theory, and it is believed that the gas flow law in coal seams is basically in line with the linear percolation law. In the international aspect, the former Soviet Union scholars used Darcy’s law for the first time to analyze and study the problem of linear gas seepage in coal seams by using the nature of gas adsorption, which opened a new chapter for the study of linear gas seepage theory. In China, experts such as Zhou and Sun [46] agreed with the linear gas seepage theory, and on this basis, they analyzed and studied the flow of gas in coal seams for the first time, and regarded coal seams as continuous and uniformly distributed porous media, which laid the theoretical foundation for the analysis and study of gas flow theory by domestic experts and scholars. The linear seepage theory believes that the flow rate of gas in coal and the pressure difference obey a linear mathematical relationship; that is, the flow rate is proportional to the pressure difference of gas. Since the 1980s, scholars in China have carried out research on the numerical equations of linear flow theory and have modified and improved them on this basis.In 1986, Tan and Yuan [47] combined the theory of seepage mechanics and thermodynamics, treated coal seam gas as real gas, proposed the seepage equation of real gas in the coal seam of the mine, and simplified the equation. In 1989, Yu et al. [48] proposed the idea that the amount of gas involved in seepage in coal seam is the partial amount of gas content in coal body, and under the assumption that the process of gas adsorption and desorption in coal body is completely reversible, the controlling equation of gas seepage in the coal seam was established; in 2010, Yang and Li [49] established a mathematical model of gas flow in the coal body around the borehole based on Darcy’s law and the assumptions of coal seam gas flow theory and applied the law of mass conservation to obtain the equation of gas flow around the borehole and its fixed solution conditions.


## 3.2.2. Nonlinear Gas Seepage Theory

Darcy’s law has a good fit for laminar flow motion that conforms to the linear drag relationship, but in practice the flow of wind flow in the gob is often very complex, with a large portion of the area being transitional seepage or turbulent flow. In these regions, fluid motion does not conform to Darcy’s law. Therefore, many scholars have proposed nonlinear seepage laws based on this problem and conducted research based on nonlinear seepage theory.In 1984, Higuchi Sumitoshi, a professor at Hokkaido University, Japan, determined the methane permeability of coal samples by varying the pressure difference and found that Darcy’s law was not quite consistent with the coal seam gas flow law. Based on a large number of experimental studies, he proposed the basic law of coal seam gas flow that is more consistent with the power law:(23)VN=−Adpdxm,where VN is the uncaused gas flow rate in the corresponding standard state; A is the uncaused gas permeability coefficient; m is the state constant; and dp/dx is the uncaused gas pressure gradient along the X-axis.Sun [50] considered that the mathematical model of gas flow proposed by Professor Higuchi was not rigorous and concluded that the gas flow in coal seams is actually a mixed unsteady flow of compressible fluid by permeation-diffusion in an anisotropic and nonhomogeneous pore-fissure twofold medium. He established the first partial differential equation for compressible gas flow in coal seams under the conditions of homogeneous and nonhomogeneous coal seams according to the generalized form of power law and made practical calculations to verify the flow of gas flow field in homogeneous coal seams based on the actual measured gas flow parameters in 23051 working face of the Jiaozuo Zhongmacun Mine, and the results showed that it was closer to the actual situation than Darcy’s law. In 1994, Liu [51] concluded that the coal seam gas flow model derived in the literature [41] on the basis of the power law and the application of the invariant had errors in the process, and the numerical results obtained by it were unreliable, and derived the correct flow model when the gas flow followed the power law in the paper.Deng et al. [52] analyzed a large amount of experimental data, on this basis generalized the universal characteristic equation of motion for seepage curves, and proposed a nonlinear equation of motion for gas seepage in low-permeability pore-fracture media. However, this type of model mainly adopts the method of phenomenology, which is only a generalization and distillation of the experimental phenomena and does not have the function of in-depth explanation, which is manifested in this study mainly by scholars fitting mathematical equations according to the physical phenomena of nonlinear seepage, and it is difficult to reflect the intrinsic causes and influencing factors that cause the nonlinear flow of gas in coal.To fill this research gap, Zhang and Cheng [53] took nonlinear seepage mechanism as the theoretical basis, combined with actual experimental research results, and established a nonlinear seepage equation describing the characteristics of gas seepage based on mechanical equilibrium equation; the physical meaning of which is as follows: the seepage resistance consists of two parts, the first part is viscous resistance, which is proportional to the primary side of the seepage velocity, and the second part is the additional seepage resistance caused by adsorption, which is proportional to the adsorbed gas content in the coal. It is precisely because of the existence of the intermolecular adsorption between coal and gas, which reduces the flow rate of gas in the coal body. The equation is the Darcy flow equation when the adsorbed gas content in the coal is zero. While most scholars currently use heat shrink tubing to seal the coal sample, using oil pressure to offset the air pressure to implement an annular seal on the sample so that the specimen is partially subjected to the oil pressure, this literature utilizes a sealing cylinder sleeve combined with a sealing reagent to implement a coal sample seal under constrained annular strain, fundamentally eliminating the influence of this factor on the experimental results.


## 3.2.3. Ground-Field Effect Flow-Solid Coupling Gas Flow Theory

With the in-depth research on the flow of gas in porous media, domestic and foreign experts and scholars also gradually realize that the effect of ground-field effect on gas flow cannot be ignored, and based on Darcy’s law, a large number of analytical studies are conducted in gas percolation, and a mathematical model of gas percolation under the conditions of gas-solid coupling action is established.By varying the pressure difference, scholars determined the permeability of gas-bearing coal under different conditions, integrated a large amount of experimental data, and found that the permeability of the coal body gradually increased with the decrease of gas pressure under low-pressure conditions, and they defined this physical phenomenon as the Klinkenberg effect. It is also suggested that the mechanism underlying this effect is that the gas slips on the solid surface under low-pressure conditions, leading to an accelerated flow of gas in the coal body, which is macroscopically manifested as an increase in the permeability of the coal body. However, the literature [54] argues that this effect does not actually exist, and it is not the slip of gas molecules on the solid surface, but the adsorption between coal and gas molecules that leads to the change of permeability. Liu et al. [55] considered the Klinkenberg effect, the effective stress, and desorption shrinkage on gas seepage and coal deformation, established a fluid-solid coupling model describing gas seepage and coal skeleton deformability, and applied it to practical engineering. The results show that the Klinkenberg effect is more rapid in the vicinity of the extraction borehole than in the absence of the Klinkenberg effect, and the farther away from the extraction borehole, the less effect of the Klinkenberg effect.Experts and scholars at home and abroad have also conducted a lot of systematic analyses and studies on the deformation laws, mechanical properties, and their rheological characteristics of coal body samples containing gas. Wang et al. [56] studied the effect of electric field on the seepage properties of gas in coal using a triaxial seepage experiment device and an electric field implementation device, refined Darcy’s seepage law under the original conditions, proposed the basic law of gas seepage under the condition of ground electric field action, and established the corresponding mathematical model of gas seepage. Based on the research results of other scholars, Sun [57] and others created a gas-solid coupling theory of gas transmissive flow in the multicoal seam system.


## 4. The Influence of Environmental Factors on the Distribution of Flow Field in the Gob

### 4.1. Ventilation Method

According to the actual mining situation and considering the geological characteristics, the selection of working face ventilation method varies. At present, U-type ventilation is most commonly used in coal mines. With the increasing maturity of coal mining theory and to meet the demand for safe mining under different geological conditions, various ventilation methods have been developed based on U-type ventilation. In recent years, scholars have made specific studies on the distribution of flow field in the gob under different ventilation methods.


#### 4.1.1. U-Type Ventilation

The U-type ventilation system is mainly used in coal mining workings of low gas mines, where there is only one air inlet and one air return lane, and it can be divided into U-back and U-forward according to the direction of workings. U-type backward ventilation is shown in Figure9. In China’s coal mines, the U-back type is commonly used in the working face. The gas in the upper corner of U-ventilated working face of high gas mine has been a difficult problem for coal mine gas management. With the development of coal science and technology, coal mining working face has realized high-production and high-efficiency comprehensive mechanized production, and the unit production of working face has been improved continuously, and the gas gushing out has been increased, so that the mines with little gas in the original working face and low gas mines also have the upper corner gas accumulation over limit [58]. Many scholars have analyzed the wind flow and gas flow law of U-type ventilation working face through on-site observation and research.Figure 9
U-type backward ventilation.In order to obtain the flow field distribution in the mining void area of U-type comprehensive mining working face, Tao [59] established a two-dimensional physical model of the mining void area according to the actual situation of 11124 working face and conducted numerical simulation studies on the constructed model using FLUENT software and adaptive grid encryption technology. The study showed that the air leakage in the gob mainly occurred in the range of 0∼20 m, and some of the air leakage into the gob returned to the working face in this area, while most of the air leakage returned to the working face in the inclination of 140∼160 m; the wind velocity contours in the gob were not symmetrical in the inclination direction; in the direction of the gob, the wind velocity of air leakage showed a gradually decreasing trend.Tian et al. [60] arranged one release point and five sampling points in the working face, used the instantaneous release method to continuously release SF6 gas at the wind source location, conducted gas sampling at the sampling points at regular intervals to analyze the SF6 gas concentration in the gas, estimated the air leakage volume and minimum air leakage velocity in the gob according to the time and concentration of SF6 in the gas sample, and analyzed the measured data to make a general inference about the location of air leakage in the working face. Through the analysis of the measured data, the location of the air leakage channel in the working face can be roughly deduced, as shown in Figure 10. After that, the numerical simulation software FLUENT was used to verify the air leakage in the U-type ventilation gob. The actual measurement results of the corresponding wind speed and air leakage velocity of each monitoring point are compared with the numerical simulation results as shown in Figure 11. The numerical simulation results are consistent with the actual measurement results on-site. The air leakage in the gob is distributed in a U-type, which flows into the gob from the inlet side of the working face and deflects to the return side after a certain depth and flows out; the air leakage in the gob of the U-type ventilation system is mainly in the lower corner area of the working face, which is more prone to spontaneous combustion.Figure 10
Schematic diagram of air leakage channel in gob of 15102 working face.Figure 11
Comparison of field measurement and numerical simulation results. (a) Monitor the corresponding location. (b) Monitor the corresponding location.
(a)
(b)With the increase of coal mining at the working face, the single U-type ventilation system cannot meet the safety production of coal mining enterprises, and in order to keep the gas concentration at the working face from exceeding the limit, new ventilation methods have emerged based on the U-type ventilation system, among which U + I-type and U + L-type are more widely used at the comprehensive working face. In U + I-type ventilation, the gas tail lane is arranged on the inner side of the return air lane, and the gas carried by the leakage is discharged through the gas tail lane; while in U + L-type ventilation, the gas tail lane is arranged on the outer side of the return air lane. Although these two ventilation methods are good for controlling the gas concentration in the corner of the working face, both of them change the air leakage in the gob compared with the U-type ventilation method, causing the flow field in the gob to change.Ma et al. [61] used the UDF interface to establish the calculation models of coal gas seepage in the mined area based on the moving coordinates of the working face and the gas dispersion and transport in the retracement tunnel and quantitatively compared the U-type and U + I-type ventilation methods under different working face advance speed and working face air distribution conditions. The paper concluded that the location of air leakage in U-type ventilation was mainly concentrated in the lower corner of the working face, while in U + I-type ventilation, there were air leakage channels in the whole working face; the size of gas concentration in the return air lane under U-type ventilation is about twice that of U + I-type ventilation, and with the increase of wind speed in the inlet air lane, U + I-type ventilation and U-type ventilation have more gas concentration. With the increase of wind speed in the inlet lane, the gas concentration in the return lane of both U + I-type and U-type ventilation methods shows a decrease in the form of decrease, and with the increase of working face advancement speed, the gas concentration in the return lane of U-type ventilation method shows a decrease in the form of increase, while the gas concentration in the return lane of U + I-type ventilation method shows a decrease in the form of decrease. Compared with the U-type ventilation method, the gas concentration in the upper corner can be reduced by 25% to 50% under the U + I-type ventilation method.Wu et al. [62] conducted numerical simulation of the flow field in the gob under U-type and U + L-type ventilation using FLUENT software and compared and analyzed the distribution of gas concentration in the gob and the distribution of spontaneous combustion hazard zone in the gob under the two ventilation methods, as shown in Figure 12, and the results showed that the gas concentration in the gob of U + L-type ventilation was lower than that of U-type ventilation and the gas explosion range in the gob was lower than that of U-type ventilation. The range of gas explosion moves to the deep part of the gob, but the width of gas explosion area becomes smaller. The gas concentration in the middle of the gob increases faster along the direction of the gob, and the spontaneous combustion hazard area in the gob moves deeper in the gob with U + L-type ventilation than U-type ventilation, and the width of the spontaneous combustion hazard area in the gob also increases significantly.Figure 12
Simulation results of flow field in U-type and U + L-type ventilation gob. (a) Gas concentration distribution in U-type ventilation gob. (b) Oxygen concentration distribution in U-type ventilation gob. (c) Gas concentration distribution in U + L-type ventilation gob. (d) Oxygen concentration distribution in U + L-type ventilation gob.
(a)
(b)
(c)
(d)


#### 4.1.2. Y-Type Ventilation

The Y-type ventilation method is a ventilation system that adds fresh air flow on the return side and flows out from the side of the gob after converging with the lacking air from the working face. The use of Y-type ventilation method in the coal face increases the air volume in the return air tunnel, which can effectively solve the problem of gas overrun in the upper corner and is conducive to gas extraction, so it has good promotion value, as shown in Figure13.Figure 13
Y-type backward ventilation.Yu et al. [63] established a two-dimensional solution model for the flow field in the gob under Y-type ventilation, applied the computational fluid dynamics software FLUENT to numerically simulate the distribution of the air leakage field and gas concentration field in the gob, conducted a comprehensive comparison of the air distribution ratios of different main and subinlet lanes, and concluded that the gas management effect was best when the air distribution ratio was 5 : 1.Yang et al. [64] established the physical model of U-type and Y-type ventilation air-gob and used FLUENT software to numerically simulate the air leakage flow field, air leakage volume, and gas distribution in the air-gob under these two ventilation methods, and the simulation results are shown in Figure 14. Most of the air flow in U-type ventilation flows into the working face through the inlet alley and finally discharges through the return alley, and a small part of the air flow leaks into the air-gob from the lower corner of the working face, and in Y-type ventilation, most of the air from the machine tunnel flows through the working face, then converges with the air from the wind tunnel, and enters the air-retention lane, finally discharges through the return air lane, and some of the air leaks into the air-gob from the lower corner of the working face to the middle of the working face, and the closer the air leaks into the air-gob, the more it blows to the bottom of the air-gob, and then it flows from the air-gob to the interface between the air-gob and the air-retention lane. The more the air leaks into the air-gob from the side of the machine tunnel, the more it can be blown to the bottom of the air-gob and then flow into the air-gob from the interface between the air-gob and the alley. After analyzing the air leakage data of the two types of ventilation, it is concluded that the air leakage at the end of the working face from 0 to 30m accounts for 50% of the air leakage at the working face with Y-type ventilation, and the air leakage at the working face to the air-gob from 30 to 90m accounts for 46% of the total air leakage, and the total air leakage is more than that of U-type ventilation, which is not suitable for the coal mining face with coal spontaneous combustion. In addition, under Y-type ventilation, most of the gas in the gob is carried by the air leakage from the working face to the gob to be discharged along the air-retention lane, and very little of the low volume fraction gas (≤0.5%) is gushed to the corner of the working face, so the volume fraction of gas in the gob is lower than that in U-type ventilation, which can avoid the accumulation of high concentration gas in the gob. The use of two-in-one Y-type ventilation can fundamentally solve the problem of gas accumulation in the upper corner and gas overload in the return airway.Figure 14
Simulation results of flow field in U-ventilated and Y-ventilated gob. (a) Flow diagram of the U-type ventilation and recovery workings in the mining void area. (b) Flow diagram of gob of Y-type ventilation and retrieval working face. (c) Gas distribution map of U-type ventilation gob. (d) Y-type ventilation gob gas distribution map.
(a)
(b)
(c)
(d)Zhao et al. [65] took 5301 working face of the Sihe Mine of the Jin Coal Group as the research object and studied the gas gushing law in the gob under partial Y-type ventilation according to the measured parameters of gas gushing, and two gas gushing were considered in the literature, which were coal left in the gob and gas gushing from the adjacent layer on the gob, set at 0.3 m from the bottom plate of the gob and the adjacent layer area of the gob (19.28 to 19.80 m from the bottom plate of the gob). Bottom plate was 19.28∼19.80 m height range, and different locations in the gob defined different gas outflow rate; according to the O-ring theory, we set the gob void rate and permeability.


#### 4.1.3. W-Type Ventilation

W-type ventilation is mainly used in long working face with high gas, where both the inlet and return airway are in the coal body. In order to solve the problem of gas overrun in the upper corner in U-type ventilation, two-in-one W-type ventilation is often used with the upper and lower flat as the inlet airway and the middle roadway as the return airway. Because of the maintenance difficulties of W-type forward ventilation tunnel, the air leakage is relatively large, and the gas gushing from the gob is also large, W-type backward ventilation is generally used, as shown in Figure15.Figure 15
W-type backward ventilation.Li and Li [66] established a finite element numerical model of the flow pattern of wind leakage in the bubble fall gob based on the equation of wind leakage and seepage in nonhomogeneous porous media, combined with field examples, and the calculation was based on MATLAB to develop a calculation program to display the results graphically. The sparse and dense transition is carried out by using the frontier generation method for the sectioning of the gob and the accuracy of the boundary sectioning of the wind leakage condition in the quarry as the basis of the regional grid density function; the spatial shape of the area is reflected by the change of the flow field height. Theoretically, the wind pressure distribution contours and flow function lines (distribution solutions) of the W-type ventilation gob are depicted, and it is pointed out that the W-type ventilation gob obviously reduces the wind pressure difference between the two ends of the working face and the effect of reducing the range of wind leakage.Cheng and Xin [67] compared the W-type ventilation method with U-type ventilation, and the simulation results showed that the leakage wind flow fields of the two are obviously different, and under U-type ventilation conditions, most of the leakage wind in the gob penetrates into the deep part of the gob, while W-type ventilation reduces the pressure difference between the two ends of the working face by about 35% compared with that of U-type due to the effect of equal pressure at both ends of the working face, causing the leakage wind in the gob to be closer to the working face under W-type ventilation. At the same time, because the air leakage will carry a lot of gas from the gob into the working face, W-type ventilation reduces the air leakage from the working face to the gob, thus suppressing the gas leakage from the gob to the working face, and fundamentally solving the problem of frequent overlimit gas accumulation in the upper corner of U-type ventilation and the return airway.


#### 4.1.4. J-Type Ventilation

This type of ventilation is a new type of ventilation system proposed only in recent years. It has a stronger and more controllable gas discharge capacity than the U-type ventilation commonly used in coal mines and is mainly suitable for high gas discharge workings in coal seams that are not prone to spontaneous combustion, which can well solve the corner gas problem on the workings, as shown in Figure16.Figure 16
J-type backward ventilation.According to the actual situation of 5201 header working face of the Wangzhuang Coal Mine in Lu’an, Wang and Wu [68] compared and analyzed the flow field and gas transport characteristics of J-type ventilation and U-type ventilation by using PHOENICS software. The simulation results of U-type ventilation and J-ventilation extraction zone are shown in Figures 17 and 18, respectively. The flow field and gas transport law of J-type ventilation gob were studied systematically under different conditions of gob size, gas discharge capacity of special alley, and extraction fan installation position. The results show that the J-type ventilation system can eliminate the concentrated air leakage in the upper corner of the gob and effectively solve the problem of gas accumulation in the upper corner of U-type ventilation.Figure 17
Simulation results of flow field andφ (CH4) in U-type ventilation gob. (a) Flow line and pressure distribution. (b) φ (CH4) distribution.
(a)
(b)Figure 18
Simulation results of flow field andφ (CH4) in J-type ventilation gob. (a) Flow line and pressure distribution. (b) φ (CH4) distribution.
(a)
(b)Gao and Sun [69] investigated the air leakage pattern in the gob, the distribution pattern of gas in the gob and the upper corner under different air volume conditions in the J-type ventilation system by means of numerical simulation. With the increase of inlet air volume, the area of low concentration of gas in the gob gradually deepens to the upper right of the gob, and the distribution range increases significantly. The magnitude of gas concentration reduction in the deep gob is larger than that in the shallow gob. The gas concentration in the upper corner also decreases gradually with the increase of the inlet air volume, and the gas concentration in the upper corner is lower than the 1% upper limit standard stipulated in the safety regulations when the air volume increases to 2350 m3/min.


### 4.2. Coal Seam Dip Angle

The geological conditions of China’s coal mines are complex and variable, and the coal seam inclination is not the same even for the same well field. The coal seam inclination is the main factor affecting the changing damage form of the overlying rock layer in the mining field. With the development of coal seam inclination from gentle to sharp inclination, the deformation and damage characteristics of the roof after coal seam mining will change, and the porous media parameters in the gob are closely related to it, which in turn affects the leakage flow field and the distribution of gas and pressure in the gob.Lu and Zhang [70] used UDEC numerical simulation to study the law of overlying rock strata spanning fall in the mining void zone of inclined extrathick coal seam and obtained the distribution of porosity in the mining void zone by the amount of overlying rock subsidence in the mining void zone. The results show that the porosity in the fallout zone is maximum along the strike at the edge of the gob and minimum at the deep part; it is maximum at the upper corner and lower corner along the tendency and minimum at the middle and upper part; in the vertical direction, the porosity in the gob decreases with increasing height from the bottom plate, and the difference between the edge and the deep part also decreases gradually.Yang and Zhang [71] studied the effect of coal seam inclination on the flow field and gas distribution in the mining void area of U + L-type ventilation face in a more systematic way, which divided the mining void area into bubble zone and fracture zone, where the bubble zone was subdivided into free accumulation zone, load-influenced zone, and compaction stability zone, set the porosity of the bubble zone to be uniformly distributed in segments, and investigated the effect of the coal seam inclination of 5.75°, 35.75°, with the increase of the coal seam inclination, the pressure difference between the two ends of the working face and the pressure difference between the working face and the tailing lane increase gradually, and the increase increases first and then decreases; the increase of the pressure difference between the two ends of the working face is greater than the pressure difference between the working face and the tailing lane; under the condition of keeping the air volume in the inlet lane, the return lane and the tailing lane unchanged, and the air volume in the working face decreases with the increase of the coal seam inclination. The air volume decreases, but the proportion of decrease is small and the trend is decreasing, and the distribution law of air leakage shows that the air leakage at the upper and lower ends is large, and the air leakage at the middle of the working face is small. As the coal seam inclination increases, the gas in the gob is affected by the uplift effect and the leakage flow, and the gas concentration gradually decreases. The gas concentration gradually decreases due to the uplift effect and leakage flow.


## 4.1. Ventilation Method

According to the actual mining situation and considering the geological characteristics, the selection of working face ventilation method varies. At present, U-type ventilation is most commonly used in coal mines. With the increasing maturity of coal mining theory and to meet the demand for safe mining under different geological conditions, various ventilation methods have been developed based on U-type ventilation. In recent years, scholars have made specific studies on the distribution of flow field in the gob under different ventilation methods.


### 4.1.1. U-Type Ventilation

The U-type ventilation system is mainly used in coal mining workings of low gas mines, where there is only one air inlet and one air return lane, and it can be divided into U-back and U-forward according to the direction of workings. U-type backward ventilation is shown in Figure9. In China’s coal mines, the U-back type is commonly used in the working face. The gas in the upper corner of U-ventilated working face of high gas mine has been a difficult problem for coal mine gas management. With the development of coal science and technology, coal mining working face has realized high-production and high-efficiency comprehensive mechanized production, and the unit production of working face has been improved continuously, and the gas gushing out has been increased, so that the mines with little gas in the original working face and low gas mines also have the upper corner gas accumulation over limit [58]. Many scholars have analyzed the wind flow and gas flow law of U-type ventilation working face through on-site observation and research.Figure 9
U-type backward ventilation.In order to obtain the flow field distribution in the mining void area of U-type comprehensive mining working face, Tao [59] established a two-dimensional physical model of the mining void area according to the actual situation of 11124 working face and conducted numerical simulation studies on the constructed model using FLUENT software and adaptive grid encryption technology. The study showed that the air leakage in the gob mainly occurred in the range of 0∼20 m, and some of the air leakage into the gob returned to the working face in this area, while most of the air leakage returned to the working face in the inclination of 140∼160 m; the wind velocity contours in the gob were not symmetrical in the inclination direction; in the direction of the gob, the wind velocity of air leakage showed a gradually decreasing trend.Tian et al. [60] arranged one release point and five sampling points in the working face, used the instantaneous release method to continuously release SF6 gas at the wind source location, conducted gas sampling at the sampling points at regular intervals to analyze the SF6 gas concentration in the gas, estimated the air leakage volume and minimum air leakage velocity in the gob according to the time and concentration of SF6 in the gas sample, and analyzed the measured data to make a general inference about the location of air leakage in the working face. Through the analysis of the measured data, the location of the air leakage channel in the working face can be roughly deduced, as shown in Figure 10. After that, the numerical simulation software FLUENT was used to verify the air leakage in the U-type ventilation gob. The actual measurement results of the corresponding wind speed and air leakage velocity of each monitoring point are compared with the numerical simulation results as shown in Figure 11. The numerical simulation results are consistent with the actual measurement results on-site. The air leakage in the gob is distributed in a U-type, which flows into the gob from the inlet side of the working face and deflects to the return side after a certain depth and flows out; the air leakage in the gob of the U-type ventilation system is mainly in the lower corner area of the working face, which is more prone to spontaneous combustion.Figure 10
Schematic diagram of air leakage channel in gob of 15102 working face.Figure 11
Comparison of field measurement and numerical simulation results. (a) Monitor the corresponding location. (b) Monitor the corresponding location.
(a)
(b)With the increase of coal mining at the working face, the single U-type ventilation system cannot meet the safety production of coal mining enterprises, and in order to keep the gas concentration at the working face from exceeding the limit, new ventilation methods have emerged based on the U-type ventilation system, among which U + I-type and U + L-type are more widely used at the comprehensive working face. In U + I-type ventilation, the gas tail lane is arranged on the inner side of the return air lane, and the gas carried by the leakage is discharged through the gas tail lane; while in U + L-type ventilation, the gas tail lane is arranged on the outer side of the return air lane. Although these two ventilation methods are good for controlling the gas concentration in the corner of the working face, both of them change the air leakage in the gob compared with the U-type ventilation method, causing the flow field in the gob to change.Ma et al. [61] used the UDF interface to establish the calculation models of coal gas seepage in the mined area based on the moving coordinates of the working face and the gas dispersion and transport in the retracement tunnel and quantitatively compared the U-type and U + I-type ventilation methods under different working face advance speed and working face air distribution conditions. The paper concluded that the location of air leakage in U-type ventilation was mainly concentrated in the lower corner of the working face, while in U + I-type ventilation, there were air leakage channels in the whole working face; the size of gas concentration in the return air lane under U-type ventilation is about twice that of U + I-type ventilation, and with the increase of wind speed in the inlet air lane, U + I-type ventilation and U-type ventilation have more gas concentration. With the increase of wind speed in the inlet lane, the gas concentration in the return lane of both U + I-type and U-type ventilation methods shows a decrease in the form of decrease, and with the increase of working face advancement speed, the gas concentration in the return lane of U-type ventilation method shows a decrease in the form of increase, while the gas concentration in the return lane of U + I-type ventilation method shows a decrease in the form of decrease. Compared with the U-type ventilation method, the gas concentration in the upper corner can be reduced by 25% to 50% under the U + I-type ventilation method.Wu et al. [62] conducted numerical simulation of the flow field in the gob under U-type and U + L-type ventilation using FLUENT software and compared and analyzed the distribution of gas concentration in the gob and the distribution of spontaneous combustion hazard zone in the gob under the two ventilation methods, as shown in Figure 12, and the results showed that the gas concentration in the gob of U + L-type ventilation was lower than that of U-type ventilation and the gas explosion range in the gob was lower than that of U-type ventilation. The range of gas explosion moves to the deep part of the gob, but the width of gas explosion area becomes smaller. The gas concentration in the middle of the gob increases faster along the direction of the gob, and the spontaneous combustion hazard area in the gob moves deeper in the gob with U + L-type ventilation than U-type ventilation, and the width of the spontaneous combustion hazard area in the gob also increases significantly.Figure 12
Simulation results of flow field in U-type and U + L-type ventilation gob. (a) Gas concentration distribution in U-type ventilation gob. (b) Oxygen concentration distribution in U-type ventilation gob. (c) Gas concentration distribution in U + L-type ventilation gob. (d) Oxygen concentration distribution in U + L-type ventilation gob.
(a)
(b)
(c)
(d)


### 4.1.2. Y-Type Ventilation

The Y-type ventilation method is a ventilation system that adds fresh air flow on the return side and flows out from the side of the gob after converging with the lacking air from the working face. The use of Y-type ventilation method in the coal face increases the air volume in the return air tunnel, which can effectively solve the problem of gas overrun in the upper corner and is conducive to gas extraction, so it has good promotion value, as shown in Figure13.Figure 13
Y-type backward ventilation.Yu et al. [63] established a two-dimensional solution model for the flow field in the gob under Y-type ventilation, applied the computational fluid dynamics software FLUENT to numerically simulate the distribution of the air leakage field and gas concentration field in the gob, conducted a comprehensive comparison of the air distribution ratios of different main and subinlet lanes, and concluded that the gas management effect was best when the air distribution ratio was 5 : 1.Yang et al. [64] established the physical model of U-type and Y-type ventilation air-gob and used FLUENT software to numerically simulate the air leakage flow field, air leakage volume, and gas distribution in the air-gob under these two ventilation methods, and the simulation results are shown in Figure 14. Most of the air flow in U-type ventilation flows into the working face through the inlet alley and finally discharges through the return alley, and a small part of the air flow leaks into the air-gob from the lower corner of the working face, and in Y-type ventilation, most of the air from the machine tunnel flows through the working face, then converges with the air from the wind tunnel, and enters the air-retention lane, finally discharges through the return air lane, and some of the air leaks into the air-gob from the lower corner of the working face to the middle of the working face, and the closer the air leaks into the air-gob, the more it blows to the bottom of the air-gob, and then it flows from the air-gob to the interface between the air-gob and the air-retention lane. The more the air leaks into the air-gob from the side of the machine tunnel, the more it can be blown to the bottom of the air-gob and then flow into the air-gob from the interface between the air-gob and the alley. After analyzing the air leakage data of the two types of ventilation, it is concluded that the air leakage at the end of the working face from 0 to 30m accounts for 50% of the air leakage at the working face with Y-type ventilation, and the air leakage at the working face to the air-gob from 30 to 90m accounts for 46% of the total air leakage, and the total air leakage is more than that of U-type ventilation, which is not suitable for the coal mining face with coal spontaneous combustion. In addition, under Y-type ventilation, most of the gas in the gob is carried by the air leakage from the working face to the gob to be discharged along the air-retention lane, and very little of the low volume fraction gas (≤0.5%) is gushed to the corner of the working face, so the volume fraction of gas in the gob is lower than that in U-type ventilation, which can avoid the accumulation of high concentration gas in the gob. The use of two-in-one Y-type ventilation can fundamentally solve the problem of gas accumulation in the upper corner and gas overload in the return airway.Figure 14
Simulation results of flow field in U-ventilated and Y-ventilated gob. (a) Flow diagram of the U-type ventilation and recovery workings in the mining void area. (b) Flow diagram of gob of Y-type ventilation and retrieval working face. (c) Gas distribution map of U-type ventilation gob. (d) Y-type ventilation gob gas distribution map.
(a)
(b)
(c)
(d)Zhao et al. [65] took 5301 working face of the Sihe Mine of the Jin Coal Group as the research object and studied the gas gushing law in the gob under partial Y-type ventilation according to the measured parameters of gas gushing, and two gas gushing were considered in the literature, which were coal left in the gob and gas gushing from the adjacent layer on the gob, set at 0.3 m from the bottom plate of the gob and the adjacent layer area of the gob (19.28 to 19.80 m from the bottom plate of the gob). Bottom plate was 19.28∼19.80 m height range, and different locations in the gob defined different gas outflow rate; according to the O-ring theory, we set the gob void rate and permeability.


### 4.1.3. W-Type Ventilation

W-type ventilation is mainly used in long working face with high gas, where both the inlet and return airway are in the coal body. In order to solve the problem of gas overrun in the upper corner in U-type ventilation, two-in-one W-type ventilation is often used with the upper and lower flat as the inlet airway and the middle roadway as the return airway. Because of the maintenance difficulties of W-type forward ventilation tunnel, the air leakage is relatively large, and the gas gushing from the gob is also large, W-type backward ventilation is generally used, as shown in Figure15.Figure 15
W-type backward ventilation.Li and Li [66] established a finite element numerical model of the flow pattern of wind leakage in the bubble fall gob based on the equation of wind leakage and seepage in nonhomogeneous porous media, combined with field examples, and the calculation was based on MATLAB to develop a calculation program to display the results graphically. The sparse and dense transition is carried out by using the frontier generation method for the sectioning of the gob and the accuracy of the boundary sectioning of the wind leakage condition in the quarry as the basis of the regional grid density function; the spatial shape of the area is reflected by the change of the flow field height. Theoretically, the wind pressure distribution contours and flow function lines (distribution solutions) of the W-type ventilation gob are depicted, and it is pointed out that the W-type ventilation gob obviously reduces the wind pressure difference between the two ends of the working face and the effect of reducing the range of wind leakage.Cheng and Xin [67] compared the W-type ventilation method with U-type ventilation, and the simulation results showed that the leakage wind flow fields of the two are obviously different, and under U-type ventilation conditions, most of the leakage wind in the gob penetrates into the deep part of the gob, while W-type ventilation reduces the pressure difference between the two ends of the working face by about 35% compared with that of U-type due to the effect of equal pressure at both ends of the working face, causing the leakage wind in the gob to be closer to the working face under W-type ventilation. At the same time, because the air leakage will carry a lot of gas from the gob into the working face, W-type ventilation reduces the air leakage from the working face to the gob, thus suppressing the gas leakage from the gob to the working face, and fundamentally solving the problem of frequent overlimit gas accumulation in the upper corner of U-type ventilation and the return airway.


### 4.1.4. J-Type Ventilation

This type of ventilation is a new type of ventilation system proposed only in recent years. It has a stronger and more controllable gas discharge capacity than the U-type ventilation commonly used in coal mines and is mainly suitable for high gas discharge workings in coal seams that are not prone to spontaneous combustion, which can well solve the corner gas problem on the workings, as shown in Figure16.Figure 16
J-type backward ventilation.According to the actual situation of 5201 header working face of the Wangzhuang Coal Mine in Lu’an, Wang and Wu [68] compared and analyzed the flow field and gas transport characteristics of J-type ventilation and U-type ventilation by using PHOENICS software. The simulation results of U-type ventilation and J-ventilation extraction zone are shown in Figures 17 and 18, respectively. The flow field and gas transport law of J-type ventilation gob were studied systematically under different conditions of gob size, gas discharge capacity of special alley, and extraction fan installation position. The results show that the J-type ventilation system can eliminate the concentrated air leakage in the upper corner of the gob and effectively solve the problem of gas accumulation in the upper corner of U-type ventilation.Figure 17
Simulation results of flow field andφ (CH4) in U-type ventilation gob. (a) Flow line and pressure distribution. (b) φ (CH4) distribution.
(a)
(b)Figure 18
Simulation results of flow field andφ (CH4) in J-type ventilation gob. (a) Flow line and pressure distribution. (b) φ (CH4) distribution.
(a)
(b)Gao and Sun [69] investigated the air leakage pattern in the gob, the distribution pattern of gas in the gob and the upper corner under different air volume conditions in the J-type ventilation system by means of numerical simulation. With the increase of inlet air volume, the area of low concentration of gas in the gob gradually deepens to the upper right of the gob, and the distribution range increases significantly. The magnitude of gas concentration reduction in the deep gob is larger than that in the shallow gob. The gas concentration in the upper corner also decreases gradually with the increase of the inlet air volume, and the gas concentration in the upper corner is lower than the 1% upper limit standard stipulated in the safety regulations when the air volume increases to 2350 m3/min.


## 4.1.1. U-Type Ventilation

The U-type ventilation system is mainly used in coal mining workings of low gas mines, where there is only one air inlet and one air return lane, and it can be divided into U-back and U-forward according to the direction of workings. U-type backward ventilation is shown in Figure9. In China’s coal mines, the U-back type is commonly used in the working face. The gas in the upper corner of U-ventilated working face of high gas mine has been a difficult problem for coal mine gas management. With the development of coal science and technology, coal mining working face has realized high-production and high-efficiency comprehensive mechanized production, and the unit production of working face has been improved continuously, and the gas gushing out has been increased, so that the mines with little gas in the original working face and low gas mines also have the upper corner gas accumulation over limit [58]. Many scholars have analyzed the wind flow and gas flow law of U-type ventilation working face through on-site observation and research.Figure 9
U-type backward ventilation.In order to obtain the flow field distribution in the mining void area of U-type comprehensive mining working face, Tao [59] established a two-dimensional physical model of the mining void area according to the actual situation of 11124 working face and conducted numerical simulation studies on the constructed model using FLUENT software and adaptive grid encryption technology. The study showed that the air leakage in the gob mainly occurred in the range of 0∼20 m, and some of the air leakage into the gob returned to the working face in this area, while most of the air leakage returned to the working face in the inclination of 140∼160 m; the wind velocity contours in the gob were not symmetrical in the inclination direction; in the direction of the gob, the wind velocity of air leakage showed a gradually decreasing trend.Tian et al. [60] arranged one release point and five sampling points in the working face, used the instantaneous release method to continuously release SF6 gas at the wind source location, conducted gas sampling at the sampling points at regular intervals to analyze the SF6 gas concentration in the gas, estimated the air leakage volume and minimum air leakage velocity in the gob according to the time and concentration of SF6 in the gas sample, and analyzed the measured data to make a general inference about the location of air leakage in the working face. Through the analysis of the measured data, the location of the air leakage channel in the working face can be roughly deduced, as shown in Figure 10. After that, the numerical simulation software FLUENT was used to verify the air leakage in the U-type ventilation gob. The actual measurement results of the corresponding wind speed and air leakage velocity of each monitoring point are compared with the numerical simulation results as shown in Figure 11. The numerical simulation results are consistent with the actual measurement results on-site. The air leakage in the gob is distributed in a U-type, which flows into the gob from the inlet side of the working face and deflects to the return side after a certain depth and flows out; the air leakage in the gob of the U-type ventilation system is mainly in the lower corner area of the working face, which is more prone to spontaneous combustion.Figure 10
Schematic diagram of air leakage channel in gob of 15102 working face.Figure 11
Comparison of field measurement and numerical simulation results. (a) Monitor the corresponding location. (b) Monitor the corresponding location.
(a)
(b)With the increase of coal mining at the working face, the single U-type ventilation system cannot meet the safety production of coal mining enterprises, and in order to keep the gas concentration at the working face from exceeding the limit, new ventilation methods have emerged based on the U-type ventilation system, among which U + I-type and U + L-type are more widely used at the comprehensive working face. In U + I-type ventilation, the gas tail lane is arranged on the inner side of the return air lane, and the gas carried by the leakage is discharged through the gas tail lane; while in U + L-type ventilation, the gas tail lane is arranged on the outer side of the return air lane. Although these two ventilation methods are good for controlling the gas concentration in the corner of the working face, both of them change the air leakage in the gob compared with the U-type ventilation method, causing the flow field in the gob to change.Ma et al. [61] used the UDF interface to establish the calculation models of coal gas seepage in the mined area based on the moving coordinates of the working face and the gas dispersion and transport in the retracement tunnel and quantitatively compared the U-type and U + I-type ventilation methods under different working face advance speed and working face air distribution conditions. The paper concluded that the location of air leakage in U-type ventilation was mainly concentrated in the lower corner of the working face, while in U + I-type ventilation, there were air leakage channels in the whole working face; the size of gas concentration in the return air lane under U-type ventilation is about twice that of U + I-type ventilation, and with the increase of wind speed in the inlet air lane, U + I-type ventilation and U-type ventilation have more gas concentration. With the increase of wind speed in the inlet lane, the gas concentration in the return lane of both U + I-type and U-type ventilation methods shows a decrease in the form of decrease, and with the increase of working face advancement speed, the gas concentration in the return lane of U-type ventilation method shows a decrease in the form of increase, while the gas concentration in the return lane of U + I-type ventilation method shows a decrease in the form of decrease. Compared with the U-type ventilation method, the gas concentration in the upper corner can be reduced by 25% to 50% under the U + I-type ventilation method.Wu et al. [62] conducted numerical simulation of the flow field in the gob under U-type and U + L-type ventilation using FLUENT software and compared and analyzed the distribution of gas concentration in the gob and the distribution of spontaneous combustion hazard zone in the gob under the two ventilation methods, as shown in Figure 12, and the results showed that the gas concentration in the gob of U + L-type ventilation was lower than that of U-type ventilation and the gas explosion range in the gob was lower than that of U-type ventilation. The range of gas explosion moves to the deep part of the gob, but the width of gas explosion area becomes smaller. The gas concentration in the middle of the gob increases faster along the direction of the gob, and the spontaneous combustion hazard area in the gob moves deeper in the gob with U + L-type ventilation than U-type ventilation, and the width of the spontaneous combustion hazard area in the gob also increases significantly.Figure 12
Simulation results of flow field in U-type and U + L-type ventilation gob. (a) Gas concentration distribution in U-type ventilation gob. (b) Oxygen concentration distribution in U-type ventilation gob. (c) Gas concentration distribution in U + L-type ventilation gob. (d) Oxygen concentration distribution in U + L-type ventilation gob.
(a)
(b)
(c)
(d)


## 4.1.2. Y-Type Ventilation

The Y-type ventilation method is a ventilation system that adds fresh air flow on the return side and flows out from the side of the gob after converging with the lacking air from the working face. The use of Y-type ventilation method in the coal face increases the air volume in the return air tunnel, which can effectively solve the problem of gas overrun in the upper corner and is conducive to gas extraction, so it has good promotion value, as shown in Figure13.Figure 13
Y-type backward ventilation.Yu et al. [63] established a two-dimensional solution model for the flow field in the gob under Y-type ventilation, applied the computational fluid dynamics software FLUENT to numerically simulate the distribution of the air leakage field and gas concentration field in the gob, conducted a comprehensive comparison of the air distribution ratios of different main and subinlet lanes, and concluded that the gas management effect was best when the air distribution ratio was 5 : 1.Yang et al. [64] established the physical model of U-type and Y-type ventilation air-gob and used FLUENT software to numerically simulate the air leakage flow field, air leakage volume, and gas distribution in the air-gob under these two ventilation methods, and the simulation results are shown in Figure 14. Most of the air flow in U-type ventilation flows into the working face through the inlet alley and finally discharges through the return alley, and a small part of the air flow leaks into the air-gob from the lower corner of the working face, and in Y-type ventilation, most of the air from the machine tunnel flows through the working face, then converges with the air from the wind tunnel, and enters the air-retention lane, finally discharges through the return air lane, and some of the air leaks into the air-gob from the lower corner of the working face to the middle of the working face, and the closer the air leaks into the air-gob, the more it blows to the bottom of the air-gob, and then it flows from the air-gob to the interface between the air-gob and the air-retention lane. The more the air leaks into the air-gob from the side of the machine tunnel, the more it can be blown to the bottom of the air-gob and then flow into the air-gob from the interface between the air-gob and the alley. After analyzing the air leakage data of the two types of ventilation, it is concluded that the air leakage at the end of the working face from 0 to 30m accounts for 50% of the air leakage at the working face with Y-type ventilation, and the air leakage at the working face to the air-gob from 30 to 90m accounts for 46% of the total air leakage, and the total air leakage is more than that of U-type ventilation, which is not suitable for the coal mining face with coal spontaneous combustion. In addition, under Y-type ventilation, most of the gas in the gob is carried by the air leakage from the working face to the gob to be discharged along the air-retention lane, and very little of the low volume fraction gas (≤0.5%) is gushed to the corner of the working face, so the volume fraction of gas in the gob is lower than that in U-type ventilation, which can avoid the accumulation of high concentration gas in the gob. The use of two-in-one Y-type ventilation can fundamentally solve the problem of gas accumulation in the upper corner and gas overload in the return airway.Figure 14
Simulation results of flow field in U-ventilated and Y-ventilated gob. (a) Flow diagram of the U-type ventilation and recovery workings in the mining void area. (b) Flow diagram of gob of Y-type ventilation and retrieval working face. (c) Gas distribution map of U-type ventilation gob. (d) Y-type ventilation gob gas distribution map.
(a)
(b)
(c)
(d)Zhao et al. [65] took 5301 working face of the Sihe Mine of the Jin Coal Group as the research object and studied the gas gushing law in the gob under partial Y-type ventilation according to the measured parameters of gas gushing, and two gas gushing were considered in the literature, which were coal left in the gob and gas gushing from the adjacent layer on the gob, set at 0.3 m from the bottom plate of the gob and the adjacent layer area of the gob (19.28 to 19.80 m from the bottom plate of the gob). Bottom plate was 19.28∼19.80 m height range, and different locations in the gob defined different gas outflow rate; according to the O-ring theory, we set the gob void rate and permeability.


## 4.1.3. W-Type Ventilation

W-type ventilation is mainly used in long working face with high gas, where both the inlet and return airway are in the coal body. In order to solve the problem of gas overrun in the upper corner in U-type ventilation, two-in-one W-type ventilation is often used with the upper and lower flat as the inlet airway and the middle roadway as the return airway. Because of the maintenance difficulties of W-type forward ventilation tunnel, the air leakage is relatively large, and the gas gushing from the gob is also large, W-type backward ventilation is generally used, as shown in Figure15.Figure 15
W-type backward ventilation.Li and Li [66] established a finite element numerical model of the flow pattern of wind leakage in the bubble fall gob based on the equation of wind leakage and seepage in nonhomogeneous porous media, combined with field examples, and the calculation was based on MATLAB to develop a calculation program to display the results graphically. The sparse and dense transition is carried out by using the frontier generation method for the sectioning of the gob and the accuracy of the boundary sectioning of the wind leakage condition in the quarry as the basis of the regional grid density function; the spatial shape of the area is reflected by the change of the flow field height. Theoretically, the wind pressure distribution contours and flow function lines (distribution solutions) of the W-type ventilation gob are depicted, and it is pointed out that the W-type ventilation gob obviously reduces the wind pressure difference between the two ends of the working face and the effect of reducing the range of wind leakage.Cheng and Xin [67] compared the W-type ventilation method with U-type ventilation, and the simulation results showed that the leakage wind flow fields of the two are obviously different, and under U-type ventilation conditions, most of the leakage wind in the gob penetrates into the deep part of the gob, while W-type ventilation reduces the pressure difference between the two ends of the working face by about 35% compared with that of U-type due to the effect of equal pressure at both ends of the working face, causing the leakage wind in the gob to be closer to the working face under W-type ventilation. At the same time, because the air leakage will carry a lot of gas from the gob into the working face, W-type ventilation reduces the air leakage from the working face to the gob, thus suppressing the gas leakage from the gob to the working face, and fundamentally solving the problem of frequent overlimit gas accumulation in the upper corner of U-type ventilation and the return airway.


## 4.1.4. J-Type Ventilation

This type of ventilation is a new type of ventilation system proposed only in recent years. It has a stronger and more controllable gas discharge capacity than the U-type ventilation commonly used in coal mines and is mainly suitable for high gas discharge workings in coal seams that are not prone to spontaneous combustion, which can well solve the corner gas problem on the workings, as shown in Figure16.Figure 16
J-type backward ventilation.According to the actual situation of 5201 header working face of the Wangzhuang Coal Mine in Lu’an, Wang and Wu [68] compared and analyzed the flow field and gas transport characteristics of J-type ventilation and U-type ventilation by using PHOENICS software. The simulation results of U-type ventilation and J-ventilation extraction zone are shown in Figures 17 and 18, respectively. The flow field and gas transport law of J-type ventilation gob were studied systematically under different conditions of gob size, gas discharge capacity of special alley, and extraction fan installation position. The results show that the J-type ventilation system can eliminate the concentrated air leakage in the upper corner of the gob and effectively solve the problem of gas accumulation in the upper corner of U-type ventilation.Figure 17
Simulation results of flow field andφ (CH4) in U-type ventilation gob. (a) Flow line and pressure distribution. (b) φ (CH4) distribution.
(a)
(b)Figure 18
Simulation results of flow field andφ (CH4) in J-type ventilation gob. (a) Flow line and pressure distribution. (b) φ (CH4) distribution.
(a)
(b)Gao and Sun [69] investigated the air leakage pattern in the gob, the distribution pattern of gas in the gob and the upper corner under different air volume conditions in the J-type ventilation system by means of numerical simulation. With the increase of inlet air volume, the area of low concentration of gas in the gob gradually deepens to the upper right of the gob, and the distribution range increases significantly. The magnitude of gas concentration reduction in the deep gob is larger than that in the shallow gob. The gas concentration in the upper corner also decreases gradually with the increase of the inlet air volume, and the gas concentration in the upper corner is lower than the 1% upper limit standard stipulated in the safety regulations when the air volume increases to 2350 m3/min.


## 4.2. Coal Seam Dip Angle

The geological conditions of China’s coal mines are complex and variable, and the coal seam inclination is not the same even for the same well field. The coal seam inclination is the main factor affecting the changing damage form of the overlying rock layer in the mining field. With the development of coal seam inclination from gentle to sharp inclination, the deformation and damage characteristics of the roof after coal seam mining will change, and the porous media parameters in the gob are closely related to it, which in turn affects the leakage flow field and the distribution of gas and pressure in the gob.Lu and Zhang [70] used UDEC numerical simulation to study the law of overlying rock strata spanning fall in the mining void zone of inclined extrathick coal seam and obtained the distribution of porosity in the mining void zone by the amount of overlying rock subsidence in the mining void zone. The results show that the porosity in the fallout zone is maximum along the strike at the edge of the gob and minimum at the deep part; it is maximum at the upper corner and lower corner along the tendency and minimum at the middle and upper part; in the vertical direction, the porosity in the gob decreases with increasing height from the bottom plate, and the difference between the edge and the deep part also decreases gradually.Yang and Zhang [71] studied the effect of coal seam inclination on the flow field and gas distribution in the mining void area of U + L-type ventilation face in a more systematic way, which divided the mining void area into bubble zone and fracture zone, where the bubble zone was subdivided into free accumulation zone, load-influenced zone, and compaction stability zone, set the porosity of the bubble zone to be uniformly distributed in segments, and investigated the effect of the coal seam inclination of 5.75°, 35.75°, with the increase of the coal seam inclination, the pressure difference between the two ends of the working face and the pressure difference between the working face and the tailing lane increase gradually, and the increase increases first and then decreases; the increase of the pressure difference between the two ends of the working face is greater than the pressure difference between the working face and the tailing lane; under the condition of keeping the air volume in the inlet lane, the return lane and the tailing lane unchanged, and the air volume in the working face decreases with the increase of the coal seam inclination. The air volume decreases, but the proportion of decrease is small and the trend is decreasing, and the distribution law of air leakage shows that the air leakage at the upper and lower ends is large, and the air leakage at the middle of the working face is small. As the coal seam inclination increases, the gas in the gob is affected by the uplift effect and the leakage flow, and the gas concentration gradually decreases. The gas concentration gradually decreases due to the uplift effect and leakage flow.


## 5. Study on the Prevention and Control of Fire and Explosion Risks in Gobs

### 5.1. Study on the Prevention and Control of Fire in the Gob

The development of mechanization and automation of coal mining technology has greatly improved the production efficiency of mines, but accordingly, it has also caused an increase in coal left in the gob and increased air leakage, resulting in spontaneous combustion accidents in the gob. According to statistics, more than 50% of natural fires occur in mines in the gob, and it is important to study the flow field in the gob to prevent coal spontaneous combustion from causing mine fires [72].The study of the flow field of the gob on the prevention and control of fire risk is mainly from the perspective of the flow field of the gob to analyze and predict the natural fire in the gob, and also to study whether the injection of inert gas, grouting, fracture filling, or pressure balancing measures in the gob are beneficial to the prevention and control of spontaneous combustion fires, and the research means are mainly numerical simulation software, fire beam monitoring system arranged on-site, etc. Combined with a series of existing theoretical theories of coal spontaneous combustion, the analysis is carried out.In order to analyze and predict the natural fire in the gob, the method of dividing the spontaneous combustion “three zones” (i.e., dispersal zone, spontaneous combustion zone, and choking zone) in the gob has been adopted, and many scholars and research institutions at home and abroad have determined the indexes of dividing the three zones in the gob after doing a lot of in-depth research on this, including the temperature rise rate of the coal left in the gob, the wind speed of the coal leakage in the gob, and the wind speed of the coal leakage in the gob. The indicators of air leakage are often the air velocity and the oxygen concentration in mine gobs [73–77]. In addition, a large number of researchers have studied the prediction and analysis of the flow field in the gob for natural fires and other conditions. For example, Zhou [72] obtained the change pattern of the temperature field in the gob during the advancement of the working face through a four-dimensional dynamic simulation study. The migration rate of the high-temperature zone is affected by factors such as porosity. At the same time, there is a dynamic change between the high-temperature zone and the high oxygen zone in the gob. Zhang et al. [78] used the discrete element modeling software Particle Flow Code (PFC) to simulate the pore evolution and porosity distribution of the overlying rock in the mined-out area and established a dynamic porosity model for the mined-out area. The flow field in the area was simulated, and the regular gas flow distribution in the mined-out area was obtained. Compared with the actual mine measurement results, it was verified that the flow field distribution in the mined-out area considering the porosity dynamics can accurately and reliably identify residual coal spontaneous combustion prone area. Zheng et al. [79] studied the competitive adsorption characteristics of multicomponent gas and its influence on coal oxidation and found that the gas adsorption selectivity order on the coal surface is CO2 > CH4 > O2 > N2. Under ambient conditions, the strong adsorption capacity of CH4 and CO2 on the coal surface hinders the spontaneous combustion of coal due to the limited adsorption sites on the coal surface. Especially, CH4 adsorbed by residual coal will affect the chemical reaction between coal and O2, while the reaction between CH4 and O2 on the surface of high-temperature coal will accelerate the spontaneous combustion of coal and increase the risk of coal spontaneous combustion. Zuo et al. [80] used a self-built experimental platform to determine the effect of different air leakage conditions on the spontaneous combustion of residual coal and obtained an exponential relationship between air supply to the working face and the spontaneous combustion of broken coal in the gob. Wang et al. [81] used the particle flow numerical simulation software PFC3D to simulate the collapse of the overlying rock layer in the void area, extracted the quantitative porosity data of the void area, imported it into FLUENT software, simulated the leakage flow field in the void area, and obtained the main leakage area of the working face. Lu and Qin [82] used numerical simulation to study the plastic deformation process of the coal column and deduced the distribution law of porosity of the coal column. On this basis, the mathematical model of oxygen consumption and seepage in the coal pillar was established using ANSYS-FLUENT software, the distribution of oxygen concentration in the coal pillar was obtained, and the potential spontaneous combustion area in the coal pillar in the gob was determined according to the oxygen concentration index.Not only does the research on the prevention of natural fires stop at prediction and analysis, but also researchers at home and abroad have studied the effect of measures to change the flow field in the gob (such as inert gas injection, slurry injection, fracture filling, or pressure equalization) on the prevention of spontaneous combustion fires. Wang and Tang [83] used the beam tube monitoring system to study that a comprehensive method combining pressure equalization, grouting, and crack filling has a significant effect on preventing spontaneous combustion of coal and complex air leakage in longwall gob. Qin et al. [84] established a theoretical and geometric model to quickly determine the high-temperature zone of fire by simulating the oxygen concentration and temperature distribution in the gob and evaluated the effect of liquid nitrogen infusion by comparing field measurements and numerical simulations. Zhang et al. [85] established a three-dimensional transient nonequilibrium thermal CFD model based on the actual situation of the long wall gob and studied the thermal evolution and active inserting scheme of the long wall gob. Wang et al. [86] used FLUENT software to numerically simulate the spontaneous combustion of residual coal in the air leakage area of the gob and studied that the slurry injection method can prevent the spontaneous combustion of residual coal by reducing the area of porosity in the gob and the area of oxidation zone behind the working face, which is a cost-effective method.


### 5.2. Study on the Prevention and Control of Explosion Risk in the Gob

With the advancement of the working face, the increase of the gob area, and the influence of the change of airflow in the roadway on the gas composition of the gob in the production process, it indirectly affects the change of the gas explosion risk area in the gob. Studying the law of gas migration in gob plays an important role in preventing gas explosion in the gob. The spontaneous combustion of coal remains in the gob fire, rock friction initiation of electricity, and other sources of spontaneous combustion as a source of ignition is the cause of ignition detonation gob gas-air mixture [87–95].Yang et al. [93, 96, 97] and Li [98] combined the theoretical analysis of coal spontaneous combustion in the gob, mine pressure, gas flow field in the gob, and gas explosion to analyze the causes of the coupled disaster of spontaneous combustion and explosion in the gob and analyzed the possible range of gas occurrence in the gob with numerical simulation, as shown in Figure 19.Figure 19
Distribution of gas explosion area in the gob under different advancing distances. (a)L = 80 m. (b) L = 130 m. (c) L = 180 m.
(a)
(b)
(c)Li [99] realized the experimental environment of coal spontaneous combustion buoyancy effect in the gob and established the corresponding numerical model, systematically studied the gas movement accumulation law of coal spontaneous combustion environment in the gob, and revealed the disaster formation process of gas explosion induced by coal spontaneous combustion in the gob. Qin et al. [100] through theoretical and experimental analysis pointed out that the main product of spontaneous combustion in the gob CO greatly increases the explosive concentration limit range of gas-air mixture, and under the action of fire and wind pressure, the full mixture of CH4, CO, and fresh air and convective exchange of heat are constantly formed between the firing and nonfiring areas, which easily lead to gas explosion accidents.The closed gas in the extraction area experiences a critical period when the gas concentration fluctuates between the upper explosive limit and the lower explosive limit [101, 102]. Since gas influx from the extraction zone cannot be detected, some scholars have analyzed gas transport patterns by scaling down experimental devices or numerical simulations [103–107]. These methods mainly focus on the changes of oxygen or methane concentration distribution in the extraction zone.By making a good estimation of the porosity and permeability of the mined area, the airflow characteristics of the mined area can be accurately assessed, and the distribution of the coal mine gas flow field can be predicted more accurately, which can provide a theoretical basis for the prevention of heavy gas explosions in the mined area [82, 104, 108]. Bell [109] considered the seepage pattern of the mined area as anisotropic three-dimensional according to the porous media theory. Gao et al. [110, 111] established an ideal physical model with spatial characteristics of the gob for the spatial characteristics of the gob and studied the effect of different rock diameters and arrangements on the gas explosion in the gob. The effect of different rock diameters and arrangements on the gas explosion characteristics of the extraction zone was studied.In Vlasin et al.’s study [112], computational simulations of the explosion of a stoichiometric air-methane mixture in an enclosed space containing obstacles were implemented by changing parameters in the ANSYS-FLUENT application model to customize these data results for a specific gas explosion domain in the form of spatial maps of flame front development and graphs representing the time evolution of characteristic parameters. Li et al. [104] investigated the effect of airflow and gas outflow on the three-dimensional distribution of oxygen and gas concentrations in longwall gobs and mapped the methane-oxygen coupled explosion hazard zone based on modified Coward’s triangle and linear coupling zone formulations. Che [1] determined the concentration distribution, pressure distribution, and temperature distribution of gas in the extraction zone through large-size experiments on gas seepage and heat transfer in the extraction zone. The spatial distribution equations of porosity and permeability of the porous media in the gob were established, and the physical model of the gob and the CFD simulation mathematical model of the gas flow in the gob were improved, and the multifield coupled change law of the three-dimensional seepage field, concentration field, and temperature field of the gas components in the gob before and after the spontaneous combustion of coal relics was derived by the solution. Once the air-gob is closed, the gas components in the closed mine area will change with time. Cheng et al. [113] developed a personal computer-based software package CCMER (Comprehensive Consulting Model for Mine Air Explosion Risk), through which CCEMR can well understand and control the flow field and its change trend in the underground air-gob, and also use the mathematical model of explosion risk assessment to estimate the change from “nonexplosive” state to “explosive” state or from “explosive” state to “nonexplosive” state change time span of the explosion trend analysis model. He [114] and Yi [115] and others analyzed the combustible gas explosion characteristics in the gob from different perspectives. He established a mathematical model of gas explosion in the gob through theoretical analysis and reasonable simplification, simulated the explosion temperature and explosion pressure distribution at different gas concentrations, and derived the effect of gas concentration on multicomponent gas explosion, in which the gas concentration of 10% of the explosion pressure and explosion temperature at 10% gas concentration. Yi used a comprehensive research method combining theoretical analysis and experimental simulation to analyze the combustible gas release law at different temperature bands; finally, according to the explosion ignition source, combined with the gas explosion rheology-sudden change theory, the combustible gas explosion characteristics and influencing factors were qualitatively studied and quantitatively determined, and the control countermeasures of gas explosion accident were proposed.


## 5.1. Study on the Prevention and Control of Fire in the Gob

The development of mechanization and automation of coal mining technology has greatly improved the production efficiency of mines, but accordingly, it has also caused an increase in coal left in the gob and increased air leakage, resulting in spontaneous combustion accidents in the gob. According to statistics, more than 50% of natural fires occur in mines in the gob, and it is important to study the flow field in the gob to prevent coal spontaneous combustion from causing mine fires [72].The study of the flow field of the gob on the prevention and control of fire risk is mainly from the perspective of the flow field of the gob to analyze and predict the natural fire in the gob, and also to study whether the injection of inert gas, grouting, fracture filling, or pressure balancing measures in the gob are beneficial to the prevention and control of spontaneous combustion fires, and the research means are mainly numerical simulation software, fire beam monitoring system arranged on-site, etc. Combined with a series of existing theoretical theories of coal spontaneous combustion, the analysis is carried out.In order to analyze and predict the natural fire in the gob, the method of dividing the spontaneous combustion “three zones” (i.e., dispersal zone, spontaneous combustion zone, and choking zone) in the gob has been adopted, and many scholars and research institutions at home and abroad have determined the indexes of dividing the three zones in the gob after doing a lot of in-depth research on this, including the temperature rise rate of the coal left in the gob, the wind speed of the coal leakage in the gob, and the wind speed of the coal leakage in the gob. The indicators of air leakage are often the air velocity and the oxygen concentration in mine gobs [73–77]. In addition, a large number of researchers have studied the prediction and analysis of the flow field in the gob for natural fires and other conditions. For example, Zhou [72] obtained the change pattern of the temperature field in the gob during the advancement of the working face through a four-dimensional dynamic simulation study. The migration rate of the high-temperature zone is affected by factors such as porosity. At the same time, there is a dynamic change between the high-temperature zone and the high oxygen zone in the gob. Zhang et al. [78] used the discrete element modeling software Particle Flow Code (PFC) to simulate the pore evolution and porosity distribution of the overlying rock in the mined-out area and established a dynamic porosity model for the mined-out area. The flow field in the area was simulated, and the regular gas flow distribution in the mined-out area was obtained. Compared with the actual mine measurement results, it was verified that the flow field distribution in the mined-out area considering the porosity dynamics can accurately and reliably identify residual coal spontaneous combustion prone area. Zheng et al. [79] studied the competitive adsorption characteristics of multicomponent gas and its influence on coal oxidation and found that the gas adsorption selectivity order on the coal surface is CO2 > CH4 > O2 > N2. Under ambient conditions, the strong adsorption capacity of CH4 and CO2 on the coal surface hinders the spontaneous combustion of coal due to the limited adsorption sites on the coal surface. Especially, CH4 adsorbed by residual coal will affect the chemical reaction between coal and O2, while the reaction between CH4 and O2 on the surface of high-temperature coal will accelerate the spontaneous combustion of coal and increase the risk of coal spontaneous combustion. Zuo et al. [80] used a self-built experimental platform to determine the effect of different air leakage conditions on the spontaneous combustion of residual coal and obtained an exponential relationship between air supply to the working face and the spontaneous combustion of broken coal in the gob. Wang et al. [81] used the particle flow numerical simulation software PFC3D to simulate the collapse of the overlying rock layer in the void area, extracted the quantitative porosity data of the void area, imported it into FLUENT software, simulated the leakage flow field in the void area, and obtained the main leakage area of the working face. Lu and Qin [82] used numerical simulation to study the plastic deformation process of the coal column and deduced the distribution law of porosity of the coal column. On this basis, the mathematical model of oxygen consumption and seepage in the coal pillar was established using ANSYS-FLUENT software, the distribution of oxygen concentration in the coal pillar was obtained, and the potential spontaneous combustion area in the coal pillar in the gob was determined according to the oxygen concentration index.Not only does the research on the prevention of natural fires stop at prediction and analysis, but also researchers at home and abroad have studied the effect of measures to change the flow field in the gob (such as inert gas injection, slurry injection, fracture filling, or pressure equalization) on the prevention of spontaneous combustion fires. Wang and Tang [83] used the beam tube monitoring system to study that a comprehensive method combining pressure equalization, grouting, and crack filling has a significant effect on preventing spontaneous combustion of coal and complex air leakage in longwall gob. Qin et al. [84] established a theoretical and geometric model to quickly determine the high-temperature zone of fire by simulating the oxygen concentration and temperature distribution in the gob and evaluated the effect of liquid nitrogen infusion by comparing field measurements and numerical simulations. Zhang et al. [85] established a three-dimensional transient nonequilibrium thermal CFD model based on the actual situation of the long wall gob and studied the thermal evolution and active inserting scheme of the long wall gob. Wang et al. [86] used FLUENT software to numerically simulate the spontaneous combustion of residual coal in the air leakage area of the gob and studied that the slurry injection method can prevent the spontaneous combustion of residual coal by reducing the area of porosity in the gob and the area of oxidation zone behind the working face, which is a cost-effective method.


## 5.2. Study on the Prevention and Control of Explosion Risk in the Gob

With the advancement of the working face, the increase of the gob area, and the influence of the change of airflow in the roadway on the gas composition of the gob in the production process, it indirectly affects the change of the gas explosion risk area in the gob. Studying the law of gas migration in gob plays an important role in preventing gas explosion in the gob. The spontaneous combustion of coal remains in the gob fire, rock friction initiation of electricity, and other sources of spontaneous combustion as a source of ignition is the cause of ignition detonation gob gas-air mixture [87–95].Yang et al. [93, 96, 97] and Li [98] combined the theoretical analysis of coal spontaneous combustion in the gob, mine pressure, gas flow field in the gob, and gas explosion to analyze the causes of the coupled disaster of spontaneous combustion and explosion in the gob and analyzed the possible range of gas occurrence in the gob with numerical simulation, as shown in Figure 19.Figure 19
Distribution of gas explosion area in the gob under different advancing distances. (a)L = 80 m. (b) L = 130 m. (c) L = 180 m.
(a)
(b)
(c)Li [99] realized the experimental environment of coal spontaneous combustion buoyancy effect in the gob and established the corresponding numerical model, systematically studied the gas movement accumulation law of coal spontaneous combustion environment in the gob, and revealed the disaster formation process of gas explosion induced by coal spontaneous combustion in the gob. Qin et al. [100] through theoretical and experimental analysis pointed out that the main product of spontaneous combustion in the gob CO greatly increases the explosive concentration limit range of gas-air mixture, and under the action of fire and wind pressure, the full mixture of CH4, CO, and fresh air and convective exchange of heat are constantly formed between the firing and nonfiring areas, which easily lead to gas explosion accidents.The closed gas in the extraction area experiences a critical period when the gas concentration fluctuates between the upper explosive limit and the lower explosive limit [101, 102]. Since gas influx from the extraction zone cannot be detected, some scholars have analyzed gas transport patterns by scaling down experimental devices or numerical simulations [103–107]. These methods mainly focus on the changes of oxygen or methane concentration distribution in the extraction zone.By making a good estimation of the porosity and permeability of the mined area, the airflow characteristics of the mined area can be accurately assessed, and the distribution of the coal mine gas flow field can be predicted more accurately, which can provide a theoretical basis for the prevention of heavy gas explosions in the mined area [82, 104, 108]. Bell [109] considered the seepage pattern of the mined area as anisotropic three-dimensional according to the porous media theory. Gao et al. [110, 111] established an ideal physical model with spatial characteristics of the gob for the spatial characteristics of the gob and studied the effect of different rock diameters and arrangements on the gas explosion in the gob. The effect of different rock diameters and arrangements on the gas explosion characteristics of the extraction zone was studied.In Vlasin et al.’s study [112], computational simulations of the explosion of a stoichiometric air-methane mixture in an enclosed space containing obstacles were implemented by changing parameters in the ANSYS-FLUENT application model to customize these data results for a specific gas explosion domain in the form of spatial maps of flame front development and graphs representing the time evolution of characteristic parameters. Li et al. [104] investigated the effect of airflow and gas outflow on the three-dimensional distribution of oxygen and gas concentrations in longwall gobs and mapped the methane-oxygen coupled explosion hazard zone based on modified Coward’s triangle and linear coupling zone formulations. Che [1] determined the concentration distribution, pressure distribution, and temperature distribution of gas in the extraction zone through large-size experiments on gas seepage and heat transfer in the extraction zone. The spatial distribution equations of porosity and permeability of the porous media in the gob were established, and the physical model of the gob and the CFD simulation mathematical model of the gas flow in the gob were improved, and the multifield coupled change law of the three-dimensional seepage field, concentration field, and temperature field of the gas components in the gob before and after the spontaneous combustion of coal relics was derived by the solution. Once the air-gob is closed, the gas components in the closed mine area will change with time. Cheng et al. [113] developed a personal computer-based software package CCMER (Comprehensive Consulting Model for Mine Air Explosion Risk), through which CCEMR can well understand and control the flow field and its change trend in the underground air-gob, and also use the mathematical model of explosion risk assessment to estimate the change from “nonexplosive” state to “explosive” state or from “explosive” state to “nonexplosive” state change time span of the explosion trend analysis model. He [114] and Yi [115] and others analyzed the combustible gas explosion characteristics in the gob from different perspectives. He established a mathematical model of gas explosion in the gob through theoretical analysis and reasonable simplification, simulated the explosion temperature and explosion pressure distribution at different gas concentrations, and derived the effect of gas concentration on multicomponent gas explosion, in which the gas concentration of 10% of the explosion pressure and explosion temperature at 10% gas concentration. Yi used a comprehensive research method combining theoretical analysis and experimental simulation to analyze the combustible gas release law at different temperature bands; finally, according to the explosion ignition source, combined with the gas explosion rheology-sudden change theory, the combustible gas explosion characteristics and influencing factors were qualitatively studied and quantitatively determined, and the control countermeasures of gas explosion accident were proposed.


## 6. Conclusion and Prospect

After nearly 40 years of development, Chinese scholars have been made with breakthrough progress in the field of flow field in the gob. The understanding of the geomechanical characteristics of the overlying rock layer, the principle of rock destruction, the law of fracture development, and the distribution of porosity and permeability in the gob has been deepened; the research methods have realized the leap from graphical method, filtered flow field, and similar simulation to CFD numerical simulation, and the solution means has transformed to computer independent operation from solely on manual metal relying on the advance of modern computer technology, which not only greatly shortens the calculation time but also improves the accuracy of the resolving results; numerical model equations are also gradually enriched with scholars’ exploration, the emergence of linear gas seepage theory, nonlinear gas seepage theory, and ground effect flow-solid coupling gas flow theory which make the model description more closely match the real situation of gas flow in the gob; in terms of research objects, the study has extended from simple U-type to various ventilation systems, and the influence of different ventilation parameters, which are wind speed, negative pressure of extraction, etc., and geological conditions of coal mines, which is coal seam inclination angle on the flow field in the gob, also be considered, and the scope of study is gradually comprehensive. At the same time, the research and development of the flow field in the gob has a pivotal role in the prevention and control of fire and explosion in the gob. China has formed a flow field research and control system with Chinese characteristics, which has significantly improved the effect and safety of gas concentration control at the working face and laid a solid foundation for the safe and efficient production of coal mines.Nevertheless, there are still many problems in China’s coal mining flow field research, which need to continue to be explored in depth:(1)
The geological conditions of the overlying rock layer in the coal gob are characterized by complexity and variability. The interior of the gob is a nonuniform porous medium, and the compaction degree of gangue and coal remains is very large in various places, and the porosity distribution is not uniform. The lithology of the overlying rock layer in the extraction area is an important factor affecting the distribution of pore size inside the extraction area. At present, there is a lack of in-depth research on the mechanical parameters of the overlying rock and the influence of different sizes and directions of the ground stress on the destructive movement of the rock layer, etc. The traditional means of measuring the ground stress destroys the original rock stress, and the measurement results have large errors, so it is necessary to develop rapid and accurate measurement technology for the ground stress. The mining stress has the characteristics of constant movement and change in space and time, and there is a lack of technology and instruments that can monitor the mining stress in a long-term, stably and effectively, and the spatial and temporal evolution of the three-dimensional mining stress has not been completely clarified; the overlying rock structure formed after the mining of working faces with different coal mining methods and parameters, and the stability study under coal mining is not sufficient; it needs to be studied from stress testing, monitoring, and simulation, and it is necessary to systematically carry out research on mining dynamics from stress testing, monitoring, simulation, theoretical modeling, and other aspects to truly reflect the process of rock-breakage movement.(2)
The study of the flow field in the gob is a flow field problem and must be based on fluid theory. The wind flow movement in the extraction area, gas transportation, etc. is the process of evolution with the participation of fluid movement, which cannot be simply simplified as the addition of gas components into and extraction out, but it should focus on its internal connection. The gas outflow and natural fire in the quarry is in the same space, and the two affect each other; however, most scholars will study the two separately, or only the gas outflow considered as a qualitative factor affecting the natural fire. Obviously, Due to the deviation between the actual situations and simulation ones, the field flow study should be integrated to consider and further clarify the internal mechanism and coupling relationship between influence factors.(3)
The numerical simulation method is the mainstream method to study the flow field, the actual gob can provide very limited field test conditions, and with the help of numerical simulation, we can understand the inner mechanism of gas gushing out of the gob; however, the model in the process of application requires a large number of field parameters, and some of which cannot be measured directly and must be determined by fitting with the field observation results to have practical significance, and how to use experimental means and simple measurement work to obtain more accurate and complete parameters of the initial value is also the focus of future research.



---

*Source: 3730606-2021-08-17.xml*