The environmental conditions due to unreasonable mining in underground stopes, the slurry diffusion mechanism in the grouting reinforcement of a stope within its influence, the causes of large-scale instability collapse, and the catastrophic stope process are analyzed, and limit upper line analysis theory and numerical analysis methods are comprehensively adopted, revealing the continuous catastrophic collapse mode of large-scale underground stopes. The method of determining the stope instability collapse boundary and the slip surface within the range based on the theory of the maximum shear strain increment is proposed, and the diffusion radius and range of the grouting slurry during the reinforcement process, which considers the multifield coupling factors, are obtained. The results show that the U-shaped hidden danger area formed after the collapse of the large-scale underground stope. The influence range reaches six adjacent stopes, which are symmetrically distributed around the collapse; the mining instability is manifested as a catastrophic chain process of stress change, energy accumulation, state change, and instability collapse. The damage mode of instability collapse is a combination method of wedge collapse, circular arc rotation, triangular translation, and strip slip. According to the multiphysics coupling numerical calculation, the diffusion radius of the grouting slurry is 12 m, exhibiting an elliptical distribution. The research results can be used to comprehensively control the underground mining environment, thus effectively solving the safety problems faced by tunnel or roadway excavations above the goaf.

1. Introduction

Mineral resources are one of the most basic sources of human production and living materials and are an important material for economic development. In addition, most mineral resources are nonrenewable [1]; thus, the rational and effective application of mineral resources is a societal responsibility. With the rapid development of China’s social economy, the consumption of mineral resources has increased, and mineral resources with high grades and easy exploitation have also been depleted [2, 3]. To meet the growing demand of economic construction for mineral resources and to maintain stable economic growth, the exploitation of mineral resources should not only continue to extend deeper but also fully recover the hidden resources left in the mining process under the existing technical and economic conditions [4, 5]. Given the unreasonable mining that took place in the early stages of mining, a series of complex engineering problems, such as the subsidence of surrounding rock, structural failure of rock masses, and instability of roadways, the results from the collapse of underground goafs, cause substantial hidden dangers to the subsequent mining and repair processes of underground roadway spaces [6]. Problems such as the determination of the slip boundary of the underground collapse zone and its influence range [710], deformation of the upper ore body of the collapsed zone and two areas of surrounding rock (filling body), evaluation and prediction of subsidence and movement [1013], and reinforcement and repair of surrounding rock in the roadway within the scope of the collapse zone [1416] all affect the selection and optimization of the hidden resource recovery plan in the entire underground collapse zone and the budget of the project investment.

Determining the influence range of mining instability in complex large-scale underground stopes and providing the corresponding reinforcement is a key link for safe mining of hidden dangerous resources and is crucial to mining safety [17]. The main research methods for the instability and collapse of underground stopes include numerical simulation analyses, the catastrophe theory, similar material simulations, and the uncertainty method [1823]. Yao, followed by Wood [24], studied the conditions and heights of overlying rock collapse caused by stoping through the finite-element method and boundary-element method [25]. Freidin et al., based on finite-element software, conducted stress distribution and stability analysis studies on goaf roofs and obtained the safe range of underground mining spans [7]. Peralta studied the stability of fluid-filled roadways under mining influences by the perturbation method, and the change pattern of impacts exerted by different dynamic waves on existing structures was analyzed [26]. Grgic studied the mechanical behavior of iron ore and marl through rock mechanics tests and FLAC 2D software, and Malan proposed an elastic-viscoplastic method to study the time dependence of the creep deformation of an ore rock fault zone during surface excavation works [27]. With respect to the stoping of resources with hidden dangers and the countermeasures of hidden dangers, developed countries studied the exploitation of resources containing hidden dangers in underground mines much earlier than China. By adopting a completely different concept, residual ore was fully and effectively recovered from the San Manuel copper mine and White Pine copper mine in the United States by using the dissolving leaching mining method, which was safe and had low mining costs [28]. Tan studied the stress distribution pattern of semi-U-shaped valley terrain [29]. In addition, research on stope/rock mass instability caused by the failure of rock discontinuities (mechanical joints, bedding, or veins) is also common in rock engineering. For example, Junlong Shang et al. studied the fracture of vein-like granite in multiaxial compression using a discrete element method model [30]. An engineering method was proposed by Grenon for the analysis and reinforcement design of open stopes in jointed rocks [31]. Shang et al. discussed the tensile strength of the initial discontinuity in rock, and a result from a laboratory test procedure was proposed to quantify this parameter [32]. For the above research, the choice of a specific method is determined by a series of trade-offs between technology and practice.

The collapse and instability mechanism of underground stopes is complex; thus, it is crucial to analyze the causes of instability and collapse of underground stopes for mining safety, to determine the catastrophic mechanism of instability and collapse area, to determine the slip surface, to carry out effective grouting reinforcement, and to control the hidden danger area. Meanwhile, there are no general approaches or methods for reference. The collapse instability mechanism of the underground stope mentioned above is very complicated, with no fixed method for a recovery plan. There is still no systematic theoretical research on the methods and means of repair, or the scope and depth of repair and reinforcement, with no standard available to be followed.

Based on this condition, which is one of the cutting-edge topics in sustainable mining development under the condition of increasingly exhausted resources as the mechanism of large-scale instability and collapse of underground stopes under complex environmental conditions, the reinforcement and restoration of the geological environment in the collapsed area are carried out to realize the safe and efficient mining of hidden danger resources in underground collapsed areas.

In this paper, a method to determine the sliding surface of a stope collapse boundary and scope based on the maximum shear strain increment theory is proposed to improve mining safety and reveal the mechanical behavior mechanism of collapse and control risks, which is combined with a metal ore mine as a case study, in view of the complex mining conditions of the collapsed area. The theory of the upper-bound theorem of limit analysis and rock pressure and rock mechanics theory are used, and the mechanical behavior mechanism of collapse and the diffusion mechanism of grouting in the collapse area are revealed to obtain the causes and catastrophic process of large-scale instability and collapse of stopes. Then, based on the Richard model, the grout diffusion trend under time series is explored to guide roadway or tunnel grouting reinforcement and modification work in the future.

2. Current Situation of Instability in the Collapse Area of an Underground Stope

Based on the unreasonable metal mining practices, the collapse of the stope forms a large-scale collapse area, whose elevation is −298.3 m to −232.0 m, and the coordinates are 313.5 to 396.5 in the north-south direction (see Figure 1). The area comprises a goaf, filling body, and ore body, and the width of the stope is 6.5–8 m. After the collapse of the stope, which is approximately 22.5 m in height at the lower part of the #0 stope that is completely explored, two bangs of fillings rush into the goaf due to the collapse, and the upper ore body sinks as a whole. According to the site survey, the mining environment conditions of the ore bodies in the area are very complicated.

3. Analysis of the Causes of Instability Collapse

The large-scale instability and collapse of the underground stope not only cause continuous movement of the upper part of the goaf and the two areas but also produce discontinuous damage. The discontinuous collapse affects the adjacent stope. The mining methods, mining time, and filling materials of different stopes differ, and the repeated mining of each stope makes the overall trend of the rock mass around the collapsed area move continuously toward the goaf. In the position with considerable uneven sinking, the “wedge-shaped” crack is likely to occur. Figure 2 shows the specific causes of instability collapse.

After the formation of the underground goaf, the collapse of the two-package body is induced, resulting in crack development at the stope interface until a local part or complete penetration. The filling body slides and moves along the surface of the fracture structure. The process of mining collapse and catastrophe also refers to the process of the gradual accumulation of energy, that is, the process from quantitative change to qualitative transformation. When the energy accumulation reaches a certain value, the filling body is destroyed, and a large crack is formed at the unstable position, finally promoting large-scale instability collapse, whose process can be expressed by the catastrophic chain process shown in Figure 3.

4. Upper-Bound Theorem of the Limit Analysis of Instability Collapse

When the method of the upper-bound theorem of limit analysis is applied in engineering practice, the stress field and ultimate load under the failure state can be determined through the equilibrium condition, yield condition, flow rule, and corresponding boundary condition with reasonable construction of the failure mechanism. This method can easily form explicit expressions of objective functions, is simple and fast in solving processes, and has been fully recognized and widely used [33].

The underground stope collapse is analyzed based on the upper-bound theorem of limit analysis (the energy method). By constructing the displacement field and velocity field, the analytical expression of the damage load is obtained based on the principle that the internal and external powers are equal [34, 35]. Figure 4 shows the expanded view of the velocity vector.

The external power calculation is as follows [36]:

In the above formula, refers to the initial velocity vector, are the respective components of the velocity vector, denote the angles, represents the gravity of the rock mass (), and is the goaf height ().

Therefore,where denotes the internal friction angle of the rock mass (°).

The calculation of the power of the bearing reaction is as follows:where refers to the support force on the unit area of the roof (), specifies the goaf span (), and is the support force provided for the unraised ore:where is the side pressure coefficient and , , , and are angles (°).

The dissipation power of the internal energy for GCk is as follows:where represents the cohesion of the rock mass () and , , and represent the energy dissipation of the shear triangles.

The following equation can be obtained according to the above simultaneous equations:where and are coefficients in the power of the work due to gravity.

The stope outside the collapse area can be explained by the strip method [37]:where Deq represents the equivalent goaf height, represents the angle between the slip surface and the horizontal plane, is the severity of the surrounding rock, is the average effective vertical stress acting on the top of the strip filling body, and represents the horizontal force of the strip-shaped filling body; when the physical and mechanical parameters are constant, is a function of (θ, π/2).

Figure 5 shows a schematic of the deformation mechanism of mining collapse. In the process of ore extraction, the supporting force of the two sides of the filling body and the stability of the surrounding rock gradually decrease. The No. 1 and S0-1 stopes of the two sides of the 0 stope and the internal backfill slag present a wedge-shaped collapse mode, and the area along D1, D2... Dn presents a circular arc-shaped rotation; that is, an arc-shaped slip surface is formed until point B1, and the ore body above the goaf and a part of the filling body of the two sides will show a triangular translational pattern. When the force q provided by the roof cannot continue to support the upper ore body, the filling bodies in each stope will show a stripe slip along the interface between them and the slip surface of C1 ∼ Cn, D1 ∼ Dn, and the distance the farther away the goaf is, the smaller the slip angle, that is, θ1 < θ2... < θn in the figure for S0-1, S1-1, S1-2, N0-1, and N1-2. The N2-3 stope slides successively, and at the same time, the upper ore body will sink or slide along with the movement of the lower filling body, which will eventually cause a large area of collapse. The mechanism of instability and collapse of underground stopes is explained from the perspective of energy dissipation and energy transfer. The system consisting of each stope in the ore body is a subsystem composed of the ore and filling body. The instability collapse is essentially the process of energy conversion, transfer, and dissipation of the subsystem during mining. The final energy conversion of the system mainly includes the work due to gravity, the work due to external forces, and the internal energy dissipation. When there is work due to gravity, energy dissipation occurs on the slip surface. According to the above derivation, the work due to gravity mainly includes Ws, , and Wn, and the external force is WT. That is, the energy released by the external body after the absorption of the external force in the process of mining and the critical value exceed e and q, respectively, and instability occurs. At the same time, the internal energy is consumed with rupture development. That is, the internal energy is lost along the ABG, FCnG, C3CnG, C2C3G, C1C2G, C0C1G, and BC0G lines until the initial equilibrium state is reached. In this process, energy is exchanged, restricted, and irreversibly fed back to the outside world. Afterward, gravity and external forces continue to work, and energy is transferred to the goaf. The internal energy is lost along the D0E, EF, FG, and GH lines and the C0I, IJ, and JK lines until the final equilibrium state is reached.

In summary, the deformation and failure mode of the filling body near the goaf of stope #0 is determined by the downwardly turning wedge. The deformation and failure mode of the surrounding rock of the stope positioned at a distance from the goaf is determined by the strip method. The failure mode of mining collapse of “wedge collapse” ⟶ “circular arc rotation” ⟶ “triangle translation” ⟶ “bar slip” is formed.

5. Searching Potential Failure Surface Based on the Maximum Shear Strain Increment

5.1. The Search Method for the Sliding Surface of the Unstable Collapse Area

The hybrid discrete method is used to simulate plastic failure and plastic flow in FLAC3D, which is more accurate and reasonable than the finite-element method. Even if simulating a static system, the dynamic equation is still used, which allows the simulation results of FLAC3D to more accurately reflect the real situation than that of other models.

According to the elastic-plastic mechanics theory and Mohr–Coulomb strength theory [38], the failure of rocks is caused by the shear stress on one side of the rock exceeding the ultimate stress, which inevitably leads to large shear deformation. The shear strain increment in FLAC3D is a physical quantity related to the node displacement; thus, choosing the position of the maximum shear strain increment as the criterion for determining the sliding surface has the advantages of providing a clear physical meaning, which is consistent with the conclusions of scholars. Therefore, the sliding surface can be considered to be connected by the position of the maximum shear strain increment in the vertical direction. On this basis, this paper presents a method to determine the potential slip surface by searching for the maximum shear strain increment.

To determine the slip surface and influence range of the mining instability area, an analysis model of the instability collapse state is constructed to obtain the mechanical state and maximum shear strain incremental cloud diagram of the failure area of the underground stope after the clarification of the slip surface search method based on the maximum shear strain increment. Then, the three-dimensional shape and influence range of the slip surface of the collapse area are determined.

The rock mass within the collapsed area is discretized into a series of vertical lines. From the definition of the above limit states, the rock mass on the slip surface is notably in a plastic state. Therefore, the unit that enters the plastic state is the position where the slip surface passes. That is, on any vertical line, several unit bodies reach a plastic state, and the maximum shear strain increment is selected here as the criterion. The tensor of the strain increment denotes the physical quantity associated with the node displacement [39, 40]. For the tetrahedral element component in the model, its calculation formula in an infinitely short time is shown in the following formula:where refers to the displacement in the direction of the first node of the unit; represents the normal vector of one of the faces of the tetrahedron corresponding to the first node of the unit; and corresponds to the area of one of the tetrahedrons corresponding to the first node. That is, during the time period , the node displacement of the unit determines the magnitude of the displacement increment of the strain increment in one time step.

The least square method is used to fit the curve, which can smooth the slip surface of the original form of fluctuation. Assuming that the three-dimensional coordinates of the above maximum shear strain increment are and , the form of the fitted curve is shown in the following formula:

5.2. Construction of the Instability and Collapse Mechanics Model of the Underground Stope

The three-dimensional numerical analysis model is constructed based on the current situation of ore mining, with a length of 200 m, a width of 100 m, a height of 100 m, a total of 187,500 units, and 199,576 nodes. The filling body in the model uses a double-yield constitutive model that simulates those with porous medium, whereas the rest is described by the Mohr–Coulomb model. The interface between the various stopes and filling bodies is connected by the interface command. Then, the material parameter assignment, meshing, boundary condition application, initial geostress generation, and numerical calculation analysis are completed. The ore and filling are generalized. The material is assumed to be isotropic and a continuous uniform medium without joints or cracks. The results of the mechanical response results within the collapse area are monitored by means of a layout of monitoring points. Table 1 shows the relevant parameters obtained by the rock mechanics experiments.

After the excavation of the underground space, the destruction of the ore body is caused, and the destruction of the ore body is represented by the size of the plastic zone. Figure 6 shows the distribution diagram of the plastic area formed after the excavation of stope #0. The figure shows the shear and tensile damage plastic zones with a large distribution range of the plastic zone. The goaf roof is represented by shear failure, and its two areas work jointly for shear and tensile damage. Moreover, the plastic area of the goaf top plate and the two areas expand and extend to the adjacent stope and form a connection, resulting in large-scale instability and goaf failure. Therefore, the two areas of filling bodies in the goaf are pulled down in the lateral direction, resulting in pieces and collapse. The roof also experiences large-scale falling, which eventually leads to large-scale instability and collapse of the entire empty area.

Figure 7 shows the maximum shear strain increment cloud diagram of the unstable collapse state. The figure shows the collapsed shape of the collapsed area. A larger maximum shear strain increment in the area results in more serious damage. However, the specific location and geometrical characteristics of the slip surface at the collapse area are still uncertain. The search method of the slip surface based on the maximum shear strain increment can be used to determine the position of the slip surface of the collapse area.

5.3. The Search Method of the Slip Surface in the Instability Collapse Area

The specific search method for the slip surface is as follows (see Figure 8):(1)The stress and strain data output by the software can be used to assess the penetration state of the plastic area.(2)A vertical line is set along the horizontal direction. The FISH language is deployed to search the coordinate position , where the maximum shear strain increment appears. n denotes the number of discrete points and is saved in a .txt format.

As the software can only obtain the maximum shear strain increment in the center of the unit and smooth the curve on the slip surface, a smooth collapse zone slip surface curve can be derived with the above search procedures.

According to the search method of the slip surface of mining instability, the influence range and three-dimensional shape of the mining collapse area can be delineated. The specific method is as follows:(1)FLAC3D can be used to calculate the stress and strain of the collapsed area in the limit state, and a series of vertical lines is set in the calculated section.(2)The FISH language programs are compiled, and the position of the maximum coordinate of the shear strain increment on the vertical line is searched to obtain the coordinates of the discrete points on the slip surface. The coordinates of the discrete points are fitted by the least squares method to obtain the shape and position of the slip surface.(3)Figure 9 shows the coordinate output of the shape and position of the slip surface, which is saved in a .txt format, and the three-dimensional slip surface is generated by SigmaPlot software.

The method of searching the potential failure surface based on the maximum shear strain increment has the advantage that the criterion for judging the sliding surface is clear in terms of its physical meaning, and the method is simple and practical.

6. The Reinforcement Repair Research of Grouting

Given the collapse of the large-scale underground stope, a large number of voids and fissures form in the area. These voids and fissures are usually filled and strengthened by grouting [41]. The key to grouting reinforcement and repair is to explore the diffusion law of the grouting slurry with time. Therefore, based on the Richard model, a three-dimensional multifield coupling model is constructed. The time step is changed, that is, the set time steps are 0 s, 60 s, 120 s, and 180 s. The study of the law of slurry diffusion is performed in the time sequence, as shown in Figure 10.

As shown in the figure below, the comparison after setting the time step shows that the diffusion proceeds upwards due to the initial pressure value. Then, the process enters into the rock mass fissure along the other pores and fills them. Thus, the slurry saturation directly below the slurry hole is the lowest. As time passes, given the effect of gravity, the slurry enters the entire area. Similarly, with the further passage of time, the concentration of the slurry reaches a certain value, the saturation diffusion diagram shows no substantial change, and the grouting time can be obtained. In actual engineering applications, the grouting numerical simulation can be performed by demand. If the project is in progress, the water-cement ratio of the material can be adjusted to change the density and viscosity of the material in the model. The optimum water-cement ratio of the material can be obtained through a simulation in the same area, effectively shortening the construction period and preventing slurry stagnation due to condensation that occurs too early and diffusion, which is beneficial to optimizing the grouting design.

Through the above simulation, the pores and cracks in the filling body that are injected into the original slurry injection are in a homogeneous state, and their diffusion radius is R0. After the slurry is filled and coagulated to become a “capsule”, the effective filling hole and fracture rate of the filling body are Nz, and the total volume of the injected slurry can be obtained fromwhere Wz denotes the total volume of the injected slurry (m3); R0 is the diffusion radius of the slurry (m); and Hz represents the total section height of the injected grouting (m). According to the grouting experience, the Nz value is 3.0% to 6.0%.

Figure 11 shows the relationship curve between the diffusion radius of the grouting and porosity rate. Remarkably, when the amount of slurry injected is constant, the range of the slurry diffusion increases with decreasing porosity. The results of the above research can be used to guide the grouting reinforcement modification work onsite.

7. Conclusion

To improve mining safety and reveal the collapse failure mode of large-scale instability in underground mining and to guide roadway or tunnel grouting reinforcement applications, this paper takes the collapse area of a metal mine as the research object, proposes a method for determining the potential slip surface by searching for the maximum shear strain increment, and obtains the radius and scope of grout diffusion in the solid treatment process.(1)The analysis and research on the mechanism of instability and collapse of underground stopes are carried out. An analysis model of the collapse mechanism is constructed for the complex mining conditions in the collapsed area. The analysis reveals the catastrophic chain process of the ore body dynamic instability as “stress change” ⟶ “energy accumulation” ⟶ “state change” ⟶ “stability collapse” and its collapse mode of “wedge collapse” ⟶ “circular rotation”⟶ “triangle parallel motion” ⟶ “bar sliding”.(2)The method of searching the failure surface based on the maximum shear strain increment is proposed. The influence range of the instability and collapse of the underground stope and the position of the slip surface are determined. A three-dimensional numerical model formed by the mechanical properties of the stope instability analysis is established to reveal the mechanical response characteristics of the collapse instability of the stope. The maximum shear strain slip surface search method is used to determine the three-dimensional appearance and spatial position of the slip surface of the collapse zone.(3)The change trend of the slurry with the time of grouting injection in the underground stope is received, and the relationship between the grouting diffusion radius and porosity rate is revealed to obtain the reinforcement and repair range of grouting in the collapsed area.

The research results can provide theoretical support for the further recovery of underground mineral resources in the collapsed area or the tunnels above the goaf.

Data Availability

All the data included in this study are available upon request by contact with the corresponding author.

Conflicts of Interest

The authors declare that they have no conflicts of interest.


This work was financially supported by the Hunan Provincial Department of Education General Project (no. 19C1744), Hunan Province Science Foundation for Youth Scholars (no. 2018JJ3510), Beijing Natural Science Foundation (no.9194027), China Postdoctoral Science Foundation (no. 2019M650750), and the National Science Foundation of China (grant nos. 51904266 and 51804270).