On the Computational Precision of Finite Element Algorithms in Slope Stability Problems

Although the finite element method (FEM) has been used extensively to analyse the slope stability problems, the computational precision and definition of failure are still two main key concepts of finite element algorithms that attract the attention of researchers. In this paper, the modified Euler algorithm and the explicit modified Euler algorithm with stress corrections are used to analyse two dimensional (2D) slope stability problems with the associated flow rule, based on the shear strength reduction method. +e rounded hyperbolic Mohr-Coulomb (M-C) yield surface is applied. Effects of the element type and various definitions of failure on the computational precision of 2D slope stability problems are evaluated. Conclusions can be drawn that the modified Euler scheme is applicable when the factor of safety (FOS) is small; however, the explicit modified Euler algorithm with stress corrections is more precise if the factor of safety is relatively large. +e fully integrated quadrilateral isoparametric element is better than the triangular element in terms of the precision. With respect to the definition of failure, the displacement mutation of the characteristic point combining with the continuums of the plastic zone can be regarded as a reliable definition of failure and can be widely used to perform and analyse numerical simulations of slope stability problems.


Introduction
e traditional limit equilibrium method [1][2][3][4] and limit analysis method [5][6][7][8] have been widely used in analysing a majority of slope stability problems over the course of previous decades of research. It should be noted that a closeform solution of elastic-plastic slope problems is only possible for basic cases where the loading and geometry are simple. e finite element (FE) approach has a number of merits, e.g., assumptions about the failure shape and location being no need to be previously determined, assumptions about the slice side forces being not necessary, information of deformations having been given at working stress levels, and the progressive failure process being able to be monitored. ese advantages make the FE method attractive in the realm of slope stability problems over the utility of limit equilibrium method or limit analysis method. It is noted that the FE method includes FE strength reduction approach [9] and FE limit analysis approach [10][11][12], and only the former approach is discussed herein. e FE approach in analyzing slope stability problems can be catalogued to the realm of elastic-plastic mechanics. Marcal and King [13] and Yamada et al. [14] developed a method to deal with the continuum elastic-plastic problems through the FE approach. e method was based on a plastic stress-strain matrix facilitating the incremental treatment of elastic-plastic problems. eir study demonstrated that the FE method is convincing and powerful. In particularly, the assumption of the perfectly plastic material or nonhardening plastic-rigid body, which was indispensable for limit analysis, was no longer a must for the FE method [15]. e FE approach in terms of incremental plasticity has been widely applied to the slope stability problems then [16][17][18][19][20][21][22][23][24][25]. Zienkiewicz et al. [9] proposed a shear strength reduction technique by which the original shear strength parameters must be divided into a number (i.e., the FOS), in order to bring the slope to the state of failure. Since then, the application of FE method in slope stability problems enjoyed a fruitful outcome. However, the computational precision and definition of failure probably are two of the main concerns that still need continuous research. e computational efficiency of the FE approach is typically affected by the integration algorithm, density of the mesh, shape of the element, and so forth. A number of studies have been conducted to improve the precision of the FE method [26][27][28], among which Sloan [10] proposed an efficiency substepping scheme (i.e., the modified Euler scheme) for integrating elastoplastic stress-strain relations. His methods were applicable to a general type of constitutive law, and the error was controlled in the integration process of elastoplastic constitutive laws by selecting the size of each substep automatically over each time interval. Abbo [29] improved the scheme proposed by Sloan by giving a method of stress corrections. eir algorithms are particularly suitable for analyzing typical boundary value problems in geotechnical engineering. e numerical simulations of slope stability problems have been performed by using the triangle element [20,23] and the isoparametric element [16][17][18][19]. ey identified that a dense mesh can result in an increase in precision, however with a sacrifice of computing time. Hence, a balance between the number of mesh and the computational precision should be considered. Different definitions of failure have been widely studied; however, a systematic study of comparing these criteria in slope stability problems has not been fully performed.
In this paper, both the modified Euler algorithm and the explicit modified Euler algorithm with stress corrections proposed by Sloan [10] and Abbo [29], respectively, are applied to integrate the elastic-plastic stress-strain relationship. e rounded hyperbolic M-C yield surface is used to smooth the vertices, thus eliminating the computational difficulties. Case studies of two types of homogeneous slopes (2D plane strain condition) are performed in an FE platform to analyse the computational precision in terms of different integration algorithms, shape of the element, and effectiveness of the three definitions of failure. e equivalent strain nephograms will be presented with respect to different magnitudes of FOS.

Finite Element Method for Slope Stability Analysis
Following Sloan [10] and Abbo [29], the explicit modified Euler, which is a family of explicit methods, is used in this study. is method is associated with the shear strength reduction scheme to present a systematic analysis on three definitions of failure in slope stability problems.

Numerical Integration Scheme.
e explicit modified Euler integration scheme requires determination of the intersection with the yield surface when the stresses experience a transition from an elastic state to plastic state (e.g., [10,29,30]). e aim of this approach is to compute the stress-strain response over each substep by integrating the elastic-plastic constitutive matrix D ep . In order to determine the portion of the stress increment that lies within the yield surface, a scaler α must be found. After that, the modified Euler scheme is accurate for very small time steps, and thus smaller substeps are required by subdividing ΔT (0 < ΔT < 1). e error is controlled in the integration process of elastoplastic constitutive laws by selecting the size of each substep automatically over each time interval. is error control can be achieved by using a local error measure. Obviously, the size of each subincrement may vary throughout the integration process instead of assuming substeps to an empirical standard and of the same size. e formulation and numerical implementation of the stressstrain relationship in the incremental form is as follows: where ∆σ e denotes a vector of elastic stresses; D denotes the stress-strain matrix; b � zg/zσ and g is the plastic potential; and ∆λ is displayed as where f represents the yield surface and ∆ε is the strain rate. Following Sloan [10], the nonlinear equation in the light of variable α is solved by the secant and Newton-Raphson method for its quick convergence. Abbo [29] argued that the drawback of this algorithm is that it may diverge in some circumstances as it does not constrain the solution. Hence, in his study, the modified regula-falsi procedure was used. As argued by Potts and Gens [31], who found that it is necessary to apply some forms of stress corrections because a cumulative effect does not satisfy the yield condition. Abbo [29] has given a stress correction method, and details of this method can be found in his publication. Both the two integration algorithms are used to perform the stress-strain relations in the present study. e flow chart of the integration algorithm is shown in Figure 1. e M-C yield criterion is applied. e rounded hyperbolic M-C method is used to solve the computational difficulties due to the gradient discontinuities which occur at the tip, and details can be found in the Appendix. As an associated flow rule is used, the plastic potential follows the same form as that of the M-C yield criterion with the dilation angle ѱ substituting the friction angle ϕ.

Shear Strength Reduction Technique.
e shear strength reduction technique is to define a number, which is normally named as the factor of safety (FOS). e original shear strength parameters, in terms of the cohesion c and the friction angle ϕ, are divided to a number of FOS, in order to bring the slope to the point of failure [17]. Hence, the technique is displayed in terms of the following equations: where F s is the FOS and c′ and ϕ′ are the reduced cohesion and friction angle, respectively. In addition, the adjustment of the Young's modulus E and Poisson's ratio υ is adopted to allow the shear strength reduction scheme to be realized in the finite element simulations [25,32,33].

Definition of Failure.
ree definitions of failure in the slope stability problems consist of (1) nonconvergence of the solution (F1); (2) continuums of the plastic zone (F2); and (3) mutation displacement of the characteristic point (F3). e nonconvergence option indicates that within a user-specified maximum number of iterations, no stress distribution can be found that is simultaneously able to satisfy both the failure criterion and global equilibrium [17,34]. e second criterion describes that the plastic zone is continuous throughout from the toe of the slope to the top of the slope [18,19], while the third indicator is the mutation displacement of the characteristic point [17,20,21,23,35]. All of the three definitions of failure will be analysed throughout this paper.
Stress correction:

Case Study
Two examples of typical FE soil slope models are established and analysed, and validation against the literature data is given where possible. Two examples of different foundation layers have been selected since Griffiths and Lane [17] argued that the foundation layer has a significant effect on the FOS. Effects of the element type and various definitions of failure on the computational efficiency of the slope stability problems are systematically studied.
3.1. Case Study 1. Following Dawson et al. [16], the slope is assumed to be homogeneous with a slope height of 10 m and slope angle of 45°( Figure 2). e foundation layer D � 3 m. e material of the soil is discretized with a quadrilateral isoparametric plane strain element (except the case of studying the effects of the element type), with a total number of 330 elements and 369 nodes. Both the left-hand and righthand side boundaries allow for vertical movement. e condition on the bottom boundary is fixed in both vertical and horizontal directions. e soil model consists of six parameters, as shown in Table 1. e associated flow rule is used. e moving direction of the characteristic point (i.e., node a in Figure 2) opposite to the positive x-axis is taken as positive.

Comparison of the Accuracy of the Definition of Failure.
e explicit modified Euler algorithm with stress corrections proposed by Abbo [29] is used to perform the integration of stress-strain relations. e convergence fails until the case of F s � 1.12. e equivalent strain nephograms with different values of F s are shown in Figure 3. Obviously, the case of F s � 1.12 indicates failure of the slope in terms of the F1 definition of failure. e plastic zone is initially identified at the toe of the slope and then develops from the toe to the top of the slope, with the increase in the magnitude of F s . When F s � 1.05, the plastic zone is continuous from the toe of the slope to the top of the slope. Hence, the case of F s � 1.05 indicates failure if the F2 definition of failure is used. e horizontal displacement of the representative node a, which is shown in Figure 2, is plotted against the FOS ( Figure 4).
As shown in Figure 4, the displacement experiences a few increases before approaching F s � 1.05. A sharp deduction of the displacement can be found beyond F s � 1.05. e mutation displacement of the characteristic point (i.e., node a) is observed at the case of F s � 1.05, which is an indicator of failure if the F3 criterion is used. is result is similar to that concluded by Dawson et al. [16], who demonstrated that F s is equal to 1.03 when failure is identified. F1 definition of failure gives a larger magnitude of F s . e safety factors of the slope predicted by F2 and F3 definitions of failure are more conservative and similar to the literature data. It is suggested that the larger value of either F2 or F3 should be selected as the FOS.

Effects of the Integration Scheme on the Computational
Precision. In order to analyse the effects of integration method on the computational precision, the modified Euler algorithm without stress corrections proposed by Sloan [10] is used to perform the integration of stress-strain relations.
e results, as presented in Figure 5, are compared with those by Abbo [29] as shown in Figure 3. is time, the convergence fails before the case of F s � 1.1. e development of the plastic zone is reasonable before F s � 1.1. However, the plastic zone is random when F s � 1.1, which is against the common sense. e reason may lay in the fact that no stress correction method is incorporated in the algorithm proposed by Sloan, and the accumulation of the global error may lead to the incorrect plastic strain. In addition, this error may result in nonconvergence when F s is even relatively very small. Hence, the integration method with stress corrections is more applicable, especially when the magnitude of F s is relatively large. When F s � 1.06, the plastic zone is continuous from the toe of the slope to the top of the slope. Hence, the case of F s � 1.06 indicates failure if the F2 definition of failure is used. In Figure 6, the mutation displacement of the characteristic point (i.e., node a) is observed when F s � 1.09, which is an indicator of failure if the F3 criterion is used. It is obvious that if the modified Euler algorithm is used, combining the failure criterions F2 and F3, the safety factor of this slope is 1.09.
Comparing the two integration algorithms, the explicit modified Euler integration with stress corrections yields a more conservative result, and its estimation of FOS is similar to the literature data.

Effects of the Element Type on the Computational
Decision. To investigate the effects of the element type, the material of the soil is discretized with a triangular plane strain element, with a total number of 1920 elements and 1027 nodes. Parametric studies have been performed on the number of the elements. It is found that if the equal numbers of the triangular elements are used as those of the quadrilateral isoparametric elements, the results are not correct when validated with the literature data. As aforementioned, the modified Euler scheme with stress corrections is used to perform the integration. e convergence fails until the magnitude of F s is equals 1.12. Hence, the indicator of failure is that F s � 1.12 with the utilization of F1 criterion. As shown in Figure 7 Figure 8. e horizontal displacement increases rapidly after the characteristic value of F s � 1.05. e conclusions are coincident with those obtained by using the quadrilateral isoparametric plane strain elements. However, the total number of elements when using the triangular elements is much larger than those when using the quadrilateral     Mathematical Problems in Engineering isoparametric elements. As a result, the computation time is larger when comparing the triangular elements to their quadrilateral isoparametric counterparts. e velocity fields obtained from case study 1 are illustrated in Figure 9, in terms of both quadrilateral isoparametric elements (Figures 9(a) and 9(b)) and triangular elements (Figures 9(c) and 9(d)). e directions of arrows represent the flow of velocity while the length of arrows represents the magnitude of displacement. It can be found in Figure 9 that the velocity patterns are coincident with the theoretical pattern predicted by the slip line method. In addition, the velocity zone indicated by the F1 criterion is larger and wider than that indicated by the F3 criterion. It can be expected that the failure zone is wider when the F1 criterion is used. Hence, the F3 criterion is more conservative, which is consistent with the conclusion drawn from the equivalent strain nephograms and horizontal displacement plots.

Case Study 2.
e soil slope model is assumed to be homogeneous. e dimensions are just the same as those studied by Su and Li [35], with a slope height of 10 m and a slope angle of 26.57°, as shown in Figure 10. e foundation layer D � 10 m. e boundary conditions and the flow rule are the same as those used in case study 1. e material of the soil is discretized with a quadrilateral isoparametric plane strain element (except the case of studying the effects of the element type), with a total number of 1000 elements and 1081 nodes. e material properties are present in Table 2.
e positive moving direction of the characteristic point (i.e., node c in Figure 10) coincident with positive x-axis is taken as positive.

Comparison of the Accuracy of the Definition of Failure.
e modified Euler scheme with stress corrections is used to perform the integration. Convergence cannot be achieved after F s reaches 1.8. Likewise, the plastic zone commences on the toe of the slope and then develops through the toe to the top of the slope. e final transfixion of the plastic zone from the toe to the top of the slope can be found at F s � 1.7 (Figure 11(c)). As illustrated in Figure 12,    the indicator of failure is found at the case of F s � 1.67 if the F3 criterion is used to define the failure, which agrees well with the conclusion of F s � 1.667 proposed by Su and Li [35]. Similar conclusions can be drawn that the combining the F2 and F3 definitions of failure, the estimation of FOS is more conservative than that of the F1 definition of failure.

Effects of the Integration Scheme on the Computational
Precision. e equivalent strain nephograms obtained by using the modified Euler scheme without stress corrections are plotted in Figure 13. When compared to Figure 11, the results are the same by using both of the two integration schemes. However, when the value of F s is larger than 1.5,   the results are random and the plastic belt is very difficult to be achieved. e convergence problem cannot be well controlled, and it takes much longer computational time. Hence, the modified Euler scheme with stress corrections is more applicable, in particular for the case with a higher value of F s .

Effects of the Element Type on the Computational
Decision.
e triangular plane strain elements are used to be compared with the quadrilateral isoparametric elements, to investigate the effects of the element type on the computation decision in slope stability problems. e total number of elements is 2000. Likewise, the explicit modified Euler scheme with stress corrections is applied. Convergence fails until F s approaches 1.8. As illustrated in Figure 14, the plastic zone is continuous throughout the toe of the slope to the top of the slope when F s � 1.67 ( Figure 14). e mutation displacement of the characteristic point is found at the case of F s � 1.67 ( Figure 15).
ough consistent conclusions can be obtained by using the triangular elements, more computational elements are required compared with the quadrilateral isoparametric elements. e accuracy of the computation cannot be well guaranteed if the mesh is coarse.

Discussion
e main objective of the present study is to evaluate the effects of the failure criterion on the slope stability, which have not been fully revealed theoretically in the literature. Apart from the above two cases, two more cases involving the conditions of heterogeneity and irregular geometrical shapes have been studied as well. Similar conclusions can be drawn that if F1 is used as the definition of failure, the calculated safety factors of the slope depend on the integration algorithms that are used. It means that different integration algorithms result in different outcomes. However, if we combine the F2 and F3 definitions of failure, the larger value obtained by these two failure criterions is selected as the FOS, and the value of FOS is rarely affected by the integration algorithms.
In addition, the effects of the parameters on the slope stability, including friction angle, cohesion, unit weight, and Poisson's ratio, are studied. e results are shown in Figures 16 and 17.
e parametric studies show that the slope stability increases with the increase in the friction angle or cohesion and, however, decreases with the increase in the unit weight, as expected. Poisson's ratio exhibits negligible effects on the FOS.         Mathematical Problems in Engineering 13

Concluding Remarks
In the framework of elastic-plastic mechanics, the shear strength reduction technique was utilized to analyse the slope stability. e effects of the finite element calculation conditions (e.g., the definition of failure, the element type, and the numerical integration algorithm) on the computational efficiency of 2D slope stability were studied. Conclusions were drawn as follows: (1) Because of the accumulative error of the computational stress, the modified Euler scheme without stress corrections was applicable only when FOS is relatively small, and it is not applicable to predict the slip surface of the slope for some cases (e.g., case 2); however, the explicit modified Euler algorithm with stress corrections was more precise even if the factor of safety was relatively large. (2) e displacement mutation of the characteristic point combined with the continuums of the plastic zone can be regarded as the most reliable definition of failure and can be widely used to perform and analyse numerical simulations of slope stability problems. (3) Compared with the fully integrated quadrilateral isoparametric element, the triangular element can yield the same FOS but with more elements.

A Rounded Hyperbolic M-C Yield Function
For the sake of implementation of geotechnical constitutive laws into finite analysis, many technical problems must be taken into consideration. Great efforts must be made for parametric control in finite element analysis allowing the newly proposed theory to run successfully in finite element codes, among which the computational difficulties due to the gradient discontinuities which occur at the tip or vertex of the yield curve are the most important. In the present study, the rounded hyperbolic M-C yield function is used to eliminate the computational difficulties [10,29]. e M-C yield surface in the three dimensional stress space is a hexagonal yield surface pyramid. e rounded hyperbolic M-C approach is to use a hyperbolic approximation in the meridional plane to eliminate the tip singularity and a trigonometric rounding in the octahedral plane to eliminate the edge singularities.
As shown in Figure 18(a), the expression of the straight line can be determined directly by the determination of the conventional M-C yield criterion and the slope of the straight line is given by sin ϕ. e straight line intercepts the p-axis at p � − c cot ϕ. Following Sloan [10] and Abbo [29], the hyperbolic approximation to the Mohr-Coulomb criterion is utilized to remove the apex singularity.
As shown in Figure 19(a), in the octahedral plane, there are six vertices. A rounded surface is used to smooth these vertices, and the complete yield surface is displayed as where A and B are parameters defined in Abbo [29].

Data Availability
No data were used to support this study.

Conflicts of Interest
e authors declare that they have no conflicts of interest.