# Numerical Simulation Study on Stability Control Technology of Large-Area Wall Caving Area (LAWCA) in Large-Inclined Face

Authors: Xun Liu; Shihao Tu; Hongsheng Tu; Long Tang; Jieyang Ma; Yan Li; Wenlong Li; Kaijun Miao; Hao Tian
Journal: Geofluids
 (2023)

---

## Abstract

The caving characteristics of the goaf significantly impact the mining stability of the working face. The coal-forming environment of the steeply dipping coal seam is complex. After mining, the caving rock mass rolls to the lower part of the working face due to gravity, and the caving compaction characteristics of the working face are different from those of the flat and gently inclined working face, taking Huainan Xinji No. 2 Mine 211112 high-angle composite roof working face large-area gang roof fall area as the engineering background. To safely pass through a large roof fall area, a FLAC3D numerical simulation method considering the goaf’s compaction effect and the excavation’s additional damage is proposed. The Weibull distribution function is introduced to modify the elastic modulus parameters of the rock mass after excavation, which more accurately reflects the damage process of the rock mass after hole. Based on the goaf compaction theory, the elastic model is developed for the second time, and the accurate simulation of the caving rock mass is realized. On this basis, the proposed numerical simulation method is used to simulate the large-area roof caving of the large-angle working face, a scheme to control the stability of the working face by using the paste filling technology is proposed, and the material parameters of the filling body are determined. The research results show that when the water-cement ratio of the backfill is 3 : 1, the simulation and field trial results are good, and the safe and efficient mining of the large-angle composite roof working face is realized, which verifies the feasibility and correctness of the coupled simulation method. At the same time, the research results can provide a reference for controlling the roof fall area of the working face with a large inclination angle.


---

## Body

## 1. Introduction

In underground coal resource mining, the compaction characteristics of the roof rock mass in the goaf after the collapse affect the roof settlement and the conduction of the cracks in the overlying rock strata [1]. Therefore, the study of the compaction characteristics of the goaf is of great significance to the complex response of the underground rock formation caused by mining [2]. The coal-forming environment of the high-dip coal seam is complex. After mining, the caving rock mass rolls to the lower part of the working face due to gravity, forming nonuniform filling in the inclined direction. The working face’s caving and compaction characteristics differ from those of the flat and gently inclined working faces [3, 4]. The rib spalling of coal wall and roof caving in front of working face in large-dip-angle working face are urgent problems to be solved in the safe and efficient mining of large dip angle working face. Currently, the primary methods are to control the support resistance, strengthen the coal wall, and reduce the pseudooblique angle [5, 6]. These methods focus on prevention and control in advance [7]. However, due to the complexity of the occurrence conditions of high-dip coal seams, large-area rib caving (LARC) accidents still occur. After the LARC occurs, grouting and filling become the key to controlling the stability of the working face—an effective method [8]. Because the roofing area of the gang is filled and compacted, part of the stress of the overlying rock layer is transferred to the backfill body, and the support pressure in front of the working face is reduced. During the mining process of the working face, the roof damage is limited due to the backfilling false roof, reducing the degree of stress concentration [9]. During the filling process, it is necessary to consider the characteristics of the filling body itself and the influence of the caving and compaction behind the working face.Due to the complexity of rock strata, the calculation or analysis of surrounding rock excavation is too complex to be analyzed by the analytical method. Numerical simulation overcomes the limitation of laboratory size and can well analyze the stability of rock mass under complex geological conditions [10–12]. Compared with experiments, numerical simulation can quantitatively control the properties and quantities of particles or elements, and numerical simulation is widely used to analyze large-scale engineering problems [13, 14]. However, it is still tricky to correctly reproduce the spatial mechanical response of the stope in the numerical simulation of finite element and finite difference methods [15, 16]. On the one hand, it is not easy to simulate the anisotropy of rock mass structure. On the other hand, due to the influence of excavation disturbance after coal seam excavation, the original defects of surrounding rock deteriorate under external load, and the initial damage is usually considered in numerical simulation.In the numerical simulation calculation, the usual treatment method is to assign the goaf as an empty model, and the progressive compaction characteristics of the caving zone are not considered [17]. Another way is to assign weak lithology to the caving rock mass in the goaf, but these do not conform to the mechanical properties of the caving zone compaction process. Some scholars choose Double-Yield for simulation according to the characteristics of the caving zone. However, due to the lack of basis for selecting mechanical parameters, it is not easy to check them [18]. In recent years, some scholars have carried out an equivalent simulation of the goaf according to the mechanical properties of the caving rock mass in the goaf, such as Bai et al. [19]. Yadav et al. [20] analyzed the sensitivity of Double-Yield, obtained the influence of Double-Yield parameters on goaf compaction, and carried out an example analysis. According to Jiang et al. [21], based on the goaf compaction theory, the coupling analysis method of mining stress and goaf compaction bearing was proposed. However, these simulation methods did not consider the influence of roof rock damage. Zhu et al. and Zhang et al. [15, 22] proposed 3D simulation methods for broken rock masses in goaves, but they mainly focused on the laboratory scale.In this study, based on the finite difference simulation software FLAC3D, a new numerical simulation method considering the additional damage to excavation and the compaction effect of the goaf is proposed. The LAWCA is the engineering background. The numerical simulation method proposed is used to analyze the control effect of the backfill on the surrounding rock of the working face, and the parameters of the backfill body and the backfill technology are determined, which have been tested by on-site practice.


## 2. Project Overview

Xinji No. 2 Mine is located in Huainan City, Anhui Province, China. The 211112 working face is located in the west wing of the second-level central mining area of Xinji No. 2 Mine. It starts from the 2111 mining area in the east, goes up to the central parking lot and cross-measure drift, and reaches the west boundary to the 1211 mining area to protect the coal pillars. The southern boundary is 7.8-27.6 m away from the 111111 goaf, the northern boundary is near the contour line of the 11-2# coal floor, and the upper distance is 56.1-95.3 m from the 13-1# coal seam, with an average of 72.5 m, and the downward distance is 11-1# coal 10.6 ~ 43.3 m, average 28.0 m, buried depth 586.3 ~ 706.3 m, average inclination angle 35°, and working face length 140 m. The 211112 working face has a composite roof with a thickness of 5.6 m, the primary lithology is sandy mudstone and mudstone, the compressive strength is generally about 35 MPa, and its overlying is siltstone and quartz sandstone; the lithology of the floor is mainly mudstone and silty sandstone and sandy mudstone. The 11-2 coal is mainly dark, followed by bright coal, which belongs to semidark briquette, with an average thickness of 5 m and a compressive strength of 0.88~8.15 MPa, as shown in Figure1.Figure 1
Working face layout.Due to the large inclination angle of the working face, the high mining height, the easy collapse of the composite roof layer, and the loose coal body, at a distance of 20 m from the mining stop line, the number of loose rock layers of the composite roof on the top roof of the working face increases, and the thickness of the single layer becomes thinner. During the usual mining process of the working face, it appears at 70# ~90# (within the range of about 30 m, the middle and upper part of the working face, the depth of the rib spalling is about 8 m in front of the working face, and the height of the roof is about 10 m above the top plate of the working face). The LAWCA will seriously affect the regular mining of the working face.The falling coal and the falling roof gangue slide down along the goaf under the action of the inclination of the coal seam, and the middle and lower gobs are filled with self-slide, and a large area of the roof is formed in front of and behind the middle and upper supports. The roof of the empty roof area is in an unsupported state, and there is a danger of sinking and instability. Under the supporting pressure of the working face, the coal wall in front of the working face and the upper rock formation in the empty roof area tend to further instability and roof collapse. The roof of the mined-out area includes wooden stack support and anchor cable support methods. The wood stacking support consumes a large amount of pit wood, the labour intensity of workers is significant, the additional transportation workload is large, and the construction speed is slow. Instability failure may affect the stability of the coal rock mass without lamellae below. In addition, according to the geological data of the working face, there is thicker silty sandstone above the empty roof area, and it is challenging to use bolts and cables to support it. Therefore, the comprehensive analysis decided to carry out grouting and filling in the roof fall area so that the following products can be carried out safely.


## 3. Methodology

### 3.1. Additional Damage Characterization of Surrounding Rock after Excavation

This study uses the Mohr-Coulomb model commonly used in coal mining as an example. According to the existing research, the five parameters that determine the mechanical behaviour of rock are elastic modulus, Poisson’s ratio, cohesion, friction angle, and tensile strength [23, 24]. Moreover, the change of the elastic modulus significantly influences the mechanical properties of the surrounding rock, so this study only characterizes the additional damage to the surrounding rock by the deterioration of the elastic modulus.According to previous studies, Weibull distribution is one of the most effective statistical methods to describe the physical properties of rock microelements. Considering the engineering complexity and the heterogeneity of rocks, the Weibull distribution function is used to describe rock parameters (such as elastic modulus, strength, stress, and strain) in the spatial distribution of the roof. Its expression [25] is
(1)fx=mx0xx0m−1e−x/x0m,where fx is the probability density function of rock parameters, x0 is the scale parameter, and m is the shape parameter of the Weibull distribution. The roof displacement damage variable D can be obtained by integrating as follows:
(2)D=0,ε≤ε01−e−ε−ε0/x0m,ε>ε0,where ε0 is the damage strain threshold, to describe the change of the elastic modulus of the rock mass under the action of external force, and the loading process is divided according to the damage to the rock unit body, as shown in Figure 2. When ε<ε0, the rock is in the elastic deformation stage, at this time, D=0, and when the axial strain is greater than the strain threshold ε0, the rock material is damaged. When ε0<ε<εp, where εp is the strain at peak strength, it is in the microdamage zone, and the rock shows the development of microscopic fractures. After the stress peak, when εp<ε<εc, where εc is the strain when the specimen is completely destroyed, it is the macroscopic crack damage zone. When ε>εc, it is the destruction stage.Figure 2
Schematic diagram of rock damage zoning.According to previous studies, it is assumed that the elastic modulus of the roof is similar to the distribution law of the strain parameter. Based on the theory of continuum damage mechanics, Kawamoto et al. [26] proposed that the elastic modulus E of the microelement is related to the damage variable D as follows:
(3)E=1−DE0,where E0 is the mean value of the elastic modulus of the rock, and the relationship between the elastic modulus of the roof and its strain can be obtained by combining the Equations (2) and (3):
(4)E=E0,ε≤ε0E0e−ε−ε0/x0m,ε>ε0.As shown in Figure3, when E0=30GPa, different m and x0 have different effects on elastic modulus. It can be seen from the figure that with the increase of m, the curve rotates clockwise around the fixed point, the rate of change of the elastic modulus increases, the strain area where the modulus of the rock mass changes is reduced, and the distribution of the elastic modulus of the rock mass is more uniform; with the increase of x0, the elastic modulus change rate decreases, and the elastic modulus curve is approximately horizontally stretched in the direction of increasing axial strain.Figure 3
Effects ofm and x0 on the relationship between elastic modulus and strain.
(a)
The influence ofm(b)
The influence ofx0m and x0 can characterize the inhomogeneity of rock materials. When appropriate m and x0 are used, the numerical analysis results of the fracturing process or mechanical behaviour of rock samples can be in good agreement with the experimental results [27, 28]. According to the classical damage theory [29], during uniaxial compression, the damage constitutive model of the Weibull distribution of rock mass can be transformed into
(5)σ=E0εe−ε/x0m.Since the increment of rock strength is 0 at the peak value of rock strength, the derivative of Equation (5) with respect to ε can be obtained:
(6)εp=x0m−1/m,where εp is the rock peak strain and m is introduced as follows:
(7)m=1lnE0εp/σp.Simultaneously with Equations (5)–(7), x0 and m can be obtained by directly selecting the parameters obtained from the rock mechanics experiment.In order to simulate the change of elastic modulus in the additional damage zone caused by rock excavation, a new numerical simulation technique is developed and applied to the Mohr-Coulomb model in FLAC3D. This method is based on formulas (5)–(7) to formulate the FISH program so that the elastic modulus changes with rock damage. Based on the existing rock mechanics test experience, it is assumed that when the rock damage reaches 70%, the bearing capacity will be lost, and its elastic modulus will no longer decrease [30]. This method overcomes the limitation of conventional simulation techniques to keep the elastic modulus constant without considering additional damage to the rock mass [31, 32].


### 3.2. Mechanical Characterization of Gob Gangue

Gob gangue is generally simulated by the Double-Yield or elastic models [33]. However, the Double-Yield model has many parameters, and the parameter checking is more complicated, so the elastic model is used for simulation. The elastic modulus of the backfill is not fixed during the compaction process. In order to accurately represent the compactional characteristics of the gangue, the elastic model needs to be redeveloped.According to Salamon’s theory of goaf compaction [34], the stress-strain relationship of gangue in goaf during compaction is
(8)σ=E1ε11−ε1/εm,where σ is the vertical stress of the goaf; E1 is the initial elastic modulus of the gangue; ε is the current vertical strain; εm is the maximum vertical strain.According to the research of Yavuz [35], εm can be expressed by the following equation:
(9)εm=b−1b,where b is the comprehensive crushing expansion coefficient of the filling body, taken as 1.3 according to the site conditions.Simultaneous (8)–(9) equations, the elastic modulus E2 of gangue is
(10)E2=bb−1E1bε1−b+12ε−E1b−1bε1−b+1.Based on Equation (10), the FISH program can be programmed to change the elastic modulus with gangue compaction to achieve an accurate simulation of gangue compaction.In order to verify the feasibility of secondary development of the elastic model, Double-Yield, Mohr-Coulomb, and elastic models are used to compare the gangue in goaf, as shown in Figure4. It can be seen from the diagram that the elastic model and Mohr-Coulomb model have little difference with the theoretical analysis stress when the gangue compression is slight. There is a significant error when the gangue compression is large. Therefore, the elastic and Mohr-Coulomb models are unsuitable for full compaction. The improved elastic and Double-Yield models have little error with the theoretical analysis. However, when the compaction degree is not large, the stress ratio of the Double-Yield model is lower than the theoretical results, so it is unsuitable when the gangue compaction is insufficient. The numerical simulation results of the improved elastic model agree with the theoretical formula. They have good consistency with the laboratory’s stress-strain curve of the gravel compaction test [36]. Therefore, it is feasible to use the method proposed in this study to simulate the compaction of the goaf.Figure 4
Compaction characteristics of models under different constitutive models.


### 3.3. Mechanical Characterization of Filling Materials

Inorganic paste filling material belongs to rock-like material, composed of two kinds of materials A and B, mixed in the same proportion and then poured with a certain proportion of water. Due to its prominent plastic characteristics after reaching compressive strength, the Strain-Softening model can simulate the mechanical characteristics. However, the parameters need to be checked before simulation to improve the accuracy of the simulation. According to the field experience, four filling material schemes are set [37], and the water-cement ratio is 2 : 1, 3 : 1, 4 : 1, and 5 : 1, respectively. The stress-strain curve obtained from the test and the checked curve are shown in Figure 5, and the checking parameters are shown in Table 1.Figure 5
Checking the parameters of paste filling material.Table 1
Simulated mechanical parameters of paste materials.
Scheme numberWater-cement ratioElastic modulus/MPaPoisson’sratioCohesion/MPaAngle of internal friction/°Tensile strength/MPaAngle of dilatancy/°12 : 11360.250.70331523 : 11320.260.32331534 : 1450.30.19301545 : 1350.30.153015


### 3.4. Numerical Simulation Method

Considering the compaction bearing capacity characteristics of the goaf caving zone and the deterioration characteristics of roof rock mass affected by mining, the coupling analysis method of mining-induced stress field and mining-induced grazing rock is proposed, as shown in Figure6.Figure 6
Flow chart of numerical simulation coupling method.Before the model is established, it is necessary to determine the gangue compaction area, the rock caving zone. From literature [38], it can be known that the calculation formula of the width of the filling zone at the lower part of the fully mechanized mining face of the high-dip angle coal seam is
(11)L1=LH1+H2bM+H1+H2b,where L is the length of the working face, L1 is the width of the filling zone, H1 is the thickness of the rock layer that is quickly caving, H2 is the thickness of the rock that is lagging behind, and M is the mining height.Combined with the 211112 working face parameters of Xinji No. 2 Mine:(i)
The working face length is 139 m(ii)
The mining height is 3.6 m(iii)
The thickness of the highly caving rock layer is the thickness of the composite roof of coal and sandy mudstone (5.6 m)(iv)
Thicker sandy mudstone is 4.7 m thickSubstituting the above parameter values into Equation (11), the width of the 211112 working face filling belt is 122 m, and the height of the caving belt is 13.9 m.According to the range of the caving zone calculated above, the acceptable numerical model of the caving rock mass is established. The FISH program based on the theory of goaf compaction is used to simulate the mechanical characteristics of the compaction process of the caving zone in the goaf. At the same time, the FISH program, considering the additional damage from the excavation, is used to simulate the degradation process of the overlying strata. After excavating to the front of the filling section, determine the paste filling area of the working face, assign the filling area to the checked backfill material mechanical parameters and constitutive model, and then continue to run the above procedure until the model is balanced.The following problems are solved by simulation with the method proposed in this article:(i)
The Weibull distribution function is used to describe the rock parameters, which more accurately reflects the damage process of the roof rock mass after excavation. The accuracy of the simulation can be improved by selecting a more reliable method for determining the Weibull distribution parameters(ii)
The numerical simulation method proposed in this paper can accurately simulate the caving and compaction characteristics of the goaf. After the height of the caving, the zone is determined, and the damage variable can be used to predict the size of fracture development quantitatively(iii)
Based on the numerical simulation method of goaf compaction and roof damage characteristics, it can also provide an analysis method for surface subsidence prediction, coal pillar retention, and roof disaster prevention and control in coal seam mining


## 3.1. Additional Damage Characterization of Surrounding Rock after Excavation

This study uses the Mohr-Coulomb model commonly used in coal mining as an example. According to the existing research, the five parameters that determine the mechanical behaviour of rock are elastic modulus, Poisson’s ratio, cohesion, friction angle, and tensile strength [23, 24]. Moreover, the change of the elastic modulus significantly influences the mechanical properties of the surrounding rock, so this study only characterizes the additional damage to the surrounding rock by the deterioration of the elastic modulus.According to previous studies, Weibull distribution is one of the most effective statistical methods to describe the physical properties of rock microelements. Considering the engineering complexity and the heterogeneity of rocks, the Weibull distribution function is used to describe rock parameters (such as elastic modulus, strength, stress, and strain) in the spatial distribution of the roof. Its expression [25] is
(1)fx=mx0xx0m−1e−x/x0m,where fx is the probability density function of rock parameters, x0 is the scale parameter, and m is the shape parameter of the Weibull distribution. The roof displacement damage variable D can be obtained by integrating as follows:
(2)D=0,ε≤ε01−e−ε−ε0/x0m,ε>ε0,where ε0 is the damage strain threshold, to describe the change of the elastic modulus of the rock mass under the action of external force, and the loading process is divided according to the damage to the rock unit body, as shown in Figure 2. When ε<ε0, the rock is in the elastic deformation stage, at this time, D=0, and when the axial strain is greater than the strain threshold ε0, the rock material is damaged. When ε0<ε<εp, where εp is the strain at peak strength, it is in the microdamage zone, and the rock shows the development of microscopic fractures. After the stress peak, when εp<ε<εc, where εc is the strain when the specimen is completely destroyed, it is the macroscopic crack damage zone. When ε>εc, it is the destruction stage.Figure 2
Schematic diagram of rock damage zoning.According to previous studies, it is assumed that the elastic modulus of the roof is similar to the distribution law of the strain parameter. Based on the theory of continuum damage mechanics, Kawamoto et al. [26] proposed that the elastic modulus E of the microelement is related to the damage variable D as follows:
(3)E=1−DE0,where E0 is the mean value of the elastic modulus of the rock, and the relationship between the elastic modulus of the roof and its strain can be obtained by combining the Equations (2) and (3):
(4)E=E0,ε≤ε0E0e−ε−ε0/x0m,ε>ε0.As shown in Figure3, when E0=30GPa, different m and x0 have different effects on elastic modulus. It can be seen from the figure that with the increase of m, the curve rotates clockwise around the fixed point, the rate of change of the elastic modulus increases, the strain area where the modulus of the rock mass changes is reduced, and the distribution of the elastic modulus of the rock mass is more uniform; with the increase of x0, the elastic modulus change rate decreases, and the elastic modulus curve is approximately horizontally stretched in the direction of increasing axial strain.Figure 3
Effects ofm and x0 on the relationship between elastic modulus and strain.
(a)
The influence ofm(b)
The influence ofx0m and x0 can characterize the inhomogeneity of rock materials. When appropriate m and x0 are used, the numerical analysis results of the fracturing process or mechanical behaviour of rock samples can be in good agreement with the experimental results [27, 28]. According to the classical damage theory [29], during uniaxial compression, the damage constitutive model of the Weibull distribution of rock mass can be transformed into
(5)σ=E0εe−ε/x0m.Since the increment of rock strength is 0 at the peak value of rock strength, the derivative of Equation (5) with respect to ε can be obtained:
(6)εp=x0m−1/m,where εp is the rock peak strain and m is introduced as follows:
(7)m=1lnE0εp/σp.Simultaneously with Equations (5)–(7), x0 and m can be obtained by directly selecting the parameters obtained from the rock mechanics experiment.In order to simulate the change of elastic modulus in the additional damage zone caused by rock excavation, a new numerical simulation technique is developed and applied to the Mohr-Coulomb model in FLAC3D. This method is based on formulas (5)–(7) to formulate the FISH program so that the elastic modulus changes with rock damage. Based on the existing rock mechanics test experience, it is assumed that when the rock damage reaches 70%, the bearing capacity will be lost, and its elastic modulus will no longer decrease [30]. This method overcomes the limitation of conventional simulation techniques to keep the elastic modulus constant without considering additional damage to the rock mass [31, 32].


## 3.2. Mechanical Characterization of Gob Gangue

Gob gangue is generally simulated by the Double-Yield or elastic models [33]. However, the Double-Yield model has many parameters, and the parameter checking is more complicated, so the elastic model is used for simulation. The elastic modulus of the backfill is not fixed during the compaction process. In order to accurately represent the compactional characteristics of the gangue, the elastic model needs to be redeveloped.According to Salamon’s theory of goaf compaction [34], the stress-strain relationship of gangue in goaf during compaction is
(8)σ=E1ε11−ε1/εm,where σ is the vertical stress of the goaf; E1 is the initial elastic modulus of the gangue; ε is the current vertical strain; εm is the maximum vertical strain.According to the research of Yavuz [35], εm can be expressed by the following equation:
(9)εm=b−1b,where b is the comprehensive crushing expansion coefficient of the filling body, taken as 1.3 according to the site conditions.Simultaneous (8)–(9) equations, the elastic modulus E2 of gangue is
(10)E2=bb−1E1bε1−b+12ε−E1b−1bε1−b+1.Based on Equation (10), the FISH program can be programmed to change the elastic modulus with gangue compaction to achieve an accurate simulation of gangue compaction.In order to verify the feasibility of secondary development of the elastic model, Double-Yield, Mohr-Coulomb, and elastic models are used to compare the gangue in goaf, as shown in Figure4. It can be seen from the diagram that the elastic model and Mohr-Coulomb model have little difference with the theoretical analysis stress when the gangue compression is slight. There is a significant error when the gangue compression is large. Therefore, the elastic and Mohr-Coulomb models are unsuitable for full compaction. The improved elastic and Double-Yield models have little error with the theoretical analysis. However, when the compaction degree is not large, the stress ratio of the Double-Yield model is lower than the theoretical results, so it is unsuitable when the gangue compaction is insufficient. The numerical simulation results of the improved elastic model agree with the theoretical formula. They have good consistency with the laboratory’s stress-strain curve of the gravel compaction test [36]. Therefore, it is feasible to use the method proposed in this study to simulate the compaction of the goaf.Figure 4
Compaction characteristics of models under different constitutive models.


## 3.3. Mechanical Characterization of Filling Materials

Inorganic paste filling material belongs to rock-like material, composed of two kinds of materials A and B, mixed in the same proportion and then poured with a certain proportion of water. Due to its prominent plastic characteristics after reaching compressive strength, the Strain-Softening model can simulate the mechanical characteristics. However, the parameters need to be checked before simulation to improve the accuracy of the simulation. According to the field experience, four filling material schemes are set [37], and the water-cement ratio is 2 : 1, 3 : 1, 4 : 1, and 5 : 1, respectively. The stress-strain curve obtained from the test and the checked curve are shown in Figure 5, and the checking parameters are shown in Table 1.Figure 5
Checking the parameters of paste filling material.Table 1
Simulated mechanical parameters of paste materials.
Scheme numberWater-cement ratioElastic modulus/MPaPoisson’sratioCohesion/MPaAngle of internal friction/°Tensile strength/MPaAngle of dilatancy/°12 : 11360.250.70331523 : 11320.260.32331534 : 1450.30.19301545 : 1350.30.153015


## 3.4. Numerical Simulation Method

Considering the compaction bearing capacity characteristics of the goaf caving zone and the deterioration characteristics of roof rock mass affected by mining, the coupling analysis method of mining-induced stress field and mining-induced grazing rock is proposed, as shown in Figure6.Figure 6
Flow chart of numerical simulation coupling method.Before the model is established, it is necessary to determine the gangue compaction area, the rock caving zone. From literature [38], it can be known that the calculation formula of the width of the filling zone at the lower part of the fully mechanized mining face of the high-dip angle coal seam is
(11)L1=LH1+H2bM+H1+H2b,where L is the length of the working face, L1 is the width of the filling zone, H1 is the thickness of the rock layer that is quickly caving, H2 is the thickness of the rock that is lagging behind, and M is the mining height.Combined with the 211112 working face parameters of Xinji No. 2 Mine:(i)
The working face length is 139 m(ii)
The mining height is 3.6 m(iii)
The thickness of the highly caving rock layer is the thickness of the composite roof of coal and sandy mudstone (5.6 m)(iv)
Thicker sandy mudstone is 4.7 m thickSubstituting the above parameter values into Equation (11), the width of the 211112 working face filling belt is 122 m, and the height of the caving belt is 13.9 m.According to the range of the caving zone calculated above, the acceptable numerical model of the caving rock mass is established. The FISH program based on the theory of goaf compaction is used to simulate the mechanical characteristics of the compaction process of the caving zone in the goaf. At the same time, the FISH program, considering the additional damage from the excavation, is used to simulate the degradation process of the overlying strata. After excavating to the front of the filling section, determine the paste filling area of the working face, assign the filling area to the checked backfill material mechanical parameters and constitutive model, and then continue to run the above procedure until the model is balanced.The following problems are solved by simulation with the method proposed in this article:(i)
The Weibull distribution function is used to describe the rock parameters, which more accurately reflects the damage process of the roof rock mass after excavation. The accuracy of the simulation can be improved by selecting a more reliable method for determining the Weibull distribution parameters(ii)
The numerical simulation method proposed in this paper can accurately simulate the caving and compaction characteristics of the goaf. After the height of the caving, the zone is determined, and the damage variable can be used to predict the size of fracture development quantitatively(iii)
Based on the numerical simulation method of goaf compaction and roof damage characteristics, it can also provide an analysis method for surface subsidence prediction, coal pillar retention, and roof disaster prevention and control in coal seam mining


## 4. Numerical Simulation Analysis

### 4.1. Model Parameters and Schemes

To verify the correctness of the above numerical simulation method and the feasibility of using the FISH language to realize the algorithm, the new approach and the traditional method in this paper are used to analyze the stope stress of the above model. The established numerical model is shown in Figure7. The simulated 211112 fully mechanized mining face has a length of 139 m, a strike length of 200 m, and a mining height of 3.6 m, advancing along the model strike. In order to reduce the numerical simulation boundary effect, 50 m boundary protection coal pillars were set around the model excavation. The calculation model adopts the elastic-plastic model and the Mohr-Coulomb failure criterion. The roof rock layer adopts the additional damage numerical simulation method, and the caving zone adopts the numerical simulation program based on the goaf compaction theory. The specific mechanical parameters of coal and rock layers can be obtained from the rock mechanical test results provided by the relevant tests. At the same time, the additional damage characterization parameters can be obtained using formulas (2)–(4), as shown in Table 2.Figure 7
Numerical model.Table 2
Rock mass mechanical properties.
LithologyDensity/kg·m-3K/GPaG/GPac/MPaφ/°σrm/MPaQuartzose sandstone28005.574.5311286.7Siltstone26803.62.78253.5Sandy mudstone26912.82.14.6252.4Coal14101.730.820.2200.2Mudstone25674.32.80.7301.68Fine sandstone280016.0412.13.5435.0Sandy mudstone26916.95.64.6352.4According to the working process of working face mining, working face excavation can be divided into three processes: (1) working face mining and excavation stage; (2) confirming the filling area and constructing the filling body stage; (3) excavating the filling body and moving the support forward. Among them, four options in Section3.3 are selected for the filling material in the filling area in front of the working face.


### 4.2. Analysis on Stress Characteristics of Working Face Excavation

The traditional and new methods proposed in this paper are used to simulate the excavation of the working face. As shown in Figure8, after the working face is pushed 180 m, a stress measurement line is arranged along the middle of the caving area of the upper roof of the coal seam, and the vertical stress distribution before and after the work is recorded. It can be seen from Figure 8 that the stress in the goaf shows a zigzag increase and decrease, which is mainly due to the step-by-step excavation and the sequence of roof stress release. Compared with the traditional method, the stress recovery of the gob in the new method is more pronounced, but the stress in the gob is challenging to recover to the original rock stress in both methods. This phenomenon is caused by the fact that the old roof is siltstone with high strength and large thickness (9.3 m), which supports the overlying rock to a certain extent. This phenomenon is also consistent with the research results of other scholars [39]. In addition, due to the consideration of the deterioration characteristics of the overlying strata, the maximum stress in front of the working face using the method in this paper is significantly lower than that using the traditional method.Figure 8
Stress distribution characteristics of goaf.In the mining process of working face, the stress of surrounding rock is redistributed due to the transfer of stress, which causes different degrees of damage to the roof. According to the additional damage characterization method of surrounding rock in Section3.1, the damage distribution of the roof is displayed by the FISH function, as shown in Figure 9.Figure 9
Roof damage distribution.It can be seen from Figure9 that in the middle and lower part of the working face, the damage amount to the roof in the rear caving area along the strike direction shows the characteristics of first increasing, then decreasing, and then increasing. The damage reaches a small value in the middle of the goaf and a maximum value near the working face. In the upper part of the working face, the damage amount does not change much, about 0.7; along the inclination direction, the damage in the middle and lower parts of the working face is relatively small, and the damage amount increases sharply in the upper part. The damage to the roof with a large inclination angle presents an asymmetric structure that is large in the upper part and small in the lower part, similar to roof cracks in general large-inclined working face. The field practice also shows that the LAWCA in the front of the working face appears in the middle and upper part of the working face, indicating that the simulation method in this paper can truly reflect the impact of the gob caving on the roof damage.


### 4.3. Optimization of the Filling Material Ratio

According to the scheme described in Section3.3, the filling simulation of the empty roof area is carried out to form a numerical model of caving at the rear of the working face and filling with paste at the front. Three measuring lines are set on the filling body from top to bottom, namely, Line 1 ~ Line 3, as shown in Figure 10. Since the filling area is a pressure relief area, the stress on the filling body of each scheme is relatively small, the vertical stress in most areas is below 1 MPa, and the stress concentration is mainly at the contact between the filling body and the coal wall. The backfill body of the four cases can stabilise the working face when it is not mined.Figure 10
Vertical stress distribution of backfill before excavation.
(a)
Case 1(b)
Case 2(c)
Case 3(d)
Case 4When the working face continues to be mined, the roof above the backfill body tends to rotate and sink to the side of the goaf again, the supporting pressure of the working face changes, and the backfill body begins to support the roof. Two survey lines are set at 2 m (in the filling body) and 6 m (in the coal wall) in front of the excavation, as shown in Figure11. It can be seen from Figure 11 that the filling body, as an artificial roof, plays a leading role in the work of the support and the roof of the working face. When the working face starts to advance, the vertical stress in the backfill body and the coal wall decreases with the increase of the water-cement ratio of the backfill body. The vertical stress in the backfill body of scheme 1 and scheme 2 is not significantly different, and the stress difference in the coal wall is also less than 1 MPa. According to the strength characteristics of the backfill in Section 3.3, most of the backfill in case 4 has entered the yield stage, and the bearing capacity is weak and unsuitable for maintaining the stability of the working face.Figure 11
Vertical stress distribution of backfill and coal wall after excavation.
(a)
Vertical stress of filling body in different schemes(b)
Vertical stress of coal wall in different schemes


## 4.1. Model Parameters and Schemes

To verify the correctness of the above numerical simulation method and the feasibility of using the FISH language to realize the algorithm, the new approach and the traditional method in this paper are used to analyze the stope stress of the above model. The established numerical model is shown in Figure7. The simulated 211112 fully mechanized mining face has a length of 139 m, a strike length of 200 m, and a mining height of 3.6 m, advancing along the model strike. In order to reduce the numerical simulation boundary effect, 50 m boundary protection coal pillars were set around the model excavation. The calculation model adopts the elastic-plastic model and the Mohr-Coulomb failure criterion. The roof rock layer adopts the additional damage numerical simulation method, and the caving zone adopts the numerical simulation program based on the goaf compaction theory. The specific mechanical parameters of coal and rock layers can be obtained from the rock mechanical test results provided by the relevant tests. At the same time, the additional damage characterization parameters can be obtained using formulas (2)–(4), as shown in Table 2.Figure 7
Numerical model.Table 2
Rock mass mechanical properties.
LithologyDensity/kg·m-3K/GPaG/GPac/MPaφ/°σrm/MPaQuartzose sandstone28005.574.5311286.7Siltstone26803.62.78253.5Sandy mudstone26912.82.14.6252.4Coal14101.730.820.2200.2Mudstone25674.32.80.7301.68Fine sandstone280016.0412.13.5435.0Sandy mudstone26916.95.64.6352.4According to the working process of working face mining, working face excavation can be divided into three processes: (1) working face mining and excavation stage; (2) confirming the filling area and constructing the filling body stage; (3) excavating the filling body and moving the support forward. Among them, four options in Section3.3 are selected for the filling material in the filling area in front of the working face.


## 4.2. Analysis on Stress Characteristics of Working Face Excavation

The traditional and new methods proposed in this paper are used to simulate the excavation of the working face. As shown in Figure8, after the working face is pushed 180 m, a stress measurement line is arranged along the middle of the caving area of the upper roof of the coal seam, and the vertical stress distribution before and after the work is recorded. It can be seen from Figure 8 that the stress in the goaf shows a zigzag increase and decrease, which is mainly due to the step-by-step excavation and the sequence of roof stress release. Compared with the traditional method, the stress recovery of the gob in the new method is more pronounced, but the stress in the gob is challenging to recover to the original rock stress in both methods. This phenomenon is caused by the fact that the old roof is siltstone with high strength and large thickness (9.3 m), which supports the overlying rock to a certain extent. This phenomenon is also consistent with the research results of other scholars [39]. In addition, due to the consideration of the deterioration characteristics of the overlying strata, the maximum stress in front of the working face using the method in this paper is significantly lower than that using the traditional method.Figure 8
Stress distribution characteristics of goaf.In the mining process of working face, the stress of surrounding rock is redistributed due to the transfer of stress, which causes different degrees of damage to the roof. According to the additional damage characterization method of surrounding rock in Section3.1, the damage distribution of the roof is displayed by the FISH function, as shown in Figure 9.Figure 9
Roof damage distribution.It can be seen from Figure9 that in the middle and lower part of the working face, the damage amount to the roof in the rear caving area along the strike direction shows the characteristics of first increasing, then decreasing, and then increasing. The damage reaches a small value in the middle of the goaf and a maximum value near the working face. In the upper part of the working face, the damage amount does not change much, about 0.7; along the inclination direction, the damage in the middle and lower parts of the working face is relatively small, and the damage amount increases sharply in the upper part. The damage to the roof with a large inclination angle presents an asymmetric structure that is large in the upper part and small in the lower part, similar to roof cracks in general large-inclined working face. The field practice also shows that the LAWCA in the front of the working face appears in the middle and upper part of the working face, indicating that the simulation method in this paper can truly reflect the impact of the gob caving on the roof damage.


## 4.3. Optimization of the Filling Material Ratio

According to the scheme described in Section3.3, the filling simulation of the empty roof area is carried out to form a numerical model of caving at the rear of the working face and filling with paste at the front. Three measuring lines are set on the filling body from top to bottom, namely, Line 1 ~ Line 3, as shown in Figure 10. Since the filling area is a pressure relief area, the stress on the filling body of each scheme is relatively small, the vertical stress in most areas is below 1 MPa, and the stress concentration is mainly at the contact between the filling body and the coal wall. The backfill body of the four cases can stabilise the working face when it is not mined.Figure 10
Vertical stress distribution of backfill before excavation.
(a)
Case 1(b)
Case 2(c)
Case 3(d)
Case 4When the working face continues to be mined, the roof above the backfill body tends to rotate and sink to the side of the goaf again, the supporting pressure of the working face changes, and the backfill body begins to support the roof. Two survey lines are set at 2 m (in the filling body) and 6 m (in the coal wall) in front of the excavation, as shown in Figure11. It can be seen from Figure 11 that the filling body, as an artificial roof, plays a leading role in the work of the support and the roof of the working face. When the working face starts to advance, the vertical stress in the backfill body and the coal wall decreases with the increase of the water-cement ratio of the backfill body. The vertical stress in the backfill body of scheme 1 and scheme 2 is not significantly different, and the stress difference in the coal wall is also less than 1 MPa. According to the strength characteristics of the backfill in Section 3.3, most of the backfill in case 4 has entered the yield stage, and the bearing capacity is weak and unsuitable for maintaining the stability of the working face.Figure 11
Vertical stress distribution of backfill and coal wall after excavation.
(a)
Vertical stress of filling body in different schemes(b)
Vertical stress of coal wall in different schemes


## 5. Case Study

The numerical simulation results selected the water-cement ratio filling material used in scheme 2. According to the experimental test, the material has a fast setting speed, high early strength, and adjustable and good later stability. When the strain reaches 10%, the compressive strength can be maintained at more than 70% of the peak strength, and the pumping distance is more than 3000 m. It can withstand 800~1000°C high temperature and has good flame retardancy. Therefore, it is an excellent choice to use paste filling material to fill areas that need to improve the integrity of coal and rock mass, such as roof crushing and side leakage.


### 5.1. Filling Process

In front of the No. 72 bracket on the working face, set up a “#”-shaped wooden stack to the coal wall. The wooden stack is constructed with a 1.4 m-long wooden road board. The wood board with a thickness of not less than 20 mm will seal the bottom of the filling area. After the board is nailed, a layer of air duct cloth will be placed on the outside of the board to prevent slurry leakage.Combined with the above numerical simulations, the water-cement ratio is determined to be 3 : 1; 8 workers are required for each shift, of which 6 are feeding materials, 1 is operating the equipment, and 1 is at the filling point. Filling speed and construction period: the actual flow rate of the grouting pump is 15 m3/h. Calculated by filling 8 h per shift, the filling volume of a tiny shift is 15×8=120m3, and the entire leaking roof area needs 9 small shifts. Setting time: the initial setting speed of the slurry is fast (5 min~20 min). The final strength of the solidified body can reach more than 2.5 MPa. The integrity of the filling body is good. The specific process flow in the filling process is as follows:
(i)
Layout of the grouting pump station. As shown in Figure 12, the pumping station is arranged in the return airway of 211112. Two 2ZBYSB300-90/5–15-55 double-liquid grouting pumps are used. Material A and material B are equipped with two JDW-1000S mixing drums. The volume of each mixing drum is 1000 L. The material field is arranged near the mixing drum.(ii)
Pipe connection. Connect the pipeline for conveying the slurry and use a matching straight-through connection between the high-pressure hoses with a diameter of 51 mm. The double-trip high-pressure hose passes through the mixer at the mixing point and becomes a one-trip mixing pipe. Connect from the 211112 return air roadway to the 72# hydraulic support.(iii)
Test system and pulping. First, use clean water to test whether the communication between the mixing tank, the grouting pump, the filling pipeline, the filling point, and the pump station is typical. After everything is normal, the slurry can be filled. When grouting, it should be noted that materials A and B cannot be mixed, and the stirring time is not less than 5 minutes. Only the slurry that is stirred evenly and meets the requirements of the stirring time can be pumped.(iv)
Pumping and cleaning equipment. The suction head of the grouting pump is put into the mixing drum with uniform mixing, and the grouting pump is started to pump. After the mixing pipe flows out of the uniform mixture slurry, the mixing pipe is placed in the filling area to fill the 72~95# bracket officially. After the slurry pump or filling is completed, it is necessary to pump clean water in time to clean the grouting pump, mixing tank, and pipeline.Figure 12
Grouting equipment configuration diagram.


### 5.2. Filling Effect Verification

As shown in Figure13, after filling the LAWCA with filling materials, when the working face moves forward through the filling area, the rib spalling and roof fall do not occur in the working face. The filling body does not appear in the phenomenon of water absorption adhesive belt conveyor and other transportation equipment and silo blocking. The coal-winning machine cuts the filling body to form a block structure, and the filling body and coal are mined together.Figure 13
Filling and excavation effect of the project site.As shown in Figure14, the support cannot fully contact the roof before filling due to the roof falling off the sheet in front of the working face. The stress of the upper and middle hydraulic support in the working face before filling body construction is generally tiny (blue area). After the filling is completed, the filling body undertakes the partial load of the upper strata to replace the interaction between the immediate roof and the hydraulic support so that the stress of the support tends to be normal, the stability of the working face is increased, and the safe and efficient mining of the working face is ensured. The working face passes through the leakage area smoothly, which proves that the filling material ratio parameters determined by the above simulation method can meet the requirements and verify the correctness of the simulation method.Figure 14
Working resistance distribution of hydraulic support.


## 5.1. Filling Process

In front of the No. 72 bracket on the working face, set up a “#”-shaped wooden stack to the coal wall. The wooden stack is constructed with a 1.4 m-long wooden road board. The wood board with a thickness of not less than 20 mm will seal the bottom of the filling area. After the board is nailed, a layer of air duct cloth will be placed on the outside of the board to prevent slurry leakage.Combined with the above numerical simulations, the water-cement ratio is determined to be 3 : 1; 8 workers are required for each shift, of which 6 are feeding materials, 1 is operating the equipment, and 1 is at the filling point. Filling speed and construction period: the actual flow rate of the grouting pump is 15 m3/h. Calculated by filling 8 h per shift, the filling volume of a tiny shift is 15×8=120m3, and the entire leaking roof area needs 9 small shifts. Setting time: the initial setting speed of the slurry is fast (5 min~20 min). The final strength of the solidified body can reach more than 2.5 MPa. The integrity of the filling body is good. The specific process flow in the filling process is as follows:
(i)
Layout of the grouting pump station. As shown in Figure 12, the pumping station is arranged in the return airway of 211112. Two 2ZBYSB300-90/5–15-55 double-liquid grouting pumps are used. Material A and material B are equipped with two JDW-1000S mixing drums. The volume of each mixing drum is 1000 L. The material field is arranged near the mixing drum.(ii)
Pipe connection. Connect the pipeline for conveying the slurry and use a matching straight-through connection between the high-pressure hoses with a diameter of 51 mm. The double-trip high-pressure hose passes through the mixer at the mixing point and becomes a one-trip mixing pipe. Connect from the 211112 return air roadway to the 72# hydraulic support.(iii)
Test system and pulping. First, use clean water to test whether the communication between the mixing tank, the grouting pump, the filling pipeline, the filling point, and the pump station is typical. After everything is normal, the slurry can be filled. When grouting, it should be noted that materials A and B cannot be mixed, and the stirring time is not less than 5 minutes. Only the slurry that is stirred evenly and meets the requirements of the stirring time can be pumped.(iv)
Pumping and cleaning equipment. The suction head of the grouting pump is put into the mixing drum with uniform mixing, and the grouting pump is started to pump. After the mixing pipe flows out of the uniform mixture slurry, the mixing pipe is placed in the filling area to fill the 72~95# bracket officially. After the slurry pump or filling is completed, it is necessary to pump clean water in time to clean the grouting pump, mixing tank, and pipeline.Figure 12
Grouting equipment configuration diagram.


## 5.2. Filling Effect Verification

As shown in Figure13, after filling the LAWCA with filling materials, when the working face moves forward through the filling area, the rib spalling and roof fall do not occur in the working face. The filling body does not appear in the phenomenon of water absorption adhesive belt conveyor and other transportation equipment and silo blocking. The coal-winning machine cuts the filling body to form a block structure, and the filling body and coal are mined together.Figure 13
Filling and excavation effect of the project site.As shown in Figure14, the support cannot fully contact the roof before filling due to the roof falling off the sheet in front of the working face. The stress of the upper and middle hydraulic support in the working face before filling body construction is generally tiny (blue area). After the filling is completed, the filling body undertakes the partial load of the upper strata to replace the interaction between the immediate roof and the hydraulic support so that the stress of the support tends to be normal, the stability of the working face is increased, and the safe and efficient mining of the working face is ensured. The working face passes through the leakage area smoothly, which proves that the filling material ratio parameters determined by the above simulation method can meet the requirements and verify the correctness of the simulation method.Figure 14
Working resistance distribution of hydraulic support.


## 6. Conclusion

(1)
A FLAC3D numerical simulation method considering the additional damage of excavation and goaf compaction theory is proposed, and the Weibull distribution function is introduced to correct the mechanical parameters of the elastic modulus of the rock mass after excavation, which more accurately reflects the excavation. The damage process of the roof rock mass after excavation, through the secondary development of the elastic model, has realized the accurate simulation of the caving rock mass(2)
Combined with the actual situation of the 211112 working face of Xinji No. 2 Mine, the numerical model was established by FLAC3D software, and the distribution of stress and roof damage was analyzed from the perspective of stress and roof damage. Compared with the current literature research results, it was verified that the additional damage of excavation and goaf pressure were considered—reliability and feasibility of real-theoretic algorithms(3)
The influence of paste filling in front of the working face on the excavation stability of the working face is studied by using the method proposed in this paper, and four kinds of filling material ratio schemes are proposed. The failure characteristics of the filling body, coal wall, and roof under different filling parameters were studied. Combined with the mechanical requirements for the backfilling of the backfill, the backfill body, and the rotation of the suspended roof, the backfilling scheme of inorganic grouting reinforcement materials with a water-solid ratio of 3 : 1 was determined(4)
According to the numerical simulation results, the relevant backfilling process was designed, and the corresponding safety guarantee technical measures were formulated. Field practice shows that the proposed backfill material ratio scheme has realized the coal and rock stability control in the LAWCA of the large-dip-angle working face



---

*Source: 5352974-2023-02-06.xml*