#### Abstract

Slope stability analysis is a core issue in geotechnical engineering. This paper proposes a method of upper bound limit stability analysis for a slope with multiparameter coordinated variation based on the comprehensive consideration of the nonuniform distribution of slope soil parameters. This method starts from the perspective of energy balance, establishes a slope failure mechanism which meets velocity separation requirements, deduces its calculation formula for external force power and internal energy dissipation power, develops a cycle program for the most dangerous slip surface searching and stability coefficient calculation through computer programming technology, and finally forms a calculation method of upper bound limit stability analysis for the soil slope with nonuniform multiparameter distribution. At the same time, this method takes a dump slope in an open-pit mine as the engineering background, considers the nonuniformity of density, cohesion, and internal friction angle of the slope soil under subsidence, applies upper bound limit analysis to analyze the slope stability, and evaluates the accuracy of analysis results by using the residual thrust method. The results show that upper bound limit analysis has highly accuracy in stability coefficient calculation; compared with the residual thrust method, the stability coefficient calculation result by upper bound limit analysis is a strict upper bound solution, and the calculation error is easy to be estimated and eliminated. Simultaneously, the most dangerous slip surface obtained by upper bound limit analysis can fully satisfy the velocity separation requirement and has a greater engineering reference value.

#### 1. Introduction

Slope stability analysis is a core issue in geotechnical engineering [1]. Commonly used slope stability analysis methods mainly include the limit equilibrium method and the limit analysis method. Limit equilibrium method analyzes slope stability through a static equilibrium equation, and this method has certain limitations because the stress-strain relationship of soil is not effectively considered [2, 3]. On the contrary, the limit analysis method analyzes slope stability from the perspective of energy balance, and the analysis results are highly reliable because of its strict theoretical foundation and rigorous derivation process [4–6].

In the recent years, with the rapid development of plastic mechanics, the limit analysis method has become an important measure in slope stability analysis [7]. Among them, upper bound limit analysis seeks the failure mechanism of the slope by constructing a kinetically admissible velocity field, and this method has the advantages of simple calculation process and wide application in engineering [8, 9]. In 1975, Chen introduced upper bound limit analysis into the geotechnical engineering field for the first time, and then he published a detailed treatise on the application of upper bound limit analysis in foundation bearing capacity, soil pressure of the retaining wall, slope stability, and the like [10]. On this basis, Michalowski proposed upper bound limit analysis of the translational kinetic field through vertically dividing the slip mass in the slope failure area [11]. Donald divided the slip mass in the slope failure area into strips diagonally and proposed upper bound limit analysis combining the diagonal strip division and the translational kinetic field, which broadened the application of upper bound limit analysis [12]. Wang proposed upper bound limit stability analysis for a slope in seismic regions based on the effects of earthquakes [13]. Yang proposed upper bound limit stability analysis for a slope under the seepage effect based on seepage effect of the drainage layer [14]. Upper bound limit analysis has given a good solution to homogeneous slope stability; however, the application of upper bound limit analysis for a heterogeneous slope is extremely limited due to the limitation of the solution process [15–17]. In natural conditions, mechanical parameters of slope soil are not uniformly distributed due to natural forces such as weathering, sedimentation, and rainfall. The nonuniformity is particularly obvious in accumulation slopes. In this regard, Kumar considered the sedimentation effect of the slope and proposed upper bound limit analysis for slope stability with nonuniform distribution of density [18]. Nian et al. proposed upper bound limit analysis for slope stability with nonuniform distribution of cohesion based on weathering and rainfall effects of the slope [19]. In the existing literatures, upper bound limit analysis of heterogeneous slopes only considers the nonuniform distribution of density or cohesion; therefore, upper bound limit analysis for slope stability needs to be further improved [20]. First of all, a large number of engineering projects show that cohesion is negatively correlated with the internal friction angle for slope soil; that is to say, when cohesion increases along the depth below the slope, the internal friction angle must be gradually reduced [21–23]. When the internal friction angle is nonuniformly distributed, the slip surface shape will also change according to associated flow rules, so the application of upper bound limit analysis is bound to be severely limited [24–26]. As mechanical parameters of soil mass change in coordination under the action of natural forces, a large number of scholars proposed that in order to reflect the distribution characteristics of soil mechanical parameters more realistic, we should explore upper bound limit analysis for the true state of soil by considering the characteristics of nonuniform multiparameter distribution [27–30].

Based on the above analysis, this paper creatively proposes an upper bound limit stability analysis for a soil slope with coordinated variation of multiparameter. The method discretizes the slope failure mechanism, and each discrete block can fully meet the associated flow rules. At the same time, the method takes the whole slip surface as the research object without establishing a complicated velocity field equation for the internal slip surface and has the advantages of simple calculation process, accurate calculation result, strong engineering applicability, and the like, thus providing a new idea for heterogeneous slope stability.

#### 2. Upper Bound Limit Analysis of Slope Stability

##### 2.1. Nonuniformity Analysis of Slope Soil

Soil is a kind of a natural product, and so its mechanical parameter distribution is extremely complicated. Among them, the soil mechanical parameters of accumulation slopes are approximately linearly distributed along the depth below the slope due to natural sedimentation [31–34]. For a typical accumulation slope, assume that the slope height is *h* and the slope surface angle is *α.* Density, cohesion, and internal friction angle are, respectively, *γ*_{0}, *c*_{0}, and *φ*_{0} for the soil at the top of the slope. In the depth below the slope, the variable quantity of density, cohesion, and internal friction angle is, respectively, *λ*_{γ}, *λ*_{c}, and *λ*_{φ} for per unit length (the soil internal friction angle decreases linearly with increasing depth below the slope due to the negative correlation between cohesion and the internal friction angle). The variation of soil mechanical parameters along the depth below the slope is shown in Figure 1.

**(a)**

**(b)**

**(c)**

##### 2.2. Establishment of Discrete Model for Failure Mechanism

According to the distribution of slope soil mechanical parameters, the slope failure mechanism is dispersed into the *n* horizontal block in depth below the slope (*n* approaches to infinity), so the discrete failure mechanism of the slope is composed of *n* logarithmic spirals according to velocity separation requirements and associated flow rules [35, 36]. The discrete failure mechanism of the slope is shown in Figure 2.

According to the uniqueness of velocity for each point of the failure mechanism, it can be determined that the discrete block has the same angular velocity and center of rotation. Obviously, each discrete block can be approximated as a homogeneous block soil when *n* approaches to infinity. Therefore, density, cohesion, and internal friction angle of no. *i* (*i* = 1, 2, …, *n*) discrete block can be expressed aswhere *γ*_{i} (kN/m^{3}), *c*_{i} (kPa), and *φ*_{i} (°) are the density, cohesion, and internal friction angle of the *i* discrete block. At the same time, the failure mechanism equation of the *i* discrete block can be expressed aswhere *r* (*m*) is the distance between the center of rotation and the logarithmic spiral, *a*_{i} is the pending parameter in the logarithmic spiral equation, and *θ* (rad) is the horizontal angle between the center of rotation and the logarithmic spiral. For convenient calculation, a rectangular plane coordinate system is established with the center of rotation as the origin of coordinates, horizontal direction as the *x*-axis (rightward is positive direction), and vertical direction as the *y*-axis (downward is positive direction). In the rectangular plane coordinate system, assume that the slope surface equation is (*y*) = cot*α*·*y* + *b* and the slope top surface equation is *y* = *y*_{0}, so the slope bottom surface equation is *y* = *y*_{0} + *h*. Regarding the *i* discrete block, the failure mechanism equation in the rectangular plane coordinate system can be expressed aswhere *x* and *y* are the abscissa and the ordinate of the failure mechanism of the *i* discrete block in the rectangular plane coordinate system. Formula (3) shows that *x*/*a*_{i} and *y*/*a*_{i} are implicit functions of *θ* onl*y*, and there is a great difficulty in the manifestation of the implicit function. This paper discretizes the value of *θ* in the domain with 0.001 as the effective step size, calculates the rectangular plane coordinates of each discrete point, and applies the least square method to fit polynomial discrete point coordinates. The fitting equation of the failure mechanism of the *i* discrete block can be expressed as *x* = *f*_{i} (*a*_{i}, *y*), so the failure area of the *i* discrete block can be expressed as

##### 2.3. Derivation of Energy Balance Equation

Upper bound limit analysis analyzes slope stability from the perspective of energy balance can effectively avoid inconsistent analysis results caused by different mapping methods. When the slope is in the limit equilibrium state, external force power should be equal to internal energy dissipation power in the slope failure area. In natural conditions, external force power is provided by gravity power only [37], so the calculation formula of external force power in the slope failure area can be expressed aswhere *E* (kW) is the external force power in the slope failure area, *S* (m^{2}) is the acreage of the slope failure area, d*S* (m^{2}) is the arbitrary microelement area in the slope failure area (m/s) is the component velocity of the microelement in the direction of gravity, and *γ* (kN/m^{3}) is the density of the microelement. According to the law of velocity distribution in the slope failure area, component velocity in the direction of gravity can be expressed aswhere *ω* (rad/s) is the angular velocity of the slope failure area. Combined with the failure area equation, external force power of the *i* discrete block can be expressed aswhere *E*_{i} (kW) is the external force power of the *i* discrete block. Therefore, external force power of the slope failure area can be expressed as

There are two assumptions about slope soil when calculating the internal energy dissipation power: the slip mass is regarded as a rigid body; the slip surface is regarded as a velocity intermittent surface [38, 39]. Based on the above assumptions, there is no internal energy dissipation power in the inner side and the outer side of the slip surface, and internal energy dissipation power is concentrated on the velocity intermittent surface, that is, logarithmic spirals [40]. At the same time, the frictional energy power and the assumed dilatant energy power can be completely offset, so internal energy dissipation is provided by cohesive power only [41]. Regarding the *i* discrete block, *r*d*θ*·(cos*φ*_{i})^{−1} in the logarithmic spiral is selected as the microelement length, and the cohesion corresponding to the microelement length is *c*_{i}·*r*d*θ*·(cos*φ*_{i})^{−1}. Internal energy dissipation power of the *i* discrete block can be expressed aswhere *D*_{i} (kW) is the internal energy dissipation power of the *i* discrete block, *θ*_{i−1} (rad) is the horizontal angles between the center of rotation and the failure mechanism starting point of the *i* discrete block, *θ*_{i} (rad) is the horizontal angles between the center of rotation and the failure mechanism end point of the *i* discrete block, *r*_{i−1} (*m*) is the distance between the center of rotation and the failure mechanism starting point of the *i* discrete block, and *r*_{i} (*m*) is the distance between the center of rotation and the failure mechanism end point of the *i* discrete block. According to the calculation formula for distance between two points in the rectangular plane coordinate system, *r*_{i-1} and *r*_{i} can be expressed as

Internal energy dissipation power of the *i* discrete block can be obtained with the calculation formula for distance between two points in the rectangular plane coordinate system, so internal energy dissipation power of the slope failure area can be expressed as

Until now, external force power and internal energy dissipation power of the slope failure area can be obtained. In order to establish an energy balance equation, constraint conditions should be determined according to the slope failure law. First of all, each discrete point of the slope failure mechanism should have consistent coordinates, and the constraint condition can be expressed as

Secondly, the failure mechanism intersects the slope surface at the slope toe, and this constraint condition can be expressed as

At last, the failure mechanism should be located inside the slope surface except the point at the slope toe, and this constraint condition can be expressed as

The slope energy balance equation can be finally established with the energy calculation method and equation constraint conditions.

##### 2.4. Strength Reduction Calculation

When the slope is in the limit equilibrium state, external force power should be equal to internal energy dissipation power in the slope failure area. A gradual transition to the limit equilibrium state through repeated strength reduction is needed for the slope under the nonlimit equilibrium state [42, 43]. This paper applies the strength reduction method to analyze slope stability, and the strength reduction formula iswhere *F* is the strength reduction coefficient, *c* (kPa) and *φ* (°) are the initial soil cohesion and the internal friction angle of slope soil and *c*′ (kPa) and *φ*′ (°) are the soil cohesion and the internal friction angle of slope soil after strength reduction. Internal friction angle needs to be reduced, which will affect the fitting results of the least square method mentioned above. Therefore, it is extremely complicated to reduce to limit equilibrium state by artificial calculation, and computer technology is required to achieve strength reduction.

##### 2.5. Compilation of Cycle Program

According to the definition of the stability coefficient, when the strength reduction coefficient *F* is equal to the stability coefficient *Fs*, external force power should be equal to internal energy dissipation power in the slope failure area; that is to say, the minimum ratio between external force power and internal energy dissipation power is equal to 1. At the same time, the minimum ratio has a strict positive correlation with the strength reduction coefficient. Based on the above characteristics, this paper applies the dichotomy method to solve the slope stability. Based on an established function *ε* = (*D*/*E*)_{min}, it is obvious that *ε*-1 has only one null point in the range of a positive real number. Based on practice demand of slope engineering, the calculation results of the slope stability coefficient are accurate to 10^{−3}. The calculation process of upper bound limit analysis for slope stability is shown in Figure 3.

In Figure 3, the interval formed by initial values *a*_{1} and *b*_{1} should include the slope stability coefficient. The strength reduction coefficient *F* can be finally output through the calculation process shown in Figure 3, and the strength reduction coefficient *F* is the slope stability coefficient. At the same time, the minimum value is obtained based on the ratio of internal energy dissipation power and external force power, and the shape of slope slip surface can be drawn based on the value of each parameter, and this slip surface is the most dangerous slip surface.

#### 3. Case Analysis

##### 3.1. Engineering Overview

This paper takes a dump slope in an open-pit mine as the engineering background. The dump slope is composed of loose sandstone, the slope height is 150 meters, and the slope surface angle is 20°. Due to soil sedimentation, density, cohesion, and internal friction angle of soil at the top of the slope are, respectively, 18 kN/m^{3}, 20 kPa, and 21°. At the bottom of the slope, density, cohesion, and internal friction angle are, respectively, 24 kN/m^{3}, 65 kPa, and 18°. In the depth below the slope, the variable quantity of density, cohesion, and internal friction angle is, respectively, 0.04 kN/m^{3}, 0.3 kPa, and 0.02° for per unit length. The slope does not have large cracks and slips, and the slope has been stable since formation from the perspective of geological characteristics and monitoring data and slope state analysis. The slope shape and mechanical parameter change law are shown in Figure 4.

According to the calculation process of upper bound limit analysis above, the slope stability with nonuniform distribution of each parameter can be obtained. The relationship between the number of discrete blocks and the slope stability coefficient is shown in Figure 5.

It can be shown from Figure 5 that with the increasing number of discrete blocks in the slope failure area, the slope stability coefficient converges gradually, and the final convergence value is the slope stability coefficient. At the same time, the minimum value is obtained based on the ratio of internal energy dissipation power and external force power, and the corresponding parameter value can be used to draw the most dangerous slip surface of the slope. Four significant digits in the calculation result of the slope stability coefficient in this paper are retained; when the number of discrete blocks is between 15 and 25, calculation results of the slope stability coefficient are accurate to 10^{−3}, which can fully meet the calculation accuracy. Upper bound limit analysis calculation results of the slope stability coefficient are shown in Table 1, and the most dangerous slip surface shape is shown in Figure 6.

##### 3.2. Accuracy Evaluation

Upper bound limit analysis analyzes slope stability from the perspective of energy balance, while the limit equilibrium method analyzes slope stability through a static equilibrium equation. Two methods are quite different in the analysis perspective, but the results should be highly consistent. Therefore, this paper selects the residual thrust method to evaluate the accuracy of results of upper bound limit analysis. Residual thrust method assumes that the thrust direction of the current block to the next block at the interface is parallel to the bottom slip surface of the block, and the values based on zero resultant force in the two directions of the parallel bottom slip surface and the vertical bottom slip surface and the zero residual thrust of the leading edge strip are iteratively obtained [44–46]. Residual thrust method calculation results of the slope are shown in Table 1, and the most dangerous slip surface shape is shown in Figure 6.

From the perspective of soil parameter influence on the slope stability coefficient, compared with a homogeneous slope (the slope is a homogeneous slope when the number of discrete block is 1), linear distribution of density and cohesion will improve the slope stability, but linear distribution of internal friction angle will reduce the slope stability. Among these three parameters, linear distribution of the internal friction angle has the most significant influence on the slope stability; therefore, the reason for reduction in accumulation slope stability is closely related to the linear distribution of the internal friction angle. In engineering practice, keeping relatively a stable internal friction angle is important for accumulation slope stability.

From the perspective of calculation results of the stability coefficient, since the residual thrust method assumes that the thrust direction of the current block to the next block at the interface is parallel to the bottom slip surface of the block, shear force at the interface may exceed shear strength of the slope when the bottom slip surface of the block is steep, and it may cause stability coefficient calculation results to be critical. To the opposite, shear force at the interface may be far less than shear strength of the slope when the bottom slip surface of the block is flat, and it can fully meet reasonable conditions but may cause stability coefficient calculation results to be too conservative. That is to say, the stability coefficient obtained by the residual thrust method is not an upper bound limit solution or a lower limit solution in strict sense; calculation errors are difficult to be estimated and eliminated although the calculation results have certain accuracy, while upper bound limit analysis takes the velocity field satisfying the velocity boundary condition and the volume invariable condition as the admissible function; the obtained functional value is an upper bound limit value, which results in the greater value shown in the calculation results of the stability coefficient by upper bound limit analysis, and calculation errors are easily estimated and eliminated. The calculation result error rate of two methods is less than 5%, which meet the requirements of engineering practice and fully verify the accuracy of the calculation results. At the same time, the calculation results show that the slope is stable, which is highly consistent with engineering analysis.

From the perspective of the most dangerous slip surface, the starting point and the end point of the most dangerous slip surface obtained by two methods are basically the same, and the slip surface shape is similar. However, the most dangerous slip surface shape obtained by the residual thrust method is a shear arc slip surface, and the most dangerous slip surface shape obtained by upper bound limit analysis is a combined logarithmic spiral slip surface. According to the basic theory of plastic mechanics, the most dangerous slip surface obtained by upper bound limit analysis can fully meet velocity separation requirement. Therefore, the most dangerous slip surface obtained by upper bound limit analysis has a greater engineering reference significance.

#### 4. Conclusion

Main conclusions of this paper are as follows:(1)A slope discrete failure mechanism that meets the velocity separation requirements is established based on nonuniformity analysis of slope soil.(2)A calculation method of external force power and internal energy dissipation rate is proposed, and calculation process of slope stability upper bound limit analysis method is compiled.(3)The engineering case is used to calculate the stability of the slope with nonuniform multiparameter distribution by upper bound limit analysis, and the residual thrust method is used to verify the accuracy of the calculation results.

#### Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

#### Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

#### Acknowledgments

This paper was supported by the Key Research Project for Higher Education in Henan Province (20B560009) and the Scientific and Technological Project in Henan Province (202102310567).