#### Abstract

For the stability of the intervening pillar of the sublevel drilling open-stope subsequent filling mining method, the multifactor stability mechanical model of the intervening pillar under two different constraint conditions (Model 1 and Model 2) was established based on the elastic thin plate theory. Then, the cusp catastrophe equation and the necessary and sufficient conditions for the instability of the intervening pillar under two different constraint conditions were obtained by using the cusp catastrophe theory. Furthermore, the minimum thickness formula for the intervening pillar without instability under two different constraint conditions was derived, and the relationships between the minimum thickness of the intervening pillar and the factors, including the depth of the stope, the inclination of the orebody, the thickness of the orebody, the height of the stage, the length of the stope, and the mechanical properties of the orebody, were analyzed. Finally, the formula was used in the design of an intervening pillar between stopes 4-1 and 4-2 in Panlong Lead-Zinc Mine. The designed thickness of the pillar was 6.01 m by calculation, its actual thickness was 6.35–7.25 m in the mining process, and its average thickness was 6.5 m. Compared with the previously designed thickness of 7-8 m at the same stage, the pillar was 0.5–1.5 m smaller, which more effectively improved the recovery rate of the ore under the premise of ensuring the stability of the intervening pillar. This example of industrial application proves that it is feasible to use the cusp catastrophe theory to analyze the stability and parameter design of the intervening pillar under different constraints.

#### 1. Introduction

The sublevel drilling open-stope subsequent filling mining method is the combination of the empty field mining method and the filling mining method. The mining method has a large production capacity and good economic benefits, and it can consume a number of tailings to reduce the risk of tailings ponds and effectively maintain the stability of a high-level goaf [1]. However, in order to achieve safe and efficient mining, the analysis of the stability of the intervening pillar and the determination of its optimal thickness are the primary technical problems that need to be solved.

At present, substantial amounts of research on the stability analysis and structural parameters optimization of pillar have been conducted. For example, Zhang et al. [2] used the ultimate strength theory to analyze the instability mechanism of coal pillars. González Nicieza et al. [3] introduced an RMR value to determine a reasonable pillar size on the basis of the empirical formula. Wu et al. [4] used numerical simulation and orthogonal tests to optimize the structural parameters of stope and pillar. Luo et al. [5] used the numerical simulation method to obtain the optimal structural parameters of an isolated pillar in a pan area. Cao et al. [6] used the elastic theory to obtain the analytical solution of the mechanical model of a filling pillar. Lan et al. [7] applied the genetic algorithm of the Pareto optimal solution set to the multiobjective optimization problem of the response surface model and obtained the optimal pillar structure parameters. Wang et al. [8] used numerical simulation and the multiobjective ideal point method to optimize the structural parameters of stope pillars. Yin et al. [9] established a simplified formula for calculating the safety factor of a rectangular pillar based on the orthogonal range analysis method and used it to optimize the structural parameters of the stope pillar. In summary, the research methods of stability analysis and structural parameters of pillars mainly include numerical simulation, elastic plastic theory, mathematical methods, and so on.

The catastrophe theory, as a branch of the nonlinear theory, can study the characteristics of a system changing with the change of its control parameters, especially when the parameters change the performance of the system to cause catastrophe under certain conditions [10–12]. At present, the catastrophe theory has been widely used in economics and physics [13, 14]. Due to failure of rock mass with the characteristics of energy accumulation and sudden release, in recent years, many achievements have been made in the research of rock mass instability based on the catastrophe theory. Bala Padmaja et al. [15] used catastrophe theory to analyze the stability of the plane landslide mathematical model. Wang et al. [16] used the cusp catastrophe model to analyze the instability mechanism of a temporary ore wall and calculate the optimal thickness of the ore wall. Zhang et al. [17] constructed the cusp catastrophe model of stope instability and deduced the sufficient and necessary conditions of stope instability to optimize the stope structural parameters. Wang et al. [18] established the cusp catastrophe model of pillar instability based on the catastrophe theory and analyzed the mechanism of pillar instability and failure. Li et al. [19] deduced the critical load and the stress of pillar instability by applying the energy principle and the catastrophe theory; furthermore, they proposed the buckling model of pillar instability. Jiang et al. [20] established the cusp catastrophe model of a pillar according to the constitutive relation of the mechanical system and reasonably optimized the structural parameters of the room and the pillar.

In summary, the catastrophe theory has been applied in pillar instability and structural parameter optimization, but current research studies are focused mainly on the 2D state of the pillar, which only considers the height and width of the pillar and ignores the length and dip angle of the pillar, causing incomplete constraints. Based on the elastic thin plate and cusp catastrophe theory, the failure mechanism of the 3D intervening pillar between stopes under two different constraints was analyzed in this paper. The factors, including the orebody dip angle, orebody thickness, stage height, mining room length, mining depth, elastic modulus of the orebody, and constraint conditions of the pillar, are introduced to establish the instability model of the intervening pillar. Finally, the minimum thickness formula for the intervening pillar was derived. This formula is more comprehensive and suitable for practical application, and the design errors caused by neglecting critical factors can be effectively avoided.

#### 2. Intervening Pillar between Stopes and Its Mechanical Model

##### 2.1. Destruction of the Intervening Pillar between Stopes

Figure 1 shows a schematic diagram of the sublevel drilling open-stope subsequent filling mining method. Orebody is divided into stopes and intervening pillars along orebody strike. Sublevel contact roadways and rock drilling roadways are excavated between stage transport roadways, and the stope is divided into subsections. The medium-deep hole blast is used for ore caving. The caved ore is concentrated in the bottom of the stope, and the ore is transported out through the loading roadways and the ore-drawing roadway. The intervening pillar divides adjacent stopes and plays a supporting role for the mining room.

**(a)**

**(b)**

**(c)**

There are two failure types of the intervening pillar between stopes. The first failure type occurs when adjacent stopes are not filled in time and both sides of the intervening pillar are empty; the intervening pillar will be subjected to concentrated stress in the vertical direction, and it may form vertical cracks or collapse towards one side. The other failure type occurs when one side has been filled and the other side is empty; the pillar will be subjected to concentrated stress in the vertical direction and lateral pressure of the filling body, and the pillar may collapse towards the side of the goaf. Because the first failure type is mainly caused by problems in the management and production arrangement, the most common failure type is mainly the second failure type. In summary, structural parameters play a pivotal role in the stability of the intervening pillar between stopes, and the thickness of the pillar is the core of the structural parameters.

##### 2.2. Mechanics Model of an Intervening Pillar between Stopes

The structure of the intervening pillar between two stopes is shown in Figure 2. Among them, *a* is the width of the intervening pillar, *b* is the height of intervening pillar, *c* is the thickness of the intervening pillar, and *φ* is the dip angle of the orebody. The intervening pillar is mainly subjected to vertical crustal stress (*σ*_{v}), horizontal tectonic stress (*σ*_{h}), and lateral pressure of filling body (*q*).

**(a)**

**(b)**

According to the structural stability theory, the stability of the intervening pillar depends not only on the above three forces but also on its structural parameters and constraint conditions. As a part of the orebody, the intervening pillar is connected with the top and bottom pillars, so it is constrained by fixed support in the *y*-direction. For the hanging and heading wall of the intervening pillar, there are two cases. The first case is that the hanging and heading wall belong to the complete orebody and are artificially divided according to ore grade, so it is considered to be fixed support constraints, as shown in Figure 3(a). The second case is that the boundary between the orebody and the surrounding rock is clear or broken, so it is considered to have simply supported constraints, as shown in Figure 3(b).

**(a)**

**(b)**

The formula for the vertical crustal stress on the intervening pillar is as follows:where *ρ* is the average density of the overlying strata (kg/m^{3}), is the acceleration of gravity (9.8 N/kg), and *h* is the thickness of the overlying strata (m).

If there is no in situ stress data, the distribution law of the in situ stress field in China proposed by Zhao et al. [21] can be used as a reference.

The pillar area bearing theory is widely used to calculate the pillar load. The theory holds that the pillar load is the weight of the overlying strata directly reaching the surface of the mining space supported by the pillar. The pillar support area is the sum of the allocated mining area and the pillar’s own area [22]. According to the theory, the force balance equation of the pillar is obtained as follows:where *σ*_{y} is average axial stress of the column (Pa) and *L* is the width of the mined-out area (m).

Therefore, the average stress on the intervening pillar in the *y*-direction can be obtained as follows:

Similarly, the average stress in the *x*-direction between pillars can be obtained as follows:

Furthermore, average pressure per unit lengths along *x*-direction and *y*-direction can be obtained as follows:where *N*_{y} is pressure per unit lengths of the *x* side of the intervening pillar (N/m) and *N*_{x} is the pressure per unit lengths of the *y*-side of the intervening pillar (N/m).

#### 3. Mechanical Characteristics of the Instability of the Intervening Pillar between Stopes

##### 3.1. Elastic Thin Plate Model of the Intervening Pillar between Stopes

###### 3.1.1. Four-Sided Clamped Elastic Thin Plate Model (Model 1)

The boundary conditions of the pillar are as follows:where *ω* is the deflection of the plate.

A deflection surface equation satisfying the above boundary conditions is constructed based on the Rayleigh–Ritz method [23] as follows:where *A* is the setting constant.

###### 3.1.2. Elastic Thin Plate Model with Two Opposite Sides Clamped and Two Opposite Simply Supported Pillars (Model 2)

The boundary conditions of the pillar are as follows:

A deflection surface equation satisfying the above boundary conditions is constructed based on the Rayleigh–Ritz method [23] as follows:

##### 3.2. Total Potential Energy of Intervening Pillar (*V*)

According to the stress condition of the intervening pillar, the total potential energy *V* is composed of the deformation potential energy (*U*), the external potential energy ( and ) caused by the combination of vertical crustal stress (*σ*_{v}) and horizontal tectonic stress (*σ*_{h}), and the external potential energy () caused by lateral pressure of filling body (*q*). It is as follows:

###### 3.2.1. Deformation Potential Energy (*U*)

The calculation formula for the deformation potential energy (*U*) is as follows:where *D* is flexural rigidity of the intervening pillar, *k*_{x} is the curvature of the intervening pillar in the *x*-direction, *k*_{y} is the curvature of the intervening pillar in the *y*-direction, *E* is the elastic modulus of the orebody, and *μ* is Poisson’s ratio of the orebody.

Combining equations (7) and (11) and expanding the Taylor series of equation (11) to five terms, the deformation potential energy of Model 1 is obtained as follows:

Similarly, the deformation potential energy of Model 2 is obtained as follows:

###### 3.2.2. External Force Potential Energy ( and )

The calculation formula for the external force potential energy ( and ) is as follows:

Combining equation (7) with (17), the external force potential energy ( and ) of Model 1 is obtained as follows:

Similarly, the external force potential energy ( and ) of Model 2 is obtained as follows:

###### 3.2.3. External Potential Energy ()

For the lateral pressure of the filling body, we can use the Rankine soil pressure theory [24] to calculate the total active lateral pressure. The calculation formula for the active lateral pressure is as follows:where *γ*_{c} is the weight of filling body and *φ*_{c} is the internal friction angle of the filling body.

The calculation formula for the external potential energy () is as follows:

Combining equation (7) with (21), the external force potential energy () of Model 1 is obtained as follows:

Similarly, the external force potential energy () of Model 2 is obtained as follows:

###### 3.2.4. Total Potential Energy of Intervening Pillar (*V*)

The total potential energy of the intervening pillar (*V*) of Model 1 is obtained by using equations (9), (15), (18), and (22) as follows:

Similarly, the total potential energy of the intervening pillar (*V*) of Model 2 is obtained by using equations (9), (16), (19), and (23) as follows:

##### 3.3. Cusp Catastrophe Model of the Intervening Pillar

For the total potential energy (*V*) of Model 1, we order

For the total potential energy (*V*) of Model 2, we order

Dimensionless parameters (*x*, *u*, ) are introduced as follows:

In addition, the dimensionless parameters of Model 1 and Model 2 are introduced, respectively, as follows:

The potential function *V*(*x*) of Model 1 and Model 2 is expressed as a standard form:where *x* is the state parameter and and are the control parameters.

It is easy to find that equation (30) is the potential function form of the cusp catastrophe model, which indicates that the stability model of the intervening pillar conforms to the cusp catastrophe model. The potential function *V*(*x*) is used to solve the first derivative of the state parameter *x* and is made to equal to 0, and the stable equilibrium curved surface equation of the intervening pillar is obtained as follows:

Furthermore, we can obtain Figure 4(a). The curved surface is divided into upper, middle, and lower lobes. The state of the intervening pillar varies with the change of the state parameter *x*. The critical state of the intervening pillar is obtained by deriving *x* again as follows:

**(a)**

**(b)**

**(c)**

The bifurcation set equation of the intervening pillar can be obtained by combining equations (31) and (32), as shown in Figure 4(c):

##### 3.4. Instability Conditions of the Intervening Pillar

###### 3.4.1. Necessary Instability Conditions of the Intervening Pillar

It can be seen from Figure 4(c) that only when *u*^{3} ≤ 0 can the equilibrium point of the intervening pillar across the bifurcation be set to lead to instability. Therefore, *u*^{3} ≤ 0 is a necessary instability condition of the intervening pillar. Furthermore, we can obtain *s* ≤ 0 as a necessary instability condition of the intervening pillar.

###### 3.4.2. Sufficient Instability Conditions of the Intervening Pillar

According to cusp catastrophe theory, when the system is unstable, the system state must satisfy the bifurcation set equation. At this time, the system is most sensitive to external disturbances and is in the critical state of instability. Therefore, the bifurcation set equation (33) is a sufficient condition for the intervening pillar. There is always , so if the bifurcation set equation (33) is established, *u*^{3} ≤ 0 must be satisfied. Therefore, *s* ≤ 0 is a necessary instability condition of the intervening pillar.

###### 3.4.3. Sufficient and Necessary Instability Conditions of the Intervening Pillar

Because the sufficient and necessary instability conditions of the intervening pillar are consistent, the sufficient and necessary instability conditions of the intervening pillar of Model 1 and Model 2 are as follows:

##### 3.5. Calculation of the Minimum Thickness of the Intervening Pillar

According to equations (34) and (35), the following can be deduced:

The minimum thickness of the intervening pillar can be obtained as follows:

Among them, *B* = *B*_{1} for Model 1 and *B* = *B*_{2} for Model 2:

##### 3.6. Designing Thickness of the Intervening Pillar

To guide the industrial scene design and construction more safely and scientifically, it is necessary to multiply the theoretical thickness of the intervening pillar by a certain safety factor. Zhao et al. [25] used fractal theory to obtain the safety factor of the pillar to ensure that the stability of the stope is 1.44; the designing thickness of intervening pillar is as follows:

#### 4. Analysis of the Influencing Factors of the Intervening Pillar Thickness

Equations (37)–(40) show that the thickness of intervening pillar mainly depends on the vertical and horizontal stress, the dip angle of the orebody, the thickness of the orebody, the stage height, the length of the mining room, the mechanical properties of the orebody, and the constraint conditions. To intuitively analyze the influence of various factors on the thickness of the intervening pillar, the control variable method can be used to analyze these factors one by one: stope depth *h* = 250 m, orebody dip angle *φ* = 70°, orebody thickness *a* = 20 m, stage height *b* = 50 m, the length of the mining room *L* = 40 m, orebody elastic modulus *E* = 2.5 GPa, and orebody Poisson’s ratio *μ* = 0.3.

##### 4.1. Thickness of the Overlying Strata *h*

Based on the above conditions, the diagram of the minimum thickness of the intervening pillar is drawn with the thickness of the overlying strata in Model 1 and Model 2 (as shown in Figure 5).

It can be seen from Figure 5 that the minimum thickness of Model 1 and Model 2 is positively correlated with the thickness of the overlying strata; however, with the increase in the thickness of the overlying strata, the minimum thickness of Model 2 increases faster than that of Model 1, and the ratio *c*_{1}/*c*_{2} is between 0.65 and 0.66.

##### 4.2. Inclination Angle of the Orebody *φ*

Based on the above conditions, the diagram of the minimum thickness is drawn with the inclination angle of the orebody in Model 1 and Model 2 (as shown in Figure 6).

From Figure 6, it can be seen that the minimum thickness both Model 1 and Model 2 first increases and then decreases with the dip angle of the orebody; it basically reaches the maximum value with the dip angle of the orebody of 60°. However, with the increase in the dip angle of the orebody, the minimum thickness of Model 2 increases faster than that of Model 1. Furthermore, the ratio *c*_{1}/*c*_{2} decreases with the dip angle in the orebody, and it is between 0.65 and 0.685.

##### 4.3. Orebody Thickness *a*

Based on the above conditions, the diagram of the minimum thickness is drawn with the orebody thickness in Model 1 and Model 2 (as shown in Figure 7).

From Figure 7, it can be seen that the minimum thickness of Model 1 and Model 2 is positively correlated with the orebody thickness, but with the increase in orebody thickness, the minimum thickness of Model 2 increases faster than that of Model 1. Furthermore, the ratio *c*_{1}/*c*_{2} increases with the increase in the orebody thickness, and the ratio is between 0.62 and 0.74.

##### 4.4. Stage Height *b*

Based on the above conditions, the diagram of the minimum thickness is drawn with the stage height in Model 1 and Model 2 (as shown in Figure 8).

It can be seen from Figure 8 that the minimum thickness of Model 2 is positively correlated with the stage height; the minimum thickness of Model 1 first increases and then stabilizes smoothly with the increase of the stage height. This finding shows that when the stage height is greater than the orebody thickness, the orebody thickness plays a decisive role in determining the minimum thickness of the intervening pillar. Furthermore, when the stage height is equal to 15 m, the minimum thickness of Model 1 and Model 2 is equal, and the gap increases with the increase in the stage height. The ratio *c*_{1}/*c*_{2} decreases with the increase in the stage height, and it is between 0.65 and 1.00.

##### 4.5. Mining Room Length *L*

Based on the above conditions, the diagram of the minimum thickness is drawn with the mining room length (as shown in Figure 9).

It can be seen from Figure 9 that the minimum thickness of Model 1 and Model 2 increases linearly with the mining room length, and the ratio *c*_{1}/*c*_{2} keeps at 0.65.

##### 4.6. Elastic Modulus of Orebody *E*

Based on the above conditions, the diagram of the minimum thickness is drawn with the elastic modulus of orebody in Model 1 and Model 2 (as shown in Figure 10).

It can be seen from Figure 10 that the minimum thickness of Model 1 and Model 2 decreases with the increase in the elastic modulus of the orebody, and the ratio *c*_{1}/*c*_{2} is maintained at 0.65.

#### 5. Engineering Example

The No. 2 orebody of Panlong Lead-Zinc Mine in Guangxi occurs in the upper interlayer fracture zone of the Lower Devonian Upper Lun Dolomite, which belongs to submarine jet-sedimentary lead-zinc deposit. The orebody is mainly stratiform, layered and lenticular, locally cystic, generally thick in the middle and thin at the edges, wedge-shaped branching or pinching, with pinching reappearance, and complex branching phenomena. The orebody is generally trending to the northeast with a length of 2000 m, a tendency of 340°, a depth of 1100 m, and an inclination of 77°∼88°. The orebody occurrence elevation is 62 to 1140 m. The thickness of the orebody is 0.89 to 56.06 m, and the average thickness is 17 m.

The surface elevation of the mine is approximately +72 m, and the mining stage is 50 m. In this design, the intervening pillar between stope 4-1 and 4-2 at the depth of 120∼70 m is selected for calculation. The average density of the overlying rock mass is 2800 kg/m^{3}, and the vertical stress is 3.976 MPa. Because the upper orebody is dominated by sedimentary lead-zinc deposits, the horizontal stress is 6.480 MPa as calculated according to [13]. The average thickness of the orebody *a* = 22 m, the average dip angle is 78°, the length of the stope 4-1 *L*_{1} = 44.27 m, and the length of the stope 4-2 *L*_{2} = 42.28 m. The stope layout is shown in Figure 11. In addition, the hanging and heading wall of the intervening pillar are obvious and fragmented, so they can be considered as simply supported constraints. Finally, the mechanical parameters of the orebody *E* = 2130 MPa, Poisson’s ratio *u* = 0.245, and the internal friction angle of filling body *φ*_{c} = 35°.

The thickness of the intervening pillar between stopes was calculated to be 6.01 m by equations (37), (39), and (40). After stope 4-1 was filled and stope 4-2 was mined, the thickness of the intervening pillar was 6.35∼7.25 m and the average thickness was 6.5 m. Compared with the previous design, the thickness of the intervening pillar was reduced by 0.5∼1.5 m, its stability is good, and there is no large area slope or instability damage based on the 3D scanning result of the goafs. The 3D scanning graphic of the goafs is shown in Figure 12. The intervening pillar effectively improves the ore recovery rate.

#### 6. Conclusion

(1)To study the stability of the intervening pillar in sublevel drilling block open-stope subsequent filling mining method, multifactor stability mechanical models suitable for two boundary conditions (four-sided clamped, two opposite sides clamped and two opposite simply supported) were built based on the elastic thin plate theory. Combined with the cusp catastrophe theory, the cusp catastrophe equation and judgment criterion suitable for two boundary conditions were obtained to analyze the instability of the intervening pillar.(2)Furthermore, the minimum thickness calculation formulas for Model 1 and Model 2 were derived based on the instability judgment criterion of the intervening pillar. Through analysis of the influence of factors on the minimum thickness, we can find (1) the minimum thickness of Model 1 is less than that of Model 2 and that explains that simply supported boundary condition has higher requirements for the thickness of the intervening pillar in order to satisfy the stability. (2) The thickness of the intervening pillar is positively correlated with some factors including mining depth, orebody thickness, the stage height, and mining room length. Among these factors, mining depth, orebody thickness, and the stage height have great influence on the thickness of the intervening pillar. (3) The thickness of the intervening pillar first increases and then decreases with the dip angle of the orebody. In addition, the thickness of the intervening pillar basically reaches the maximum value at 60° dip angle of orebody. (4) The thickness of the intervening pillar is negatively correlated with an elastic modulus of the orebody, and the elastic modulus of the orebody has great influence on the thickness of the intervening pillar. (5) The smaller stage height and orebody thickness play a decisive role in its influence. That is, when the stage height is larger than orebody thickness, orebody thickness will have a decisive effect on the thickness of the intervening pillar and vice versa.(3)Finally, the theory is applied to the calculation and design of the intervening pillar between stopes 4-1 and 4-2 in Panlong Lead-Zinc Mine, and its theoretical calculation result is 6.01 m. After stope 4-1 was filled and stope 4-2 was mined, the thickness of the intervening pillar was 6.35∼7.25 m and the average thickness was 6.5 m. Compared with the previous design, the thickness of the intervening pillar was reduced by 0.5∼1.5 m, its stability is good, and there is no large area slope or instability damage based on the goaf monitoring data. Therefore, the intervening pillar effectively improves the ore recovery rate. Meanwhile, this example proves the feasibility of this method.

#### Data Availability

The orebody and stope data used to support the findings of this study are included within the article.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This work was supported by the National Key R&D Program of China during the Thirteenth Five Year Plan Period: The Continuous Mining Theory and Technology on Spatiotemporal Synergism of Multi-mining Areas within a Large Ore Block for Deep Metal Deposit (2017YFC0602901).