Abstract
The approximate assumptions of limit equilibrium methods are the fundamental cause of the deviation between their calculation results and actual situation. This study proposes a new finite-element evaluation method to reflect the progressive failure characteristics of rock slopes. By comparing the results of limit equilibrium and finite-element methods, the influence factors of stability for planar landslides are systematically studied. The factors include the plastic parameters of sliding surfaces in progressive failure, the elastic parameters of sliding mass, the elastic deformation of sliding beds, and excavation stress release. Moreover, the stress distribution rules on sliding surfaces and the diversity of safety factors are explored. Finally, the error source and calculation accuracy of the limit equilibrium method in slope analysis are obtained. The study provides scientific references for analyzing and evaluating the stability of such slopes.
1. Introduction
With the development of computers, numerical techniques are widely used in the analysis of landslides and engineered slopes, especially the application of the finite-element (FE) method. However, conventional limit equilibrium (LE) approaches still play important roles. These methods have clear mechanical concepts and a wide variety of software available for different failure modes (planar, wedge, toppling, etc.). The calculation results of LE methods are exact analytical solutions under the conditions that satisfy their basic assumptions. These unavoidable approximate assumptions and the lack of balance equations can lead to some errors between the methods and actual situations [1–3].
In contrast, numerical methods are able to provide more accurate solutions to problems, which would not be solvable using LE techniques alone. The problems of rock slope stability involve complexities related to geometry, material anisotropy, nonlinear behavior, in situ stresses, and coupling processes (pore pressures, seismic loads, etc.) [4, 5]. Stead et al. [6] compared in detail the advantages and disadvantages of LE and numerical methods to rock slope analysis.
Many researchers have studied the inherent defects of LE methods and the calculation errors of the methods in practical application. For methods that satisfied all equilibrium conditions (e.g., Janbu, Morgenstern–Spencer, Spencer), Duncan [7] believed that the differences between calculated results would not exceed 12% and the error should be less than 6%, which was considered accurate. Liu et al. [8] contrasted the LE methods with two FE methods (strength reduction and overloading). The results show that the FE methods have a good consistency, and the safety factors of LE methods were slightly lower than those of the FE methods. Yu et al. [9] compared the results of the Bishop method with rigorous upper- and lower-bound solutions, resulting in that LE methods were reliable for homogeneous slopes. However, the methods tended to underestimate the stability of inhomogeneous slopes with a low slope angle. Huang et al. [10] discussed the influence of the lack of some simplified and static equilibrium conditions of LE methods on calculation results. Safety factors of elastic-plastic FE methods could be close to those of LE methods. Zheng et al. [11] pointed out that the premise was that Poisson’s ratio must be adjusted with the reduction of the cohesion and friction angle. Lane and Griffiths [12] claimed that traditional analysis methods could not include all boundary conditions, as well as fully reflect the heterogeneity of slopes and the complexity of forces. Hammouri et al. [13] suggested that critical slopes should be analyzed by both FE and LE methods.
Most of the existing studies are compared with the final calculation results of FE and LE methods, and only the comprehensive errors under various assumptions are obtained. However, fundamental causes of the generated error under each assumption and the impact of relevant factors have not been analyzed. Moreover, the mechanical mechanism and error range of LE methods are rarely involved.
In this study, we regard the basic assumptions that cause deviations between LE methods and actual situations as a starting point. By comparative analysis of FE and LE methods, the assumptions are separated to study the effect of each assumption on calculation results. Therefore, the error source, mechanism, and error range of LE methods are further explored.
2. Analysis and Evaluation Methods for Planar Landslides
Planar sliding, which has a universal characteristic of existing planar discontinuities, is very common in natural and engineered slopes [14]. The 2D case of such slopes should not be ignored, as many valuable lessons can be learned considering the mechanics of this simple failure mode. Planar failure is particularly useful for demonstrating the sensitivity of the slope to changes in shear strength and ground water conditions. These changes are not obvious when dealing with 3D slope failure [15].
2.1. Limit Equilibrium Analysis
For all shear failure slopes, the rock can be assumed to be the Mohr–Coulomb material, and its shear strength is represented by the cohesion and friction angle . If the effective normal stress on a sliding surface is , the shear stress developed on this surface is expressed as
The effective normal stress is the normal component of the vertical stress due to the weight of the sliding bed. Figure 1 shows a slope containing a continuous joint dipping out of the face and forming a sliding body. Calculation of the safety factor for the body involves resolution of the force acting on the sliding surface into two components that are perpendicular and parallel to this surface. If the dip of the sliding surface is , its area is , and the weight of the sliding body is , then the normal and shear stresses on the sliding surface areand Equation (1) can be expressed asor
In Equation (4), the expression () defines the resultant force acting down the sliding surface and is termed the sliding force (). The expression () defines the total shear strength forces acting up the surface and is termed the antislide force (). The safety factor of the slope can be quantified by the ratio of the antislide and sliding forces, which is expressed as
If the sliding surface is smooth and contains no infilling, the cohesion is zero and the above equation is simplified as Equation (6). When the slope is in a limit equilibrium state (), the safety factor :
2.2. A Finite-Element Evaluation Method
In numerical modelling of slopes, strength reduction and overloading techniques are commonly used methods for calculating safety factors [16, 17]. For comparison with the LE method and manifesting the actual stress states of sliding surfaces, a stability evaluation method is proposed to reflect the progressive failure characteristics of rock slopes based on FE calculation results. The method still defines safety factors as a ratio of antislide forces to sliding forces on a whole sliding surface and uses the concept of material safety reservation. In the calculation of safety factors, peak strength is used for unbroken materials and residual strength is selected for damaged materials [18]. For planar-type slopes, the antislide and sliding forces of a sliding surface are on a same line; thus, vector superposition can be performed. In summary, the global safety factor for the slope is calculated aswhere is the normal stress of the unbroken sliding surface element, is the normal stress of the damaged sliding surface element, is the shear stress of the unbroken sliding surface element, is the shear stress of the damaged sliding surface element, is the peak cohesion of the unbroken sliding surface element, is the friction angle of the unbroken sliding surface element, is the residual cohesion of the damaged sliding surface element, is the residual friction angle of the damaged sliding surface element, and is the length of sliding surface elements.
In general, the cohesion of slip-zone rocks decays relatively fast during the shear-sliding process of rock slopes, but the change of the friction angle is not very severe [19]. Therefore, residual shear strength is mainly provided by the friction angle. For this study, the parameters are assumed to be and in error analysis experiments. The selection of residual strength values in actual slope calculation should be based on laboratory tests.
3. Introduction of Interface Elements
The traditional joint model Goodman elements [20] is used to simulate faults, shear fracture zones, and discontinuities in rock mass. However, because the virtual normal stiffness and tangential stiffness of fractures are not easily determined, the practical application of the elements is greatly restricted [21]. Katona [22] proposed a friction-contact interface element without stiffness coefficients to simulate the sliding friction and opening and closing processes between two contact surfaces. However, the element uses a two-node simple element with constant contact forces, which makes it difficult to adapt to complex interface problems.
According to Katona’s work, Swoboda and Marence [23, 24] developed a new contact-friction interface element (COJO). The element can directly select the normal and tangential stresses on a contact surface as additional unknowns and simulate the complex geometry of the contact surface with a six-node triangular isoparametric element, which overcomes the shortcomings of traditional fracture elements. This study uses COJO elements to simulate the potential sliding surface of slopes. Geometric models are shown in Figure 2, where and are the normal and tangential stiffness of contact surfaces, and are the displacements of nodes, and and are the coupling forces of nodes.
(a)
(b)
(c)
Contact conditions of COJO elements are divided into three categories: fixed, sliding, and opening. Fixed means that the contact surfaces of rock mass are closed without relative displacements. Sliding indicates that the overall shear stress on contact surfaces exceeds their shear strength, and there is a relative displacement between the interfaces. Opening indicates that the tensile stress of an interface exceeds its tensile strength.
Frictional contact problems require repeated iterations to obtain correct solutions [25]. In calculation, first assume that elements are in a certain contact state and thus calculate the equivalent element stiffness matrix and load vector. After solving the FE equation, a set of test solutions can be obtained to check the contact state. If the state is the same as the original assumption, it is proven to be correct and the calculation is complete. If they are different, the test solution is selected as a new hypothesis state and the load vector is modified to perform a new iteration until convergence. Table 1 shows the contact states of COJO elements. In the table, and are the normal and tangential allowable stresses, where . and are the cohesion and friction coefficient of contact surfaces, the positive value of is specified as tensile stress, and is the initial displacement in the normal direction.
4. Numerical Tests for Error Analysis
Slope instability is a result of the interaction between the elastic-plastic deformation of sliding mass (or sliding beds) and the plastic sliding of sliding surfaces. In this study, we select the FE method to analyze slope stability under the basic assumptions of LE methods. Safety factors calculated by the two methods are used to calibrate the error caused by each assumption [26].
The planar sliding method has the following basic assumptions: (1) The sliding body is rigid. (2) Points of the sliding surface are destroyed at the same time. (3) The sliding surface follows the Mohr–Coulomb strength theory.
As shown in Figure 3, this research is carried out based on the following five aspects:(i)The FE method is used to simulate a slope under the above assumptions, and the systematic errors of a simulation platform are compared and verified.(ii)Removing assumption (2) makes a progressive failure of the slope, and the errors caused by plastic sliding of the sliding surface are compared and analyzed.(iii)We continue to remove assumption (1) and regard the sliding body as an elastic-plastic mass, and the errors caused by elastic-plastic deformation of the sliding body are comparatively analyzed. Normally, sliding surfaces are the weakest positions of planar landslides, and slopes are likely to be unstable before the sliding body has not entered a plastic stage. Therefore, we only discuss the effect of elastic deformation of the sliding body on calculation results.(iv)The errors caused by changes in sliding bed stiffness are further studied.(v)Based on the above research, the influence of stress release due to excavation on slope stability is analyzed.
4.1. Comparative Analysis Models
FE models are shown in Figure 4, which can be regarded as 2D plane-strain problems. Geometric parameters of the slope are shown in Table 2. Boundary conditions of the models are as follows: the slope surface is a free boundary, and the bottom and left sides are zero-displacement boundaries. The sliding body and sliding bed are regarded as elastic-plastic isotropic materials. The sliding surface is simulated by the interface elements, and the Mohr–Coulomb criterion is used to judge the friction states. There is no interface element between the excavation parts and the bed rock.
(a)
(b)
4.2. Systematic Error of the Simulation Platform
Although the LE method cannot reflect the actual states of slopes to a certain extent, the analytical solution under its basic assumption is an accurate value.
4.2.1. FE Simulation Approaches
(i)Sliding body is assumed to be rigid: in FE analysis, we select a large value of the elastic modulus and a small value of Poisson’s ratio for rocks; thus, they are approximated as rigid bodies. The sliding bed is reflected by the form of forces in LE methods; therefore, it is considered as a rigid support to simulate the assumption in finite elements. Linear-elastic analysis is used for the sliding body and sliding bed in solutions, which can avoid the plastic effect of c and φ.(ii)Points of the sliding surface are destroyed simultaneously: to make the points of the sliding surface reach the limit equilibrium (or instability) state at the same time, we set the friction angle of interface elements equal to the dip angle of the sliding surface and select the cohesion of zero. The interface elements are analyzed by the elastic-plastic model in FE solution. Calculation parameters of rock mass and the sliding surface are shown in Table 3.
4.2.2. Results Analysis
It can be seen from Figure 5 that the normal stress and shear stress on the sliding surface obtained by FE calculation increase linearly from top to bottom of the slope. The shear stress τ on the sliding surface is , which accords with the Mohr–Coulomb criterion and shows that the points reach the limit state at the same time. Thus, the setting that the friction angle of the sliding surface is equal to its dip angle (35°) and the cohesion is zero is correct.
(a)
(b)
As shown in Table 4, the safety factor calculated by the FE method is 0.99 and the error with the LE method is only 1%, which can prove the reliability of the numerical simulation platform.
4.3. Error Analysis of “Simultaneous Failure and Progressive Failure”
Removing the assumption of “simultaneous failure” for the planar sliding method, the influence of progressive failure on slope stability is studied. The plastic slip of the sliding surface is closely related to the cohesion and friction angle of interface elements, so the above problem is reflected by changing the value of c and φ of the sliding surface.
4.3.1. Analysis Schemes
(i)Assuming that the value of φ is equal to the dip angle of the sliding surface, the effect to calculation errors is studied by changing the value of c, as shown in Schemes 1–6(ii)When the value of φ is less than the dip angle of the sliding surface, calculation errors are affected by a combination of c and φ, as shown in Schemes 7–9(iii)When the value of c is a constant, the effect to calculation errors is explored by changing the value of φ, as shown in Schemes 10–12
Physical and mechanical parameters of the sliding body and sliding bed are the same as in Table 3, and parameters of the sliding surface and calculation schemes are shown in Table 5.
4.3.2. Results Analysis
When the sliding body and sliding bed are rigid, normal stress increases linearly along with the sliding surface from top to bottom of the slope. The magnitude of normal stress on the sliding surface is irrelevant to plastic parameters but only relates to the weight of the sliding body.
Figures 6(a)–6(c) show the shear stress on the sliding surface of Schemes 1–6, 7–9, and 10–12, respectively. The values of c and φ of the sliding surface directly affect the magnitude and distribution of the shear stress. According to the Mohr–Coulomb criterion, the sliding surface is in a steady state when shear stress is horizontal and linearly distributed. When the shear stress distribution appears in a partial oblique line, it indicates that the corresponding points of the oblique line reach limit instability states and local failure occurs in the slope. A complete oblique line illustrates that the whole sliding surface reaches the limit state and the slope is entirely unstable.
(a)
(b)
(c)
The shear stress distribution of Schemes 1–6 () can obviously reflect the progressive failure characteristics of the slope. As the cohesion decreases, sliding zones gradually expand from top to bottom of the sliding surface. Meanwhile, the slope presents a state from overall stability to local instability and eventually evolves into complete failure.
Table 6 shows the safety factors and errors of each scheme. Safety factors obtained by the FE method are less than those by the planar sliding method. In Schemes 5, 6, 9, and 12, slopes are steady and calculation errors are basically controlled within 10%. The reason of the errors is that the planar sliding method regards the sliding body as a particle without considering the inhomogeneity of forces. However, the errors of safety factors for local and overall failure slopes are relatively large. Due to the different roles played by the cohesion and friction angle, the calculation errors vary between 5% and 88%.
When the value of φ is equal to the dip angle of the sliding surface, the error range retains within 10% (Schemes 1–6 and Scheme 12). When the value of φ is less than the dip angle of the sliding surface, the influence of plastic parameters of the sliding surface on calculation errors increases significantly (Schemes 7–11).
In summary, the progressive failure effect of slopes on calculation errors between the two methods is that the ultimate failure modes of slopes depend on a combined strength of c and φ of the sliding surface. A basic reason is that safety factors solved by the FE method use peak and residual strengths to reflect different states of sliding surface elements. However, LE methods can only select a unified strength parameter, and the stable (or unstable) state of a whole sliding surface is consistent.
4.4. Error Analysis of “Rigid and Elastic Assumptions”
Removing the “rigid assumption” of the LE method, calculation results of the two methods based on progressive failure and elastic assumptions are analyzed. Elastic parameters of the sliding body are reflected by the elastic modulus E and Poisson’s ratio μ.
4.4.1. Analysis Schemes
The E and μ values of the sliding body are changed from rigid to general parameters of rock mass. At this time, the sliding body is gradually adjusted from rigid to elastic, and the sliding bed is still considered as a rigid support.
The bulk density of rock mass is . Parameters of the sliding surface for stable slopes are and ; for local failure slopes, they are and . Calculation schemes and results are shown in Tables 7 and 8, respectively.
4.4.2. Influence of the Elastic Modulus on Slope Stability
Figure 7 shows the influence of the sliding body elastic modulus on the stress distribution of the sliding surface. The distributions of the normal and shear stresses are roughly consistent. When the elastic modulus of the sliding body is large (close to the sliding bed), stress distribution behaves linear and the stress value of various parts of the sliding surface has little difference. As the value of E continues to decrease, stress distribution changes from linear to “valleys” and the inhomogeneity gradually increases. Besides, the maximum normal stress and shear stress are basically increasing. When the elastic modulus of the sliding body is equal to or less than that of actual rock mass, it has almost no effect on stress distribution.
(a)
(b)
The effect of the sliding body elastic modulus on safety factors and errors is shown in Table 7. As the value of E decreases, the error of two methods decreases first and then increases. For stable slopes, the error reduces from 15% to 1% and then increases to 7%. For local failure slopes, the error drops from 49% to 2% and then increases to 17%. As the elastic modulus gradually approaches a general value of rocks and continues to soften, the error eventually drops to 4%.
The main reason for such a change is that when the sliding body begins to transform from rigid to elastic, its slight deformation properly relieves the plastic slip of each point on the sliding surface. However, as the sliding body continues to soften, the increasing deformation exacerbates the degree of plastic slip, and the expansion of damage zones in the slope increases the calculation error. When the elastic modulus of the sliding body becomes much smaller than that of the sliding bed, the plastic slip of the sliding surface is restrained.
4.4.3. Influence of Poisson’s Ratio on Slope Stability
As mentioned above, stress distribution properties of the sliding surface are determined by the elastic modulus of the sliding mass. However, when the elastic modulus is a constant, the position and magnitude of the maximum shear stress on the sliding surface are decided by Poisson’s ratio. With the increase of Poisson’s ratio, the shear stress on the upper part (slope crest) of the sliding surface decreases and on the lower part (slope toe) increases (Figure 8). The distribution of the maximum shear stress tends to shift toward the crest.
(a)
(b)
For stable slopes, the shear stress on the sliding surface shows a regular and equal proportion change. Failure zones of local instability slopes appear in the middle and upper parts of the sliding surface, and damage scopes gradually enlarge with the increase of Poisson’s ratio.
Table 8 shows the influence of the sliding body Poisson’s ratio on safety factors and errors. The FE results of all the schemes are less than the results of the planar sliding method. With the increase of μ, the error of stable slopes changes little from 10% to 6%. The change in local failure slopes is not monotonic but first reduces from 19% to 10% and then rapidly increases to 31%. The reason for the phenomenon is the same as that of the elastic modulus, which is a result of the interaction between elastic deformation of the sliding body and elastic-plastic deformation of the sliding surface. Relatively, the elastic modulus exhibits more sensitivity than Poisson’s ratio in FE analysis; hence, the elastic modulus has a greater influence on coupling action.
4.5. Influence of Sliding Bed
4.5.1. Analysis Schemes
According to Section 4.4, the effect of coordinate deformation between the sliding bed and sliding body on calculation results is further studied. The problem can be reflected by changing the elastic modulus and Poisson’s ratio of the sliding bed so that it transforms from a rigid support to an elastic support. Parameters of the sliding body and sliding surface are shown in Table 8. Sliding bed parameters and schemes are shown in Table 9.
4.5.2. Results Analysis
Table 10 shows the effect of the sliding bed stiffness on safety factors and errors. For both stable and local failure slopes, calculation errors are gradually increasing and safety factors of the FE method are still less than those of the LE method. Moreover, the change in the sliding bed stiffness has great influence on the error. When the sliding bed is a rigid or an approximate rigid support, the error has little change between 7% and 11%. However, when the elastic parameter of the sliding bed is gradually close to that of the sliding body, the safety factor calculated by the FE is greatly reduced and the error rapidly increases to over 50%. As the sliding bed is further softened, the error continues to increase, but the growth rate slows down significantly.
4.6. Influence of Excavation Stress Release
In this section, we study the effect of excavation stress release on calculation results when the geometric parameters (slope shape and dip angle of sliding surfaces) of slopes are determined.
4.6.1. Analysis Schemes
Most of the excavated slopes have no large tectonic stress; thus, an initial stress field is dominated by the gravity stress. A main factor that influences the magnitude of the gravity stress is the bulk density γ of rock mass. The larger value of γ indicates that the more tectonic stress the rock mass gathered was before excavation. In this study, the initial stress changes are simulated by altering the bulk density of the sliding body and sliding bed. Calculation errors are explored under the cases considering and without considering in situ stress release, and the relevant schemes are shown in Table 11.
4.6.2. Results Analysis
As shown in Figure 9, after excavation and in situ stress release, the normal stress and shear stress on the sliding surface tend to transfer to the bottom of the slope, especially the shear stress changes more significantly. Under the same bulk density, stresses on the sliding surface obtained by considering excavation schemes are larger than those without considering. In the case of the largest amount of in situ stress release (Scheme 1), local instability occurs in the slope.
(a)
(b)
Stability analysis of actual slopes should consider the stress release caused by excavation, and the resulting error is difficult to be neglected. Table 12 shows the safety factors and error range for slopes with and without considering in situ stress release. It is obvious that the in situ stress release is very unfavorable to slope stability. Calculation errors of the slopes without considering stress release are between 4% and 13%. However, safety factors reduce rapidly with the increase of in situ stress release, and the calculation error reaches more than 30%. In particular, the error of local damage slopes is as high as 54%.
It should be noted that the comparative analysis in this section is conducted after slope excavation. This is different from the concept of “cutting slope and reducing load” to achieve stabilization because the latter is to study the impact of excavation on slope stability, that is, to compare the slope states before and after excavation.
5. Discussion
This study only discusses the error source and error range of LE methods when the sliding surface is planar. Further research is needed for biplanar and circular sliding slopes. The following are some preliminary understandings of these two types of landslides:(i)The stability of a biplanar landslide depends not only on the occurrence and mechanical index of sliding surfaces but also on the interaction of the two sliding surfaces. For example, retrogressive and thrust-type landslides are different in sliding forms and failure mechanisms. Safety factors of the biplanar landslide cannot be simply calculated by scalar quantity of the upper and lower sliding surface results, and an overall stability of the slope needs to be evaluated.(ii)In numerical modelling, we need to determine the critical sliding surface of a circular sliding slope by LE methods. Combined with the FE evaluation method proposed in this study, the safety factor is defined as the ratio of the total antislide moment and sliding moment on the surface, and it is solved by means of moment balance.
6. Conclusions
(i)For a planar landslide, safety factors obtained by the FE method without considering residual strength are close to those obtained by the LE method, and the error is within 10%. The FE evaluation method considering residual strength proposed in this study is more practical. For slope stability, the results are generally more dangerous than the results of the LE method.(ii)The influence of progressive failure makes the calculation error of stable slopes within 10%; the error of local failure slopes is between 6% and 35%, and the error of overall failure is above 38%. Error values are closely related to the damage range of sliding surfaces and the selection of strength parameters. Because LE methods assume that the strength of sliding surfaces is reduced simultaneously, an unrealistic stress distribution is obtained.(iii)As the elastic deformation of the sliding body increases, the calculation error of stable slopes is between 5% and 15%. For local failure slopes, the error is first reduced and then increased. The reason is that the energy of the sliding body can be firstly sustained by its coordinated deformation, thereby alleviating the plastic slip of the sliding surface, but the continuous softening of the sliding body further aggravates the sliding state. When the sliding body parameters are close to those of general rock mass, the error grows to 10%–30%. The elastic modulus of the sliding body has a greater effect on calculated results than Poisson’s ratio.(iv)The sliding bed rigidity has the same influence on stable and local failure slopes. When the sliding bed is a rigid support, the error is between 7% and 10%. When elastic parameters of the sliding bed are close to those of the sliding body, the error rapidly grows to 50% and beyond. Further softening of the sliding bed continues to increase the error, but the growth rate slows down.(v)In situ stress release has a significant impact on the stability of excavated slopes, and calculation errors are proportional to the magnitude of the release. For a rock slope with a height of 80 m under self-weight stress, the release of in situ stress caused by stepped excavation can reduce safety factors by more than 30%.
Data Availability
The finite-element program used to support the findings of this study was supplied by Prof. Ning Li under license and so cannot be made freely available. Requests for access to these data should be made to the corresponding author.
Conflicts of Interest
The authors declare no conflicts of interest.
Acknowledgments
This research was funded by the National Natural Science Foundation of China (Nos. 11572246, 51779207, and 41372295). Besides, we would like to thank Dr. Gao Lv for his guidance on research ideas and Miss Qi’s help in English writing.