*4.3. Determination of Dynamic Excavation Parameters for the Numerical Model*

In the discrete element numerical simulation software, the calculation of the model is mainly controlled by two different methods [37]. The first method is by setting the number of the software's operation cycle steps, and stopping the operation before moving on to the next stage. Moreover, the operation will stop or move on to the next stage if the maximum unbalance force ratio (the ratio of the maximum unbalanced force to the initial unbalanced force) is less than the default value even if the number of calculation steps has not reached the set cycle steps. The second method is by setting the maximum unbalanced force ratio (MUFR). When the maximum unbalanced force ratio calculated by the software is less than the aforementioned set value, the operation will stop or move on to the next stage. The calculation principle is shown in Figure 8.

**Figure 8.** The procedure of continuous excavation in the numerical modeling.

For the dynamic representation of coal seam excavation, a certain excavation step is generally set, and a certain cycle step is given, or the value of MUFR is changed. When the calculation is completed, the model enters the next excavation operation. However, in the continuous excavation, the deformation and failure of the overlying strata are a constantly changing process. It is necessary to set appropriate excavation steps and cycle steps (or MUFR) to accurately simulate an actual mining operation. Therefore, it is necessary to check the parameters used in the numerical simulation through relevant data measured on site, in order that the numerical simulation can be used as a reliable means to study the movement of the surface and overlying strata.

#### (1) Determination of excavation step

Kbj-60 III rock pressure observation system is installed in Peigou Coal Mine to accurately grasp the ground stress and other data of roof strata. During the mining, the monitoring system continuously logs the stress of the hydraulic support for the underground working face in the support extension. The extension data are gathered by specialized personnel and uploaded to a specialized computer software. Based on the statistical data of the mine pressure monitoring system during the mining, it is known that the average periodic stress step distance of the panel is 18.5 m, and the numerical simulation excavation step distance is determined to be 20 m.

(2) Setting of MUFR

According to the schematic diagram of the dynamic calculation scheme in Figure 9, an initial value is first set, and the numerical simulation calculation is compared with the measured surface subsidence curve. Then, the maximum unbalanced force ratio is continuously adjusted in order that the numerical simulation results are consistent with the field-measured results. It is found that when the excavation step is 20 m, the maximum unbalanced force ratio is 4.8 e−5. After the numerical simulation of dynamic excavation, the dynamic subsidence data of the surface during the mining process are obtained, and the subsidence curve of survey station *Z*8 and the corresponding numerical model positions are drawn, as shown in Figure 10. The subsidence curves of the survey stations on the surface and the corresponding points in the numerical model are very similar. The relationship curve between the survey station subsidence and mining distance presents an asymmetric

"s" shape with a small front end and a large rear end. The two curves can show three stages of dynamic development of the surface.

**Figure 9.** The dynamic calculation scheme.

**Figure 10.** The subsidence velocity between the numerical modeling and the survey data.

Figure 11 illustrates a comparison of measured data and numerical simulation results of subsidence velocity. It can be noticed that when two adjacent survey stations with large measurement error are excluded with a mining distance of about 400 m, the curve of subsidence velocity obtained by measurement is similar to the one obtained from the numerical model in curve shape, the MSV on the curve are 124.32 mm/d and 132.78 mm/d, and the relative error rate is 6.4%. This demonstrates that the numerical model can accurately calculate the surface dynamic movement and can therefore be used to study the dynamic movement of the surface and overlying strata.

**Figure 11.** Surface subsidence velocity between the numerical modeling and the survey data.

#### *4.4. Dynamic Movement Characteristics of Overlying Strata*

(1) Variation of overlying strata MSV

Considering the surface subsidence curve of the numerical model as an example, the MSV of the surface is fitted according to the analysis method of the surface dynamic movement, as shown in Figure 4. The functional relationship between the surface MSV and the mining distance is obtained:

$$V\_{\text{max}} = 137.96(1 - e^{-0.0083d}) \tag{3}$$

The subsidence data of rock stratus are then selected to study the dynamic movement, which are 250 m, 210 m, 170 m, 130 m, 90 m, 50 m, and 10 m away from the coal seam. Due to its discreteness, the MSV of rock strata under different mining distances is fitted and defined by referral to the form of the relationship function between the MSV of the surface and the mining distance (Equation (3)):

$$V\_{\text{max}}^i = V'(1 - e^{-hd}) \tag{4}$$

In Equation (4), *Vi* max is the MSV of the strata with the *i* depth, mm/d; *V* 'is the stable value of the MSV, mm/d; and *b* is the growth coefficient.

By substituting the MSV of different depths of strata in the mining process into Equation (4), the parameter of the dynamic change function of the MSV of strata is obtained as shown in Table 2.


**Table 2.** The MSV parameter of the overlying strata.

The MSV of different depths of strata gradually increases with the mining distance. When the working face advances to a certain distance, the MSV of strata tends to have a fixed value. It can be seen from Table 2 that the MSV of the rock strata gradually increases

with the distance, which is between the rock strata and the coal seam. The MSV of the surface increases from 137.96 to 261.58 mm/d when the distance from the coal seam is 30 m. To express the relationship between the MSV after the rock strata is stabilized and the distance between the rock strata and the coal seam, the MSV of the rock strata in Table 2 is nonlinearly fitted. In addition, the fitting formula (Equation (5)) is obtained, as well as the variation curve of the MSV at various distances between the rock stratus and the coal seam (Figure 5):

$$V'\_{\text{max}}(z) = V\_0(\frac{H\_0 - \mathbf{z}}{H\_0})^{-0.30} \tag{5}$$

In Equation (5), *z* is the depth of the rock strata, m; *H*<sup>0</sup> is the depth of the coal seam, m; and *V*<sup>0</sup> is the MSV of the surface, mm/d, i.e., 132.75 mm/d.

It can be seen from Table 2 that the growth parameter *B* in the fitting parameter of the MSV of the rock strata increases with the decrease in the distance from the coal seam, indicating that when the MSV of the rock strata close to the coal seam is stable, the mining distance is small. According to the curve fitting of the above data, the relationship between the growth parameter *b* and the depth of rock strata *z* is obtained:

$$b(z) = 0.008(\frac{H\_0 - z}{H\_0})^{-0.12} \tag{6}$$

Substituting Equations (5) and (6) into Equation (4), the relationship between the MSV of rock strata and the mining distance is obtained:

$$V^i = V\_0(\frac{H\_0 - z}{H\_0})^{-0.30} \left[1 - e^{-0.008(\frac{H\_0 - z}{H\_0})}\right]^{-0.12} d\tag{7}$$

#### (2) Variation of overlying strata LDMSV

Due to the different distances between the rock strata and the coal seam, the mininginduced impact and the dynamic movement of the strata are also different. This paper selects the numerical calculation data with strata depths of 250 m, 210 m, 170 m, 130 m, 90 m, 50 m, and 10 m in order to study the variation of the LDMSV of rock strata in the mining process.

Similar to the form of Equation (4), *L'* is defined as the LDMSV of rock strata, and the parameter *c* represents the growth trend of the LDMSV with the mining distance, as shown in Table 3.

**Table 3.** The LDMSV parameter of the overlying strata.


When the distance between the rock strata and the coal seam gradually decreases, *L'* gradually decreases; namely, from 79.51 to 40.36 m when the distance between the rock strata and the coal seam is 30 m. However, *L'* of the rock strata has little change with the growth parameter *c*, which is between 0.011 and 0.015, and its average value can be considered as 0.013. Similarly, the stable value *L'*max of the lag distance from the MSV of the rock strata can be determined as shown in Equation (8). Therefore, the calculation formula of *L'* with different mining depths can be obtained (Equation (9)):

$$L'\_{\text{max}}(z) = 41.84 \frac{H\_0 - z}{H\_0} + 36.32 \tag{8}$$

$$L^i = (41.84 \frac{H\_0 - z}{H\_0} + 36.32)(1 - e^{-0.013d})\tag{9}$$

#### **5. Calculation Method of Subsidence Velocity**

#### *5.1. Calculation Method of Subsidence Velocity in Strike Major Section of Surface Subsidence Basin*

According to the variation of the subsidence velocity of the survey stations, the study in [13] shows that the distribution curve of the surface subsidence velocity is similar to the "quadratic curve" distribution. If the projection position on the surface at a certain time is considered as the coordinate origin, the direction of the working face on the strike is the positive direction of the *x* axis, and the subsidence velocity of the station is the *y* axis. Therefore, the relative relationship between the subsidence velocity curve of the ground surface and the working face position can be expressed. The subsidence velocity formula at one position on the strike major section is:

$$V(x) = \frac{V\_{\text{max}}}{1 + \left(\frac{x+L}{a}\right)^2} \tag{10}$$

In Equation (10), *x* is the distance to the working face, *a* is the morphological parameter, which indicates the steepness of the curve.

For the solution of the shape parameter *a*, this paper assumes that the *a* calculation point on the surface starts when it is significantly distant from the working face, and its subsidence tends to have a stable value when the working face pushes over significantly away from the point. The following results are obtained:

$$
\Delta W = \int\_{+\infty}^{-\infty} \frac{V\_{\text{max}}}{1 + \left(\frac{x+L}{a}\right)^2} d\left(-\frac{x}{c}\right) \tag{11}
$$

Then:

$$a = \frac{W\_{\text{max}}}{V\_{\text{max}}} \cdot \frac{c}{\pi} \tag{12}$$

In Equation (12), *W*max is the maximum surface subsidence and can be obtained by *<sup>W</sup>*max <sup>=</sup> *qm*√<sup>3</sup> *<sup>n</sup>*1*n*<sup>2</sup> cos *<sup>α</sup>*; *<sup>m</sup>* is the thinness of coal seam, *<sup>n</sup>*<sup>1</sup> and *<sup>n</sup>*<sup>2</sup> are respectively equal to the ratio of mining size and coal seam burial depth on the strike and incline direction, and c is the average mining velocity of the working face, m/d. Moreover, *V*max can be calculated according to Equation (1).

From Equation (12), it can be seen that when the mining velocity of the working face is constant, the value of the shape parameter *a* changes with the change in the maximum subsidence and maximum subsidence velocity values. Moreover, it can be judged that before reaching the supercritical mining, the shape parameter *a* and parameter *L* of the subsidence velocity curve will change with the increase in the goaf area, and the change law can be obtained from the relationship between the MSV, LDMSV, and mining distance. After a mining distance of 400 m, the maximum subsidence *W*max no longer increases due to supercritical mining in the strike direction of the working face, while the MSV of Equation (1) increases, but the increasing range gradually decreases. Therefore, it can be considered that the parameter *a* of the subsidence velocity curve increases slightly after supercritical mining, and the shape of subsidence velocity curve becomes steeper.

In summary, the variation of the subsidence velocity curve with the mining process can be described as follows: With the increase in the goaf area, the MSV gradually increases, the subsidence velocity curve shape gradually becomes steep, and the LDMSV gradually increases. When the working face reaches supercritical mining; namely, the working face mining distance exceeds 400 m, the subsidence velocity curve of the strike major section maintains a certain lag distance with the working face in a fixed form and moves forward with the mining. Figure 12 shows the relative position of the subsidence and subsidence velocity curves of the strike major section of the surface, when the working face mining distance is 200m, 400m, 600m, and 800 m, respectively.

**Figure 12.** Surface subsidence and subsidence velocity relative locations during mining.

The mining and geology parameters of the panel of Peigou Coal Mine are brought into Equations (1), (2), and (12), and the formula for calculating the subsidence velocity at any position on the strike major section at any time of the working face is obtained by bringing in Equation (10):

$$V(\mathbf{x},d) = \frac{50.40(1 - e^{-0.0028d})}{1 + (\frac{\mathbf{x} + 95.60(1 - e^{-0.11d})}{\sqrt{n\_1 n\_2} \cos \alpha \cdot c} \cdot 50.40(1 - e^{-0.0028d}) \cdot \pi)^2} \tag{13}$$

When the working face mining distance is equal to 360 m, the subsidence velocity calculation points on the strike major section of subsidence basin at different mining distances are obtained according to Equation (13) and compared with the measured values, as shown in Figure 13.

**Figure 13.** Surface subsidence velocity curve while *d* = 360 m.

It can be seen from Figure 13 that when the mining distance of the working face is 360 m, the deviation between the measured and calculated values of the subsidence velocity at each point is small, and the average error is 1.57 mm/d, indicating that the predicted results of dynamic subsidence velocity can meet engineering needs.

#### *5.2. Prediction Model of Overlying Strata Dynamic Movement*

For the correction of the LDMSV of rock strata, a correction factor *aL* can be added before Equation (9), and its value is equal to *L*s/*L*´(z)|*z* = *H*0. Therefore, the calculation formula of LDMSV of overlying strata based on the measured LDMSV of the surface is obtained:

$$L^i = a\_L (41.84 \frac{H\_0 - z}{H\_0} + 36.32)(1 - e^{-0.013d}) \tag{14}$$

Assuming that the subsidence velocity distribution curve on the major section of the surface and overlying strata is similar to the form of a quadratic curve during the working face mining, referring to the subsidence velocity function of the strike major section of the surface (Equation (10)), the subsidence velocity formula of the overlying strata is:

$$V^{i}(\mathbf{x}) = \frac{V^{i}\_{\text{max}}}{1 + \left(\frac{\mathbf{x} + L^{i}}{a^{i}}\right)^{2}} \tag{15}$$

According to Equation (11), the parameter *a* of the subsidence velocity curve will change with the increase in mining distance. When the working face moves from subcritical mining to supercritical mining, the shape of the subsidence velocity curve becomes steeper. The shape parameter *ai* of the distribution function of rock strata subsidence velocity is mainly calculated according to the maximum subsidence value *W*<sup>i</sup> MAX, the *V*<sup>i</sup> max (MSV), and the working face mining velocity *c* of the rock strike under this level. The maximum subsidence velocity *V*<sup>i</sup> max and the rock strata can be referred to Equation (5), and the maximum subsidence *W*<sup>i</sup> MAX of the rock strata can be determined according to the mining distance and the distance between the rock strata and the coal seam; namely:

$$\mathcal{W}\_{\text{max}}^{i} = M q\_i \sqrt[\theta]{n\_1^i n\_2^i} \cos \alpha \tag{16}$$

In Equation (16), for rock strata *i*, the mining coefficients *ni <sup>1</sup>* and *n<sup>i</sup> <sup>2</sup>* of the working face are equal to *ηD*/(*H*0-*z*), where *η* is a coefficient generally considered as 0.7~0.9, *D* is the penal length of incline or strike (m), and *z* is the depth of rock strata *i* (m).

Therefore, Equations (7), (14), and (16) are combined with Equation (12) and brought into Equation (15) to obtain the calculation formula of subsidence velocity of the surface and overlying strata on the strike major section at different mining distances:

$$V(\mathbf{x},d,z) = \frac{V\_{\mathbf{S}}(\frac{H\_{0}-z}{H\_{0}})^{-0.30}[1-e^{-0.008(\frac{H\_{0}-z}{H\_{0}})^{-0.12}}]}{1+\left\{\frac{\mathbf{x}+\mathbf{a}\_{1}(41.84\frac{H\_{0}-z}{H\_{0}}+36.32)(1-e^{-0.015})}{M\_{\mathbf{0}}\sqrt[3]{\mathbf{u}\_{1}^{\prime}\mathbf{u}\_{2}^{\prime}\cos\mathbf{a}\cdot\mathbf{c}}}\cdot V\_{\mathbf{S}}(\frac{H\_{0}-z}{H\_{0}})^{-0.30}[1-e^{-0.008(\frac{H\_{0}-z}{H\_{0}})^{-0.12}}]\cdot\pi\right\}^{2}}\tag{17}$$

According to the above analysis of the MSV and the LDMSV on the strike major section of the surface in Peigou Coal Mine, it is known that in Equation (17), *V*<sup>0</sup> is equal to 50.40 mm/d, and the correction coefficient *aL* of the LDMSV is equal to 1.22. Bringing the above parameters into Equation (17) results in the calculation formula of the subsidence velocity of the surface and overlying strata. When the mining velocity is 2.3 m/d and the mining time is selected as 20d, 40 d, 60 d, 100 d, and 140 d, the subsidence velocity distribution curves of the rock strata with the depths of 100 and 200 m, respectively, are shown in Figures 5–13.

It can be seen from Figure 14 that the subsidence velocity of rock strata shows the following laws: First, the peak value on the subsidence velocity curve gradually increases with the increase in mining distance, and when the working face mining reaches supercritical mining, the increase amplitude of the MSV gradually decreases, finally stabilizing at a certain value. In addition, the shape of the subsidence velocity curve of the rock strata becomes steeper with the increase in mining distance, and finally tends to have a stable shape. Finally, the LDMSV of the subsidence velocity curve increases gradually with the increase in mining distance. In general, the peak value of the rock subsidence velocity curve gradually increases, and the curve shape gradually becomes steep during the mining. When the mining distance reaches a certain value, the peak value and shape of the rock strata subsidence velocity curve tend to be stable and maintain a certain lag distance from the working face, moving forward with the mining of the working face. By comparing Figure 14a,b, it can be found that with the same mining distance, the mining-induced influence of the strata close to the coal seam is greater, resulting in a larger MSV of the rock strata, a smaller LDMSV, and a steeper shape of the subsidence velocity curve.

**Figure 14.** The time-dependent subsidence velocity of strata in different buried depths with mining.

The dynamic impact of mining on structures (shafts, roadways, and other structures) in overburden can be predicted and evaluated according to the variation rule of overlying strata dynamic subsidence and movement. During the mining process, the position and time of the most severe rock strata movement can be determined to reinforce, prevent in advance, and ensure the safe operation of underground mining.

#### **6. Conclusions and Discussion**

In this paper, the dynamic movement law and prediction method of the surface and overburden are studied with the stratum structure in Peigou Coal Mine. Since the dynamic movement in the overlying strata is not easy to be observed, this paper innovatively put forward a numerical simulation method to study the dynamic movement, which is used to establish the relationship between surface and overlying strata dynamic movement parameters. Therefore, the dynamic movement law and prediction method of overlying strata and surface are obtained. The research results can be used to predict the dynamic subsidence impact on the structures on the surface and in the overlying strata during the mining process, in order to obtain the mining time or position of working face, which induces the strongest impact on the these structures, and provides a reference for the protection of the structures. The key conclusions from this study are as follows:


Due to the complexity and diversity of stratum structure, the strata with different structures show different movement laws induced by underground mining, but entire or partial overburden strata movement after mining is controlled by strata structures of "key stratum". Due to the different locations of the key strata in the overlying strata, the propagation law of the movement and deformation in the rock stratum is different. Another important factor influencing the dynamic movement of overlying strata is the bulking factor of the roof strata in the goaf [38–43]. The stratum immediately above the seam bends downward and later breaks into different size blocks depending on their lithological nature. This phenomenon continues upward until the lowest uncaved stratum receives the support of the caved and bulked rock pile. During this process, the strata are broken into blocks by vertical and/or sub-vertical fractures and horizontal cracks due to the vertical bed separation. Vertical and/or sub-vertical fractures andhorizontal crackswill

be re-compacted by the stress induced by the weight of overlying strata, and the surface and rock stratum dynamic movement is closely related to the re-compaction of goaf. The above two aspects are the focus of further research.

**Author Contributions:** Writing—original draft preparation, G.X.; project administration, D.L.; resources, Y.Z.; data curation, H.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the regional first-class discipline of Ecology in Guizhou province (XKTJ [2020]22), National Natural Science Foundation of China (Grant No. 51974105), and Fund Project of Bijie Science and Technology Bureau (G [2019]1).

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The authors would like to acknowledge the support of the regional first-class discipline of Ecology in Guizhou province (XKTJ [2020] 22), National Natural Science Foundation of China (Grant No. 51974105), and Fund Project of Bijie Science and Technology Bureau (G [2019]1).

**Conflicts of Interest:** The authors declare no conflict of interest.
