#### Abstract

Strip mining with subsequent filling is one of the main mining methods for gently inclined thin ore bodies. The production process of excavating alternate strips is beneficial to the safety of mining. Reasonable stope structural parameters are key to achieving safe and efficient mining. The Tiaoshuihe phosphate mine of Sanning Mining is taken as an example in this study. Based on precision finite element modeling and simulation, a reasonable width range and the interval value of the strip are determined. A reliable and efficient strip width is calculated by using the probability model of the improved Mathews stability graph method. The results show the following. Firstly, under the geological conditions and backfill environment of the Tiaoshuihe phosphate mine, the reasonable and safe strip width interval is 7–9 m. Secondly, the probabilities of open stope stability with strip widths of 7, 8, and 9 m are 88.55%, 86.76%, and 84.94%, respectively. The reasonable probabilities of stope stability with strip widths of 7 and 8 m are higher than 85%. Thirdly, combining this with the drilling equipment operation parameters, it is suggested that the best strip width is 7 m without increasing the strength of the backfill.

#### 1. Introduction

At present, the filling mining method is becoming increasingly widely used. It is a safe and efficient method to apply strip mining with subsequent filling for gently inclined thin ore bodies [1]. The structural parameters of the stope are closely related to its productivity and recovery rate, as well as the stability of the goaf. Meanwhile, the mining equipment level and technology are required to consist with the structural parameters of the stope [2]. Gao-hui et al. [3] combined an engineering geological evaluation with a numerical simulation to determine the stope height, chamber width, and spike diameter in Baiyang ore block. Ning and Hu [4]. used Mathews stability graph method and an orthogonal numerical simulation to determine the limit exposure area of zinc-copper ore body of Dachang Tongkeng Mine. Li and Li [5] applied the 3D-sigma numerical simulation method to analyze the surrounding rock stability in different stope structure parameters and designed the reasonable size of room and pillar for a phosphate mine. Guo et al. [6] adopted the fuzzy comprehensive evaluation method to prove that the optimal stope structural parameter is 10 m × 10 m for safety and high-efficiency mine operation in Sanshandao Gold Mine. Chen et al. [7] optimized the structural parameters of Yongping Copper Mine with FLAC3D and analyzed plastic zone, displacement, and stress distribution after excavation. Li et al. [8] performed a series of three-dimensional simulation and studied the relationship between mining safety-rock movement parameters and the design parameters and then obtained the reliable structural parameters of the stope.

However, the traditional single computational analysis method has some limitations. For uncertain rock mechanics parameters, the probability of its stability has an important influence on the optimization of engineering parameters. In this study, with finite element numerical simulation determining the width range of the strip, the overall stability of the goaf is calculated by using the Mathews binary linear regression model, and the stability probability of goaf is quantized for different strip widths. Taking the parameters of underground equipment into consideration, a reasonable strip width for the Tiaoshuihe phosphate mining was determined by comparing the stability probability finally.

#### 2. Project Profile

The phosphate ore bodies of the Tiaoshuihe phosphate mine of Sanning Mining are divided into the middle phosphate layer (Ph^{2}) and the lower phosphate layer (Ph_{1}^{3}), and the average interval between the two layers is 10 m. The middle phosphate layer currently being mined is 1.63–6.79 m in thickness, with an average thickness of 3.22 m and an inclined angle of 15°. Its burial depth is 81.67–614.75 m, with an average burial depth of 254.27 m. The folds in the mine area are not well developed, while the fault structure in the area is relatively developed. Considering the limited space in the thin ore body, the DD211L low-type drilling rig manufactured by Sandvik was purchased and used in the mine, and the maximum drilling width of this rig is 7.6 m for a single-drilling application.

Full-thickness panel strip mining is adopted in this mine, with the sequence of excavating alternate strips. The waste rock is cemented with tailings filling into the open stope. At present, the method of exploiting a strip width of 5 m and leaving 5.5 m is used in the mine. To improve the production capacity of the stope and reduce the construction involved in the mining, the strip width of the stope is proposed to be enlarged, so as to enhance the recovery rate of the stope. The mechanical parameters of the rock body in the Tiaoshuihe phosphate mine are obtained by a standard reduction of the rock mechanical parameters tested in the laboratory [9], as shown in Table 1.

#### 3. Determination of Strip Width Range

##### 3.1. Accurate Finite Element Model and Simulation Schemes

FLAC3D is a type of simulation software developed by the ITASCA of America and is widely used in rock and soil mechanics calculation [10–12]. Since it comprises 11 kinds of original constitutive models of elastoplastic materials, such as seepage, static, creep, dynamic, and temperature calculation models, it is appropriate for many scientific research fields like geotechnical engineering, geological engineering, and mining engineering [13–15]. Furthermore, it can get the accurate data and graph of stress, displacement, and plastic zone after caving and filling. Hence, FLAC3D can meet the requirements of this paper.

To ensure the accuracy and authenticity of the finite element analysis, the ore body model has the same thickness as the average three-dimensional ore body extracted from the Tiaoshuihe phosphate mine in DIMINE software. Four model groups with strip widths of 6, 7, 8, and 9 m and the length of 50 m were established, in which the roadway widths were 4.5 m and the left and right sides of the roadway each contained five strips. Figure 1 shows the model with strip width of 9 m as an example, and ten strips are numbered 1–10. The intercepted ore body was meshed using a DIMINE-CAD-MIDAS/GTS coupling operation and imported into the numerical simulation software FLAC3D for calculation [16, 17].

The calculation steps of exploiting-filling-exploiting are designed according to the four sets of models established, and each model corresponds to one scheme. The simulation consists of three main steps in each scheme: the first step (I) is the exploitation of strips 3, 4, 7, and 8. The second step (II) is the filling of strips 3, 4, 7, and 8. The third step (III) is the exploitation of strips 1, 2, 5, 6, 9, and 10. The four simulation schemes of different strip widths are shown in Table 2.

##### 3.2. Simulation Results

###### 3.2.1. Displacement and Stress Analysis of the Roof

In the simulation process, it is found that the stability of the stope after the third step is the worst. Therefore, the displacement and stress of the roof after the third step in the four schemes are analyzed. The specific values are shown in Table 3.

Comparing these parameters in Table 3 with the maximum bearing capacity of the roof in Table 1, it can be concluded that the roof in the stope of the four schemes will not be damaged.

###### 3.2.2. Plastic Failure Analysis of the Backfill

After the third step, both sides of the backfill are goaf. When the pressure exerted by the roof on the backfill exceeds the compressive strength of the backfill, plastic failure will lead to overall collapse of the backfill [18]. Table 4 shows the plastic failure volumes of the backfill after the third step of the four schemes. Shear failure is mainly caused by the pressure exerted on the backfill by the top plate exceeding the compressive strength of the backfill. When the shear failure area is too large, the backfill will perforate. And the tensile failure is mainly located in the center of the roof, which is caused by the self-gravity stress of the overlying strata.

It can be concluded from Table 4 that the tensile failure volume is sufficiently small that it does not affect the stability of the backfill. Therefore, the shear failure is mainly analyzed as follows. The change curve of shear failure volume with strip width is shown in Figure 2.

As shown in Figure 2, the shear failure volume of the backfills in schemes 2, 3, and 4 exhibits a significant change compared with that of scheme 1, with a decrease of more than 1000 m^{3}. The overall stabilities of the backfills of these three schemes are better than that of scheme 1. As can be seen from Figure 3, the shear failure plastic zone passes through the backfill; thus, the backfill will lose its stability. “shear-n” in Figure 3 means the backfill is in plastic state now, which is the main factor of influencing the stability of the backfill. “shear-*n*” in schemes 2, 3, and 4 is far less than that in scheme 1.

In the two-step stopping process, there are two main forms of stope instability damage [19]. The first is when the deformation of the roof exceeds the allowable range, leading to collapse failure of the roof, or when the tensile stress produced by the roof exceeds its ultimate tensile strength, resulting in tensile failure. The second is when plastic penetration failure occurs when the pressure of the backfill exceeds its compressive strength, and the whole backfill collapses.

It can be concluded that when the width of the strip is 6 m, the overall strength of the backfill after the third step is smaller than that of the 7, 8, and 9 m strips, and plastic zone penetration failure occurs. Therefore, the analysis shows that 7, 8, and 9 m may be stable strip widths.

#### 4. Calculation of Mathews Stability Graph Probability Model

According to the previous finite element analysis, 7, 8, and 9 m are the strip widths that are likely to be applicable. The probability model of Mathews stability graph is used to further optimize the three strip widths.

##### 4.1. Mathews Stability Graph Probability Model

The original Mathews stability graph was divided into stable, potentially unstable, and potential caving zones according to the scatter of the stability data in Figure 4. The initial stability zones and graph devised by Mathews et al. were based on 50 cases [20].

After 1980, researchers like Potvin, Stewart, Forsyth, and Tureman redrew the Mathews stability graph by increasing the number of projects evaluated to 500 [21]. Mawdesley applied mathematical statistics to explain the Mathews stability graph in 2001 and 2002 [22, 23]. In 2004, Mawdesley adopted the method of logarithmic regression analysis, redefining the stability area and great destruction zone. Parallel lines in the equal probability graph are drawn to divide the stable area, unstable area, and caving zone [24].

The Mathews method is based on a stability graph relating two calculated factors: the Mathews stability number (), which represents the competency of the rock mass for a given stress condition, and the shape factor (*S*), or hydraulic radius, which accounts for the geometry of the surface.

The logistic regression line defining the stability boundaries is defined by Equations (1) and (2), where *f(z)* is the logit value. The logit value is analogous to the response variable in a linear regression model and is determined for each data point based on the shape factor and the stability parameter [25, 26]:where is the predicted log odds value, is a constant, are numerical coefficients, and is the predicted logit probability value.

In the logit model of the stability data, the probability of stability is expressed as a linear function of the shape factor, the Mathews stability number, and a constant, represented bywhere and are the shape factor and stability parameter, respectively.

##### 4.2. Parameters Determination

###### 4.2.1. Calculation of Stability Parameter

The calculation formula of stability index is as follows [27]:where is the modified -value, is the stress factor, is the joint orientation factor, and is the surface orientation factor. The modified -value is calculated using the following equation:where is the rock mass rating system and is the tunneling quality index system [28].

The indicators were graded by consulting the relevant data of the geological report of the Tiaoshuihe phosphate mine: (1) rock compressive strength is 205.70 MPa, scoring 12 points; (2) rock mass quality indicator RQD is 79.18%, scoring 17 points; (3) joint spacing scores 10 points; (4) joint state scores 10 points; (5) state of groundwater on average is 15 L/min, scoring 7 points; and (6) joint bearing scores −5 points:

In summary, the RMR points of the surrounding rocks of the up panel can be obtained by Equation (6), and *Q* is calculated from the following equation:

*Q*′ approximates the value of *Q* as 2.18.

Figure 5 shows how to calculate the parameters of *A*, *B*, and *C* [23]. *A* is determined from the ratio of the intact rock strength (unconfined compressive strength) to the induced stress at the center-line of the stope surface. The intact rock strength is 205.70 MPa, and the induced stress at the center-line of the stope surface is less than 20.5 MPa. Thus, the value of *A* is 1.

**(a)**

**(b)**

**(c)**

According to the actual project, the value of *B* is 0.9.

*C* is the surface orientation factor:where is the angle of dip.

For the ore body in the Tiaoshuihe phosphate mine, *α* is 15°. Thus, the value of *C* is 1.24.

Substituting these results into Equation (4) and taking the value of *Q*′ as 2.18 determine the stability number () of the surrounding rock of the upper panel. is calculated by the following equation:

###### 4.2.2. Calculation of Shape Parameter *S*

The shape factor *S* (or hydraulic radius *R*) reflects the shape and size of the stope. It can be calculated by the following equation:where and are the width and length of the stope, respectively.

Since the length of the ore chamber is set as 50 m, the widths of the ore chamber are 7, 8, and 9 m, and the obtained shape coefficients *S* are 3.07, 3.45, and 3.81, respectively.

##### 4.3. Analysis of Calculation Results

According to the previous Mathews stability probability model of Equations (2) and (3), the lines of any probabilistic model can be represented in the Mathews stability diagram. Set any probability value ƒ(*z*), the corresponding *z* can be got. Then, Equation (3) will be a linear function between and *S*. In Figure 6, different lines can represent different stable probabilities. The model is straight line, which is different from the traditional curve graph. The three points on the graph correspond to the stable probabilities when the widths of the strip are 7, 8, and 9 m.

Comparison of the coordinate points corresponding to the strip widths of 7, 8, and 9 m in Figure 6 shows that when the strip width is 9 m, the probability of stope stability is lower than 85%. Moreover, when the strip width is 7 or 8 m, the stability probability of the stope is above 85%, and the stability of the strip width of 7 m is better than that of 8 m; thus, the stability of the strip width of 9 m is the worst.

The stability probabilities can be quantified according to Equations (2) and (3), using a strip width of 7 m as an example:

It is calculated in Equations (11) and (12) that when the strip width is 7 m, the probabilities of stope stability and failure are 88.55% and 11.45%, respectively. Similarly, when the strip width is 8 m, the probabilities of stope stability and failure are 86.76% and 13.24%, respectively. When the strip width is 9 m, the probabilities of stope stability and failure are 84.94% and 15.06%, respectively.

The Mathews stability probability model is used to study the structural parameters of the mining field of the Tiaoshuihe phosphate mine. Moreover, the stability probability of the strip width of 9 m is lower than 85%, and the possibility of stope failure is more than 15%. Based on previous engineering experience, the unstable probability more than 15% is dangerous. Therefore, the strip width of 9 m is not suitable. Considering that the maximum width of a single drilling application is 7.6 m for the DD211L low-type drilling rig in the underground, it is suggested that the best strip width is 7 m. It is being applied well in Tiaoshuihe phosphate mine now, and the engineering practice proves that this parameter is suitable.

#### 5. Conclusions

(1)An accurate finite element model was established for the Tiaoshuihe phosphate mine, and the rock mechanics parameters were based on engineering specifications. It is concluded that when the width of the strip is 6 m, the overall strength of the backfill is the worst after the third step, and plastic penetration failure occurs easily. Strip widths of 7, 8, and 9 m may be applicable.(2)The study of the Mathews stability probability model, generation of the logarithmic coordinate graph, and calculation of the quantitative probability indicate that the stability of the goaf decreases gradually when the strip width ranged from 7, 8, to 9 m. When the width of the strip is 9 m, the probability of goaf instability is greater than 15%, which is dangerous.(3)Comprehensively considering that the maximum width of a single drilling application is 7.6 m for the DD211L low-type drilling rig in the Tiaoshuihe phosphate mine, and that the stability of the 7 m strip width is higher than that of the 8 m according to the results calculated previously, a strip width of 7 m is recommended.

#### Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

The support from National Natural Science Foundation of China (Grant no. 41672298) is sincerely appreciated.