Abstract
In order to explore the stability of multi-step high steep slope, this paper considers layer distribution of soil and builds a slope failure mechanism which meets requirements of velocity separation based on plastic mechanics principles, then derives the calculation method of external force power and internal energy dissipation power in failure area using numerical analysis theory, and proposes a loop algorithm for slope stability analysis combining upper bound limit analysis theorem and strength reduction principle. At the same time, in order to enhance the applicability of the upper bound limit analysis method, this paper develops a slope stability calculation system of upper bound limit analysis combined with computer programming technology and evaluates the accuracy of the calculation system through slope engineering practice in open-pit mine. The results show that stability factor by calculation system of upper bound limit analysis is an upper limit solution compared with limit equilibrium method, but the error rate of two methods is less than 5%, which indicates that the slope stability calculation system established in this paper has high accuracy in stability analysis of multi-step high steep slope in the practice of open-pit mine. On the other hand, the most dangerous sliding surface searched by upper bound limit analysis can fully solve velocity problems of multi-step high steep slope failure, so it has stronger reference values in engineering practice.
1. Introduction
In recent years, with the vigorous development of civil engineering construction, multi-step high steep slopes are becoming more common in highways, water conservancies, mines, constructions, and other fields [1–3]. This kind of slope is characterized by wide distribution, poor stability, and difficult landslide management, and its stability has become an important issue to be urgently solved in slope practice [4, 5].
The common methods of slope stability analysis mainly include the limit equilibrium method and limit analysis method [6, 7]. In the 1920s, Fellenius proposed the Swedish circular arc method, which has established the beginning of the limit equilibrium method [8–10]. In the following decades, the simplified Bishop method, Morgenstern–Price method, Sarma method, simplified Janbu method, residual thrust method, and other limit equilibrium methods were successively proposed [11–16]. Over nearly a century’s development, the limit equilibrium method has developed from a simplified empirical algorithm to a mature method with complete theoretical systems, which is widely used in slope engineering practice [17]. These abovementioned methods only met parts of static equilibrium conditions and moment equilibrium conditions, but the relationship between stress and strain of soil was not considered effectively. A series of assumptions must be introduced to make the problem solvable, and these assumptions will lead to errors in the calculation [18, 19]. The limit analysis method based on plastic principles has gradually become a research hotspot in slope engineering practice in order to overcome deficiencies of the limit equilibrium method. Furthermore, upper bound limit analysis seeks the slope failure mechanism by establishing movement permissible velocity field, which has a strict theoretical foundation and a rigorous derivation process and can fully solve the velocity problems of slope failure [20–22]. In the 1970s, Chen led upper bound limit analysis into geotechnical engineering practices and proposed upper bound limit analysis for slope stability calculation [23]. In the 1990s, Michalowski divided sliding body within slope failure areas vertically and put forward upper bound limit analysis for translational failure forms [24]. On this basis, Donald divided sliding body within slope failure areas slantwise and put forward upper bound limit analysis for combined translational failure forms and inclined failure forms [25].
Obviously, upper bound limit analysis has its unique advantages in slope stability analysis [26]. It provides a good solution to the problem of homogeneous slope stability [27], but the applicability of upper bound limit analysis is greatly limited for complex terrain conditions due to the limitation of solving process [28, 29]. Aiming at this problem, on the basis of considering the layered distribution of soil, this paper researches on the stability of multi-step high steep slope in engineering practice, proposes upper bound limit analysis algorithm of multi-step high steep slope stability, and develops a widely applicable analysis system based on computer programming technology. The results will enhance the applicability of the upper bound limit analysis method in slope engineering practice, which has important significance and popularization values.
2. Stability Analysis of Multi-Step High Steep Slope
2.1. Slope Morphology and Soil Parameters
For the high steep slope composed by n steps, the morphological parameters of each step are defined as follows. Height of each step is h1, h2, …, hn, slope angle of each step is α1, α2, …, αn, and plate width of each step is b1, b2, …, bn−1. At the same time, considering the layered distribution characteristics of slope soil in engineering practice, the slope consists of m layers soil, and the mechanical parameters of soil in each layer are defined as follows from top to bottom. Thickness of each soil layer is h1′, h2′, …, hm′, density of each soil layer is γ1, γ2, …, γm, cohesion of each soil layer is c1, c2, …, cm, and internal friction angle of each soil layer is φ1, φ2, …, φm. The slope morphology and soil parameters are shown in Figure 1.

Obviously, for the high steep slope consisting of n steps, when the morphological parameters of each step are known, the slope height and the slope angle can be calculated through plane geometric relationship, and the slope morphology can be uniquely obtained.
2.2. Calculation Model Establishment of High Steep Slope
2.2.1. Mathematical Model Establishment of Failure Mechanism
The plastic theory shows that when shear failure occurs, the slip surface morphology is logarithmic spiral [30, 31]. For the high steep slope in Figure 1, shear failure may occur at toe point of each step A1, A2, …, An. This paper assumes that shear failure occurs at the point of An to illustrate the calculation process, and shear failure occurs at other points which can be calculated in the same way. The failure mechanism is formed by intersecting parts of m logarithmic spirals when shear failure occurs at point of An. At the same time, the failure mechanism starts at the top point of slope B0 and passes through the toe point of slope An. The failure mechanism is shown in Figure 2.

In Figure 2, r (m) is the line segment distance between the rotation center and the logarithmic spiral, and θ (rad) is the horizontal angle for the line between the rotation center and the logarithmic spiral. According to unique velocity of the slope failure mechanism at soil layer interface, it can be obtained that the combined logarithmic spiral has the same angular velocity ω and the same rotation center O, so the equation of the slope failure mechanism can be expressed as
Establish a searching area of rotation center above the shoulder point of the slope; for any point in the searching area, polar coordinate of the toe point An (rm, θm) is definite. Then, for any definite rm and θm, the parameter value am in the mth logarithmic spiral can be calculated. According to coordinate consistency of the failure mechanism at the soil layer interface, the parameters am-1, am−2, …, a2, a1 can be calculated layer by layer from bottom to top. On this basis, the slope failure mechanism equation can be obtained. So, the slope failure mechanism equation can be uniquely controlled by the position of rotation center.
2.2.2. Mathematical Model Establishment of Slope Morphology
According to the above introductions, for any point in the searching area, polar coordinate of the toe point An (rm, θm) is definite, so the plane rectangular coordinate of the toe point An can be expressed as . Therefore, slope equation of the nth step can be expressed as formula (2) according to slope angle of the nth step which is αn.
According to the step distribution morphology, the plane rectangular coordinate of the toe point of the n−1th step An−1 can be expressed as (rm · cos θm + hn · cot αn + bn-1, −rm · sin θm + hn). With the same reasoning as above, slope equation of the (n−1)th step can be expressed as
In the same way, slope equation can be obtained according to the plane rectangular coordinate of toe point and slope angle of each step, so the slope morphology equation can be obtained only by the position of rotation center.
Meanwhile, when the polar coordinate of the toe point An is (rm, θm), the floor equation of the mth layer soil can be expressed as y = −rm · sin θm. So, roof equation of the mth layer soil can be expressed as formula (4) according to the mth layer soil thickness which is .
With the same reasoning as above, according to the (m−1)th layer, soil thickness is , and roof equation of the (m−1)th layer soil can be expressed as
Similarly, floor equation and roof equation can be obtained according to each layer soil thickness, so floor equation and roof equation of slope layer soil can be determined only by the position of rotation center.
To sum up, failure mechanism, slope morphology, floor equation, and roof equation can be determined only by the position of rotation center. Then, the mathematical model of failure mechanism and slope morphology can be established through the position of rotation center.
2.3. Upper Bound Limit Analysis for Slope Stability
The upper bound limit analysis method analyzes slope stability from the way of energy and velocity. The basic theory of this method is that for any slope failure mechanism, when external force power is equal to internal energy dissipation power, the unsafe upper bound of failure loading can be obtained, and then the slope stability can be analyzed. Therefore, this paper takes the position of rotation center as initial parameter and proposes calculation methods of external force power and internal energy dissipation power of the slope. Furthermore, the upper bound limit analysis method for slope stability can be established by changing the position of the rotation center.
2.3.1. Calculation Method of External Force Power
In natural state, external force power of the failure area is all provided by gravity power [32, 33]. For the failure mechanism shown in Figure 2, the step plates are lengthened in reverse and intersect the failure mechanism (assuming that thickness of the mth soil layer is higher than height of the nth step). In the vertical direction, the step plates and soil layer interfaces divide the slope failure area into several horizontal strips. For the lowest horizontal strip AnAn−1C, its failure region S in the planar rectangular coordinate system can be expressed aswhere y0 is the ordinate of the toe point, which can be solved by the formula is the slope equation of the nth step, it is shown in formula (2), and is the logarithmic spiral equation of the mth soil layer. On the premise that the position of rotation center is known, and can be uniquely determined. Therefore, external force power E′ of the lowest horizontal strip can be expressed aswhere E′ (kW) is the external force power of the lowest horizontal strip in the slope failure area, S (m2) is the failure region of the lowest horizontal strip, dS (m2) is the arbitrary area element, γ (kN/m3) is the soil density of the area element, and is the vertical direction velocity of the area element. According to the plane geometry, vertical direction velocity of the area element can be expressed as
Therefore, formula (7) can be converted intowhere is a linear function about y, so the definite integral of can be solved. However, because the function expression of fn(y) cannot be determined, the definite integral of cannot be solved directly. In this paper, the complex Simpson method is applied to approximate calculation. Divide the lowest horizontal strip into 2k equal parts along the vertical direction and extract the coordinates of each equal point, denoted as (x1, y1), (x2, y2), …, (x2k−1, y2k−1). Then, the definite integral of fn2(y) can be approximated aswhere xAn and xC represent the corresponding abscissa of point An and point C, x1, x2, x3, …, x2k−1 represent the corresponding abscissa of the 1st, 2nd, 3rd… 2k−1th equal points, and the abscissa of each equal point can be expressed as
The abscissa of each equal point is known because the failure mechanism is unique, so the definite integral of fn2(y) can be solved and external force power of the lowest horizontal strip in the failure region can be calculated. Similarly, external force power of each horizontal strip can be calculated from bottom to top, and external force power E of the failure area of high steep slope can be finally determined by adding external force power of each strip.
2.3.2. Calculation Method of Internal Energy Dissipation Power
Upper bound limit analysis assumes that the slope soil is rigid material and the velocity of soil inside and outside the slope slip surface is continuous. Obviously, the velocity discontinuity only occurs on the slip surface, the same with combined logarithmic spiral [34–36]. Therefore, internal energy dissipation power in failure area of high steep slope concentrates on intersecting parts of m logarithmic spirals. For the soil in the mth layer, rdθ · (cos φm)−1 on the logarithmic spiral is selected as element length, and the corresponding cohesion force of the element length can be expressed as cmrdθ · (cos φm)−1, and then internal energy dissipation power of the mth soil layer can be expressed aswhere Dm (kW) is the internal energy dissipation power of the mth soil layer in the slope, θm−1 (rad) is the horizontal angle for the line between rotation center and starting point of failure mechanism in the mth soil layer, and rm−1 (m) is the length for the corresponding line. Obviously, on the premise that the position of rotation center is known, both rm−1 and θm−1 can be obtained. Therefore, internal energy dissipation power of the mth soil layer in the slope can be solved according to formula (12). Similarly, internal energy dissipation power of each soil layer can be calculated from bottom to top, and internal energy dissipation power of the high steep slope D can be determined by adding internal energy dissipation power of each soil layer.
2.3.3. Strength Reduction Algorithm for Slope Stability
For any imaginary failure mechanism, when external force power is greater than internal energy dissipation power, the slope soil cannot bear the applied load value, and then the failure will occur. The upper bound limit analysis method assumes that when external force power is just equal to internal energy dissipation power, the slope is in a limit equilibrium state, and the slope stability analysis can be calculated according to this limit equilibrium state [37–39]. For the slope in nonlimiting equilibrium state, repeated strength reduction is required to gradually transition the slope to the limit equilibrium state. The reduction formula can be expressed aswhere F is the strength reduction coefficient, c and φ are cohesion and internal friction angle of the slope soil before strength reduction, and and are cohesion and internal friction angle of the slope soil after strength reduction. This paper applied computer programming technology to realize strength reduction calculation of slope stability.
2.3.4. Slope Stability Analysis System Development
The strength reduction principle shows that when the strength reduction coefficient F is equal to the stability factor Fs, it indicates that slope is in a limit equilibrium state, and external force power E is exactly equal to internal energy dissipation power D, which means (E−D)max = 0, so build the following function:
It is clear that Δ has only one zero point in the range of arithmetic numbers. Meanwhile, F and Δ are negatively correlated. According to the above two characteristics, a slope stability analysis system is developed by using the dichotomy method. The calculation process of slope stability analysis system is shown in Figure 3.

In Figure 3, the closed interval formed by ai and bi should strictly contain the slope stability factor at the condition that i = 1. Through Figure 3, the strength reduction coefficient F can be output, which is the slope stability factor Fs. Furthermore, the stability factor is accurately 10−3 and can fully meet the requirement of engineering practice. At the same time, the polar coordinate of the toe point can be output through Figure 3, and the slope slip surface can be drawn according to the toe point, which is the most dangerous slip surface of the slope.
3. Engineering Case Analysis
Due to requirements of mining and transportation, open-pit mine slope is usually designed in the form of multi-stage steps. In this paper, the complex slope in open-pit mine composed of dump slope and mining slope is selected as engineering background, and the accuracy of the stability analysis system is evaluated.
3.1. Slope Geology
The complex slope is composed of dump slope and mining slope. Dump slope is composed of two dump steps, height of each step is 24 m, slope angle of each step is 35°, and plate width of each step is 50 m. Mining slope is composed of six mining steps, height of each step is 12 m, slope angle of each step is 65°, and plate width of each step is 35 m. The distance between dump slope and mining slope is 35 m, and soil layer of the complex slope is composed of abandoned material, quaternary, tertiary, and mud stone from top to bottom. The slope morphology is shown in Figure 4, and the mechanical parameters of the slope soil are shown in Table 1.

3.2. Stability Analysis Results
Apply the analysis system described above to analyze and evaluate the slope stability. At the same time, in order to ensure stability of the slope, toe points of all steps A1, A2, …, A8 are selected as the slope toe point, and analyze the slope stability of the above eight conditions. The stability factor distribution in the searching area is shown in Figure 5.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)
It can be seen from Figure 5 that when A8 is selected as slope toe point, the corresponding stability factor Fs = 1.294, and the stability factor reaches the minimum value, so the slope stability factor is 1.294. At the same time, the most dangerous sliding surface of the slope can be drawn according to the polar coordinates of the toe point, and the most dangerous sliding surface is shown in Figure 6.

3.3. Calculation Result Evaluation
The limit analysis method and limit equilibrium method analyze slope stability from different directions, but their calculation results should be highly consistent [40, 41]. This paper selects simplified Bishop method in the limit equilibrium method to evaluate accuracy of calculation results. The slope stability analysis results are shown in Table 2, and the most dangerous sliding surface is shown in Figure 6.
From the perspective of slope stability factor, the error rate of calculation results of two methods is less than 5%, which can fully meet the requirements of slope practice. It indicates that the slope stability calculation system established in this paper has high accuracy in stability analysis of multi-step high steep slopes in the practice of open-pit mines.
On the other hand, from the perspective of the most dangerous sliding surface, the position of starting point and ending point of landslide in the calculation results of two methods is approximately the same. But the most dangerous sliding surface searched by upper bound limit analysis can fully solve velocity problems of multi-step high steep slope failure [42], so it has stronger reference values in engineering practices.
4. Conclusion
This paper proposes a calculation method of upper bound limit analysis for multi-step high steep slope stability, develops a cycle program of stability analysis system, and then evaluates the accuracy of analysis system through engineering practice. The main conclusions are as follows:(1)The layered distribution law of slope soil in engineering practice is analyzed, and according to the basic principle of plastic mechanics, the failure mechanism of multi-step high steep slope which can meet the requirement of velocity separation is constructed.(2)The calculation method of external force power and internal energy dissipation power in the failure area is derived, and combined with upper bound limit analysis theorem and strength reduction principle, upper bound limit analysis cycle algorithm suitable for the stability of multi-step high steep slope is proposed.(3)The slope stability analysis system is developed, and the accuracy of the analysis system is evaluated according to the engineering practice, which fully verifies the applicability of the analysis system in slope engineering practice.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The author declares that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This study was supported by the Key Research Project for Higher Education in Henan Province (22B560005) and the Scientific and Technological Project in Henan Province (202102310567).