Numerical and theoretical analysis of multi-pillar instability under elastic beams

Abstract

With the increase of mining scale and depth, the cascading pillar failure (CPF) disaster has gradually become one of the core technical challenges for safe mining. This paper studies the CPF disaster of multi-pillar (P_1-2, P_2-3, P_3-4 and P_4-5) at 848 m level of Alhada Lead-Zinc Mine based on the RFPA^2D numerical simulation software and stiffness theory. The numerical simulation results indicate that the load transfer effect induces the domino instabilities of double pillar (P_2-3 and P_3-4) and double pillar (P_1-2 and P_4-5), which is characterized by the ‘2 + 2’ compound failure mode. The physical essence of load transfer effect is revealed through numerical simulation, which is the elastic rebound of surrounding rockmass (such as roof-floor). In addition, the theoretical model of rebound overload mechanism of roof-multi-pillar-floor system is established based on the catastrophe theory, and the instability criterion, sudden jump and energy release of roof-multi-pillar-floor system are derived. Finally, the main influencing factors of load transfer effect such as pillar spacing and damage fracture zone are quantitatively analyzed, and the load transfer law shows that the load transfer effect decreases with the increasing pillar spacing and is hindered by the structural plane in the surrounding rockmass.

Keywords:

• Multi-pillar
• load transfer
• instability
• numerical simulation
• stiffness theory

1. Introduction

The open stope mining has left a large number of pillars to support the goaf. The load transfer effect between pillars may induce the CPF, resulting in the large-scale collapse of goaf group (Sai Vinay et al. 2022). The CPF of multi-pillar can be summarized into three failure modes: domino failure mode, compound failure mode and successive failure mode, and has the characteristics of catastrophic, suddenness, and unpredictability (Dong et al. 2023). For example, the Coalbrooke disaster in South Africa killed 437 miners and created a 2 km² subsidence area within a few minutes on January 21, 1960 (Zipf 2001). Such catastrophic disasters occur from time to time and have caused tremendous casualties and property losses, where the CPF of multi-pillar has become one of the core technical challenges that restrict the safe mining (Wang et al. 2008; Dehghan et al. 2013; Li 2016; Vardar et al. 2017; Zhou et al. 2019).

Due to the energy storage characteristics of rockmass, the roof-multi-pillar-floor system stores a large amount of elastic energy during its whole deformation and failure process. When pillar occurs instability, the unsupported roof-floor will rebound and compress adjacent pillars, inducing a load transfer effect among multi-pillar (Dong et al. 2022b). In addition, Wang et al. (2021) verifies the feasibility of the elastic foundation hypothesis and studies the instability mechanism of roof-pillar system. Therefore, the load transfer effect is not only related to the mechanical properties of pillars, but also influenced by the stiffness of roof-floor (Kaiser and Tang 1998; Wang et al. 2011; Gao et al. 2019). The local mine stiffness (LMS) is commonly used to evaluate the stability and impact tendency of engineering rockmass structures (Cook 1965; Starfield and Fairhurst 1968; Salamon 1970). When LMS of roof-floor in the pillar support area is less than the post-peak stiffness of pillar, the pillar will experience unstable failure. Kan et al. (2022) estimates the rockmass stiffness of roof-floor and analyzes the energy release mode of roof-single-pillar-floor system by numerical simulation. However, due to the complex geological conditions of surrounding rockmass, accurately estimating the LMS of rockmass has become the main challenge. Some researchers estimate the LMS by establishing a true 3D geological model that fully investigate the geological conditions and combining the numerical simulation methods (Zipf 1992; Jaiswal and Shrivastva 2012; Kias 2013; Gu and Ozbay 2015).

Numerical simulation is an important approach for studying large-scale engineering problems (Zhao et al. 2020), especially the load transfer effect of roof-multi-pillar-floor system. Some studies have used numerical simulation method to analyze the instability of multi-pillar or goaf group (Kaiser and Tang 1998; Wang et al. 2011; Dehghan et al. 2013; Li et al. 2013; Hauquin et al. 2016; Zhang et al. 2017; Zhou et al. 2018). Guggari et al. (2023) employs a finite element approach to investigate the stress concentration and relaxation process that may cause the dynamic failure of multi-pillar. Sai Vinay et al. (2022) study the influence of different extraction lines (such as steep-diagonal line and straight line) on the stability of overlying strata through numerical simulation. Poulsen and Shen (2013) make quantitative assumptions about load transfer and use probability statistical methods to assess the risk of CPF disasters. However, with the increasing scale and depth of mining, the physical essence of load transfer undergoes a fundamental transformation from stress redistribution to rebound overload of roof-floor under elastic foundation. Thus, the researches on the load transfer law of roof-multi-pillar-floor system is insufficient, and the disaster process analysis of multi-pillar instability is still lack of theoretical basis.

In this study, the CPF disaster of multi-pillar (e.g. P_1-2, P_2-3, P_3-4 and P_4-5) at 848 m level of Alhada Lead-Zinc Mine is taken as the research object. Firstly, a numerical model is established to reproduce the CPF disaster of multi-pillar, and the CPF modes of multi-pillar are analyzed based on the stiffness theory in Section 2. Secondly, the catastrophe theory model is established to derive the analytical expressions of sudden jump and energy release of roof-multi-pillar-floor system in Section 3. Finally, the numerical simulation and theoretical derivation results are compared and the load transfer behavior during CPF is quantitatively analyzed to reveal the load transfer laws in Section 4. This case study will provide technical reference and theoretical support for safe mining.

2. Numerical simulation of multi-pillar instability

In this Section, a multi-pillar numerical model is established based on the engineering case. The basic mechanical parameters of rockmass are obtained by field sampling and laboratory test, and the numerical simulation parameters are obtained according to the transformation relationship between macro and mesoscopic parameters. The instability behavior of multi-pillar is studied based on the local mine stiffness (LMS) theory and RPPA^2D numerical simulation software.

2.1. Numerical model and results

2.1.1. Engineering background

The Alhada Lead-Zinc Mine is located in Inner Mongolia at an altitude of 922 to 1045 m. It is a polymetallic skarn deposit (Lead-zinc-silver) formed by the interaction between hypabyssal and ultrahypabyssal intrusive rocks and carbonate rocks. The length of the orebody is 2100 m, with an average thickness of 1 to 10 m and an average inclination angle of 75°. The rockmass quality is classified as Grade III, and the rockmass structure is relatively fragmented. Before the year 2014, the mining filling system had not yet developed, and the open stope method was used at 808, 848 and 888 m level, resulting in more than 120 goafs totaling approximately 601000 m³. Especially the goaf group at 848 m level has the characteristics of large scale and interconnectivity, which can induce domino instability of multi-pillar and large-scale surface subsidence, making it the focus of this study. Figure 1 shows the model of multi-pillar and goaf group located at 848 m level with a burial depth of 150 m, consisting of 4 pillars (P_1-2, P_2-3, P_3-4 and P_4-5) and 5 goafs (1#, 2#, 3#, 4# and 5#). The mining sequence first involves extracting the 1# mining room, followed by the sequential extraction of 2#, 3#, 4# and 5# mining rooms, ultimately forming a single-layer horizontal goaf group and multi-pillar. The mining room has a span of 35 m, with pillars 30 m high, 10 m wide, and 6 m thick.

Figure 1. Model of multi-pillar and goaf group.

2.1.2. Model setup and test plan

Figure 2 shows the RFPA numerical model of multi-pillar and goaf group at 848 m level. The mesh size is taken as 500 × 300 with 150000 elements being used to discretize the whole numerical model 500 m wide × 300 m high. There are 5 goafs in the middle of the numerical model along the horizontal direction, represented by 1#, 2#, 3#, 4# and 5# respectively. The pillar between goafs 1# and 2# is represented by P_1-2, similarly, the other pillars are represented by P_2-3, P_3-4 and P_4-5, respectively. Horizontal displacement constraint is applied to the boundary on both sides of the numerical model, normal displacement constraint is applied to the bottom boundary, and the upper surface is a free boundary.

Figure 2. Numerical test model of multi-pillar and goaf group at 848 m level.

The method of applying incremental loads at the top of the model is adopted to approximate the gravity increase method. The loading rate is determined by referring to the field monitoring data and rock mechanics experiments. Thus, the vertical constant loading rate (5 mm/step) is applied at upper surface of the numerical model. The stress of goaf roof reaches 12.51 MPa at step 32, which is used to simulate the initial ground stress at 150 m deep. And then continuous loading is applied to simulate the creep loading behavior from surrounding rockmass until multi-pillar collapse. The rockmass is regarded as a non-uniform material, and its mechanical properties follow the Weibull distribution. According to Equation (1)(1) $$\{\begin{array}{c}\frac{f_{\mathrm{\text{macro}}}}{f_{\mathrm{\text{meso}}}}=0.2062\mathit{\text{lnm}}+0.0233(1.2\le m\le 50)\\ \frac{E_{\mathrm{\text{macro}}}}{E_{\mathrm{\text{meso}}}}=0.1412\mathit{\text{lnm}}+0.6476(1.2\le m\le 50)\end{array}$$(1) (Zhu 2008), the macroscopic mechanical parameters of the rockmass in Table 1 can be transformed into mesoscopic mechanical parameters in Table 2. (1) $$\{\begin{array}{c}\frac{f_{\mathrm{\text{macro}}}}{f_{\mathrm{\text{meso}}}}=0.2062\mathit{\text{lnm}}+0.0233(1.2\le m\le 50)\\ \frac{E_{\mathrm{\text{macro}}}}{E_{\mathrm{\text{meso}}}}=0.1412\mathit{\text{lnm}}+0.6476(1.2\le m\le 50)\end{array}$$(1)

Table 1. Basic mechanical parameters of rockmass (macro parameters).

Rock type	Cohesion (MPa)	Internal friction angle (°)	Uniaxial compressive strength σc (MPa)	Tensile strength σt (MPa)	Elastic modulus E (GPa)	Density (g/cm^3)	Poisson’s ratio
Skarn rock	2.07	33.89	40.15	2.48	11.09	2.78	0.21

Table 2. Input Weibull parameters of rockmass (mesoscopic parameters).

Rock type	Homogeneity	Internal friction angle (°)	Uniaxial compressive strength σc (MPa)	σc/σt	Elastic modulus E (GPa)	Density (g/cm^3)	Poisson’s ratio
Skarn rock	3	33.89	160.71	16.19	13.82	2.78	0.21

$f_{\mathrm{\text{meso}}}$ and $E_{\mathrm{\text{meso}}}$ represent the strength and elastic modulus of Weibull distribution, while $f_{\mathrm{\text{macro}}}$ and $E_{\mathrm{\text{macro}}}$ represent the macroscopic strength and elastic modulus of rock specimen.

2.1.3. Numerical results

Figure 3 shows the numerical simulation results of multi-pillar instability. Representative stress cloud diagrams at different damage evolution stages are selected, including step 32 of initial stress state, step 90 of critical state of multi-pillar, step 96 of CPF of P_2-3 and P_3-4, step 104 of roof-floor tensile failure, step 106 of CPF of P_1-2 and P_4-5, step 120 of goaf group collapse. The detailed instability process of roof-multi-pillar-floor system is described as follows:

Figure 3. Numerical simulation results of CPF process of multi-pillar (or goaf group).

1. Step 32 shows the initial stress state (12.51 MPa) of roof-multi-pillar-floor system at 150 m deep. The multi-pillar (P_1-2, P_2-3, P_3-4 and P_4-5) is in a stable state and forms a stable support system with roof-floor, whose maximum stress is 10.65 MPa.

2. Step 90 shows the critical state of multi-pillar, and the maximum stress inside the multi-pillar (P_1-2, P_2-3, P_3-4 and P_4-5) is 34.7 MPa, which is close to the uniaxial compressive strength 40.15 MPa. Therefore, obvious damage elements generate in pillars P_1-2, P_2-3 and P_3-4, and tensile cracks begin to form in the roof-floor of goafs 1#, 2# and 4#.

3. Step 96 shows the CPF of P_2-3 and P_3-4, where P_2-3 and P_3-4 fail at step-in-step 96-6 and 96-11, respectively. The CPF of P_2-3 and P_3-4 leads to the connection of goafs 2#, 3# and 4#. The load transfer effect induced by CPF of P_2-3 and P_3-4 causes significant damage elements in P_1-2 and P_4-5, and accelerates the crack propagation in the floor of P_3-4 and P_4-5. The crack in the floor of P_3-4 further expands to the deep under tensile stress, while the crack in the floor of P_4-5 stops expanding due to the tensile stress release of crack in the floor of P_3-4.

4. At step 104, the connection of goafs 2#, 3# and 4# results in the significant increase of hanging arch length and tensile stress of roof-floor. Thus, two main tensile cracks generate in the middle of roof-floor.

5. Step 106 shows the CPF of P_1-2 and P_4-5, where P_1-2 and P_4-5 fail simultaneously from step-in-step 106-1 to 106-7. The CPF of P_1-2 and P_4-5 results in the connection of goaf group 1#, 2#, 3#, 4# and 5# and significant increase of hanging arch length and tensile stress of roof-floor. The two main tensile cracks in the middle of roof-floor further expand. In addition, the load transfer effect induced by CPF of P_1-2 and P_4-5 causes significant stress concentration, which mainly distributes on the edges of goafs 1# and 5# with a maximum principal stress of 34.7 MPa.

6. The two main tensile cracks in the middle of the roof-floor extend deeper, ultimately causing the goaf group collapse at step 120.

Additionally, the stability of goaf group at 848 m level of Alhada Lead-Zinc Mine are also evaluated based on the modified stability graph method, and it is considered that all the tested stopes are in the caving zone (Jia et al. 2020). Furthermore, it is found that all goaf group have collapsed in order from middle to both sides according to the on-site investigation, further verifying the numerical simulation results. The concept of ‘step-in-step’ refers to the step-by-step loading and analysis technique used by software to simulate the rock failure process (Tang 1997). For example, when some parts of the rock model are damaged, the stress will be redistributed within the model, and this stress redistribution process is realized by the step-by-step technique.

2.2. Instability analysis of multi-pillar

2.2.1. Rockmass LMS calculation

Salamon (1970) proposed the concept of Local Mine Stiffness (LMS), which is defined as the equivalent stiffness of roof-floor within specific area. And K_1-2、K_2-3、K_3-4 and K_4-5 are the LMS of roof-floor where pillars P_1-2、P_2-3、P_3-4 and P_4-5 are located. According to the method of Jaiswal and Shrivastva (2012), the calculation process of K_1-2、K_2-3、K_3-4 and K_4-5 is as follows:

1. As shown the Case I of Table 3, initial stress is applied to the numerical model with pillar support. The deformation (i.e. d_1-2, d_2-3, d_3-4 and d_4-5) and load (i.e. f_1-2, f_2-3, f_3-4和f_4-5) of each pillar are recorded. The monitoring points of deformation and load are represented by four red dots in Table 3, which are the intersection points between the central axis of pillar and the boundary of roof.

2. As shown the Case II of Table 3, initial stress is applied to the numerical model without pillar support. The deformation (i.e. d’_1-2, d’_2-3, d’_3-4 and d’_4-5) and load (i.e. f’_1-2, f’_2-3, f’_3-4 and f’_4-5) of each imaginary pillar are recorded.

3. The LMS K_1-2、K_2-3、K_3-4 and K_4-5 of roof-floor where P_1-2、P_2-3、P_3-4 and P_4-5 located are calculated according to Equation (2)(2) $$\{\begin{array}{c}{K}_{1-2}=\frac{f_{1-2}}{d_{1-2}^{\prime }-d_{1-2}}\\ K_{2-3}=\frac{f_{2-3}}{d_{2-3}^{\prime }-d_{2-3}}\\ K_{3-4}=\frac{f_{3-4}}{d_{3-4}^{\prime }-d_{3-4}}\\ K_{4-5}=\frac{f_{4-5}}{d_{4-5}^{\prime }-d_{4-5}}\end{array}$$(2) , e.g. K_1-2 = 37.39 GN/m, K_2-3 = 32.69 GN/m, K_3-4 = 31.85 GN/m and K_4-5 = 38.20 GN/m, respectively.

(2) $$\{\begin{array}{c}{K}_{1-2}=\frac{f_{1-2}}{d_{1-2}^{\prime }-d_{1-2}}\\ K_{2-3}=\frac{f_{2-3}}{d_{2-3}^{\prime }-d_{2-3}}\\ K_{3-4}=\frac{f_{3-4}}{d_{3-4}^{\prime }-d_{3-4}}\\ K_{4-5}=\frac{f_{4-5}}{d_{4-5}^{\prime }-d_{4-5}}\end{array}$$(2)

Table 3. Stress-strain State of multi-pillar system under initial in-situ stress.

	Case I				Case II
Numerical model
Pillar No.	P1-2	P2-3	P3-4	P4-5	P1-2	P2-3	P3-4	P4-5
Convergence displacement/mm	d1-2	d2-3	d3-4	d4-5	d’1-2	d’2-3	d’3-4	d’4-5
	167.75	170.45	168.75	169.70	241.93	260.94	261.00	242.15
Load/MN	f1-2	f2-3	f3-4	f4-5	f’1-2	f’2-3	f’3-4	f’4-5
	2761.37	2957.85	2933.05	2761.33	0	0	0	0

2.2.2. Instability properties of multi-pillar

The rapid rebound of elastic foundation (or roof-floor) was the physical essence of instability of roof-multi-pillar-floor system, which will have a strong impact on multi-pillar (Dong et al. 2022a, 2022b). Thus, the instability properties of multi-pillar are analyzed as follows:

1. Stability analysis of multi-pillar failure process

   When evaluating the stability of roof-multi-pillar-floor system, the LMS of roof-floor is compared with the post-peak stiffness of multi-pillar (Zipf and Mark 1997). If the LMS of roof-floor is less than the absolute value of post-peak stiffness of multi-pillar, it is considered that the pillar will occur unstable failure (Yavuz 2001). Therefore, by comparing the LMS of roof-floor (i.e. K_1-2、K_2-3、K_3-4 and K_4-5) with the post-peak stiffness (i.e. k_p1-2、k_p2-3、k_p3-4 and k_p4-5) of pillars P_1-2、P_2-3、P_3-4 and P_4-5, the stability of the pillar failure process can be obtained. As shown in Figure 4(a), the post-peak stiffness of P_1-2, P_2-3, P_3-4 and P_4-5 is k_p1-2 = 46.38 GN/m, k_p2-3 = 47.73 GN/m, k_p3-4 = 44.08 GN/m and k_p4-5 = 47.95 GN/m, respectively. The stiffness ratios (K_1-2/k_p1-2, K_2-3/k_p2-3, K_3-4/k_p3-4 and K_4-5/k_p4-5) are all less than 1, so the multi-pillar will occur unstable failure.

2. Load transfer and sudden jump of roof-multi-pillar-floor system

   As shown in Figure 4(a), the CPF of P_2-3 and P_3-4 induces the sudden load increases (or load transfer) in P_1-2 and P_4-5, which is 714.81 MN and 480.03 MN, respectively. The transferred load does not exceed the bearing capacity of P_1-2 and P_4-5. As shown in Figure 4(b), each pillar instability is accompanied by a sudden jump displacement, which is essentially caused by the elastic rebound of roof-floor due to the loss of pillar support. Thus, the sudden jump displacement (or elastic rebound of local roof-floor) of P_1-2, P_2-3, P_3-4 and P_4-5 is 105.26, 63.91, 106.19 and 69.36 mm, respectively. Furthermore, the detailed analysis of load transfer law will be introduced in Section 4.2.

3. Energy release of roof-floor (or elastic foundation)

   The roof-floor that loses pillar support will rebound and release elastic energy, and the calculation formula is $W=\frac{1}{2}K{\Delta d}^2.$ K and Δd is the LMS and elastic rebound of roof-floor. Thus, the elastic energies of roof-floor where P_1-2, P_2-3, P_3-4 and P_4-5 are located are calculated as 102.87, 133.84, 137.59 and 100.26 MJ, respectively. And the total elastic energy accumulated inside the roof-floor is 474.56 MJ.

Figure 4. Load transfer and sudden jump displacement of multi-pillar. (a) Load-displacement curves of pillars, (b) Displacement-step curves of pillars.

2.2.3. CPF mode of multi-pillar

By drawing the load-displacement curve of multi-pillar, the CPF mode of multi-pillar at 848 m level of Alhada Lead-Zinc Mine can be obtained according to Dong et al. (2022b) and Dong et al. (2023). As shown in Figure 5, the multi-pillar experienced two times of instability, namely a-b and c-d. The first instability is the CPF of pillars P_2-3 and P_3-4, which is accompanied by the elastic rebound (86.66 mm) and energy release (254 MJ) of roof-floor. The second instability is the CPF of pillars P_1-2 and P_4-5, which is accompanied by the elastic rebound (119.82 mm) and energy release (486 MJ) of roof-floor. Thus, the multi-pillar instability consists of two instability events of double-pillar (i.e. the CPF a-b of double-pillar P_2-3-P_3-4, and the CPF c-d of double-pillar P_1-2-P_4-5). And there is also a loading process b-c between two instability events, so the two instabilities occur under different load levels. Therefore, the CPF mode of multi-pillar at 848 m level of the Alhada Lead-Zinc Mine can be summarized into the ‘2 + 2’ compound failure mode according to the theory of Dong et al. (2022b).

Figure 5. Load-displacement curves and stability analysis of multi-pillar system.

3. Theoretical analysis of multi-pillar instability

In this section, the theoretical model of rebound overload mechanism and potential function that can characterize the evolution characteristics of roof-multi-pillar-floor system are established based on the catastrophe theory. The analytical expressions of sudden jump and energy release of roof-multi-pillar-floor system are derived in this section.

3.1. Model setup

3.1.1. Simplified mechanical model

Figure 6 shows the mechanical model of roof-multi-pillar-floor system. The stiffness of roof-floor is K. The loading displacement of multi-pillar is u. The loading displacement from surrounding rockmass is a. Thus, the difference between a and u is the deformation of roof-floor. The constitutive equations of single pillars P_1-2, P_2-3, P_3-4 and P_4-5 in roof-multi-pillar-floor system can be expressed as follows (Wang et al. 2006): (3) $$f_i\left(u\right)=k_{i}u{e}^{-{\left(\frac{u}{u_i}\right)}^{m_i}}$$(3) where i = 1-2, 2-3, 3-4 and 4-5. Thus, k_1-2, k_2-3, k_3-4 and k_4-5 are the initial stiffnesses of P_1-2, P_2-3, P_3-4 and P_4-5; u_1-2, u_2-3, u_3-4 and u_4-5 are the load displacement corresponding to the peak strength of P_1-2, P_2-3, P_3-4 and P_4-5; m_1-2, m_2-3, m_3-4 and m_4-5 are the homogeneity coefficients of P_1-2, P_2-3, P_3-4 and P_4-5. These parameters can be determined by the uniaxial compression test and parameter estimation; and it is more convenient to determine the parameters by drawing according to the test curves.

Figure 6. Mechanical model of roof-multi-pillar-floor system.

3.1.2. Instability evolution process

Figure 7 shows the schematic diagram of instability evolution process of multi-pillar (P_1-2, P_2-3, P_3-4 and P_4-5) and elastic rebound of roof-floor according to the field observation and numerical simulation results in Section 2. The numerical simulation results conclude that the CPF mode of multi-pillar is the ‘2 + 2’ compound failure mode, which include (1) CPF of P_2-3-P_3-4: P_2-3 and P_3-4 fail simultaneously at step 2 with an elastic rebound of Δu_2-4, because the load transfer from P_2-3 exceeds the bearing capacity of P_3-4; and (2) CPF of P_1-2-P_4-5: With the continued loading of surrounding rockmass from step 2 to step 3, P_1-2 and P_4-5 fail simultaneously at step 4 with an elastic rebound of Δu_1-5.

Figure 7. Schematic diagram of instability evolution process of roof-multi-pillar-floor system.

3.2. Catastrophe analysis of multi-pillar

3.2.1. CPF of P_2-3-P_3-4

For CPF of P_2-3-P_3-4 from steps 0 to 2 shown in Figure 7, its potential function is the total strain energy expression of the roof-multi-pillar-floor system, which is calculated as (4) $$W_{23-34}\left(u\right)={\int }_0^u{F}_{23-34}\left(u\right)\mathit{\text{du}}+\frac{1}{2}\left(k_{1-2}+k_{4-5}\right)u^2+\frac{1}{2}K{\left(a_{2-3}-u\right)}^2$$(4) where $F_{23-34}\left(u\right)=f_{2-3}\left(u\right)+f_{3-4}\left(u\right).$ Let $\frac{d{W}_{23-34}}{\mathit{\text{du}}}=0$ and the equilibrium surface equation is given by (5) $$W_{23-34}^{\prime }\left(u\right)=F_{23-34}\left(u\right)+\left(k_{1-2}+k_{4-5}\right)u-K\left(a_{2-3}-u\right)=0$$(5)

Let $\frac{d^2{W}_{23-34}}{d{u}^2}=0$ and the equation of bifurcation set can be expressed as (6) $$W_{23-34}^{″}\left(u\right)=F_{23-34}^{\prime }\left(u\right)+k_{1-2}+k_{4-5}+K=0$$(6) where $F_{23-34}^{\prime }\left(u\right)=k_{2-3}\left[1-m_{2-3}{\left(\frac{u}{u_{2-3}}\right)}^{m_{2-3}}\right]e^{-{\left(\frac{u}{u_{2-3}}\right)}^{m_{2-3}}}+k_{3-4}\left[1-m_{3-4}{\left(\frac{u}{u_{3-4}}\right)}^{m_{3-4}}\right]e^{-{\left(\frac{u}{u_{3-4}}\right)}^{m_{3-4}}}.$ According to the smoothness property of equilibrium surface, let $\frac{d^3{W}_{23-34}}{d{u}^3}=0$ and the cusp point can be obtained as (7) $$W_{23-34}^{″\prime }\left(u\right)=F_{23-34}^{″}\left(u\right)=k_{2-3}{e}^{-{\left(\frac{u}{u_{2-3}}\right)}^{m_{2-3}}}\left(\frac{m_{2-3}^2{u}^{2m_{2-3}-1}}{u_{2-3}^{2m_{2-3}}}-\frac{m_{2-3}{u}^{m_{2-3}-1}}{u_{2-3}^{m_{2-3}}}-\frac{m_{2-3}^2{u}^{m_{2-3}-1}}{u_{2-3}^{m_{2-3}}}\right)+k_{3-4}{e}^{-{\left(\frac{u}{u_{3-4}}\right)}^{m_{3-4}}}\left(\frac{m_{3-4}^2{u}^{2m_{3-4}-1}}{u_{3-4}^{2m_{3-4}}}-\frac{m_{3-4}{u}^{m_{3-4}-1}}{u_{3-4}^{m_{3-4}}}-\frac{m_{3-4}^2{u}^{m_{3-4}-1}}{u_{3-4}^{m_{3-4}}}\right)=0$$(7) where the absolute value of slope at cusp point $u_{t_{23-34}}$ of multi-pillar is $k_{t_{23-34}}=-F_{23-34}^{\prime }\left(u_{t_{23-34}}\right)-(k_{1-2}+k_{4-5}).$ It should be noted that the analytical solution of cusp point $u_{t_{23-34}}$ cannot be obtained due to the complexity of Equation (7)(7) $$W_{23-34}^{″\prime }\left(u\right)=F_{23-34}^{″}\left(u\right)=k_{2-3}{e}^{-{\left(\frac{u}{u_{2-3}}\right)}^{m_{2-3}}}\left(\frac{m_{2-3}^2{u}^{2m_{2-3}-1}}{u_{2-3}^{2m_{2-3}}}-\frac{m_{2-3}{u}^{m_{2-3}-1}}{u_{2-3}^{m_{2-3}}}-\frac{m_{2-3}^2{u}^{m_{2-3}-1}}{u_{2-3}^{m_{2-3}}}\right)+k_{3-4}{e}^{-{\left(\frac{u}{u_{3-4}}\right)}^{m_{3-4}}}\left(\frac{m_{3-4}^2{u}^{2m_{3-4}-1}}{u_{3-4}^{2m_{3-4}}}-\frac{m_{3-4}{u}^{m_{3-4}-1}}{u_{3-4}^{m_{3-4}}}-\frac{m_{3-4}^2{u}^{m_{3-4}-1}}{u_{3-4}^{m_{3-4}}}\right)=0$$(7) , the numerical solution ($u_{t_{23-34}}$ = 342.8 mm) is used.

Making Taylor series expansion with respect to cusp point $u_{t_{23-34}}$ for Equation (5)(5) $$W_{23-34}^{\prime }\left(u\right)=F_{23-34}\left(u\right)+\left(k_{1-2}+k_{4-5}\right)u-K\left(a_{2-3}-u\right)=0$$(5) , and introduction of the dimensionless variable $x=\frac{u-u_{t_{23-34}}}{u_{t_{23-34}}},$ Equation (5)(5) $$W_{23-34}^{\prime }\left(u\right)=F_{23-34}\left(u\right)+\left(k_{1-2}+k_{4-5}\right)u-K\left(a_{2-3}-u\right)=0$$(5) can be transformed into the standard form of equilibrium surface equation as (8) $$x^3+p_{23-34}x+q_{23-34}=0$$(8)

The parameters $p_{23-34}$ and $q_{23-34}$ in Equation (8)(8) $$x^3+p_{23-34}x+q_{23-34}=0$$(8) are derived as (9) $$\{\begin{array}{c}{p}_{23-34}=\frac{6\left[F_{23-34}^{\prime }\left(u_{t_{23-34}}\right)+k_{1-2}+k_{4-5}+K\right]}{F_{23-34}^{″\prime }\left(u_{t_{23-34}}\right)u_{t_{23-34}}^2} \\ q_{23-34}=\frac{6\left[F_{23-34}\left(u_{t_{23-34}}\right)-K\left(a_{2-3}-u_{t_{23-34}}\right)\right]}{F_{23-34}^{″\prime }\left(u_{t_{23-34}}\right)u_{t_{23-34}}^3}\end{array}$$(9)

According to the catastrophe theory, only when $p_{23-34}\le 0$ can the bifurcation set be crossed and then cause a catastrophic action. Thus, the instability criterion of roof-multi-pillar-floor system can be deduced as (10) $$K\le k_{t_{23-34}}\mathit{\text{ or }}{k}_{1-2}+k_{4-5}+K\le -F_{23-34}^{\prime }\left(u_{t_{23-34}}\right)$$(10)

When $p_{23-34}\le 0,$ the three real roots of Equation (8)(8) $$x^3+p_{23-34}x+q_{23-34}=0$$(8) are obtained as (11) $$\{\begin{array}{c}{x}_1={2\left(\frac{-p_{23-34}}{3}\right)}^{\frac{1}{2}} \\ x_2=x_3=-{\left(\frac{-p_{23-34}}{3}\right)}^{\frac{1}{2}}\end{array}$$(11)

The dimensionless sudden jump Δx_23-34 and the dimensionless elastic energy release ΔU_23-34 of multi-pillar is obtained as (12) $$\{\begin{array}{c}\Delta x_{23-34}=x_1-x_2={\left(-3p_{23-34}\right)}^{\frac{1}{2}}\\ \Delta U_{23-34}={\int }_{x_2}^{x_1}\left(x^3+p_{23-34}x+q_{23-34}\right)\mathit{\text{dx}}=\frac{7}{12}{p}_{23-34}^2\end{array}$$(12)

Thus, the theoretical sudden jump Δu_23-34 and the theoretical elastic energy release ΔW_23-34 of the roof-multi-pillar-floor system (caused by the collapse of P_2-3 and the elastic rebound of roof-floor) are calculated as (13) $$\{\begin{array}{c}\Delta u_{23-34}=u_{t_{23-34}}\Delta x_{23-34}=u_{t_{23-34}}{\left(-3p_{23-34}\right)}^{\frac{1}{2}}\\ \Delta W_{23-34}=\frac{7}{288}{F}_{23-34}^{″\prime }\left(u_{t_{23-34}}\right)u_{t_{23-34}}^3{p}_{23-34}^2\end{array}$$(13)

3.2.2. CPF of P_1-2-P_4-5

For CPF of P_1-2-P_4-5 from steps 2 to 4 shown in Figure 7, its potential function is the total strain energy expression of the roof-multi-pillar-floor system, which is calculated as (14) $$W_{12-45}\left(u\right)={\int }_0^u{F}_{12-45}\left(u\right)\mathit{\text{du}}+\frac{1}{2}K{\left(a_{1-2}-u\right)}^2$$(14) where $F_{12-45}\left(u\right)=f_{1-2}\left(u\right)+f_{4-5}\left(u\right).$ Let $\frac{d{W}_{12-45}}{\mathit{\text{du}}}=0$ and the equilibrium surface equation is given by (15) $$W_{12-45}^{\prime }\left(u\right)=F_{12-45}\left(u\right)-K\left(a_{1-2}-u\right)=0$$(15)

Let $\frac{d^2{W}_{12-45}}{d{u}^2}=0$ and the equation of bifurcation set can be expressed as (16) $$W_{12-45}^{″}\left(u\right)=F_{12-45}^{\prime }\left(u\right)+K=0$$(16) where $F_{12-45}^{\prime }\left(u\right)=k_{1-2}\left[1-m_{1-2}{\left(\frac{u}{u_{1-2}}\right)}^{m_{1-2}}\right]e^{-{\left(\frac{u}{u_{1-2}}\right)}^{m_{1-2}}}+k_{4-5}\left[1-m_{4-5}{\left(\frac{u}{u_{4-5}}\right)}^{m_{4-5}}\right]e^{-{\left(\frac{u}{u_{4-5}}\right)}^{m_{4-5}}}.$ According to the smoothness property of equilibrium surface, let $\frac{d^3{W}_{12-45}}{d{u}^3}=0$ and the cusp point can be obtained as (17) $$W_{12-45}^{″\prime }\left(u\right)=F_{12-45}^{″}\left(u\right)=k_{1-2}{e}^{-{\left(\frac{u}{u_{1-2}}\right)}^{m_{1-2}}}\left(\frac{m_{1-2}^2{u}^{2m_{1-2}-1}}{u_{1-2}^{2m_{1-2}}}-\frac{m_{1-2}{u}^{m_{1-2}-1}}{u_{1-2}^{m_{1-2}}}-\frac{m_{1-2}^2{u}^{m_{1-2}-1}}{u_{1-2}^{m_{1-2}}}\right)+k_{4-5}{e}^{-{\left(\frac{u}{u_{4-5}}\right)}^{m_{4-5}}}\left(\frac{m_{4-5}^2{u}^{2m_{4-5}-1}}{u_{4-5}^{2m_{4-5}}}-\frac{m_{4-5}{u}^{m_{4-5}-1}}{u_{4-5}^{m_{4-5}}}-\frac{m_{4-5}^2{u}^{m_{4-5}-1}}{u_{4-5}^{m_{4-5}}}\right)=0$$(17) where the absolute value of slope at cusp point $u_{t_{12-45}}$ of multi-pillar is $k_{t_{12-45}}=-F_{12-45}^{\prime }\left(u_{t_{12-45}}\right).$ It should be noted that the analytical solution of cusp point $u_{t_{12-45}}$ cannot be obtained due to the complexity of Equation (17)(17) $$W_{12-45}^{″\prime }\left(u\right)=F_{12-45}^{″}\left(u\right)=k_{1-2}{e}^{-{\left(\frac{u}{u_{1-2}}\right)}^{m_{1-2}}}\left(\frac{m_{1-2}^2{u}^{2m_{1-2}-1}}{u_{1-2}^{2m_{1-2}}}-\frac{m_{1-2}{u}^{m_{1-2}-1}}{u_{1-2}^{m_{1-2}}}-\frac{m_{1-2}^2{u}^{m_{1-2}-1}}{u_{1-2}^{m_{1-2}}}\right)+k_{4-5}{e}^{-{\left(\frac{u}{u_{4-5}}\right)}^{m_{4-5}}}\left(\frac{m_{4-5}^2{u}^{2m_{4-5}-1}}{u_{4-5}^{2m_{4-5}}}-\frac{m_{4-5}{u}^{m_{4-5}-1}}{u_{4-5}^{m_{4-5}}}-\frac{m_{4-5}^2{u}^{m_{4-5}-1}}{u_{4-5}^{m_{4-5}}}\right)=0$$(17) , the numerical solution ($u_{t_{12-45}}$ = 389.9 mm) is used.

Making Taylor series expansion with respect to cusp point $u_{t_{12-45}}$ for Equation (15)(15) $$W_{12-45}^{\prime }\left(u\right)=F_{12-45}\left(u\right)-K\left(a_{1-2}-u\right)=0$$(15) , and introduction of the dimensionless variable $x=\frac{u-u_{t_{12-45}}}{u_{t_{12-45}}},$ Equation (15)(15) $$W_{12-45}^{\prime }\left(u\right)=F_{12-45}\left(u\right)-K\left(a_{1-2}-u\right)=0$$(15) can be transformed into the standard form of equilibrium surface equation as (18) $$x^3+p_{12-45}x+q_{12-45}=0$$(18)

The parameters $p_{12-45}$ and $q_{12-45}$ in Equation (18)(18) $$x^3+p_{12-45}x+q_{12-45}=0$$(18) are derived as (19) $$\{\begin{array}{c}{p}_{12-45}=\frac{6\left[F_{12-45}^{\prime }\left(u_{t_{12-45}}\right)+K\right]}{F_{12-45}^{″\prime }\left(u_{t_{12-45}}\right)u_{t_{12-45}}^2}\\ q_{12-45}=\frac{6\left[F_{12-45}\left(u_{t_{12-45}}\right)-K\left(a_{1-2}-u_{t_{12-45}}\right)\right]}{F_{12-45}^{″\prime }\left(u_{t_{12-45}}\right)u_{t_{12-45}}^3}\end{array}$$(19)

According to the catastrophe theory, only when $p_{12-45}\le 0$ can the bifurcation set be crossed and then cause a catastrophic action. Thus, the instability criterion of roof-multi-pillar-floor system can be deduced as (20) $$K\le k_{t_{12-45}}$$(20)

Similarly, the theoretical sudden jump Δu_12-45 (Equation (21)(21) $$\Delta u_{12-45}=u_{t_{12-45}}{\left(-3p_{12-45}\right)}^{\frac{1}{2}}$$(21) ) and the theoretical elastic energy release ΔW_12-45 (Equation (22)(22) $$\Delta W_{12-45}=\frac{7}{288}{F}_{12-45}^{″\prime }\left(u_{t_{12-45}}\right)u_{t_{12-45}}^3{p}_{12-45}^2$$(22) ) within the roof-multi-pillar-floor system can be deduced based on the established standard form of equilibrium surface equation. (21) $$\Delta u_{12-45}=u_{t_{12-45}}{\left(-3p_{12-45}\right)}^{\frac{1}{2}}$$(21) (22) $$\Delta W_{12-45}=\frac{7}{288}{F}_{12-45}^{″\prime }\left(u_{t_{12-45}}\right)u_{t_{12-45}}^3{p}_{12-45}^2$$(22)

4 Analysis of results

4.1. Comparison of results

Based on the above theoretical analysis, the theoretical values of sudden jump displacement and energy release of roof-multi-pillar-floor system are compared with the numerical simulation results, as shown in Figure 8. Figure 8(a) shows the comparison results of sudden jump of roof-multi-pillar-floor system, and it should be noted that the sudden jump numerical simulation result of CPF of P_2-3 and P_3-4 is larger than the theoretical value, because the tensile failure of rockmass leads to a greater rebound of roof-floor. Figure 8(b) shows the comparison results of energy release of roof-multi-pillar-floor system. In this respect, both numerical results are in good agreement with theoretical results. In summary, the numerical simulation and theoretical results can well reproduce the CPF characteristics of multi-pillar system, and accurately describe the sudden jump, energy release and load transfer behavior in multi-pillar compound failure mode.

Figure 8. Comparison between numerical and theoretical results. (a) Sudden jump of roof-multi-pillar-floor system, (b) energy release of roof-multi-pillar-floor system.

4.2. Load transfer effect analysis

Figure 9(a) and (b) show the detailed cascading failure and load transfer behavior of P_2-3-P_3-4 and P_1-2-P_4-5, respectively. Figure 9(a) shows domino failure mode of P_2-3 and P_3-4 from step-in-step 96-1 to 96-11 by amplifying the cascading failure process of step 96. It can be seen from Figure 9(a) that the P_2-3 instability at step-in-step 96-6 causes the elastic rebound of roof-floor, inducing the load increases in adjacent pillars P_1-2, P_3-4 and P_4-5 (increasing by 6.2%, 4.7% and 3.8%, respectively); Moreover, the induced load transfer at step-in-step 96-6 exceeded the bearing capacity of P_3-4, resulting in the domino failure behavior of P_2-3 and P_3-4. Subsequently, the P_3-4 instability at step-in-step 96-11 causes the elastic rebound of roof-floor, resulting in the load increases in adjacent pillars P_1-2 and P_4-5 (increasing by 7.8% and 7.1%, respectively). Figure 9(b) shows the simultaneous failure of P_1-2 and P_4-5 from step-in-step 106-1 to 106-7 by amplifying the cascading failure process of step 106. It can be seen from Figure 9(b) that P_1-2 and P_4-5 fail simultaneously from step-in-step 106-1 to 106-7 without causing a significant load increases in adjacent pillars.

Figure 9. Cascading pillar failure and load transfer effect analysis.

Above all, it is concluded that: (1) The load transfer effect induced by pillar instability may cause sudden load increases in adjacent pillars under elastic foundation. (2) The influence range of elastic rebound of roof-floor (or load transfer effect) is limited and inversely proportional to pillar spacing. For example, when P_2-3 fails at step-in-step 96-6, the load increase of P_3-4 is 4.7%, which is 0.9% more than that of P_4-5 farther away from P_2-3; and the load increase of P_4-5 caused by the P_3-4 instability at step-in-step 96-11 is 7.8%, which is 0.7% more than that of P_1-2 farther away from P_3-4, as shown in Figure 9(a). (3) The load transfer effect is also affected by the damage fracture zone. For example, the damage fracture zone caused by the tensile stress in the roof-floor, as well as the long load transfer distance or pillar spacing, hinder the load transfer behavior between P_1-2 and P_4-5, as shown in Figure 9(b).

5 Conclusions

This paper studies the instability and load transfer process of roof-multi-pillar-floor system at 848 m level of Alhada Lead-Zinc Mine based on the RFPA^2D numerical simulation software. And the CPF mode of multi-pillar system is analyzed based on the stiffness theory. The conclusions can be summarized as follows:

1. The dynamic instability evolution process and CPF characteristics (e.g. elastic rebound and energy release) of multi-pillar and goaf group are simulated and analyzed based on the RFPA^2D numerical simulation software. The simulation results show that the multi-pillar at 848m level of Alhada Lead-Zinc Mine experiences a ‘2 + 2’ compound failure mode including two CPF events. The first CPF of P_2-3-P_3-4 is characterized by the elastic rebound of 86.66 mm and energy release of 254 MJ of local roof-floor; while the second CPF of P_1-2-P_4-5 is characterized by the elastic rebound of 119.82 mm and energy release of 486 MJ of local roof-floor, which verify the field observation results.

2. The theoretical model of rebound overload mechanism of roof-multi-pillar-floor system is established to derive the analytical expressions of stress criterion, elastic rebound and energy release based on catastrophe theory. When the LMS of roof-floor is greater than the post-peak stiffness at cusp point, the pillar will fail progressively; conversely, the pillar will occur unstable failure. The theoretical results are in good agreement with numerical results, and both can well reproduce the CPF characteristics of multi-pillar system and accurately describe the sudden jump, energy release and load transfer behavior of roof-multi-pillar-floor system.

3. The influencing factors of load transfer effect such as pillar spacing and damage fracture zone are quantitatively analyzed. The analytical results show that the influence range of load transfer effect is inversely proportional to the distance between pillars. When P_2-3 fails at step 96-6, the load increase of P_3-4 is 4.7%, which is 0.9% more than that of P_4-5 farther away from P_2-3; and the load increase of P_4-5 caused by P_3-4 instability is 7.8%, which is 0.7% more than that of P_1-2 farther away from P_3-4. The load transfer effect is also affected by the damage fracture zone. The tension fractures generated within the roof-floor hinder the load transfer behavior between P_1-2 and P_4-5.

List of symbols
P1-2	=	Pillar between goafs 1# and 2#
P2-3	=	Pillar between goafs 2# and 3#
P3-4	=	Pillar between goafs 3# and 4#
P4-5	=	Pillar between goafs 4# and 5#
fmeso	=	Mesoscopic parameters of strength
Emeso	=	Mesoscopic parameters of elastic modulus
fmacro	=	Macroscopic parameters of strength
Emacro	=	Macroscopic parameters of elastic modulus
m	=	Homogeneity coefficient
CPF	=	Cascading pillar failure
LMS	=	Local mine stiffness
K1-2, K2-3, K3-4, K4-5	=	LMS of P1-2, P2-3, P3-4 and P4-5
d1-2, d2-3, d3-4, d4-5	=	Deformation of P1-2, P2-3, P3-4 and P4-5
f1-2, f2-3, f3-4, f4-5	=	Load of P1-2, P2-3, P3-4 and P4-5
d’1-2, d’2-3, d’3-4, d’4-5	=	Deformation of imaginary P1-2, P2-3, P3-4 and P4-5
f’1-2, f’2-3, f’3-4, f’4-5	=	Load of imaginary P1-2, P2-3, P3-4 and P4-5
kp1-2, kp2-3, kp3-4, kp4-5	=	Post-peak stiffness of P1-2, P2-3, P3-4 and P4-5
W	=	Elastic energy release of roof-floor
K	=	LMS of roof-floor
Δd	=	Elastic rebound of roof-floor
u	=	Loading displacement of multi-pillar
a	=	Loading displacement of surrounding rockmass
k1-2, k2-3, k3-4, k4-5	=	Initial stiffnesses of P1-2, P2-3, P3-4 and P4-5
u1-2, u2-3, u3-4, u4-5	=	Load displacement of peak strength of P1-2, P2-3, P3-4 and P4-5
m1-2, m2-3, m3-4, m4-5	=	Homogeneity coefficients of P1-2, P2-3, P3-4 and P4-5
ut12-45, ut23-34	=	Cusp point of roof-multi-pillar-floor system
Δu12-45, Δu23-34	=	Theoretical sudden jump value of roof-multi-pillar-floor system
ΔW12-45, ΔW23-34	=	Theoretical elastic energy release of roof-multi-pillar-floor system

Acknowledgements

Great appreciation goes to the editorial board and the reviewers of this paper.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The datasets analyzed during the current study are available from the corresponding author on reasonable request.

Additional information

Funding

This work was financially supported by the Enlisting and Leading Project of the Key Scientific and Technological Innovation in Heilongjiang Province, China (Grant No. 2021ZXJ02A03, 2021ZXJ02A04), the Youth Fund of the Natural Science Foundation of Henan Province (Grant No. 242300421457), the National Natural Science Foundation of China (Grant No. 42277478), and the North China University of Water Resources and Electric Power Launch Fund for High-level Talents Research (Grant No. 40937).