# Numerical Study on the Movement Laws of Overlying Strata in Shallow-Buried Stope Based on the Goaf Compaction Effect

Authors: Guangchao Zhang; Guangzhe Tao; Miao Chen; You Li; Pan Li; Yonghui Lai
Journal: Shock and Vibration
 (2021)

---

## Abstract

This study presents an integrated approach including the theoretical analysis and numerical modelling to investigate the failure characteristics of the overlying strata in the shallow-buried stope. The mechanical characteristics of the caving zone are first revealed and then calibrated by using the double-yield model. The theoretical results show that the mechanical properties of the collapsed rock mass are closely related to its crushing expansion coefficient and uniaxial compressive strength. The vertical stress of the collapsed rock mass increases slowly with the strain and then increases exponentially after a certain critical strain. The simulation indicates that the fracture zone volume is 1.7-1.8 times that of the caving zone in the 31108 working face, and the failure volume of the overlying strata is 9-10 times that of the stope. The simulated height of the caving zone and fracture zone is 9 m and 20 m, respectively. The comparison between the numerical and field measurement results demonstrates that the new evaluation method using shear-tensile strain behaviors can accurately predicate the height of the two zones. The proposed numerical method could be a viable alternative approach to two zones height calculation.


---

## Body

## 1. Introduction

Western China is rich in coal resources, accounting for 65% of the country’s total reserves, but water resources reserves only account for 3.9% [1, 2]. The lack of water resources seriously limits the sustainable development of energy in Western China. Considering the serious contradiction between coal resources development and water resources protection in the western region, Qian et al. [3] presented the concept of green coal mining and water conservation mining technology in the western mining area. Accordingly, underground reservoir technology of coal mine was first proposed by Shenhua Group. In this way, water storage in goaf formed after coal mining has become an important technical way to handle the contradiction between coal mining and water resources protection in Western China. In the process of coal mining, the caving zone, fracture zone, and bending subsidence zone are formed from bottom to top with the movement of the overlying strata [4–6]. The water storage capacity of the underground reservoir is directly determined by the failure mode and height of the caving zone and fracture zone [7]. Therefore, it is urgent to understand the failure mode and height of the overlying strata during coal seam mining.To date, many studies have been conducted to investigate the failure law of the overlying strata during underground reservoir construction. For example, Guo et al. [8] investigated the process of overburden failure transmission and further proposed a method for predicting overburden failure height. Xu et al. [9] presented the development law of a water flowing fractured zone with the breaking of the key stratum based on the key stratum theory. Huang et al. [10] proposed a method to characterize the distribution of the fractured zone by using the tensile deformation rate of the overlying strata. Yang et al. [11] performed a field observation to reveal the development law of the water flowing fractured zone of the soft overburden above a thick coal seam by monitoring borehole water injection leakage. Zhang et al. [12] revealed the influence of mining height and the panel layout on the height of the two zones in the soft overlying strata using discrete element numerical simulation. The above research studies provide a better understanding of the movement laws of overlying strata.Currently, theoretical analysis, field measurement, and numerical modelling are the most used method to predict the development height of the water flowing fractured zone. There into, theoretical calculation methods are mostly based on various assumptions with low accuracy, while field measurement improves the accuracy of the measurement results, but the test process is time-consuming and difficult. Therefore, numerical simulation has become a simple and effective method for analyzing the failure characteristics of overlying strata [13, 14]. In addition, numerical modelling can involve various influence factors in the analysis. Relevant studies have shown that the goaf compaction effect can result in characteristics changes of the overlying strata movement by providing an additional support resistance to the roof strata [15, 16]. Nevertheless, existing studies have rarely taken this influence into account. Therefore, it is urgent to conduct numerical research on the movement laws of overlying strata considering the compaction effect of goaf.In this study, a numerical method by using FLAC3D is presented to investigate the range of caving and fracture zones by discerning the plastic strain value of the zones in the model. In this modelling method, the compaction effect of the goaf is achieved by using the double-yield constitution model. The test site is located in the Lijiahao coal mine, Ordos City, Inner Mongolia Autonomous Region, China. The results presented in this study further improve the accuracy of numerical simulation, which is of great help to calculate the underground reservoir capacity.


## 2. Project Overview

The Lijiahao coal mine is located in Dongsheng District, Ordos City, Inner Mongolia Autonomous Region, China, with a designed production capacity of six million tons yearly. Currently, the mining seam is 3–1# coal seam. The average overburden depth and thickness of the coal seam are 201–254 m and 2.5–6.3 m, respectively. The coal seam is a relatively simple structure, and the immediate roof is sandy mudstone, with an average thickness of 3.72 m. The main roof is fine sandstone with an average thickness of 5.00 m, and the floor is mainly sandstone. The test site is located in the 31108 working face. Figure1 shows the coal seam histogram at this position.Figure 1
Generalized stratigraphy column of 31108 working face.


## 3. Goaf Compaction Theory and Its Mechanical Parameter Deduction Method

### 3.1. Mechanical Properties of the Collapsed and Compacted Rock Mass

A large number of field observation data on strata movement has demonstrated that overlying strata can be divided into caving, fracture, and bending subsidence zone from bottom to top after the coal seam is excavated. The collapsed rock mass is transformed from a loose body to a supporting body under the compaction action induced by the gravity and overlying strata movement [17–19]. The fractures development and stress distribution of the overlying strata is significantly affected by the strain hardening behavior of the collapsed rock mass. Salamon [20] proposed the stress-strain relationship of rock mass in the caving zone by taking the caving block as a granular material. The expression is shown in the following equation:(1)σcap=E0ε1−ε/εmax,where σcap, ε, εmax, and E0 are the vertical stress of the collapsed rock mass, the strain of collapsed rock mass under the action of vertical stress, the maximum vertical strain, and the initial elastic modulus of the collapsed rock mass, respectively. εmax depends on the crushing expansion coefficient Kp of the collapsed rock mass, and the expression is as follows [21]:(2)εmax=Kp−1Kp.The relationship between the initial modulusE0, crushing expansion coefficient Kp, and compressive strength of the collapsed rock mass σc is expressed as follows [22]:(3)E0=10.39σc1.042Kp7.7.Figure2 illustrates the effect of the crushing expansion coefficient and uniaxial compressive strength on the initial modulus of a collapsed rock mass. The initial modulus of the collapsed rock mass decreases exponentially with the increase in the crushing expansion coefficient but decreases gently with the uniaxial compressive strength decreasing. The initial modulus of the collapsed rock mass decreases greatly when the crushing expansion coefficient is less than 1.6; conversely, it decreases gently when the crushing expansion coefficient is larger than 1.6.Figure 2
Effect of the crushing expansion coefficient and uniaxial compressive strength on the initial modulus of the collapsed rock mass.From equations (1)–(3), the expression of the mechanical properties of the collapsed rock is deduced, as shown in the following equation.(4)σcap=10.39σc1.042Kp7.7ε1−εKp/Kp−1.The relationship of the stress-strain of the collapsed rock mass with different crushing expansion coefficients and the compressive strength is illustrated in Figure3. With the strain increasing, the vertical stress of the collapsed rock mass increases slowly; and then, it increases rapidly at an exponential function when the strain reaches a certain critical value. When the crushing expansion coefficient is equal, the compressive stiffness of the collapsed rock mass increases with the increase of its uniaxial compressive strength. In other words, the strain value of the low-strength rock mass is larger than that of the high-strength rock mass with the same stress value [23, 24].Figure 3
Effect of the bulking coefficient and strength of the rock mass on the mechanical property of collapsed rock masses.


### 3.2. Determination of Mechanical Parameters of the Caving Zone

In this study, the double-yield constitutive model in FLAC3D is employed to depict the stress recovery behavior of collapsed rock mass, which considers the permanent volumetric strain under isotropic pressure [25, 26]. Two kinds of parameters of cap pressure and material properties are needed in the double-yield constitutive model. And the material properties include density, bulk modulus, shear modulus, internal friction angle, and dilatancy angle. The cap pressure can be calculated by a theoretical equation, while the material parameters can be determined by the “trial-and-error” method, as shown in Figure 4.Figure 4
Back analysis of the parameters of the double-yield constitutive model.For the 31108 working face, the bulking coefficient is 1.3, and the uniaxial compression strength is 41 MPa. Based on equations (2)-(3), the maximum strain and the initial modulus of the gob materials can be estimated as 0.23 and 66 MPa, respectively. Furthermore, the cap pressures for the double-yield model are obtained and plotted in Figure 5.Figure 5
Numerical simulation inversion of mechanical properties of the double-yield model.A trial-and-error method was employed to determine the gob materials parameters by matching; the stress-strain curve is obtained by numerical modelling to that found by equation (1) [16]. For this purpose, a single-element submodel with dimensions 1 m × 1 m × 1 m was generated. A constant velocity was applied to the top of the model in the negative z-direction to generate vertical loading on the model. The velocity magnitude was set at 10–5 m/s. The displacement of the four vertical planes of the model was restricted in the normal direction, and a zero vertical displacement condition was set at the base of the model. The input parameters were calibrated by an iterative change in the bulk modulus, shear modulus, the angle of dilation, and the angle of friction of the gob materials. The final material parameters are listed in Table 1.Table 1
Material parameters for the double-yield model.
Density (kg/m3)Bulk modulus (GPa)Shear modulus (GPa)Internal friction angle (degree)Dilatancy angle (degree)10009.585.323015


## 3.1. Mechanical Properties of the Collapsed and Compacted Rock Mass

A large number of field observation data on strata movement has demonstrated that overlying strata can be divided into caving, fracture, and bending subsidence zone from bottom to top after the coal seam is excavated. The collapsed rock mass is transformed from a loose body to a supporting body under the compaction action induced by the gravity and overlying strata movement [17–19]. The fractures development and stress distribution of the overlying strata is significantly affected by the strain hardening behavior of the collapsed rock mass. Salamon [20] proposed the stress-strain relationship of rock mass in the caving zone by taking the caving block as a granular material. The expression is shown in the following equation:(1)σcap=E0ε1−ε/εmax,where σcap, ε, εmax, and E0 are the vertical stress of the collapsed rock mass, the strain of collapsed rock mass under the action of vertical stress, the maximum vertical strain, and the initial elastic modulus of the collapsed rock mass, respectively. εmax depends on the crushing expansion coefficient Kp of the collapsed rock mass, and the expression is as follows [21]:(2)εmax=Kp−1Kp.The relationship between the initial modulusE0, crushing expansion coefficient Kp, and compressive strength of the collapsed rock mass σc is expressed as follows [22]:(3)E0=10.39σc1.042Kp7.7.Figure2 illustrates the effect of the crushing expansion coefficient and uniaxial compressive strength on the initial modulus of a collapsed rock mass. The initial modulus of the collapsed rock mass decreases exponentially with the increase in the crushing expansion coefficient but decreases gently with the uniaxial compressive strength decreasing. The initial modulus of the collapsed rock mass decreases greatly when the crushing expansion coefficient is less than 1.6; conversely, it decreases gently when the crushing expansion coefficient is larger than 1.6.Figure 2
Effect of the crushing expansion coefficient and uniaxial compressive strength on the initial modulus of the collapsed rock mass.From equations (1)–(3), the expression of the mechanical properties of the collapsed rock is deduced, as shown in the following equation.(4)σcap=10.39σc1.042Kp7.7ε1−εKp/Kp−1.The relationship of the stress-strain of the collapsed rock mass with different crushing expansion coefficients and the compressive strength is illustrated in Figure3. With the strain increasing, the vertical stress of the collapsed rock mass increases slowly; and then, it increases rapidly at an exponential function when the strain reaches a certain critical value. When the crushing expansion coefficient is equal, the compressive stiffness of the collapsed rock mass increases with the increase of its uniaxial compressive strength. In other words, the strain value of the low-strength rock mass is larger than that of the high-strength rock mass with the same stress value [23, 24].Figure 3
Effect of the bulking coefficient and strength of the rock mass on the mechanical property of collapsed rock masses.


## 3.2. Determination of Mechanical Parameters of the Caving Zone

In this study, the double-yield constitutive model in FLAC3D is employed to depict the stress recovery behavior of collapsed rock mass, which considers the permanent volumetric strain under isotropic pressure [25, 26]. Two kinds of parameters of cap pressure and material properties are needed in the double-yield constitutive model. And the material properties include density, bulk modulus, shear modulus, internal friction angle, and dilatancy angle. The cap pressure can be calculated by a theoretical equation, while the material parameters can be determined by the “trial-and-error” method, as shown in Figure 4.Figure 4
Back analysis of the parameters of the double-yield constitutive model.For the 31108 working face, the bulking coefficient is 1.3, and the uniaxial compression strength is 41 MPa. Based on equations (2)-(3), the maximum strain and the initial modulus of the gob materials can be estimated as 0.23 and 66 MPa, respectively. Furthermore, the cap pressures for the double-yield model are obtained and plotted in Figure 5.Figure 5
Numerical simulation inversion of mechanical properties of the double-yield model.A trial-and-error method was employed to determine the gob materials parameters by matching; the stress-strain curve is obtained by numerical modelling to that found by equation (1) [16]. For this purpose, a single-element submodel with dimensions 1 m × 1 m × 1 m was generated. A constant velocity was applied to the top of the model in the negative z-direction to generate vertical loading on the model. The velocity magnitude was set at 10–5 m/s. The displacement of the four vertical planes of the model was restricted in the normal direction, and a zero vertical displacement condition was set at the base of the model. The input parameters were calibrated by an iterative change in the bulk modulus, shear modulus, the angle of dilation, and the angle of friction of the gob materials. The final material parameters are listed in Table 1.Table 1
Material parameters for the double-yield model.
Density (kg/m3)Bulk modulus (GPa)Shear modulus (GPa)Internal friction angle (degree)Dilatancy angle (degree)10009.585.323015


## 4. Development Law Analysis of Water Flowing Fractured Zone in Overlying Strata

As we all know, there are many bedding, discontinuities, and joints in coal or rock mass, which affect their failure characteristics. Although the finite difference method using FLAC3D is incapable of truly capturing their characteristics, coal or rock mass properties that are upscaled from intact cores properties using the Hoek–Brown criterion can be performed in the model. In other words, the coal and rock properties in our model are effective properties that accounted for the rock discontinuities. Therefore, in this section, we use FLAC3D software to analyze the height of the caving zone and fracture zone in the 31108 working face [2, 27].


### 4.1. Numerical Simulation Scheme

#### 4.1.1. Establishment of the Numerical Model

Based on the geological conditions of 31108 working face, a three-dimensional numerical calculation model was established, as shown in Figure6. The length in the x, y, and z-axis of the model is 200 m, 100 m, and 100 m respectively. The thicknesses of the coal seam floor, coal seam, and overlying strata are 30 m, 3 m, and 67 m, respectively. The horizontal displacements of the four vertical planes of the model were restricted in the normal direction, and the vertical displacement at the base of the model was set to zero. Vertical stress of 3.25 MPa was applied at the top model boundary to simulate an overburden pressure.Figure 6
Three-dimension simulation model.The double-yield model was used for the gob modelling, and the strain-softening Mohr–Coulomb model was used for the coal and rock mass modelling. The mechanical properties for the double-yield model and strain-softening model are listed in Tables1 and 2, respectively. Note that the mechanical properties of the coal and rock mass are calculated by the RockLab software, which is based on the generalized Hoek–Brown strength criterion [27–32].Table 2
Rock mass mechanical parameters of the numerical model.
Rock stratumDensity (kg·m3)Bulk modulus (GPa)Shear modulus (GPa)Cohesion (MPa)Internal friction angle (degree)Tensile strength (MPa)Residual cohesion (MPa)Residual friction angle (degree)Softening plastic strain (%)Floor26005.564.172.00351.200.20250.01Coal seam14001.1901.000.80200.500.08250.01Mudstone22003.472.641.50280.690.15250.01Sandy mudstone22002.031.561.20250.800.12200.01Fine sandstone26004.563.171.22260.900.12250.01


#### 4.1.2. Simulation Scheme and Process

Previous studies have shown that tensile failure mainly occurs in the caving zone, while shear failure mainly occurs in the fracture zone [33]. Accordingly, during the numerical simulation process, we distinguish the caving and fracture zone by the tensile and shear strain value by a self-developed FISH language. If the plastic tensile strain value of the zone is larger than 0.05, it is judged as a caving zone, while it is judged as a rupture zone if the shear strain value is larger than 0.05.The modelling process is as follows. After the initial balance of the numerical model, preliminary excavation is conducted; then, the range of the caving and fracture zones is determined by the plastic tensile and shear strain values; meanwhile, the caving zone is defined as a double-yield model and calculated to equilibrium. The above modelling process is conducted circularly until the working face is mined out.


#### 4.1.3. Analysis of Simulation Results

Figure7 shows a three-dimensional evolution diagram of the failure zone during the working face retreating, and Figure 8 shows the volume statistical analysis diagram of failed zones in two zones.Figure 7
Three-dimensional evolution diagram at various retreating distances. (a) 25 m. (b) 50 m. (c) 75 m. (d) 100 m.
(a)(b)(c)(d)Figure 8
Statistical analysis chart of two zones volume.It can be inferred in Figures7 and 8 that when the panel retreated 25 m, 50 m, 75 m, and 100 m, the volume of failed zones is 6.6 times, 7.59 times, 9.11 times, and 9.89 times that of the stope. Meanwhile, with the working face retreating, the volume of the failed zones increases quickly first and then tends to be stable. And the ultimate volume ratio of the failed zones and the stope maintains at 9-10. When the panel retreated 25 m, 50 m, 75 m, and 100 m, the volume ratio of the fracture zone to the caving zone is 1.14, 1.74, 1.72, and 1.80, respectively. It can be observed clearly that the volume ratio changes tend to be stable with the working face retreating and finally maintains 1.7-1.8 when the working face retreats to the “square” stage (the panel retreat length equates to the dip length).Figure9 is the height evolution law of the two zones at various retreating distances. It can be observed that the height evolution law of the caving and fracture zones kept a similar tendency, that is, it first increases rapidly and then tends to be stable with the working face retreating. When the working face retreats 25 m, the heights of the caving and fracture zones are 6 m and 10 m, respectively. And when the working face retreats 50 m (square stage), the height of the caving and fracture zones are 8 m and 16 m, increased by 33% and 60%, respectively. And furthermore, the increase rate slows down when the working face retreats to 75 m and 100 m, and the ultimate height of the caving and fracture zones are 9 m and 20 m, increased by 12.5% and 25%, respectively.Figure 9
Height evolution law of the two zones at various retreating distances. (a) 25 m. (b) 50 m. (c) 75 m. (d) 100 m.
(a)(b)(c)(d)


### 4.2. Theoretical Calculation and Field Measurement

In the longwall mining process, the height of the two zones reflects the movement characteristics of the overlying strata. Therefore, it is helpful to master the height of the two zones for the understanding of the movement law of overlying strata.Based on the field measured data, the calculation method of the heights of two zones under different geological conditions is deduced by Bai et al. [34]. The obtained formula is as follows:(5)Hc=100hc1h+c2,Hf=100hc3h+c4,where σ1 is the uniaxial compressive strength; Hc is the height of the caving zone; Hf is the height of the fracture zone; c1, c2, c3, and c4 are the rock strength coefficients, as given in Table 3.Table 3
Strength coefficient of strata in caving and fractured zones.
Rock typeσ1c1c2c3c4Hard＞402.1161.22Medium20–404.7191.63.6Weak＜206.2323.15The mining height of the 31108 working face is 2.9 m. According to the uniaxial compressive test of the roof strata, the average strength of strata in the roof is 15 MPa, which belongs to the weak rock layer. Substituting the parameters into equation (5), the heights of the caving and fracture zones are calculated to 8.5 m and 20.7 m, respectively.Studies have shown that field observation of borehole water injection leakage can obtain the height of two zones accurately [11]. In the field, a test of borehole water injection leakage was performed in a nearby working face, which has similar geological conditions to 31108. The observation results show that the height of the caving and fracture zone is 11 m and 21.5 m, respectively.


### 4.3. Result from Analysis

The heights of the caving zone and fracture zone obtained by numerical simulation, theoretical calculation, and field measurement are given in Table4. It can be observed clearly that the results obtained by simulation are less than that by the field measurement. This may be attributed to the great simplification in terms of geological conditions in the numerical model. However, within the allowable error range, the measured results and simulation results are almost the same.Table 4
Comparison of numerical simulation and theoretical calculation results.
Caving zone (m)Fracture zone (m)Total (m)Simulation results92029Theoretical results8.520.729.2Field measurement results1121.532.5It is well known that the failure mode and height of the strata overlying the stope are affected by various factors in the field, such as mining height, retreating rate, and buried depth. This study is only based on a specific geological condition of Lijiahao coal mine, and more studies are needed to prove the reliability of the numerical calculation method. However, the modelling procedure and evaluation method presented in this study are helpful to the other coal mines.


## 4.1. Numerical Simulation Scheme

### 4.1.1. Establishment of the Numerical Model

Based on the geological conditions of 31108 working face, a three-dimensional numerical calculation model was established, as shown in Figure6. The length in the x, y, and z-axis of the model is 200 m, 100 m, and 100 m respectively. The thicknesses of the coal seam floor, coal seam, and overlying strata are 30 m, 3 m, and 67 m, respectively. The horizontal displacements of the four vertical planes of the model were restricted in the normal direction, and the vertical displacement at the base of the model was set to zero. Vertical stress of 3.25 MPa was applied at the top model boundary to simulate an overburden pressure.Figure 6
Three-dimension simulation model.The double-yield model was used for the gob modelling, and the strain-softening Mohr–Coulomb model was used for the coal and rock mass modelling. The mechanical properties for the double-yield model and strain-softening model are listed in Tables1 and 2, respectively. Note that the mechanical properties of the coal and rock mass are calculated by the RockLab software, which is based on the generalized Hoek–Brown strength criterion [27–32].Table 2
Rock mass mechanical parameters of the numerical model.
Rock stratumDensity (kg·m3)Bulk modulus (GPa)Shear modulus (GPa)Cohesion (MPa)Internal friction angle (degree)Tensile strength (MPa)Residual cohesion (MPa)Residual friction angle (degree)Softening plastic strain (%)Floor26005.564.172.00351.200.20250.01Coal seam14001.1901.000.80200.500.08250.01Mudstone22003.472.641.50280.690.15250.01Sandy mudstone22002.031.561.20250.800.12200.01Fine sandstone26004.563.171.22260.900.12250.01


### 4.1.2. Simulation Scheme and Process

Previous studies have shown that tensile failure mainly occurs in the caving zone, while shear failure mainly occurs in the fracture zone [33]. Accordingly, during the numerical simulation process, we distinguish the caving and fracture zone by the tensile and shear strain value by a self-developed FISH language. If the plastic tensile strain value of the zone is larger than 0.05, it is judged as a caving zone, while it is judged as a rupture zone if the shear strain value is larger than 0.05.The modelling process is as follows. After the initial balance of the numerical model, preliminary excavation is conducted; then, the range of the caving and fracture zones is determined by the plastic tensile and shear strain values; meanwhile, the caving zone is defined as a double-yield model and calculated to equilibrium. The above modelling process is conducted circularly until the working face is mined out.


### 4.1.3. Analysis of Simulation Results

Figure7 shows a three-dimensional evolution diagram of the failure zone during the working face retreating, and Figure 8 shows the volume statistical analysis diagram of failed zones in two zones.Figure 7
Three-dimensional evolution diagram at various retreating distances. (a) 25 m. (b) 50 m. (c) 75 m. (d) 100 m.
(a)(b)(c)(d)Figure 8
Statistical analysis chart of two zones volume.It can be inferred in Figures7 and 8 that when the panel retreated 25 m, 50 m, 75 m, and 100 m, the volume of failed zones is 6.6 times, 7.59 times, 9.11 times, and 9.89 times that of the stope. Meanwhile, with the working face retreating, the volume of the failed zones increases quickly first and then tends to be stable. And the ultimate volume ratio of the failed zones and the stope maintains at 9-10. When the panel retreated 25 m, 50 m, 75 m, and 100 m, the volume ratio of the fracture zone to the caving zone is 1.14, 1.74, 1.72, and 1.80, respectively. It can be observed clearly that the volume ratio changes tend to be stable with the working face retreating and finally maintains 1.7-1.8 when the working face retreats to the “square” stage (the panel retreat length equates to the dip length).Figure9 is the height evolution law of the two zones at various retreating distances. It can be observed that the height evolution law of the caving and fracture zones kept a similar tendency, that is, it first increases rapidly and then tends to be stable with the working face retreating. When the working face retreats 25 m, the heights of the caving and fracture zones are 6 m and 10 m, respectively. And when the working face retreats 50 m (square stage), the height of the caving and fracture zones are 8 m and 16 m, increased by 33% and 60%, respectively. And furthermore, the increase rate slows down when the working face retreats to 75 m and 100 m, and the ultimate height of the caving and fracture zones are 9 m and 20 m, increased by 12.5% and 25%, respectively.Figure 9
Height evolution law of the two zones at various retreating distances. (a) 25 m. (b) 50 m. (c) 75 m. (d) 100 m.
(a)(b)(c)(d)


## 4.1.1. Establishment of the Numerical Model

Based on the geological conditions of 31108 working face, a three-dimensional numerical calculation model was established, as shown in Figure6. The length in the x, y, and z-axis of the model is 200 m, 100 m, and 100 m respectively. The thicknesses of the coal seam floor, coal seam, and overlying strata are 30 m, 3 m, and 67 m, respectively. The horizontal displacements of the four vertical planes of the model were restricted in the normal direction, and the vertical displacement at the base of the model was set to zero. Vertical stress of 3.25 MPa was applied at the top model boundary to simulate an overburden pressure.Figure 6
Three-dimension simulation model.The double-yield model was used for the gob modelling, and the strain-softening Mohr–Coulomb model was used for the coal and rock mass modelling. The mechanical properties for the double-yield model and strain-softening model are listed in Tables1 and 2, respectively. Note that the mechanical properties of the coal and rock mass are calculated by the RockLab software, which is based on the generalized Hoek–Brown strength criterion [27–32].Table 2
Rock mass mechanical parameters of the numerical model.
Rock stratumDensity (kg·m3)Bulk modulus (GPa)Shear modulus (GPa)Cohesion (MPa)Internal friction angle (degree)Tensile strength (MPa)Residual cohesion (MPa)Residual friction angle (degree)Softening plastic strain (%)Floor26005.564.172.00351.200.20250.01Coal seam14001.1901.000.80200.500.08250.01Mudstone22003.472.641.50280.690.15250.01Sandy mudstone22002.031.561.20250.800.12200.01Fine sandstone26004.563.171.22260.900.12250.01


## 4.1.2. Simulation Scheme and Process

Previous studies have shown that tensile failure mainly occurs in the caving zone, while shear failure mainly occurs in the fracture zone [33]. Accordingly, during the numerical simulation process, we distinguish the caving and fracture zone by the tensile and shear strain value by a self-developed FISH language. If the plastic tensile strain value of the zone is larger than 0.05, it is judged as a caving zone, while it is judged as a rupture zone if the shear strain value is larger than 0.05.The modelling process is as follows. After the initial balance of the numerical model, preliminary excavation is conducted; then, the range of the caving and fracture zones is determined by the plastic tensile and shear strain values; meanwhile, the caving zone is defined as a double-yield model and calculated to equilibrium. The above modelling process is conducted circularly until the working face is mined out.


## 4.1.3. Analysis of Simulation Results

Figure7 shows a three-dimensional evolution diagram of the failure zone during the working face retreating, and Figure 8 shows the volume statistical analysis diagram of failed zones in two zones.Figure 7
Three-dimensional evolution diagram at various retreating distances. (a) 25 m. (b) 50 m. (c) 75 m. (d) 100 m.
(a)(b)(c)(d)Figure 8
Statistical analysis chart of two zones volume.It can be inferred in Figures7 and 8 that when the panel retreated 25 m, 50 m, 75 m, and 100 m, the volume of failed zones is 6.6 times, 7.59 times, 9.11 times, and 9.89 times that of the stope. Meanwhile, with the working face retreating, the volume of the failed zones increases quickly first and then tends to be stable. And the ultimate volume ratio of the failed zones and the stope maintains at 9-10. When the panel retreated 25 m, 50 m, 75 m, and 100 m, the volume ratio of the fracture zone to the caving zone is 1.14, 1.74, 1.72, and 1.80, respectively. It can be observed clearly that the volume ratio changes tend to be stable with the working face retreating and finally maintains 1.7-1.8 when the working face retreats to the “square” stage (the panel retreat length equates to the dip length).Figure9 is the height evolution law of the two zones at various retreating distances. It can be observed that the height evolution law of the caving and fracture zones kept a similar tendency, that is, it first increases rapidly and then tends to be stable with the working face retreating. When the working face retreats 25 m, the heights of the caving and fracture zones are 6 m and 10 m, respectively. And when the working face retreats 50 m (square stage), the height of the caving and fracture zones are 8 m and 16 m, increased by 33% and 60%, respectively. And furthermore, the increase rate slows down when the working face retreats to 75 m and 100 m, and the ultimate height of the caving and fracture zones are 9 m and 20 m, increased by 12.5% and 25%, respectively.Figure 9
Height evolution law of the two zones at various retreating distances. (a) 25 m. (b) 50 m. (c) 75 m. (d) 100 m.
(a)(b)(c)(d)


## 4.2. Theoretical Calculation and Field Measurement

In the longwall mining process, the height of the two zones reflects the movement characteristics of the overlying strata. Therefore, it is helpful to master the height of the two zones for the understanding of the movement law of overlying strata.Based on the field measured data, the calculation method of the heights of two zones under different geological conditions is deduced by Bai et al. [34]. The obtained formula is as follows:(5)Hc=100hc1h+c2,Hf=100hc3h+c4,where σ1 is the uniaxial compressive strength; Hc is the height of the caving zone; Hf is the height of the fracture zone; c1, c2, c3, and c4 are the rock strength coefficients, as given in Table 3.Table 3
Strength coefficient of strata in caving and fractured zones.
Rock typeσ1c1c2c3c4Hard＞402.1161.22Medium20–404.7191.63.6Weak＜206.2323.15The mining height of the 31108 working face is 2.9 m. According to the uniaxial compressive test of the roof strata, the average strength of strata in the roof is 15 MPa, which belongs to the weak rock layer. Substituting the parameters into equation (5), the heights of the caving and fracture zones are calculated to 8.5 m and 20.7 m, respectively.Studies have shown that field observation of borehole water injection leakage can obtain the height of two zones accurately [11]. In the field, a test of borehole water injection leakage was performed in a nearby working face, which has similar geological conditions to 31108. The observation results show that the height of the caving and fracture zone is 11 m and 21.5 m, respectively.


## 4.3. Result from Analysis

The heights of the caving zone and fracture zone obtained by numerical simulation, theoretical calculation, and field measurement are given in Table4. It can be observed clearly that the results obtained by simulation are less than that by the field measurement. This may be attributed to the great simplification in terms of geological conditions in the numerical model. However, within the allowable error range, the measured results and simulation results are almost the same.Table 4
Comparison of numerical simulation and theoretical calculation results.
Caving zone (m)Fracture zone (m)Total (m)Simulation results92029Theoretical results8.520.729.2Field measurement results1121.532.5It is well known that the failure mode and height of the strata overlying the stope are affected by various factors in the field, such as mining height, retreating rate, and buried depth. This study is only based on a specific geological condition of Lijiahao coal mine, and more studies are needed to prove the reliability of the numerical calculation method. However, the modelling procedure and evaluation method presented in this study are helpful to the other coal mines.


## 5. Conclusions

(1)
Considering the goaf compaction effect, a new numerical simulation method by using FLAC3D is proposed to judge the fracture zone and caving zone of the overlying strata. The method can accurately predicate the height of two zones.(2)
With the strain increasing, the vertical stress of the collapsed rock mass increases slowly; and then, it increases rapidly at an exponential function when the strain reaches a certain critical value. The critical strain increases with the increase in the crushing expansion coefficient.(3)
The simulation results show that the volume of the fracture zone in the 31108 working face is 1.7-1.8 times that of the caving zone, while the rock failure is 9-10 times that of the stope. The simulated height of the caving zone and fracture zone is 9 m and 20 m, respectively.



---

*Source: 6071957-2021-09-03.xml*