Theoretical unloading fracture mechanism and stability analysis of the slope rock masses in open-pit mines

Abstract

Instability of the rock slopes in open pit mines during the excavation unloading is an increasingly serious problem in the field of slope engineering. To explain the mechanical mechanism of slope unloading damage from a theoretical point of view, a simplified mechanical model of slope excavation based on the theory of fracture mechanics and the rock strength damage criterion is established. The damage process of the slope under excavation disturbance is dynamically analyzed by combining the interstructural characteristics of the slope. The solution equations for the extent of the plastic zone at the end of the crack of the excavated slope and its propagation length are derived. Calculation method of unloaded slope stability coefficient is proposed based on the mechanical model of crack propagation. The results show that (1) the stress intensity factor (SIF) at crack end in the slope under the action of unloading was larger than that under the original condition. (2) The range of the plastic zone at crack end in the slope rock mass can be attributed to the slope height, inverse logarithmic function to the slope angle, positive proportional function to the crack length, and the periodic fluctuation function of the crack angle. (3) The slope safety factor (SF) was found to be negatively related with the slope angle, slope height, crack angle, and the crack length unloading factor and positively related with the friction factor. Finally, the reasonableness of the theoretical derivation is verified by an engineering case study.

Keywords:

• Excavation unloading
• branch crack
• plastic zone at the crack end
• superposition principle
• safety factor (SF)

1. Introduction

Open pit mining is the main method of metal ore mining. With the depletion of resources in the shallow part of open pit mines, the mining of open pit mines is constantly moving to the deep part (Xie et al. 2017) and (Li et al. 2023). Deep mining has the characteristic of ‘strong disturbance’, which leads to high energy and large capacity engineering disasters (Niu et al. 2021). In the excavation process of the high slopes or ultra-high slopes, mining disturbance causes an instantaneous release of high stress and high energy that had been hidden in the rock mass, which induces crack propagation within the rock mass and coalescence with adjacent cracks of consistent orientation, resulting in serious landslide disasters (Wu et al. 2022). Therefore, it is of great engineering and theoretical significance to study the crack propagation law and rock stability state of open pit mine slopes under the action of excavation unloading.

The mechanical behavior of the unloaded rock masses, although, has been studied, it is not completely understood. For instance, Ha (1998) introduced the concept of the rock mechanics of the nonlinear unloading. Xu et al. (2019) analyzed the effects of the initial stress and the rate of unloading on the bearing capacity, and the deformation damage characteristics of Jinping marble. Liu and Deng (2019) used a combined finite-element-discrete-element modeling method (FDEM), and discussed the effect of the addition and removal rate on the cracking process of the surrounding rock. Duan et al. (2019) studied granite deformation characteristics, damage modes, and strain energy evolution patterns. Li and Shao (2016) derived a new stress-strain theoretical model and verified it through the unloading experiment of marble. Acoustic emission (AE) monitoring technology has been widely used in mining and rock mechanics research. Huang and Huang (2014) and Huang and Li (2014) and He et al. (2010) obtained the internal rupture laws of the specimens in the process of unloading by using the AE technique, which is utilized to study the internal break during the rupture of the rock masses. Du et al. (2016) and Li et al. (2018) studied the rock damage mechanism by collecting AE data during the unloading damage of the different test sample, and the energy evolution the rock samples under the different unloading paths were also studied from the energy perspective. Zhao et al. (2015) and Yang et al. (2014) proposed a new analysis method from the perspective of rock damage and failure. Zhao et al. (2014) explored the rock burst characteristics and energy evolution law of Beishan granite under different strain rates. Nie et al. (2022) studied the reinforcement effect of NPR anchor rods in slope through numerical simulation. Tao et al. (2019) explained the mechanical mechanism of the instability damage of an inverted anti-inclined slope by approximating the anti-inclined slope damage units as a cantilever beams. Liang et al. (2023) studied the fracture pattern of sandstone with different confining pressures and water soaking durations. Similar studies include, Tang et al. (2022) and Zhu et al. (2022). The above research shows that there are many studies on rock unloading. Through taking rocks of different lithology at the project site, the physical characteristics, damage type, and energy change law of rock specimens are explored through triaxial unloading test. With the continuous application of large-scale simulation calculations in the field of geotechnical engineering, we can quickly and conveniently study the physical characteristics, damage type, and energy change law of rock specimens under action of unloading by simulation calculations. The research conclusion effectively support the construction of the project site.

The study of the rock unloading mechanics is mainly based on the method of triaxial unloading test system to simulate the unloading mechanics behavior of the engineering rock masses to achieve a qualitative description or a semi-quantitative analysis, but this does not reflect the real unloading mechanics state of the slope rock masses well, and the relevant results also differ from the actual engineering monitoring results to some extent. On the other hand, the main studies about rock damage of unloading were exploring the damage characteristics and mechanical properties of cracked specimens under unloading conditions by making flaws with different numbers and distribution rules on the rock samples, and the methods included indoor mechanical experiments and numerical simulations. Despite this, there is an overall lack of studies, and understanding of the theory of unloading. Zhu and Huang (2019) simulated the shear damage pattern of sandstone induced by excavation unloading using the method of normal stress unloading. Zhong et al. (2020) studied the damage characteristics of rock samples containing parallel cracks under normal stress unloading. Huang and Huang (2014) and Huang and Li (2014) and Liu et al. (2020) studied the damage pattern of concrete specimens containing different crack distributions. Yang and Zhang (2020) simulated the fracture damage process of a specimen under unloading conditions, which had a single crack. Li et al. (2021) compared and analyzed the mechanical properties and damage pattern of cracked granite under cyclic loading and unloading and static loading. Yang et al. (2021, 2022) studied the relationship between crack propagation pattern and specimen strength of prefabricated crack specimens during loading process through DIC and numerical calculations. Similar studies include, Duan et al. (2019), Feng et al. (2019), Zhang et al. (2022), Ren et al. (2023), Yin et al. (2023), and Li et al. (2019). With the theory of unloading, Zhou and Chen (2019) derived expressions for the tip SIF for a pair of parallel offset cracks and co-linear cracks from the mechanical perspective. A similar study was done by Wang et al. (2020) and Zhu et al. (2021).

The above-mentioned studies on unloading all explored rock specimens as research objects, and some studies used simulation calculation and model tests to research the damage mechanism and the destructive patterns of the rock masses during slope unloading. Kong et al. (2011) and Li et al. (2021) studied the excavation failure patterns of slopes through a finite element model. Qu and Dang (2022) and Bao et al. (2023) investigated the slope instability mechanism under the combined effect of excavation and rainfall using finite element models. Zhu et al. (2020) and Zhao and Zhang (2020) found that the anti inclined slope exhibits a phased and regional failure pattern under excavation conditions. Fang et al. (2022) conducted many sets of physical model tests on slopes. Liu et al. (2022) reproduced the process of slope landslides triggered by underground tunnel excavation through model tests. Zhang et al. (2016) provided a method for calculating the safety factor of different types of dangerous rocks. Chang et al. (2019) explained the propagation mechanism of cracks in slope rock mass under loading and unloading conditions. At present, the study on unloading theory is mainly based on the laboratory scale to study the fracture of rocks containing flaws under the action of unloading, but not on the engineering scale.

In summary, studies on the unloading aspect initially started with rock specimens. However, with the strong progress of numerical calculations, studies are conducted to resolve the unloading damage characteristics of real slopes by means of numerical simulation calculations and model tests. Yet, there are relatively few studies regarding the theoretical exploration of the failure characteristics of unloading rock mass. Lack of theoretical study about rock mass damage under the action of unloading based on actual engineering characteristics.

To explain the mechanical mechanism of slope unloading damage from a theoretical point of view, this manuscript establishes a simplified mechanical model of slope excavation based on the theory of fracture mechanics and the rock strength damage criterion. The damage process of the slope under excavation disturbance is dynamically analyzed by combining the interstructural characteristics of the slope. The solution equations for the extent of the plastic zone at the end of the crack of the excavated slope and its propagation length are derived. Calculation method of unloaded slope stability coefficient is proposed based on the mechanical model of crack propagation. Finally, the reasonableness of the theoretical derivation is verified by an engineering case study. The study conclusions can provide a new way to explain the mechanism of slope instability in open-pit mines under the action of excavation.

2. A mechanical model of crack propagation in the rock mass for excavation-unloading slope

In Figure 1, In the process of forming a new mechanical equilibrium state after slop excavation, the rock in the unloading area will generate tensile stresses perpendicular to the excavation direction (Huang and Huang 2014; Huang and Li 2014; Huang and Huang 2010). Based on the rules of rock deformation, the slope rock mass is divided into three regions, which are ① unloading failure area (I), ② unloading influence area (II), and ③ unaffected area (III), as shown in Figure 1 (Wang et al. 2023). In Figure 1, σ₀, σ_z, Δσ, ${{\sigma }^{\prime }}_x$ and σ_x represent the original rock stress, vertical stress, horizontal tensile stress generated by rock deformation, horizontal tensile stress of unloading failure area, and horizontal stress, respectively.

Figure 1. Stress state of slope before and after excavation.

The area (I) is the closest area to the slope, this area is most affected by excavation unloading, the rock is more broken, no catastrophic landslide will occur, usually a small landslide between local steps that can be detected in time. The unloading influence area (II) rock unloading is not sufficient, discontinuous cracks in this area rock will continue to propagation under the action of unloading tensile stress. With the continuous excavation unloading, discontinuous cracks will break through the rock bridge area and penetrate, and then a catastrophic landslide will occur. The crack propagation in the unloading influence area is hidden and difficult to detect in time, so it is important to strengthen the research in this area. In this paper, we mainly focus on the fracture mechanism of rock mass in the unloading influence area.

With the continuous mining of open pit mines, the slope is subjected to the combined effect of unloading stress and self-weight stress, the mechanical state of unloading slope is relatively complex, and the SIF at the crack end cannot be directly solved based on the current theory. Based on the stress superposition theory, the force state of the rock mass under multiple loads can be equated to the sum under a single load. Therefore, the stress state of the slope in the unloading influence area (II) can be decomposed into the stress state before being disturbed by unloading (initial stress) and the stress state only by the disturbed stress (mining stress), as shown in Figure 2. σ_z and σ_x have been introduced before, and Δσ is the stress unloaded instantaneously by the slope excavation, which is collectively referred to as the unloading stress in this study.

Figure 2. Mechanical model. (a) The mechanical model of weak unloading area, (b) stress state before unloading disturbance, and (c) stress state under unloading disturbance.

2.1. Mechanical analysis of slope rock mass under initial stress

Open pit slopes are stable when not disturbed by external loads, and the initial stress is formed by the joint action of the slope self-weight stress and the tectonic stress within the rock mass. When the slope is excavated according to different design slope angles, the initial stresses within the rock mass are expressed in different forms, and the specific expressions are (Wang et al. 2015): (1) $$\begin{array}{l}{\sigma }_{\mathrm{x}}\text{\hspace{0.17em}=\hspace{0.17em}}{K}_0\gamma \mathrm{z}\left(1-\hspace{0.17em}\sin \hspace{0.17em}\beta \right)\\ {\sigma }_{\mathrm{z}}\text{\hspace{0.17em}=\hspace{0.17em}}\gamma \mathrm{z}\end{array}\}\left(\mathrm{z}\in \left[0\mathrm{, }H\right]\right)$$(1) where H is the slope height, β is the slope angle, K₀ is the lateral pressure coefficient, γ is the slope rock density, z is the vertical distance between the slope excavation and the top of the slope, and the meanings of σ_z and σ_x are the same as before.

Figure 3(a), shows a more detailed force diagram in Figure 2(b), where the stress unit is in a stable state under the action of vertical and horizontal stresses(σ_z and σ_x), and the angle between the crack in the rock mass and the horizontal direction is α. Further decomposition of σ_z and σ_x along the crack surface, the normal stress, σ_n, perpendicular to the crack surface and the tangential stress, τ_n, along the crack surface can be obtained as (the direction of tensile stress specifies the positive direction): (2) $$\begin{array}{l}{\sigma }_n=\frac{{\sigma }_{\mathrm{z}}+{\sigma }_x}{2}+\frac{{\sigma }_{\mathrm{z}}-{\sigma }_x}{2}\cos \hspace{0.17em}2\alpha \\ {\tau }_n=\frac{{\sigma }_{\mathrm{z}}-{\sigma }_x}{2}\sin \hspace{0.17em}2\alpha \end{array}\}$$(2)

Figure 3. Mechanical model of crack propagation.

The cracks in the rock mass will be closed under the action of biaxial stress, and the relative sliding between the cracks will occur along the crack surface, and the effective shear stress, τ_eff, between the cracks drives the relative sliding of the cracks, and the type II SIF, $K_{\text{II}}^{},$ at the tip of the closed cracks can be expressed as follows: (3) $$K_{\text{II}}\text{\hspace{0.17em}=\hspace{0.17em}}{\tau }_{\text{eff}}\sqrt{\pi a}\text{=}\left({\tau }_{\mathrm{n}}-f{\sigma }_{\mathrm{n}}\right)\sqrt{\pi a}$$(3)

Where f and a are the friction coefficient and crack half-length, respectively.

2.2. Mechanical analysis of slope rock mass under mining stress

Figure 3(b) shows a more detailed force diagram in Figure 2(c), the stress unit in the slope rock mass under the action of the mining stress is subjected to horizontal stress △σ, the angle between the fissure in the rock body and the horizontal direction is α. Further decomposition of △σ along the crack surface, the normal stress, ${\sigma }_n^{\prime },$ perpendicular to the crack surface and the tangential stress, ${\tau }_n^{\prime },$ along the crack surface can be expressed as: (4) $$\begin{array}{l}{\sigma }_n^{\prime }=\frac{1}{2}\Delta \sigma \left(1-\hspace{0.17em}\cos \hspace{0.17em}2\alpha \right)\\ {\tau }_n^{\prime }=\frac{1}{2}\Delta \sigma \hspace{0.17em}\sin \hspace{0.17em}2\alpha \end{array}\}$$(4)

Figure 3(b) shows that the rock force unit is subjected to horizontal tensile stress, the tensile stress is decomposed along the crack face, the force acting on the crack face includes shear stress along the crack face and tensile stress outward from the vertical crack face, at this time the crack will occur tensile-shear composite damage, the crack end SIF can be expressed as: (5) $$\begin{array}{l}{K}_{{\rm I}\mathrm{M}}=\frac{1}{2}\Delta \sigma \left(1-\hspace{0.17em}\cos \hspace{0.17em}2\alpha \right)\sqrt{\pi a}\\ K_{\text{IIM}}=\frac{1}{2}\Delta \sigma \hspace{0.17em}\sin \hspace{0.17em}2\alpha \sqrt{\pi a}\end{array}\}$$(5) where $K_{\text{IM}}$ is the type I SIF at the crack end, and $K_{\text{IIM}}$ is the type II SIF at the crack end (mining stress).

Based on the SIF superposition theory, the SIF at crack end in the rock mass of slope excavation and unloading is equal to the sum of the SIF at crack end in the initial stress state and the crack end SIF in the mining stress state, therefore, the expression of SIF at crack end in the slope rock mass can be expressed as (joint vertical (1)–(5)): (6) $$\begin{array}{l}{K}_{{\rm I}X}=\frac{1}{2}\mathrm{\Delta }\sigma \left(1-\hspace{0.17em}\cos \hspace{0.17em}2\alpha \right)\sqrt{\pi a}\\ K_{\text{IIX}}=\frac{1}{2}\sqrt{\pi a}[{\sigma }_z\left(\hspace{0.17em}\sin \hspace{0.17em}2\alpha -f-f\hspace{0.17em}\cos \hspace{0.17em}2\alpha \right)+{\sigma }_x\left(f\hspace{0.17em}\cos \hspace{0.17em}2\alpha -f-\hspace{0.17em}\sin \hspace{0.17em}2\alpha \right)]+\frac{1}{2}\mathrm{\Delta }\sigma \hspace{0.17em}\sin \hspace{0.17em}2\alpha \sqrt{\pi a}\end{array}\}$$(6) where $K_{\text{IX}}$ is the type I SIF at the crack end, and $K_{\text{IIX}}$ is the type II SIF at the crack end .(unloading stress)

2.3. The influence of excavation unloading on the slope of the rock mass damage

Before excavation, the slope rock mass is in equilibrium, the crack in the rock mass is in a closed state under the original rock stress, and the stress unit at a height of H from the horizontal plane is subject to the vertical stress, ${\sigma }_z^0$ = γH, horizontal stress, ${\sigma }_x^0$ = K₀γH, and at this point, the expression for the SIF of the crack end type II is: (7) $$K_{\text{II}}\text{=}\frac{1}{2}[{\sigma }_z^0\left(\hspace{0.17em}\sin \hspace{0.17em}2\alpha -f-f\hspace{0.17em}\cos \hspace{0.17em}2\alpha \right)+{\sigma }_x^0\left(f\hspace{0.17em}\cos \hspace{0.17em}2\alpha -f-\hspace{0.17em}\sin \hspace{0.17em}2\alpha \right)]\sqrt{\pi a}$$(7) where $K_{\text{II}}$ is the SIF of the crack end of the slope rock mass before unloading, and other parameters have the same physical meaning as before.

Combining Equations (6)(6) $$\begin{array}{l}{K}_{{\rm I}X}=\frac{1}{2}\mathrm{\Delta }\sigma \left(1-\hspace{0.17em}\cos \hspace{0.17em}2\alpha \right)\sqrt{\pi a}\\ K_{\text{IIX}}=\frac{1}{2}\sqrt{\pi a}[{\sigma }_z\left(\hspace{0.17em}\sin \hspace{0.17em}2\alpha -f-f\hspace{0.17em}\cos \hspace{0.17em}2\alpha \right)+{\sigma }_x\left(f\hspace{0.17em}\cos \hspace{0.17em}2\alpha -f-\hspace{0.17em}\sin \hspace{0.17em}2\alpha \right)]+\frac{1}{2}\mathrm{\Delta }\sigma \hspace{0.17em}\sin \hspace{0.17em}2\alpha \sqrt{\pi a}\end{array}\}$$(6) and (7)(7) $$K_{\text{II}}\text{=}\frac{1}{2}[{\sigma }_z^0\left(\hspace{0.17em}\sin \hspace{0.17em}2\alpha -f-f\hspace{0.17em}\cos \hspace{0.17em}2\alpha \right)+{\sigma }_x^0\left(f\hspace{0.17em}\cos \hspace{0.17em}2\alpha -f-\hspace{0.17em}\sin \hspace{0.17em}2\alpha \right)]\sqrt{\pi a}$$(7) , it can be obtained as follows: (8) $$K_{\text{IIX}}-K_{\text{II}}\text{=}\frac{1}{2}\sqrt{\pi a}{K}_0\gamma H\hspace{0.17em}\sin \hspace{0.17em}\beta (f+\hspace{0.17em}\sin \hspace{0.17em}2\alpha -f\hspace{0.17em}\cos \hspace{0.17em}2\alpha )+\frac{1}{2}\mathrm{\Delta }\sigma \hspace{0.17em}\sin \hspace{0.17em}2\alpha \sqrt{\pi a}$$(8)

In Equation (8)(8) $$K_{\text{IIX}}-K_{\text{II}}\text{=}\frac{1}{2}\sqrt{\pi a}{K}_0\gamma H\hspace{0.17em}\sin \hspace{0.17em}\beta (f+\hspace{0.17em}\sin \hspace{0.17em}2\alpha -f\hspace{0.17em}\cos \hspace{0.17em}2\alpha )+\frac{1}{2}\mathrm{\Delta }\sigma \hspace{0.17em}\sin \hspace{0.17em}2\alpha \sqrt{\pi a}$$(8) , when α∈[0, 90°], $(f+\hspace{0.17em}\sin \hspace{0.17em}2\alpha -f\hspace{0.17em}\cos \hspace{0.17em}2\alpha )$ and $\mathrm{\Delta }\sigma \hspace{0.17em}\sin \hspace{0.17em}2\alpha \sqrt{\pi a}$ are positive, so $K_{\text{II}X}-K_{\text{II}}$ is positive, which leads to the type II SIF at the crack end in the slope rock mass after unloading being greater than before unloading, and it is known from Equation (6)(6) $$\begin{array}{l}{K}_{{\rm I}X}=\frac{1}{2}\mathrm{\Delta }\sigma \left(1-\hspace{0.17em}\cos \hspace{0.17em}2\alpha \right)\sqrt{\pi a}\\ K_{\text{IIX}}=\frac{1}{2}\sqrt{\pi a}[{\sigma }_z\left(\hspace{0.17em}\sin \hspace{0.17em}2\alpha -f-f\hspace{0.17em}\cos \hspace{0.17em}2\alpha \right)+{\sigma }_x\left(f\hspace{0.17em}\cos \hspace{0.17em}2\alpha -f-\hspace{0.17em}\sin \hspace{0.17em}2\alpha \right)]+\frac{1}{2}\mathrm{\Delta }\sigma \hspace{0.17em}\sin \hspace{0.17em}2\alpha \sqrt{\pi a}\end{array}\}$$(6) that the fracture tip is also subject to the type I tensile damage (as defined here) after unloading. Therefore, the cracks in the slope of the rock mass before excavation do not undergo fracture damage under the original rock stress. After the excavation of the slope, the rock mass is in the state of pressure relief within a certain range near the slope surface, and the SIF at the crack end increases, which ultimately leads to the destruction of the rock mass, which in all leads to landslide disasters. The crack end SIF increases as a function of increased slope height, which indicates that the slope is more unstable with the greater slope height. Increased crack end SIF with increased unloading indicates that the rock mass is more unstable with greater unloading disturbance.

3. The solution to the plastic zone range at the crack end and crack propagation length of the slope rock mass under excavation unloading

In Figure 3(a) of Section 2, compression-shear propagation of cracks within the slope rock occurs under the initial stress, and branch crack propagation occurs as the excavation continues, and the branch crack end SIF can be expressed as (Liu et al. 2012): (9) $$K_{\mathrm{I}}\left(l\right)\text{\hspace{0.17em}=\hspace{0.17em}}{\left(K_{\mathrm{I}}\right)}_1+{\left(K_{\mathrm{I}}\right)}_2\text{=}\frac{{\tau }_{\mathit{eff}}2a\hspace{0.17em}\sin \hspace{0.17em}\alpha }{\sqrt{\pi l}}-{\sigma }_x\sqrt{\pi l}$$(9) where ${(K_{\mathrm{I}})}_1$ is the SIF generated by τ_eff, ${(K_{\mathrm{I}})}_2$ is the SIF generated by σ_x, and l is the branch crack length.

In Figure 3(b) of Section 2, tensile-shear expansion of cracks within the slope rock occurs under mining stress, and branch crack propagation occurs as excavation continues, and the SIF at branch crack end can be expressed as (Kemeny 1991): (10) $$K_{\mathrm{I}M}\left(l\right)=\frac{5.18a{\tau }_n^{\prime }\hspace{0.17em}\cos \hspace{0.17em}\alpha }{\sqrt{\pi l}}$$(10)

Based on the SIF superposition theory, the SIF at branch crack end by excavation unloading of the slope is the sum of the SIF of branch crack end in the initial stress state and mining stress state, and the specific expression is as follows: (11) $$K_{\mathrm{I}X}\left(l\right)\text{=}{K}_{\mathrm{I}}\left(l\right)+K_{\mathrm{I}M}\left(l\right)\text{= }\frac{{\tau }_{\mathit{eff}}2a\hspace{0.17em}\sin \hspace{0.17em}\alpha }{\sqrt{\pi l}}+\frac{5.18a{\tau }_n^{\prime }\hspace{0.17em}\cos \hspace{0.17em}\alpha }{\sqrt{\pi l}}-{\sigma }_x\sqrt{\pi l}$$(11)

3.1. The solution to the plastic zone range at the crack end under excavation unloading

There is always a plastic zone at the crack end area, which is under the effect of excavation unloading disturbance. The basic theory of the elastic-plastic fracture mechanics quantitatively describes the stress distribution at the crack end, and the plastic zone at the crack end area of the unloaded rock on the slope can be calculated by combining it with the rock-strength damage criterion. With the elastodynamics, the asymptotic stress field of the crack end type I can be expressed as: (12) $$\begin{array}{l}{\sigma }_x\left(r,\theta \right)\text{=}\frac{K_{{\rm I}}}{\sqrt{2\pi r}}\cos \hspace{0.17em}\frac{\theta }{2}\left(1-\hspace{0.17em}\sin \hspace{0.17em}\frac{\theta }{2}\sin \hspace{0.17em}\frac{3\theta }{2}\right)\\ {\sigma }_y\left(r,\theta \right)\text{=}\frac{K_{{\rm I}}}{\sqrt{2\pi r}}\cos \hspace{0.17em}\frac{\theta }{2}\left(1+\hspace{0.17em}\sin \hspace{0.17em}\frac{\theta }{2}\sin \hspace{0.17em}\frac{3\theta }{2}\right)\\ {\tau }_{xy}\left(r,\theta \right)\text{=}\frac{K_{{\rm I}}}{\sqrt{2\pi r}}\sin \hspace{0.17em}\frac{\theta }{2}\cos \hspace{0.17em}\frac{\theta }{2}\cos \hspace{0.17em}\frac{3\theta }{2}\end{array}\}$$(12) where $K_{\mathrm{I}}$ is the SIF of type I at crack end; ${\sigma }_x(r,\theta ),$ ${\sigma }_y(r,\theta ),$ and ${\tau }_{xy}(r,\theta )$ are radial stress, axial stress, and shear stress, respectively; θ is the crack end propagation direction angle, and r is the distance from the crack end at any point on the crack.

The expression for the calculation of the principal stress in elastic mechanics can be expressed as: (13) $$\begin{array}{l}{\sigma }_1,{\sigma }_2\text{=}\frac{{\sigma }_x+{\sigma }_y}{2}\pm \sqrt{{\left(\frac{{\sigma }_x-{\sigma }_y}{2}\right)}^2+{\tau }_{xy}^2}\\ {\sigma }_3\text{=}\mu \left({\sigma }_x+{\sigma }_y\right)\end{array}\}$$(13)

Combining Equations (12)(12) $$\begin{array}{l}{\sigma }_x\left(r,\theta \right)\text{=}\frac{K_{{\rm I}}}{\sqrt{2\pi r}}\cos \hspace{0.17em}\frac{\theta }{2}\left(1-\hspace{0.17em}\sin \hspace{0.17em}\frac{\theta }{2}\sin \hspace{0.17em}\frac{3\theta }{2}\right)\\ {\sigma }_y\left(r,\theta \right)\text{=}\frac{K_{{\rm I}}}{\sqrt{2\pi r}}\cos \hspace{0.17em}\frac{\theta }{2}\left(1+\hspace{0.17em}\sin \hspace{0.17em}\frac{\theta }{2}\sin \hspace{0.17em}\frac{3\theta }{2}\right)\\ {\tau }_{xy}\left(r,\theta \right)\text{=}\frac{K_{{\rm I}}}{\sqrt{2\pi r}}\sin \hspace{0.17em}\frac{\theta }{2}\cos \hspace{0.17em}\frac{\theta }{2}\cos \hspace{0.17em}\frac{3\theta }{2}\end{array}\}$$(12) and (13)(13) $$\begin{array}{l}{\sigma }_1,{\sigma }_2\text{=}\frac{{\sigma }_x+{\sigma }_y}{2}\pm \sqrt{{\left(\frac{{\sigma }_x-{\sigma }_y}{2}\right)}^2+{\tau }_{xy}^2}\\ {\sigma }_3\text{=}\mu \left({\sigma }_x+{\sigma }_y\right)\end{array}\}$$(13) , the type I stress component at the tip of the branch crack in polar coordinates can be expressed as: (14) $$\begin{array}{l}{\sigma }_1,{\sigma }_2=\frac{K_{{\rm I}}}{\sqrt{2\pi r}}\left(\hspace{0.17em}\cos \hspace{0.17em}\frac{\theta }{2}\pm \sqrt{{\left(\hspace{0.17em}\cos \hspace{0.17em}\frac{\theta }{2}\sin \hspace{0.17em}\frac{\theta }{2}\sin \hspace{0.17em}\frac{3\theta }{2}\right)}^2+{\left(\hspace{0.17em}\sin \hspace{0.17em}\frac{\theta }{2}\cos \hspace{0.17em}\frac{\theta }{2}\cos \hspace{0.17em}\frac{3\theta }{2}\right)}^2}\right)\\ {\sigma }_3\text{=}\frac{2\mu K_{{\rm I}}}{\sqrt{2\pi r}}\cos \hspace{0.17em}\frac{\theta }{2}\end{array}\}$$(14)

The Mohr–Coulomb (M–C) strength criterion was used to solve the plastic zone range at the branch crack end, and the specific calculation expression is shown below: (15) $${\sigma }_1=\frac{1-\hspace{0.17em}\sin \hspace{0.17em}\phi }{1+\hspace{0.17em}\sin \hspace{0.17em}\phi }{\sigma }_3+\frac{2c\hspace{0.17em}\cos \hspace{0.17em}\phi }{1-\hspace{0.17em}\sin \hspace{0.17em}\phi }$$(15) where ${\sigma }_1$ is the maximum principal stress at the time of rock failure, ${\sigma }_3$ is the minimum principal stress, $c$ is the cohesive force of the rock, and ${\sigma }_3\text{=}{\sigma }_0-{\sigma }_{3t}$ is the internal friction angle of the rock mass.

By combining Equations (11)–(15), the equation for the plastic zone at the branch crack end at $\theta$=0 (on the crack propagation) is shown below: (16) $$r_{\mathrm{t}}=\frac{K_{{\rm I}X}^2\left(l\right)}{2\pi {\left[\frac{2c\hspace{0.17em}\cos \hspace{0.17em}\phi \left(1\hspace{0.17em}+\hspace{0.17em}\hspace{0.17em}\sin \hspace{0.17em}\phi \right)}{1-{\hspace{0.17em}\sin \hspace{0.17em}}^2\phi \hspace{0.17em}-\hspace{0.17em}\mu {\left(1-\hspace{0.17em}\sin \hspace{0.17em}\phi \right)}^2}\right]}^2}$$(16) where r_t is the range of the plastic zone at the branch crack end (m).

3.2. The solution to the crack propagation length of the slope rock mass under excavation unloading

The number and length of cracks in the rock mass of the slope have a controlling effect on the slope stability, and the formation of the sliding surface under the effect of excavation disturbance is caused by the propagation and coalescence of the dominant crack group in the rock mass. In continuous slope excavation unloading, the long cracks in the rock mass first initiate propagation along the excavation surface at an approximately parallel angle, and with further slope excavation unloading, connection occurs between crack groups with basically the same direction, resulting in the occurrence of slope landslide disasters. With further crack propagation, the rock bridge between the cracks affected by the interaction between the cracks is destroyed, and a slip surface is formed within a certain distance from the unloading surface. At present, the study method for the rock bridge breaking problem is to establish the equation between the rock fracture toughness and the SIF at crack end . In this study, the function of the extremum method to solve for the range of the plastic zone at branch crack end, and crack propagation length when the rock mass reaches the damage threshold, was used. The rock bridge mechanics model of the interaction between branch cracks is shown in Figure 4.

Figure 4. The diagram of the unloading crack interaction of slope.

The length, S, between the centers of the two adjacent dominant cracks, and the length, T, of the rock bridge between the branch cracks can be expressed as follows: (17) $$\begin{array}{l}S\text{=}\frac{1}{\sqrt{N_A}}\\ T=S-2\left[l\text{\hspace{0.17em}+\hspace{0.17em}}r\text{\hspace{0.17em}+\hspace{0.17em}}a\hspace{0.17em}\cos \hspace{0.17em}\left(\omega -\alpha \right)\right]\end{array}\}$$(17) where N_A is the number of cracks per unit area, $\omega$ is the angle between the slope surface and the horizontal direction.

In the mechanical model of the cracked rock bridge shown in Figure 5, the force, $F\text{=}{T}_{\mathrm{e}}^{\prime }\hspace{0.17em}\sin \hspace{0.17em}(\omega -\alpha ),$ on the crack is balanced by the stress, ${\sigma }_{\mathrm{n}},$ between the rock bridge, where $T_{\mathrm{e}}^{\prime }\text{=2}a({\tau }_{\mathit{eff}}+{\tau }_{\mathrm{n}}^{\prime }),$ and the expression for the calculation of the mechanical equilibrium can be expressed as: (18) $${\sigma }_n=\pi A^2\cdot \frac{T_{\mathrm{e}}^{\prime }\hspace{0.17em}\sin \hspace{0.17em}\left(\omega -\alpha \right)}{\pi A^{2}S-\left[2\pi A^{2}l+{\left(\frac{{\tau }_{\mathit{eff}}2a\hspace{0.17em}\sin \hspace{0.17em}\alpha }{\sqrt{\pi l}}+\frac{5.18a{\tau }_n^{\prime }\hspace{0.17em}\cos \hspace{0.17em}\alpha }{\sqrt{\pi l}}-{\sigma }_x\sqrt{\pi l}\right)}^2+2\pi A^{2}a\hspace{0.17em}\cos \hspace{0.17em}\left(\omega -\alpha \right)\right]}$$(18)

Figure 5. Analysis of factors affecting the size of the plastic zone at the crack end under unloading conditions.

In Equation (18)(18) $${\sigma }_n=\pi A^2\cdot \frac{T_{\mathrm{e}}^{\prime }\hspace{0.17em}\sin \hspace{0.17em}\left(\omega -\alpha \right)}{\pi A^{2}S-\left[2\pi A^{2}l+{\left(\frac{{\tau }_{\mathit{eff}}2a\hspace{0.17em}\sin \hspace{0.17em}\alpha }{\sqrt{\pi l}}+\frac{5.18a{\tau }_n^{\prime }\hspace{0.17em}\cos \hspace{0.17em}\alpha }{\sqrt{\pi l}}-{\sigma }_x\sqrt{\pi l}\right)}^2+2\pi A^{2}a\hspace{0.17em}\cos \hspace{0.17em}\left(\omega -\alpha \right)\right]}$$(18) , $A=\left[2c\hspace{0.17em}\cos \hspace{0.17em}\phi (1+\hspace{0.17em}\sin \hspace{0.17em}\phi )\right]/\left[1-{\hspace{0.17em}\sin \hspace{0.17em}}^2\phi -\mu {(1-\hspace{0.17em}\sin \hspace{0.17em}\phi )}^2\right].$

Let ${\sigma }_{\mathrm{n}}$ = ${\sigma }_{\mathrm{t}}$ (${\sigma }_{\mathrm{t}}$ is the uniaxial tensile strength of the rock mass), and Equation (18)(18) $${\sigma }_n=\pi A^2\cdot \frac{T_{\mathrm{e}}^{\prime }\hspace{0.17em}\sin \hspace{0.17em}\left(\omega -\alpha \right)}{\pi A^{2}S-\left[2\pi A^{2}l+{\left(\frac{{\tau }_{\mathit{eff}}2a\hspace{0.17em}\sin \hspace{0.17em}\alpha }{\sqrt{\pi l}}+\frac{5.18a{\tau }_n^{\prime }\hspace{0.17em}\cos \hspace{0.17em}\alpha }{\sqrt{\pi l}}-{\sigma }_x\sqrt{\pi l}\right)}^2+2\pi A^{2}a\hspace{0.17em}\cos \hspace{0.17em}\left(\omega -\alpha \right)\right]}$$(18) can be reorganized into a form about the function after the variable substitution: (19) $$\varphi \left(l\right)\text{=}B{l}^2+Cl+D$$(19)

In Equation (19), $B=2{\pi }^2{A}^2{\sigma }_{\mathrm{t}}+{\pi }^2{\sigma }_x^2{\sigma }_{\mathrm{t}},$ $$C={\pi }^2{A}^2{T}_{\mathrm{e}}^{\prime }\hspace{0.17em}\sin \hspace{0.17em}\left(\omega -\alpha \right)-{\pi }^2{A}^{2}S{\sigma }_{\mathrm{t}}-4\pi a{\tau }_{\mathit{eff}}{\sigma }_x{\sigma }_{\mathrm{t}}-10.36\pi a{\tau }_{\mathrm{n}}^{\prime }{\sigma }_x{\sigma }_{\mathrm{t}}\hspace{0.17em}\cos \hspace{0.17em}\alpha -2a{\pi }^2{A}^2{\sigma }_{\mathrm{t}}\hspace{0.17em}\cos \hspace{0.17em}\left(\omega -\alpha \right),$$ $$D=4a^2{\tau }_{\mathit{eff}}^2{\sigma }_{\mathrm{t}}{\hspace{0.17em}\sin \hspace{0.17em}}^2\alpha +26.83a^2{\left({\tau }_{\mathrm{n}}^{\prime }\right)}^2{\sigma }_{\mathrm{t}}{\hspace{0.17em}\cos \hspace{0.17em}}^2\alpha +10.36a^2{\tau }_{\mathit{eff}}{\tau }_{\mathrm{n}}^{\prime }{\sigma }_{\mathrm{t}}\hspace{0.17em}\sin \hspace{0.17em}2\alpha .$$

The second derivation of Equation (19) gives: (20) $$\begin{array}{l}\frac{\partial \varphi \left(l\right)}{\partial l}\text{=2}Bl+C=0\\ \frac{{\partial }^2\varphi \left(l\right)}{{\partial }^{2}l}=2B>0\end{array}\}$$(20)

Therefore, when the unloaded rock reaches the damage threshold, the branch crack propagation limit length, $l_{\mathrm{t}},$ can be expressed as follows: (21) $$l\text{\hspace{0.17em}=\hspace{0.17em}}{l}_{\mathrm{t}}\text{\hspace{0.17em}=\hspace{0.17em}}-\frac{C}{2B}\text{=}\frac{\begin{array}{c}{\pi }^2{A}^{2}S{\sigma }_{\mathrm{t}}+2\pi a{\sigma }_z{\sigma }_x{\sigma }_{\mathrm{t}}\left(\hspace{0.17em}\sin \hspace{0.17em}2\alpha -f-f\hspace{0.17em}\cos \hspace{0.17em}2\alpha \right)+2\pi a{\sigma }_z{\sigma }_x{\sigma }_{\mathrm{t}}\left(f\hspace{0.17em}\cos \hspace{0.17em}2\alpha -f-\hspace{0.17em}\sin \hspace{0.17em}2\alpha \right)-\\ 10.36\pi a\eta {\sigma }_x\hspace{0.17em}\sin \hspace{0.17em}2\alpha {\sigma }_x{\sigma }_{\mathrm{t}}\hspace{0.17em}\cos \hspace{0.17em}\alpha +2a{\pi }^2{A}^2{\sigma }_{\mathrm{t}}\hspace{0.17em}\cos \hspace{0.17em}\left(\omega -\alpha \right)-\\ \left[{\pi }^2{A}^{2}a{\sigma }_z\left(\hspace{0.17em}\sin \hspace{0.17em}2\alpha -f-f\hspace{0.17em}\cos \hspace{0.17em}2\alpha \right)+{\pi }^2{A}^{2}a{\sigma }_x\left(f\hspace{0.17em}\cos \hspace{0.17em}2\alpha -f-\hspace{0.17em}\sin \hspace{0.17em}2\alpha \right)-{\pi }^2{A}^{2}a\eta {\sigma }_x\hspace{0.17em}\sin \hspace{0.17em}2\alpha \right]\hspace{0.17em}\sin \hspace{0.17em}\left(\omega -\alpha \right)\end{array}}{4{\pi }^2{A}^2{\sigma }_{\mathrm{t}}+2{\pi }^2{\sigma }_x^2{\sigma }_{\mathrm{t}}}$$(21)

By substituting Equation (21)(21) $$l\text{\hspace{0.17em}=\hspace{0.17em}}{l}_{\mathrm{t}}\text{\hspace{0.17em}=\hspace{0.17em}}-\frac{C}{2B}\text{=}\frac{\begin{array}{c}{\pi }^2{A}^{2}S{\sigma }_{\mathrm{t}}+2\pi a{\sigma }_z{\sigma }_x{\sigma }_{\mathrm{t}}\left(\hspace{0.17em}\sin \hspace{0.17em}2\alpha -f-f\hspace{0.17em}\cos \hspace{0.17em}2\alpha \right)+2\pi a{\sigma }_z{\sigma }_x{\sigma }_{\mathrm{t}}\left(f\hspace{0.17em}\cos \hspace{0.17em}2\alpha -f-\hspace{0.17em}\sin \hspace{0.17em}2\alpha \right)-\\ 10.36\pi a\eta {\sigma }_x\hspace{0.17em}\sin \hspace{0.17em}2\alpha {\sigma }_x{\sigma }_{\mathrm{t}}\hspace{0.17em}\cos \hspace{0.17em}\alpha +2a{\pi }^2{A}^2{\sigma }_{\mathrm{t}}\hspace{0.17em}\cos \hspace{0.17em}\left(\omega -\alpha \right)-\\ \left[{\pi }^2{A}^{2}a{\sigma }_z\left(\hspace{0.17em}\sin \hspace{0.17em}2\alpha -f-f\hspace{0.17em}\cos \hspace{0.17em}2\alpha \right)+{\pi }^2{A}^{2}a{\sigma }_x\left(f\hspace{0.17em}\cos \hspace{0.17em}2\alpha -f-\hspace{0.17em}\sin \hspace{0.17em}2\alpha \right)-{\pi }^2{A}^{2}a\eta {\sigma }_x\hspace{0.17em}\sin \hspace{0.17em}2\alpha \right]\hspace{0.17em}\sin \hspace{0.17em}\left(\omega -\alpha \right)\end{array}}{4{\pi }^2{A}^2{\sigma }_{\mathrm{t}}+2{\pi }^2{\sigma }_x^2{\sigma }_{\mathrm{t}}}$$(21) into Equation (16)(16) $$r_{\mathrm{t}}=\frac{K_{{\rm I}X}^2\left(l\right)}{2\pi {\left[\frac{2c\hspace{0.17em}\cos \hspace{0.17em}\phi \left(1\hspace{0.17em}+\hspace{0.17em}\hspace{0.17em}\sin \hspace{0.17em}\phi \right)}{1-{\hspace{0.17em}\sin \hspace{0.17em}}^2\phi \hspace{0.17em}-\hspace{0.17em}\mu {\left(1-\hspace{0.17em}\sin \hspace{0.17em}\phi \right)}^2}\right]}^2}$$(16) , the tip plasticity zone, $r_{\mathrm{t}},$ can be obtained as follows: (22) $$r_{\mathrm{t}}=\frac{{\left({\tau }_{\mathit{eff}}2a\hspace{0.17em}\sin \hspace{0.17em}\alpha +5.18a{\tau }_n^{\prime }\hspace{0.17em}\cos \hspace{0.17em}\alpha -{\sigma }_x\pi l_t\right)}^2}{2{\pi }^2{A}^2{l}_t}$$(22) where $l_{\mathrm{t}}$ and $r_{\mathrm{t}}$ are the branch crack propagation length and plastic zone range, respectively.

Hence, the main factors affecting the size of the plastic zone at the crack end within the slope of the rock mass, under the effect of excavation unloading are the crack angle, slope angle, crack length, and slope height. To further analyze the relationship between the plastic zone at the crack end and the influencing factors, the values of the internal friction angle, $\phi ,$ range from 16° to 40°, the values of cohesion, c, range from 1.0 to 10 MPa, and the values of $\phi$ and, in turn, c characterize the different stability states of the sloping rock.

The relationship between the range of the plastic zone at the crack end and the slope height is approximately exponential (Figure 5(a)) under the unloading, that is, the plastic zone width increased with the increased of the slope height until it tended to infinity, and its curvature increased gradually until it reached infinity. When the slope height is between 0 and 300 m, the change in the plastic zone width was very small, which indicated that the change in the slope height had less influence on the plastic zone at this point, and the influence of the unloading disturbance on the slope damage was smaller. When the slope height approached 800 m, the plastic zone width increased sharply. The above conclusion is also consistent with the observation that in engineering practices, in general, when the slope is overrun from a high slope to an ultra-high slope, excavation and unloading under high ground stress will result in the mechanical response to the ultra-high slope being completely different from that of the low and medium slopes.

The relationship between the range of the plastic zone at the crack end and the slope angle shows an inverse logarithmic (Figure 5(b)) function under unloading, that is, the width of the plastic zone decreased with the increased slope angle until it tended to be stable, and its curvature gradually decreased until it tended to 0. The slope angle mainly affected the size of the horizontal stress, and the range of the plastic zone is basically not affected by the inclination angle ($\omega$ > 60°), and the horizontal stress was smaller. At this point, the corresponding horizontal unloading, with the same unloading factor, was smaller, and, therefore, the unloading disturbance to the slope, the SIF at the crack end, and the plastic zone at the crack end was smaller. The influence of slope angle on the fracture range of unloading slope rock mass becomes larger ($\omega$ < 60°), which resulted in a sharp increase in the width of the plastic zone. The trend in Figure 6(b) also indirectly verifies the influence law of the unloading disturbance on the slope stability under the different stress states considered here.

Figure 6. Schematic diagram of through stress of slope rock bridge.

The relationship between the range of the plastic zone at the crack end and the crack length shows a positive proportional function (Figure 5(c)) under the unloading, that is, the plastic zone width increased with the increased slope height, and its curvature remained constant. The hierarchy between the different surfaces also shows that the denser was the distribution of slope joints, the larger was the plastic zone width.

The relationship between the range of the plastic zone at the crack end and the crack angle shows a periodic fluctuation function (Figure 5(d)) under the unloading. This was mainly caused by the fact that when the fracture angle was closer to the angle of the slip surface, which was composed of the dominant crack group, the greater was the influence of the unloading disturbance on the crack propagation. Combined with the mechanical model in Figure 5, when the crack angle gradually increased close to the slip surface angle, the force applied by the sliding mass of the side slope to decompose the load along the crack surface increased, that is, the shear load along the crack surface increased, the SIF at the rack end increased, in the end, but whose role cannot be evaluated at the moment. Therefore, the plastic zone width was that of a fluctuating change rather than it being constant.

4. The stability analysis of the slope excavation unloading

4.1. The solution to the slope safety factor

This section further evaluates the SF of the slope. In Figure 6, the dominant crack group propagated along the MN direction, the slope height is H, the slope angle is β, the slip surface inclination is ε, the unconnected length of the rock bridge is l_d, and the weight of the slip body is G. The dominant crack on the slip surface was simplified to an equivalent straight crack, and the crack connectivity rate, k, and the continuity of the rock bridge, 1-k, in the MN direction of the slip surface within the slope of the rock mass are defined as: (23) $$k\text{\hspace{0.17em}=\hspace{0.17em}}\frac{\sum a}{L_{MN}}\prime 1-k\text{\hspace{0.17em}=\hspace{0.17em}}\frac{\sum b}{L_{MN}}$$(23)

The crack propagated in the MN direction with the length of l_t+r_t, and the rock bridge penetration rate can be expressed as: (24) $$q\text{\hspace{0.17em}=\hspace{0.17em}}\frac{2\sum \left(l_t+r_t\right)}{2\sum \left(l_t+r_t\right)+\sum l_d}\text{\hspace{0.17em}=\hspace{0.17em}}\frac{2\sum \left(l_t+r_t\right)}{\sum b}\text{\hspace{0.17em}=\hspace{0.17em}}\frac{2\sum \left(l_t+r_t\right)}{\left(1-k\right)L_{MN}}$$(24)

From the geometric relationship in Figure 6, the expression for the weight of the sliding body on the side slope can be obtained as follows: (25) $$G\text{=}\frac{\left(H\cdot \mathrm{cot\hspace{0.17em}}\epsilon -H\cdot \mathrm{cot\hspace{0.17em}}\beta \right)\gamma H}{2}\text{\hspace{0.17em}=\hspace{0.17em}}\frac{\gamma H^2}{2}\cdot \frac{\hspace{0.17em}\sin \hspace{0.17em}\left(\beta -\epsilon \right)}{\hspace{0.17em}\sin \hspace{0.17em}\epsilon \hspace{0.17em}\sin \hspace{0.17em}\beta }$$(25)

The component of the slide weight, G, along the slip surface and the vertical slip surface are calculated using the following equation: (26) $$G_T\text{\hspace{0.17em}=\hspace{0.17em}}G\hspace{0.17em}\sin \hspace{0.17em}\beta \prime G_N\text{\hspace{0.17em}=\hspace{0.17em}}G\hspace{0.17em}\cos \hspace{0.17em}\beta$$(26)

The average normal stress, MN, on the slip surface is: (27) $${\sigma }_N\text{\hspace{0.17em}=\hspace{0.17em}}\frac{G\hspace{0.17em}\cos \hspace{0.17em}\beta }{L_{MN}}$$(27)

Based on the stress state of the slope in Figure 6, the maximum shear force that the slip surface, MN, can withstand, can be expressed as follows: (28) $$Q_T\text{\hspace{0.17em}=\hspace{0.17em}}\left(1-q\right)\left(1-k\right)\left(c+{\sigma }_N\hspace{0.17em}\tan \hspace{0.17em}\phi \right)L_{MN}+k\left(c_i+{\sigma }_N\hspace{0.17em}\tan \hspace{0.17em}{\phi }_i\right)L_{MN}+q\left(1-k\right)\left(c_i+{\sigma }_N\hspace{0.17em}\tan \hspace{0.17em}{\phi }_i\right)L_{MN}$$(28) where $c$ and $\phi$ are the cohesion and internal friction angle of the rock mass, respectively; $c_i$ and ${\phi }_i$ are the cohesion and internal friction angle of the original crack surface, respectively; $c_j$ and ${\phi }_j$ are the cohesion and internal friction angle of the propagation crack surface, respectively.

The rock bridge penetration rate of the slope safety factor, SF, can be expressed as follows: (29) $$SF=\frac{Q_T}{G_T}=\frac{\begin{array}{c}\left(1-q\right)\left(1-k\right)\left(c+G\hspace{0.17em}\cos \hspace{0.17em}\beta \hspace{0.17em}\tan \hspace{0.17em}\phi /L_{MN}\right)L_{MN}+k\left(c_i+G\hspace{0.17em}\cos \hspace{0.17em}\beta \hspace{0.17em}\tan \hspace{0.17em}{\phi }_i/L_{MN}\right)L_{MN}\hspace{0.17em}+\hspace{0.17em}\\ q\left(1-k\right)\left(c_i+G\hspace{0.17em}\cos \hspace{0.17em}\beta \hspace{0.17em}\tan \hspace{0.17em}{\phi }_i/L_{MN}\right)L_{MN}\end{array}}{G\hspace{0.17em}\sin \hspace{0.17em}\beta }$$(29)

When SF > l, it means that the maximum shear force that the rock bridge can withstand is greater than the sliding force acting on the crack surface, and the slope is in a stable state. When SF < l, it means that the maximum shear force that the rock bridge can withstand is less than the sliding force acting on the crack surface, and the slope is in an unstable state. When SF = l, it means that the maximum shear force that the rock bridge can withstand is equal to the sliding force acting on the crack surface, and the slope is in a critical instability state.

4.2. The analysis of the factors affecting the slope stability

The SF of the slope considering the influence of plastic zone at the end of branch cracks was relatively small compared to not considering it (Figure 7), and the increased plastic zone at the crack end made the sliding resistance decrease on the slope, the risk of landslide on the slope increased, which meant that the calculation of the slope safety factor was more reasonable by considering the influence of the plastic zone at the fissure tip. The relationship between the slope safety factor and different parameters is analyzed in Figure 7.

Figure 7. The influence of the plastic fractured region of the branch crack end on the SF.

Figure 7(a) shows the relationship between the SF of the slope and slope height, the SF decreased with the increased slope height, which indicated that the risk of the slope destabilization landslide increased with the increased slope height. Figure 7(b) shows the relationship between the SF of the slope and slope angle, the SF decreased with the increased slope angle. Figure 7(c) shows the relationship between crack angle and the SF of the slope, the SF first decreased and then tended to stabilize with the increased crack angle, which was because the fraction of the slip surface increased with increased crack angle, which led to the slope being more prone to instability damage. Figure 7(d) shows the relationship between the SF of the slope and crack length, and the slope safety factor decreased with the increased crack length. The increased crack length indicated that the degree of damage within the rock mass increased, which, thus, Accelerated the occurrence of slope landslides. Figure 7(e) shows the relationship between the the SF of the slope and the friction coefficient on the crack surface, the SF increased with the increased friction coefficient, which is due to the increase of the friction coefficient, and which made the fracture surface resist the crack propagation damage and the safety factor increased. Figure 7(f) shows the relationship between the the SF of the slope and the unloading factors, the SF decreased with the increased unloading factor, which indicated that the stronger the unloading disturbance, the worse the slope stability.

5. The engineering case analysis

The Daguoshan open-pit iron ore mine (Figure 8), located in Anshan City, Liaoning Province, has been mined for more than 60 years, forming the ‘Daguoshan Pit’, which is nearly 400 m deep. At present, the final bottom elevation of the second phase design realm is 450 m, the surface elevation is 70–190 m levels, and the maximum height of the design slope is about 650 m. The effect of the multiple geological and tectonic movements makes the joints and cracks more developed and shows a fractured structure, and the mechanical properties of the rock mass are closely related to the distribution of joints, which constitutes a more typical metal open pit joints slope. The slope of the lower plate of the southwest of the quarry is located between the elevation of the top of the slope of the ore body at +66 m level, with a step height of 12 m, an overall side slope angle of 43°, and the steep slope angle of 50°. The relevant physical and mechanical parameters of the slope rock mass are shown in Table 1.

Figure 8. Slope position and the layout crack measurement.

Table 1. Rock mechanics test results of rock samples in the study area.

Title	Rock density (γ)	Rock cohesion (c)	Crack cohesion (cj)	Crack internal friction angle (φj)	Poisson’s ratio (μ)	Lateral pressure coefficient
Value	28kN/m^3	3.5 MPa	0.1 MPa	20°	0.25	0.5

Under the excavation and unloading disturbances, the dominant crack group in the slope rock mass propagation and coalescence with each other, leading to rock mass collapse. The cracks in the rock mass of unloading slope are relatively developed. The crack length and inclination was obtained through the non-contact 3D photogrammetry, which was about 1.87 m long and with 50.2° inclination. The size and orientation of the dominant crack group were further determined through borehole photography (Figure 8). It is difficult to measure the small unloading tensile stress in the engineering site, generally through numerical simulation or indoor triaxial unloading experiments to obtain. The specific method can be divided into four steps. First, through the field experiment to measure the value of the ground stress at the side slope excavation and take rock samples. Second, in the laboratory through the real triaxial testing machine loaded crack specimens to the measured value of ground stress. Third, to keep the axial stress unchanged, the lateral began to slowly unload in stages, while the complete destruction process of the fissure specimens through high-speed cameras. Fourth, through the public ${\sigma }_3\text{–=–}{\sigma }_0-{\sigma }_{3t}$ can be found in the size of the unloading tensile stress σ₃, σ₀ is the value of the horizontal stress loaded to the beginning of the experiment, and σ_3t is the horizontal stress value recorded when the precast fracture initiation within the specimen is captured by the high-speed camera. These parameters obtained in the field and in the laboratory can be combined with Equation (29)(29) $$SF=\frac{Q_T}{G_T}=\frac{\begin{array}{c}\left(1-q\right)\left(1-k\right)\left(c+G\hspace{0.17em}\cos \hspace{0.17em}\beta \hspace{0.17em}\tan \hspace{0.17em}\phi /L_{MN}\right)L_{MN}+k\left(c_i+G\hspace{0.17em}\cos \hspace{0.17em}\beta \hspace{0.17em}\tan \hspace{0.17em}{\phi }_i/L_{MN}\right)L_{MN}\hspace{0.17em}+\hspace{0.17em}\\ q\left(1-k\right)\left(c_i+G\hspace{0.17em}\cos \hspace{0.17em}\beta \hspace{0.17em}\tan \hspace{0.17em}{\phi }_i/L_{MN}\right)L_{MN}\end{array}}{G\hspace{0.17em}\sin \hspace{0.17em}\beta }$$(29) to find the safety factor under the unloading effect of slope excavation, which can reflect the stability of the slope at all times.

The SIF (Section 2.4) at crack end after excavation unloading was greater than that before unloading, which made the initial crack propagation direction is almost the same as the direction of the slope surface. With the continuous mining of open pit mines, the stress environment within the slope rock mass continued to deteriorate, which resulted in the continued destabilization and penetration of the dominant cracks. The relationship between the number of excavation unloading and the accumulated propagation length of the crack are shown in Figure 9(a). In the early stage of slope excavation (with smaller slope height), the cumulative propagation length was small because the initial cracks were small, but when the cumulative length of newly formed dominant cracks increased, the cracks started to accelerate their propagation, with the result that if the sloping rock was not reinforced at this time, it would have significantly increased the possibility of landslides.

Figure 9. (a) The calculated crack propagation length and (b) the calculationed SF of the slope.

The safety factor of the slope decreases with the decrease of excavation unloading times. In the early stage of excavation, the decrease of SF is smaller, while in the later stage of excavation, the decrease of SF is larger, which is consistent with the results of crack propagation length (Figure 9(b)). The increased length of the new cracks significantly increased the probability of larger cracks coalescing the rock mass in a single propagation, which made the SF of slope decrease rapidly. Therefore, to be able to accurately detect the scale and size of the dominant initial crack group within the slope of the rock mass is the basis and key to slope safety control at the project site. Figure 9(b) shows that the (calculated) results are smaller when the plastic zone range is considered than otherwise, which verifies the accuracy of the theoretical derivation results in this study.

Further verification of the result in this study was done by on-site borehole camera and continuous displacement monitoring. Drilling cameras can monitor the propagation of internal cracks in slope rock masses. The continuous monitoring of displacement clarified whether the slope of the rock mass underwent a sudden change of displacement. The relevant monitoring settings and the displacement monitoring results are shown in Figure 10(a,b), respectively.

Figure 10. (a) Site-survey point layout and borehole peeping results and (b) displacement variation law of measuring points.

On-site borehole monitoring results indicate a clearly visible slide zone following the 14th excavation, and the displacement monitoring results show that the vertical displacement changed abruptly after the 14th excavation (Figure 10(b)), which indicated that the slope had already slipped at this point; with the most conservative calculations, the SF is less than 1 after 14 excavations, and the above results show that derivation process and conclusion is correct.

6. Conclusion

1. A mechanical model for unloading crack propagation was established based on fracture mechanics theory. The process of cracks in the unloading slope rock mass from initiation to the formation of the slope sliding surface was explained by mechanical analysis. The calculation formula for the SIF at the crack end of rock mass under unloading state has been obtained, and the SIF at the crack end in the unloaded state was larger than that in the unloaded state, which was explained from a quantitative perspective.

2. Combining fracture mechanics theory and the M-C rock damage criterion, the equations of branch crack propagation length and the range of plastic zone of crack end in the rock mass of unloading slope were derived, and the plastic zone range at the crack end exponentially functions with slope height, positively proportioned with crack length, inversely logarithmically functioned with slope angle, and periodically fluctuating with crack angle.

3. Dynamic process of crack propagation in slope rock mass based on excavation disturbance, the formula for calculating the SF of slope under the unloading conditions was derived, and the SF of slope was found to be more accurate when the range of the plastic zone at the branch crack end was considered than otherwise. Finally, the correctness of the theoretical calculation model was verified by engineering calculation examples.

7. Discussion

Although some useful conclusions are obtained in this study on the unloading of slopes, due to the extremely complex fracture mechanism of the unloading of the slope rock mass, there are still many problems to be further studied and explored:

1. How to apply the unloading theory of crack-containing rocks and the results of experimental research on fine-scale damage to the analysis of mining fracture damage of macroscopic slope rock mass will be a frontier research direction with high theoretical value.

2. The study of the fracture mechanism of slope rock mass under the coupling of multiple factors such as excavation unloading, rainfall infiltration and blasting vibration is also an important frontier development direction.

3. As the mining of mineral resources continues to deepen, the slope will gradually evolve into a high slope or even a super-high slope, and the mining disturbance will be more intense, so the research on the fracture mechanism of super-high slope will be the future development trend.

Author contributions

Conceptualization, J.W. (Jianming Wang) and Z.Z. (Zihan Zhou); funding acquisition, Z.C. (Zhonghui Chen); methodology, J.W. (Jianming Wang); writing—original draft, J.W. (Jianming Wang) and Z.Z. (Zihan Zhou); writing—review and editing, W.D. (Wei Dou).

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This research was funded by the National Natural Science Foundation of China (No. 52308316), China Postdoctoral Science Foundation (No. 2022M721885) and was supported by Key Laboratory of Rock Mechanics and Geohazards of Zhejiang Province (No. ZJRMG-2022-01).