Location of Tension Cracks at Slope Crests in Stability Analysis of Slopes

Abstract

Over the conventional limit equilibrium method and limit analysis method, the finite element method is advantageous, especially for slopes involving complex failure mechanisms where the critical slip surfaces cannot be represented by log spirals and other similarities. In the presence of tension cracks at slope crests, however, the finite element method encounters difficulties in convergence while handling Mohr–Coulomb’s yielding surfaces with tensile strength cut-off. Meanwhile, the commonly used load-controlled method for the system of nonlinear equilibrium equations is hard to bring the slope into the limit equilibrium state. The two drawbacks drag down the finite element method in more extensive applications. By reducing the constitutive integration of plasticity with non-smooth yielding surfaces to the mixed complementarity problem, the convergence in numerical constitutive integration is established for arbitrarily large incremental strains. In order to bring the slope to the limit equilibrium state, a new displacement-controlled algorithm is designed for the system of nonlinear equilibrium equations, which is far more efficient than the load-controlled method. A procedure is proposed to locate tension cracks. Corresponding to the Mohr–Coulomb failure criterion with and without tensile strength cut-off, the failure mechanisms differ significantly, while the difference in the factor of safety might be ignorable.

1. Introduction

Over the conventional limit equilibrium method and the limit analysis, the finite element method (FEM) is advantageous, especially for complex slopes where the critical slip surfaces cannot be approximated by circles, log spirals, and other similarities. The failure mechanism is automatically identified as long as the slope is brought into the limit equilibrium state from its natural state. Therefore, FEM and other similarities together with the strength reduction technique have been advocated in the past decades by many researchers, such as Matsui and San [1], Ugai and Leshchinsky [2], Griffiths and Lane [3], Dawson et al. [4] and Yang et al. [5]. Based on field surveys and monitoring data, Xuan Zhang et al. analyzed the deformation characteristics and disaster prediction model of a granite rockslide with a discrete element method [6]. Somaie Jolfaei and Ali Lakirouhani analyzed the sensitivity of effective parameters in borehole failure in rocky environments with neural networks [7].

Tension cracks are ubiquitous at the slope crests. Ignoring the presence of these cracks in the analysis of slope stability may reportedly overestimate the stability factor up to 70% [8]. Duncan and Wright [9] have suggested that if tensile stresses are encountered in the process of slope analysis, either a tension crack should be set or the failure envelope should be modified to eliminate tension.

In the conventional limit equilibrium method, the location and depth of the tension crack are determined by a test-and-error process, see Spencer [10], where many artificial assumptions are introduced, however. With the rapid development of computer software and hardware, various effective optimization methods have been adopted to locate the critical failure surface, such as the improved wolf swarm algorithm [11] and harmony search algorithm [12].

In the upper bound method of the limit analysis, the tension crack is a portion of the critical slip surface, which is usually represented by a log-spiral and determined by minimizing the critical slope height. Tensile failure and shear failure are unified in one yielding surface, and the most notable is the Mohr–Coulomb yielding surface with the tensile strength cut-off. The upper bound method is limited by the simple failure mechanism. Since Michalowski [13], the upper bound method appears to be a research hot point in the stability analysis of slopes with tension cracks. Utili [14] proposed a set of solutions for the stability of homogeneous slopes with open cracks by the kinematic method of limit analysis. In addition to a pre-existing crack, a crack formation is also considered as part of the failure mechanism in Michalowski [8] with reduced or no tensile strength. Such a method is extended to three-dimensional slopes [15] and infinite slopes against shallow slides [16]. However, the transformation of the yielding surface from the principal stress space to the Mohr stress plane applied in the kinematic method is somewhat intricate.

Using FEM, in principle, the tension crack should be captured as long as the tension portion is cut off the Mohr–Coulomb yielding surface and the slope is brought into the limit equilibrium state. However, this is not the case. Thus far, the slope stability analysis based on the finite element strength reduction has been limited to the yielding surfaces without tensile strength cut-off. The reason lies in the singularity created by tension cut-off even if the yielding surface itself is smooth, such as the Drucker–Prager yielding surface.

As is well known, the Mohr–Coulomb yielding surface in the principal stress space is composed of six facets with six edges and one apex, as shown in Figure 1a. At a point on the edge of the yielding surface, the normal is not unique. Huge efforts have been paid to the constitutive integration of plasticity with non-smooth yielding surfaces for either implicit algorithms such as Clausen et al. [17] or explicit ones such as Sloan et al. [18]. Tensile strength cut-off creates additional edges and vertices near the head of the yielding surface, as shown in Figure 1b, which further exacerbates the burden of constitutive integration. Since Zienkiewicz and Pande [19], followed by Menetrey and Willam [20], as a result, attempts to smooth the Mohr–Coulomb yielding surfaces have never ceased, see Abbo et al. [21] and Andy et al. [22]. Nevertheless, smoothing compromises precision in order to conserve the convexity of the yielding surface. Further, smoothing leads to the sacrifice of computational efficiency because the incremental strain has to be small enough to secure the convergence of constitutive integration when the stress point is close to the portion of the big curvature of the smoothed yielding surface.

This study directly deals with the Mohr–Coulomb yielding surface and its cut-off version instead of their smoothed versions, and thus no error is introduced. Meanwhile, a displacement-controlled method (DCM) is designed, which, compared to the load-controlled method (LCM) familiar to all of us, is easy to bring the slope into the limit equilibrium state. DCM has a very similar program structure to LCM. One of the outcomes from DCM is the load multiplier representing the load level the slope can undertake. Then a secant method is designed for the factor of safety according to the relation between the load multiplier and the strength reduction factor. Finally, a procedure is proposed for the location and depth of tension cracks.

Typical slopes are analyzed using the Mohr–Coulomb yielding surfaces without and with tensile strength cut-off, represented by intact M–C surfaces and cut-off M–C surfaces subsequently, respectively. Comparisons are made of the results from the two kinds of yielding surfaces, suggesting that significant differences exist in the failure mechanism while the factor of safety might differ little. Using the intact M–C surface will result in a notorious abnormal phenomenon, i.e., the tensile failure zone extends to the model boundary, severely overestimated. Using the cut-off M–C surface, the abnormal is eliminated. A procedure is proposed for the location and depth of tension cracks.

2. The Mixed Complementarity Problem from Constitutive Integration of Plasticity

In this section, we recapitulate the implicit constitutive integration of non-associative plasticity with a non-smooth yielding surface. Then, the mathematical model, represented by a mixed complementarity problem, is established for the numerical constitutive integration of plasticity. In order to avoid unnecessary distraction, the material is assumed to be elastic and perfectly plastic, which is also the most applied material model in geotechnical engineering. Extension to hardening materials is straight forward. Here, we only touch upon the basic formulation of finite-dimensional variational inequalities, which can be referred to in the monograph [23] for details. The mixed complementarity problem (MiCP) to be built in this study is a special case of finite-dimensional variational inequalities.

2.1. Elasto-Perfectly Plasticity in the Rate Form

As is well known, the 6-dimensional stress rate vector $\dot{\sigma }$ is related to the 6- dimensional strain rate vector $\dot{\epsilon }$ by the constitutive relationship

$$\dot{\sigma }=D\left(\dot{\epsilon }-{\dot{\epsilon }}_p\right)$$

(1)

Here, $D$ is the 6 × 6 elasticity matrix, which is symmetric and positive; see Zienkiewicz and Taylor [23]. The strain rate vector $\dot{\epsilon }$ is provided and has the six components of ${\dot{\epsilon }}_x$,${\dot{\epsilon }}_y$, ${\dot{\epsilon }}_z$, ${\dot{\gamma }}_{yz}$, ${\dot{\gamma }}_{zx}$, and ${\dot{\gamma }}_{xy}$. The 6-dimensional stress rate vector $\dot{\sigma }$ is unknown and has the same subscripts as $\dot{\epsilon }$. The 6-dimensional plastic strain rate vector ${\dot{\epsilon }}_p$ is unknown, having the same subscripts as $\dot{\epsilon }$ but following the flow rule [20]

$${\dot{\epsilon }}_p={\displaystyle \sum _{i=1}^m{\dot{\lambda }}_i\nabla g_i}$$

(2)

Here, the m scalars, ${\dot{\lambda }}_1$,$\cdots$,${\dot{\lambda }}_m$, are nonnegative plastic multiplier rates.

Except $\dot{\epsilon }$ and $D$, all the quantities in Equations (1) and (2) and the subsequent discussions are dependent nonlinearly on the current stress point $\sigma$ and plastic multiplier vector $\lambda$, written as $\dot{\sigma }\left(\sigma ,\lambda \right)$, ${\dot{\epsilon }}_p\left(\sigma ,\lambda \right)$, ${\dot{\lambda }}_i\left(\sigma ,\lambda \right)$ and $g_i\left(\sigma ,\lambda \right)$. In order to shorten the formulae, the arguments $\left(\sigma ,\lambda \right)$ attached to $\dot{\sigma }$, etc., are totally omitted.

Equation (2) deserves a bit longer explanation.

m in Equation (2) is the number of boundary patches encompassing the elastic domain $E_{\sigma }$ defined by

$$E_{\sigma }\triangleq \left\{\tau \in R^6|y\left(\tau \right)\le 0\right\}$$

(3)

Here, $y:R^6\to R^m$ is a vector-valued function with m component functions, $y_1$, $\cdots$, $y_m$, and $y_i\left(\tau \right)$ = 0 represents the ith boundary patch of $E_{\sigma }$ that is smooth and convex. $y\left(\tau \right)\le$ 0 in Equation (3) means that all the m component functions of $y$ are not positive, namely, $y_i\left(\tau \right)\le$ 0, for $i=1,\cdots ,m$. In this study, what we will concern is elastic and perfectly plastic deformation. Therefore, $E_{\sigma }$ keeps invariant in the stress space.

Take the intact M–C surface as an instance. In the principal stress space of ${\tau }_1$-${\tau }_2$-${\tau }_3$, $E_{\sigma }$ is a domain illustrated in Figure 1a and enclosed by the six planar patches as follows,

$$\begin{array}{c}{y}_1\left({\tau }_1,{\tau }_2,{\tau }_3\right)=\left({\tau }_1-{\tau }_3\right)+\left({\tau }_1+{\tau }_3\right)\sin \varphi -2c\cos \varphi \\ y_2\left({\tau }_1,{\tau }_2,{\tau }_3\right)=\left({\tau }_2-{\tau }_3\right)+\left({\tau }_2+{\tau }_3\right)\sin \varphi -2c\cos \varphi \\ y_3\left({\tau }_1,{\tau }_2,{\tau }_3\right)=\left({\tau }_2-{\tau }_1\right)+\left({\tau }_2+{\tau }_1\right)\sin \varphi -2c\cos \varphi \\ y_4\left({\tau }_1,{\tau }_2,{\tau }_3\right)=\left({\tau }_3-{\tau }_1\right)+\left({\tau }_3+{\tau }_1\right)\sin \varphi -2c\cos \varphi \\ y_5\left({\tau }_1,{\tau }_2,{\tau }_3\right)=\left({\tau }_3-{\tau }_2\right)+\left({\tau }_3+{\tau }_2\right)\sin \varphi -2c\cos \varphi \\ y_6\left({\tau }_1,{\tau }_2,{\tau }_3\right)=\left({\tau }_1-{\tau }_2\right)+\left({\tau }_1+{\tau }_2\right)\sin \varphi -2c\cos \varphi \end{array}$$

(4)

Namely, m = 6; $E_{\sigma }$ has six edges and one apex. At any point on the edges, the normal is not unique. $\varphi$= frictional angle; $c$= cohesion.

If the tension portion is cut off the intact M–C surface, we have the cut-off M–C surface as shown in Figure 1b, and thus $E_{\sigma }$ has another three planar patches as its boundary

$$\begin{array}{c}{y}_7\left({\tau }_1,{\tau }_2,{\tau }_3\right)={\tau }_1\\ y_8\left({\tau }_1,{\tau }_2,{\tau }_3\right)={\tau }_2\\ y_9\left({\tau }_1,{\tau }_2,{\tau }_3\right)={\tau }_3\end{array}$$

(5)

In this case, m = 9, and more edges and vertices are on the boundary of $E_{\sigma }$.

Corresponding to the i-th patch of the yielding surface, i.e., $y_i\left(\tau \right)$ = 0, is the potential function $g_i\left(\tau \right)$ in Equation (2). Replacement of $\varphi$ in $y_i\left(\tau \right)$ in Equation (4) by the dilatancy angle $\psi$ results in $g_i\left(\tau \right)$. Usually, 0$\le \psi \le \varphi$. Taking $\psi =\varphi$, namely, $g_i\left(\tau \right)$ = $y_i\left(\tau \right)$, Equation (2) represents the associative plasticity; otherwise, non-associative plasticity. In any case, Equation (2) tells us that ${\dot{\epsilon }}_p$ is a linear nonnegative combination of all the m gradient vectors of the potential functions, $\nabla g_i\left(\sigma \right)$, $i=1,\cdots ,m$.

The stress point $\sigma$ evolves following Equation (1) and, at the same time, cannot be outside $E_{\sigma }$; in other words, $y\left(\sigma \right)\le$ 0.

${\dot{\lambda }}_i$ in Equation (2) and $-y_i$ constitute the complementarity relationship as follows

$${\dot{\lambda }}_i\ge 0,-y_i\ge 0, {\dot{\lambda }}_i\left(-y_i\right)=0$$

(6)

for $i=1,\cdots ,m$; or the equivalent vector form

$$0\le \dot{\lambda }\perp \left(-y\right)\ge 0$$

(7)

with $\dot{\lambda }$ the m-dimensional vector having m components of ${\dot{\lambda }}_i$. Here, notation “$\perp$” represents “perpendicular to” and “$a\perp b$” is equivalent to $a^{T}b$ = 0, with $a$ and $b$ being two vectors of the same dimension.

Substituting Equation (2) into Equation (1), we have

$$\dot{\sigma }={\dot{\sigma }}_e-{\dot{\sigma }}_p$$

(8)

Here,

$${\dot{\sigma }}_e\triangleq D\dot{\epsilon }$$

(9)

$${\dot{\sigma }}_p\triangleq DJ\dot{\lambda }$$

(10)

$J$ = 6 × m Jacobean matrix of vector-valued function $g$:$R^6\to R^m$, defined as

$$J\left(\sigma \right)\triangleq \left[\nabla g_1\left(\sigma \right),\cdots ,\nabla g_m\left(\sigma \right)\right]$$

(11)

To this point, Equations (7) and (8) prescribe the constitutive relationship of elasto-perfect plasticity in the rate form. The number of Equations (7) and (8) is six and m, respectively, equal exactly to the number of the unknowns, namely, six components of $\sigma$ and m components of $\lambda$. Hence, once the initial values of $\sigma$ and $\lambda$ are provided, in principle, their values at any time or corresponding to any incremental strain $\Delta \epsilon$ can be determined.

Solving Equations (7) and (8) for $\dot{\sigma }$ and $\dot{\lambda }$ is easy for smooth yielding surfaces and has been expounded upon in the literature; for non-smooth yielding surfaces, however, that becomes quite intricate in the classic treatments, such as the return-mapping algorithm, while the incremental strain has to be small enough so as to reach convergence. Fortunately, the treatment to be proposed turns out to be very easy, no matter how many the boundary patches of $E_{\sigma }$ are.

2.2. The Numerical Constitutive Integration

Supposing that the current load step starts at the stress point ${\sigma }_0$ and experiences the incremental strain vector $\Delta \epsilon$, we find the stress point $\sigma$ at the end of the load step. This constitutes constitutive integration.

By the Euler backward integration to Equation (8), we have the nonlinear equations in $\left(\sigma ,\lambda \right)$,

$$\sigma ={\sigma }_e-{\sigma }_p$$

(12)

called the stress decomposition subsequently, where

$${\sigma }_e\triangleq {\sigma }_0+D\Delta \epsilon$$

(13)

is referred to as the elastic trial stress in the literature [23]; and

$${\sigma }_p\triangleq DJ\lambda$$

(14)

Integration of Equation (7) leads to the nonlinear complementarity relationship

$$0\le \lambda \perp \left(-y\right)\ge 0$$

(15)

Compared with the Euler backward approximation (12), the complementarity relationship (15) is accurate with no error introduced.

Among the most commonly used to solve for $\left(\sigma ,\lambda \right)$ from Equations (12) and (15) is the return-mapping algorithm. If ${\sigma }_e$ is outside the yielding surface, the algorithm needs to find the intersection ${\sigma }_I$ of the elastic path ${\sigma }_0$-${\sigma }_e$ and the yielding surface, as illustrated in Figure 2, which needs some computations.

Suppose the fourth and fifth patches of the Mohr–Coulomb yielding surface are active at the current iteration. By stating the i-th patch is active at the current iteration, we mean ${\lambda }_i$ > 0. Then, Equations (12) and (15) reduce to the system of nonlinear equations in $\sigma$, ${\lambda }_4$, and ${\lambda }_5$,

$$\{\begin{array}{l}\sigma ={\sigma }_e-D\left[{\lambda }_4\nabla g_4\left(\sigma \right)+{\lambda }_5\nabla g_5\left(\sigma \right)\right]\\ f_4\left(\sigma \right)=0\\ f_5\left(\sigma \right)=0\end{array}$$

(16)

If both ${\lambda }_4$ and ${\lambda }_5$ are positive at the end of the iteration, the guess of the active patches is right; otherwise, the active patches have to be adjusted at the next iteration, see Simo and Hughes [24]. As a result, the return-mapping algorithm is actually a test and error process.

If the above system is solved using the Newton method, the Hesse matrices of functions $g_4$ and $g_5$ have to be calculated, denoted by

$$G_4\left(\sigma \right)={\left(\frac{{\partial }^2{g}_4}{\partial {\tau }_i\partial {\tau }_j}\right)}_{\tau =\sigma } \mathrm{and} G_5\left(\sigma \right)={\left(\frac{{\partial }^2{g}_5}{\partial {\tau }_i\partial {\tau }_j}\right)}_{\tau =\sigma }$$

(17)

where the subscript i or j being one index of the six components of ${\tau }_x$, ${\tau }_y$, ${\tau }_z$, ${\tau }_{yz}$, ${\tau }_{zx}$, and ${\tau }_{xy}$; for example, $\frac{{\partial }^2{g}_4}{\partial {\tau }_2\partial {\tau }_4}$=$\frac{{\partial }^2{g}_4}{\partial {\tau }_y\partial {\tau }_{yz}}$.

Obviously, the more the non-smooth points on the yielding surface are, the more complicated the return-mapping algorithm becomes, and the more adjustments to the active patches are needed. The tensile strength cut-off creates more non-smooth points on the yielding surface and further burdens the return-mapping algorithm.

Fortunately, the test-and-error adjustments in the return-mapping algorithm are totally avoided in the procedure to be proposed.

2.3. The Mixed Complementarity Problem

Equations (12) and (15) are equivalent to the mixed complementarity problem, stated as follows. Provided the incremental strain vector $\Delta \epsilon$, find $\left(\sigma ,\lambda \right)$$\in R^6\times R_{+}^m$ such that

$$\{\begin{array}{l}0\le \lambda \perp f_I\left(\sigma ,\lambda \right)\ge 0\\ f_E\left(\sigma ,\lambda \right)=0,\hspace{0.17em}\sigma \hspace{0.17em}\mathrm{free}\end{array}$$

(18)

Here, $R_{+}^m$ is the set of all the m-dimensional vectors with all the m components being nonnegative; $f_I:$$R^6\times R_{+}^m\to R_{+}^m$ is called the n-yielding function, defined as

$$f_I\left(\sigma ,\lambda \right)\triangleq -y\left(\sigma \right)$$

(19)

and $f_E:$$R^6\times R_{+}^m\to R^6$ is referred to as the constitutive function, reading

$$f_E\left(\sigma ,\lambda \right)\triangleq \sigma +{\sigma }_p-{\sigma }_e$$

(20)

The above-mixed complementarity problem is abbreviated as MiCP, algorithm for MiCP can be seen in Appendix A.

Suppose the set of solutions to MiCP $\left(f_I,f_E\right)$ is not empty [25]; based on the Gauss–Seidel iteration, we designed a projection–contraction algorithm, abbreviated by GSPC. GSPC is unnecessary to (i) test if the elastic trial stress point ${\sigma }_e$, defined in Equation (13), is outside the elastic domain $E_{\sigma }$; (ii) find the intersection between the elastic path ${\sigma }_0$-${\sigma }_e$ and the yielding surface; (iii) guess which yielding patches are active; and (iv) form the Hesse matrix of any potential function $g_i$.

3. Displacement-Controlled Method Tailored for Plastic Deformation

It is well known that based on a provided nodal incremental load vector, the conventional finite element load-controlled analysis solves a system of nonlinear equations to obtain the nodal incremental displacement vector. Of many existing algorithms for systems of nonlinear equations, the full Newton–Raphson method and the modified Newton–Raphson method are the two most widely used. However, when the load level is near the collapsed load, namely, it is about to reach the limit equilibrium state, the convergence rate of the Newton–Raphson method is then rather poor. A slope in the limit equilibrium state is actually a mechanism at which the tangential stiffness matrix of the finite elements is ranked one deficient [26], causing an invalid load-controlled method (LCM) for solving the system of equilibrium equations. In order to traverse stably limited points of the slope on an equilibrium path, in this section, we will design the displacement-controlled method (no generally accepted names; here we refer to this method as DCM in reference to LCM), which is tailored for elastic-plastic deformation and easy to bring the slope into the limit equilibrium state.

Once the problem domain $\Omega$ is discretized with finite elements, we have the system of nonlinear equilibrium equations in the discrete form

$${\displaystyle \sum _e{\displaystyle {\int }_{{\Omega }^e}{B}^T\sigma d\Omega }}=\rho q+q_0$$

(21)

Here, ${\Omega }^e$ is the domain of a typical element in the mesh; $\sigma$$\in R^6$ is the total stress vector at the end of this load step; $B$ is the 6 × $n_e$ matrix that is related to the incremental strain vector $\Delta \epsilon$$\in R^6$ by the relationship

$$\Delta \epsilon =B{p}_e$$

(22)

where $p_e\in R^{n_e}$ is the incremental displacement degrees of freedom vector of element ${\Omega }^e$, with $n_e$ the number of displacement degrees of freedom of element ${\Omega }^e$.

In Equation (21), $q_0$ is the load vector at the beginning of this load step; $q$ the reference load vector; $\rho q$ the incremental load vector at this step; and $\rho$ the load multiplier, in the conventional LCM, $\rho$ is taken as a series of known values, now, let $\rho$ be unknown and specify some component of $p_e$ as a determined value.

Matrices of $B$, $q_0$ and $q$ can be referred to in the literature on FEM, such as Zienkiewicz and Taylor [24].

Substituting the stress decomposition (12) together with Equations (13) and (22) into Equation (21) leads to the system of nonlinear equations

$$Kp=\rho q+r$$

(23)

Here, $p$ is the incremental displacement degrees of freedom vector of the mesh; $K$ is the elastic stiffness matrix

$$K={\displaystyle \sum _e{\displaystyle {\int }_{{\Omega }^e}{B}^{T}DBd\Omega }}$$

(24)

which keeps variant during the iteration in solving Equation (23); and vector $r$

$$r={\overline{q}}_0+{\displaystyle \sum _e{\displaystyle {\int }_{{\Omega }^e}{B}^T{\sigma }_{p}d\Omega }}$$

(25)

is nonlinearly dependent on $p$, with

$${\overline{q}}_0=q_0-{\displaystyle \sum _e{\displaystyle {\int }_{{\Omega }^e}{B}^T{\sigma }_{0}d\Omega }}$$

(26)

also invariant during iteration.

Here is the iteration scheme for the solution of system (23), where vector $r$ is deemed a known vector and dependent on the previous iterate of $p$, whose initial value will be provided subsequently. As a result, the iteration version of system (23) reads

$$K{p}^{k+1}=\rho q+r^k$$

(27)

where the superscript k is the iteration number. In this way, whether plasticity is associative or non-associative, we always work with the symmetric and positive definite matrix $K$ that keeps invariant within the load step, avoiding a non-symmetric stiffness matrix in the classical Newton–Raphson method and enjoying much higher efficiency.

Multiplying by $K^{-1}$ the two sides of system (27), we have

$$p^{k+1}=\rho p_q+p_r^k$$

(28)

where $p_q$ is the elastic displacement vector due to the reference load vector $q$, namely, the solution to system

$$K{p}_q=q$$

(29)

$p_q$ keeps invariant with iteration; and $p_r^k$ in Equation (28) is the solution to system

$$K{p}_r^k=r^k$$

(30)

and changes with iteration because vector $r^k$ depends on ${\sigma }_p^k$.

The introduction of $\rho$ creates one more unknown than the equilibrium equations in system (21). If $\rho$ is kept invariant, say $\rho$ = 1, the solution scheme is the load-controlled method (LCM), but fitted only to the situation prior to collapse. By letting $\rho$ be unknown but the ith component $p_i^{k+1}$ of $p^{k+1}$ be provided, namely,

$$p_i^{k+1}={\overline{p}}_i$$

(31)

the solution scheme is the displacement-controlled method (DCM) and is fitted to the whole deformation history.

Substituting Equation (31) into the ith equation of Equation (28) and solving for $\rho$ yields

$$\rho =\frac{{\overline{p}}_i-p_{ri}^k}{p_{qi}}$$

(32)

with $p_{qi}$ and $p_{ri}^k$ the ith component of vectors $p_q$ and $p_r^k$. Substituting $\rho$ back to Equation (28), we have the vector $p^{k+1}$. Here, the provided component ${\overline{p}}_i$ of $p^{k+1}$ in Equations (31) and (32) is specified as

$${\overline{p}}_i=p_{qi}$$

(33)

where $\left|p_{qi}\right|$ is the maximum of all the absolute values of components of $p_q$, see Equation (29).

At each step, the equilibrium iteration is convergent if

$${\Vert p^{k+1}-p^k\Vert }_{\infty }<e_p{\Vert p_q\Vert }_{\infty }$$

(34)

with $e_p$ = tolerance of the displacement relative error.

DCM starts at $k$ = 0, and the initial iterate $p^0$ of $p$ is specified as the elastic displacement.

$$p^0=p_q$$

(35)

where $p_q$ is defined in Equation (29).

DCM is proved far more efficient in the nonlinear regime [27]. Since always convergent, DCM will be applied to bring the slope to the limit equilibrium state.

4. The Secant Method for the Factor of Safety

The strength reduction technique for the factor of safety of a slope, $F_s$, is simply to conduct the elastic-perfectly plastic analysis of the slope with the reduced strength parameters

$$c_F=\frac{c}{F} \mathrm{and} \tan {\varphi }_F=\frac{\tan \varphi }{F}$$

(36)

until the slope reaches the limit equilibrium state, at which the reduction factor $F$ is taken as $F_s$.

As for the adjustment to the deformation parameters: Young’s modulus E and Poission’s ratio v, we adopt the strategy proposed by Zheng et al. [28], the adjusted E and v are denoted by $E_F$ and $v_F$, respectively.

Let the gravity generate the reference load vector $q$ in Equation (21), and no other external load exists in the slope, i.e., $q_0$ = 0. With the reduced strength parameters of $c_F$ and ${\varphi }_F$, and the adjusted deformation parameters of $E_F$ and $v_F$, we apply DCM to bring the slope to the limit equilibrium state.

Let ${\rho }_{t}q$ represent the load level under which the slope is at the limit equilibrium state. ${\rho }_t$ is obtained by accumulating the load multipliers ${\rho }_i$ of all the steps before the slope reaches the limit equilibrium state, reading

$${\rho }_t={\displaystyle \sum _{i=1}^N{\rho }_i}$$

(37)

Here, ${\rho }_i$ is the load multiplier of the i-th step; N is the number of all the steps needed for the slope to reach the limit equilibrium state. At the limit equilibrium state, ${\rho }_N$ is 0 theoretically but $\left|{\rho }_N\right|<0.001$ numerically in this study.

Obviously, ${\rho }_t$ is a strictly decreasing function of the strength reduction factor $F$, denoted by ${\rho }_t(F)$, and the factor $F_s$ of safety is the root of the equation, ${\rho }_t(F)=1$, as depicted in Figure 3. Obviously, ${\rho }_t\left(1\right)\le 1$ implies $F_s$ $\le$ 1.

Let $\left(F_1,{\rho }_t^1\right)$ and $\left(F_2,{\rho }_t^2\right)$ be two points on the curve of $F$ vs. ${\rho }_t$, assuming without loss of generality that ${\rho }_t^1$ ≥ 1 and ${\rho }_t^2<1$. Then, execute the following operations.

Step 1. The interpolation of $\left(F_1,{\rho }_t^1\right)$ and $\left(F_2,{\rho }_t^2\right)$ gives rise to a guessed factor of safety

$$F_3=\frac{1-{\rho }_t^1}{{\rho }_t^2-{\rho }_t^1}\left(F_2-F_1\right)+F_1$$

(38)

Step 2. Calculate ${\rho }_t^3$ by taking $F_3$ as the new reduction factor and repeating DCM.

Step 3. If $\left|{\rho }_t^3-1\right|$ < 0.001, set $F_s=F_3$ and terminate.

Step 4. If ${\rho }_t^3$ < 1, then let $F_2=F_3$, ${\rho }_t^2={\rho }_t^3$; otherwise, let $F_1=F_3$, ${\rho }_t^1$ = ${\rho }_t^3$; return to Step 1.

The above solution is actually the secant method. It should be noted that the secant method needs to estimate two initial values; this may require some trial work. Moreover, the secant method is likely to have difficulty when the horizontal axis is tangent to the graph of the curve at a certain point, and there is no guarantee for computational convergence. Fortunately, ${\rho }_t$ is a strictly decreasing function of the strength reduction factor $F$ in this paper.

5. Illustrative Examples

Three typical examples are presented to illustrate the capability of the proposed method in the stability analysis of slopes by accounting for the tension failure. All the examples are assumed in the plain strain condition, with associative plasticity.

All the meshes consist of isoparametric quadrilateral elements, with each having four nodes and deployment of quadrature points as 2 × 2. The displacement tolerance $e_p$ = 1% in Equation (34).

The left and right boundaries are roller boundaries, while the nodes on the bottom boundary are fixed horizontally and vertically.

Two cases are considered for all the examples with the intact M–C surfaces and the cut-off M–C surfaces. $F_{MC}^s$ denotes the factor of safety provided by the intact M–C surface; $F_{Cut}^s$ the factor of safety by the cut-off M–C surface.

5.1. Homogeneous Soil Slope

Figure 4 shows a homogeneous soil slope with a height of 20 m and a slope angle of 45°, designed by Cheng et al. [29]. The soil has a unit weight $\gamma$ = 25 kN/m³, cohesion $c$ = 42 kPa, internal friction angle $\varphi$ = 30°, Young’s modulus $E$ = 30 MPa, Poisson’s ratio $v$ = 0.3.

The outcomes of $F_{MC}^s$ = 1.537 and ′$F_{Cut}^s$ = 1.512 imply that ignoring tension failure overestimates the factor of safety only by 1.65%.

Figure 5 displays the distribution of the equivalent plastic strain ${\overline{\epsilon }}_p$ at the moment when the slope reaches the limit equilibrium state, defined by,

$${\overline{\epsilon }}_p={\displaystyle \sum \sqrt{\frac{2}{3}\Delta {\epsilon }_p^T\Delta {\epsilon }_p}}$$

(39)

where the summation is with regards to all the incremental steps, and $\Delta {\epsilon }_p$ the incremental plastic strain vector within a typical step.

In Figure 5a, the tensile failure zone at the top of the slope, denoted by ${\Omega }_T$, is drawn by the condition

$${\sigma }_1>0$$

(40)

for the intact M–C surface, with ${\sigma }_1$ the major principal stress; while in Figure 5b, ${\Omega }_T$ is marked by

$$\lambda \triangleq \max \left({\lambda }_7,{\lambda }_8,{\lambda }_9\right)>0$$

(41)

for the cut-off M–C surface; ${\lambda }_7$, ${\lambda }_8$ and ${\lambda }_9$ are plastic multipliers associated with the three tensile strength surfaces $f_7$ = 0, $f_8$ = 0 and $f_9$ = 0.

Figure 5 clearly tells us that ignoring the tensile failure brings about the tension failure zone ${\Omega }_T$ extending to the left boundary, an abnormal phenomenon familiar to those who have some experience in using the finite element strength reduction. Considering the tensile failure causes a smaller ${\Omega }_T$, seemingly more reasonable.

Further, included in Figure 5 are the tension cracks. Each tension crack is represented by a vertical line segment that is within ${\Omega }_T$ and passes the quadrature point in ${\Omega }_T$ at which ${\overline{\epsilon }}_p$ assumes the maximum. In this way, the location and depth of the tension crack are automatically determined.

The influence of mesh density is shown in Figure 6, suggesting that with the increase in mesh density,

(1) the failure zone becomes narrower and narrower;

(2) the factor of safety becomes smaller and smaller. This is because the FEM we adopt is based on compatible elements of displacement-type and always approaches from below the exact displacement. Certainly, a larger deformation field results in a more dangerous slope under the same load level and the same boundary displacement constraints;

(3) the location and depth of the tension crack converge prior to the failure zone; for example, at least at the mesh density of 1500 elements, the location and depth of the tension crack have become convergent, yet the failure zone has not.

The depth of tension cracks of slopes in the limit equilibrium state is still an open issue in theory. There have been many works on this. For example, Michalowski [13] proposed an empirical interval of the depth of tension cracks for homogeneous slopes with no water

$$h=\frac{\kappa c_m}{\gamma }\tan \left(\frac{\pi }{4}+\frac{{\varphi }_m}{2}\right)$$

(42)

where $\kappa$ varies between 2 and 3.83; $c_m$ and ${\varphi }_m$ are the mobilized strength parameters, namely,

$$c_m=\frac{c}{F_s} \mathrm{and} \tan {\varphi }_m=\frac{\tan \varphi }{F_s}$$

(43)

By taking $F_s$ = $F_{Cut}^s$ = 1.512, the h-interval evaluated by Equation (42) is (3.23 m, 6.18 m).

The convergent depth estimated by the proposed analysis is 2.80 m, even below the lower bound of 3.23 m evaluated by Equation (42).

With the same mechanical parameters as the above example, let the slope angle $\beta$ be 50°, 60°, 70°, 80°, and 90°.

Table 1. Outcomes of different slope angles.

Angle	${\mathit{F}}_{\mathit{M}\mathit{C}}^{\mathit{s}}$	${\mathit{F}}_{\mathit{C}\mathit{u}\mathit{t}}^{\mathit{s}}$	Difference
45°	1.512	1.537	1.63%
50°	1.386	1.412	1.84%
60°	1.18	1.218	3.12%
70°	1.012	1.06	4.53%
80°	0.851	0.928	8.3%
90°	0.652	0.814	19.9%

lists the results, suggesting that before $\beta$ exceeds 70° the difference in the factor of safety can be ignored using the intact M–C surface or the cut-off M–C surface. If $\beta$ > 80°, the difference increases sharply and reaches 20% at $\beta$ = 90°. Consequently, the stability of deep foundation pits must be analyzed under the condition of tensile strength cut-off.

By taking $F_s$ = $F_{Cut}^s$ = 1.06 corresponding to $\beta$ = 70°, the interval of the tension crack depth is (17.02 m, 32.60 m) according to the evaluation (42), while the slope height is only 20 m, and thus the overestimation is concluded. The depth evaluated by the proposed procedure is 4.8 m, as we can see in Figure 7, a seemingly more practical value. Using the intact M–C surface will result in the tensile failure zone extending to the model boundary and is severely overestimated. Using the cut-off M–C surface, the abnormal is eliminated.

5.2. A Layered Slope with Weak Interlayer

This example concerns a layered slope where a layer of low resistance is interposed between two layers of higher strength, see Greco [30] and Arai and Tagyo [31]. Geometrical features of the slope and mechanical parameters are in Figure 8 and

Table 2. Mechanical Parameters of Section 5.2.

Soil	$\mathit{\gamma } (\mathbf{k}\mathbf{N}/{\mathbf{m}}^{\mathbf{3}})$	$\mathit{c} (\mathbf{k}\mathbf{N}/{\mathbf{m}}^{\mathbf{2}})$	$\mathit{\phi } (\mathbf{°})$	$\mathit{E} (\mathbf{M}\mathbf{P}\mathbf{a})$	$\mathit{v}$
S1	18.82	29.4	12	10	0.35
S2	18.82	9.8	5	10	0.35
S3	18.82	29.4	40	10	0.35

.

Calculated by Greco and Arai and Tagyo are $F_s$ = 0.388 and $F_s$ = 0.405 without tension crack, respectively. Using the proposed method, $F_{MC}^s$ = 0.417 and $F_{MC}^s$ = 0.391. The distributions of the equivalent strain are shown in Figure 9; it can be seen that the plastic failure characteristics at the back edge of the slope for the two kinds of yielding surfaces are quite different. As the tension cut-off is considered, due to the tensile failure of the soil at the top of the slope, a tensile crack will run through the first soil layer, the sliding body moves forward at the top of the slope, and the safety factor of the slope is reduced by about 6.2%. While only considering shear failure, a slide will occur along the weak interlayer at the rear edge of the slope. Because the factor of safety is much less than 1.0 and is impossible to sustain, neither the tensile failure zone nor the tension crack is drawn.

5.3. Multistage Non-Homogeneous Soil Slope

Figure 10a shows the geometry of a slope with a strong interlayer. The first grade of the slope is 8 m with a rate of 1:1, and the second is 10 m with a rate of 1:0.75. A thin layer of hard soil (A3) with a thickness of 0.8 m is embedded in the layered soil with slightly lower strength. Figure 10b is the mesh of 943 elements. The soil parameters are in

Table 3. Mechanical parameters of Section 5.3.

Soil	$\mathit{\gamma } (\mathbf{k}\mathbf{N}/{\mathbf{m}}^{\mathbf{3}})$	$\mathit{c} (\mathbf{k}\mathbf{N}/{\mathbf{m}}^{\mathbf{2}})$	$\mathit{\phi } (°)$	$\mathit{E} (\mathbf{M}\mathbf{P}\mathbf{a})$	$\mathit{v}$
S1	18.8	19.8	15.7	10	0.4
S2	19	22.7	16.9	12	0.4
S3	19.5	39.3	24.2	14	0.35
S4	19	27.6	21.3	12	0.35
S5	19	35.7	27.4	15	0.3

.

The proposed factors of safety are $F_{MC}^s$ = 1.142 and $F_{Cut}^s$ = 1.126, differing merely by 1.4%.

Still, significant differences exist in the failure mechanism of the slope in the limit equilibrium state; as shown in Figure 11, it is obvious that the finite strength reduction with tensile–shear composite yield surface can obtain almost the same shear sliding surface at the front edge of the slope compared with the intact M–C shear yield criterion, the difference lies in the plastic zone obtained by the tension-shear failure extends from the top of the slope to a certain depth in a nearly vertical state, which is consistent with the vertical tension cracks in the actual project. The positions of tensile cracks obtained by the two methods are also different. Again, the intact M–C surface results in extending the tensile failure zone to the left model boundary.

For a multilayered slope, such as this example, Equation (42) is not applicable to evaluate the depth of the tension crack because we are even unclear what layers are selected while calculating $c_m$ and ${\varphi }_m$ in this equation.

6. Conclusions

For complex slopes such as layered slopes, the failure mechanism is recognizable only by numerical methods such as the finite element method with tension strength cut-off. With meshes finer and finer, the factor of safety approximates above its accurate value. For deeper slopes of angle exceeding 70, the factor of safety differs significantly using the intact Mohr–Coulomb yielding surface and its cut-off version; for gentle slopes, the difference in the factor of safety is ignorable, but an essential difference exists in the failure mechanism, especially near the slope top. The empirical estimate (42) is suitable only for homogenous slopes, and its lower bounds of the depth of the tension cracks match better with the values from the proposed procedure.

The Gauss–Seidel iteration-based projection–contraction algorithm, abbreviated by GSPC, is very well qualified to calculate numerical constitutive integration in the strength reduction technique for slope stability analysis using the Mohr–Coulomb yielding surface and its cut-off version. The displacement-controlled method (DCM) tailored for the elastic-perfectly plastic deformation analysis is easy to bring the slope into the limit equilibrium state, only at which is the failure mechanism able to be accurately recognized.