Stability analysis of slopes with planar failure using variational calculus and numerical methods

Abstract

This study investigates the technique of variational calculus applied to estimate the slope stability considering the mechanism of planar failure. The critical plane failure surface should be determined because it theoretically indicates the most unfavorable plane to be considered when stabilizing a slope to rectify the instability generated by several statistically possible planes. This generates integrals that can be solved by numerical methods, such as the Newton Cotes and the finite differences methods. Additionally, a system of nonlinear equations is obtained and solved. The surface of the critical planar failure is determined by applying the condition of transversality in mobile boundaries, for which various examples are provided. The number of slices is varied in one of the examples, while the surface of the critical planar failure is determined in the others. Results are compared using analytical methods through axis rotations. All the results obtained by considering normal stress, safety factors, and critical planar failure are nearly the same; however, in this research, a study is carried out for “n” number of slices using programming methods. Sub-routines are important because they can be applied in slopes with different geometry, surcharge, interstitial pressure, and pseudo-static load.

Keywords

slopes stability / planar failure / variational calculus / numerical methods

Introduction

Plane failure occurs due to the existence of layers or strata of different composition Although plane failure in soils is practically non-existent, it can occur in soil coverings on rocks or in rocky slopes owing to the presence of parallel discontinuities in the slope. However, the critical plane failure surface should be determined to theoretically indicate the most unfavorable plane when stabilizing a slope to rectify the instability generated from several statistically possible planes.

Slope stability in soils and rock mass has been of interest to researchers of different fields. Consequently, several techniques and methods for evaluation and safety factor (SF) calculation have been developed. Limit equilibrium is one of those methods, wherein it compares the resisting forces with the mobilized forces to estimate the SF [¹,²]. The other methods comprise numerical models that show stress distribution and displacements [³], the upper bound analysis [⁴–⁶], and the slope stability analysis considering seismic effects [⁷–¹⁰].

A method to perform the three-dimensional slope stability analysis has been developed using independent cover based numerical manifold (ICMM3D) and vector sum methods (VSM) [¹¹], i.e., here, a genetic algorithm is employed to search the critical slip surface by a simple and robust three-dimensional cracking-particle method, without modification. This describes a new robust and efficient approach for modeling discrete cracks via mesh-free methods, wherein the crack is modeled by a set of cracked segments [¹²]. Additionally, a dual-horizon peridynamics (DHPD) formulation that naturally includes varying horizon sizes and completely solves the “ghost force” issue was proposed [¹³]. Ref. [¹⁴] reported that the traditional peridynamics can be derived as a special case of the present DHPD. Furthermore, a simple adaptive refinement procedure was proposed to reduce the computational cost. To demonstrate the capability of this method, tests were conducted for two- and three- dimensional examples including the Kalthoff-Winkler experiment. Tests were performed on a plate with branching cracks to demonstrate the capability of the method. Therefore, a DHPD formulation was developed that allows for simulations with dual horizon and minimal spurious wave reflection, based on which the balance of momentum and angular momentum in PD are naturally satisfied.

A phase field model that can simulate compressive-shear fractures in rock-like materials was investigated [¹⁵]. A new driving force was introduced in the evolution equation of phase field; thereafter, a hybrid formulation was established for the phase field modeling. A phase field model for fluid-driven dynamic crack propagation in poroelastic media was developed [¹⁶]. Ref. [¹⁷] employed a staggered scheme and implemented our approach into a software named COMSOL multiphysics. Finally, complex crack patterns of a plate subjected to increasing internal pressure, the (3D) Pertersson beam, and a 3D NSCB test, were presented.

This study applies the technique of variational calculus to estimate slope stability, which was first proposed by Ref. [¹⁸]. This use of this technique in homogeneous soils excludes both the interstitial pressures and external load It is used to obtain the critical combination between failure line y(x) and normal stress s(x) in order to satisfy all the equilibrium equations [¹⁹,²⁰]. Furthermore, a method based on global equilibrium that minimizes the weight of the sliding mass is proposed [⁵]. Variational calculus is applied to estimate the slope stability by generalizing Euler’s equations, natural boundary conditions and transversality to obtain the critical line [²¹]. All of the variational methods highlight its limitations [²²]. All of these research papers failed in arriving at a solution to solve the nonlinear equation system.

A general and extensive analysis [²⁴] using the technique of variational calculus to determine sliding surface and normal stress distribution [²³] is proposed. Once the minimum SF and the normal stress function has been determined, Ref. [²⁴] proposed a numerical method to obtain the normal stresses (s_n) by applying the third degree Lagrange’s polynomial, because the method proposed in Ref. [²⁴] results in an indeterminate normal stresses equation. Here, the solution to the problem lies in solving a set of nonlinear equations using the finite differences technique. The principles of variational calculus are applied to determine the two functions, s(x) and y(x), representing the normal stress along the y(x) line, i.e., the failure line of the slope. These functions relate to each other at the point of contact of the acting forces to establish three equilibrium equations. Applying Euler’s equations, a solution of one of the equations is obtained. A change of variable is performed to obtain the equation of the failure line or sliding in polar coordinates [²⁵].

However, the mass of soil or rock does not require changes in terms of its characteristics and properties, and only the equilibrium between the horizontal and vertical forces are considered for two main reasons: the failure surface is plane and it uses the condition of transversality. Assuming that the failure is produced due to the sliding action, wherein the moments generated on the failure plane are considered to have negligible effect, it is suitable to consider that all the forces pass through the center of gravity of the potential slip surface [²⁶].

Finally, a criterion of nonlinear failure for the variational analysis on the stability of the slope is established. The main advantage of the traditional slope stability estimation method is that it does not introduce any static or geometric hypothesis [²⁷]. In particular, Ref. [²⁷] demonstrated that the critical sliding surface has a spiral shape. In addition, the critical normal forces are distributed based on convex functions having a single maximum on its definition range. The critical sliding surfaces are singular in the points where the normal forces approach the tensile strength of the material. The solution of the equations product of the application of the technique has not yet been developed. However, Ref. [²⁷] showed that the logarithmic spiral surface was highly critical but did not solve the equations obtained nor determined the SF used in practical engineering cases.

The main purpose of this study is to extend the work of Ref. [²⁶], wherein the technique of variational calculus was applied to estimate the slope stability; however, this study had a limited number of slices. The advantage of the variation calculus technique over the strength reduction and limit equilibrium methods is that it considers all equations without simplifying them. Therefore, the method incorporating variational calculus is more appropriate when highly complicated geometries are used. Here, we perform programming on an infinitenumber of slices. Moreover, a new advanced modification is developed with respect to the equation of the straight line of the planar failure mechanism, which is not fixed and indefinite on the crest or on the surface of the slope. It is only considered known at the base of the slope that coincides with the origin of coordinates (0,0). Simultaneously, iteration graphs are generated where both the minimum SF and the slope of the straight line of the planar failure are determined for comparison with results obtained from both the minimum equilibrium method (Slide software) and the analytic method of axis rotation [²⁸].

Methodology

Basic assumptions

The technique of variational calculus provides mathematical procedures to obtain the shape of the curve as an extremal, i.e., the curve that maximizes or minimizes the value of the integral along the curve determined under boundary conditions [²⁹]. Therefore, by applying the aforementioned mathematical method along with the conditions of static equilibrium, potential line of failure, and formulated in the form of integrals with certain boundary conditions, it is possible to minimize those integrals [³⁰]. This enables us to determine the extreme surface that results in the minimum SF of the slope. The method is accurate because all equilibrium equations are satisfied.

When the technique of variational calculus is used in the stability of the slope, we consider the following aspects.

1) The Lagrange multipliers to determine the functional G.

2) The Euler equations used to determine a combined analytical expression of normal stresses and the curve of the failure mechanism.

3) The transversality condition for the mobile boundaries in the slope surface.

4) A numerical method to solve integrals.

5) Programming and subroutines to solve the system of nonlinear equations.

Mechanism of planar failure

The equations required to analyze stability are presented below. Applying variational calculus, where three equilibrium equations are considered, the summation of horizontal and vertical forces and moments are obtained. A system of nonlinear equations is derived considering that the SF and the normal stresses acting on the potential surface of failure must be obtained through numerical methods. Stability is analyzed by applying the Mohr-Coulomb linear failure criterion in non-homogeneous soil and rock, indicating that the shear parameters within the analyzed domain vary depending on the type of soil or rock mass. In addition, the pseudo-static load, interstitial pressure, and surcharge are considered.

The SF equation, considering the limit equilibrium method and the Mohr-Coulomb criteria, can be expressed through the following equation:

(1)FS={\displaystyle \frac{\left[C+\left(\sigma -u\right)\tan \phi \right]}{\tau }},

where C represents the cohesion, tanf is the friction coefficient, u is the pore pressure and, t is the shear stress in the potential failure surface.

The equilibrium equations are obtained considering the forces acting in Fig. 1, and SF in Eq. (1). By summing the horizontal and vertical forces, and moments, we obtain:

\mathrm{\sum }{F}_{\mathrm{h}}=0,

{\int }_{x_0}^{x_{\mathrm{m}}}[\tau -\sigma y\prime \left(x\right)-K_{\mathrm{h}}\gamma \left(g(x)-y\left(x\right)\right)+\frac{q_x}{\cos {\beta }_{\mathrm{t}}}]\mathrm{d}x=0,

(2){\int }_{x_0}^{x_{\mathrm{m}}}\left[\begin{array}{c}C+(\sigma -u)\tan \phi \hfill \\ -FS[\sigma y\prime \left(x\right)+k_{\mathrm{h}}\gamma \left(g(x)-y\left(x\right)\right)+{\displaystyle \frac{q_x}{\cos {\beta }_{\mathrm{t}}}}]\hfill \end{array}\right]\mathrm{d}x=0,

\mathrm{\sum }{F}_{\mathrm{v}}=0,

{\int }_{x_0}^{x_{\mathrm{m}}}[\tau y\prime \left(x\right)+\sigma +{\displaystyle \frac{q_y}{\cos {\beta }_{\mathrm{t}}}}-\gamma (1+K_{\mathrm{v}})(g(x)-y\left(x\right))]\mathrm{d}x=0,

(3){\int }_{x_0}^{x_{\mathrm{m}}}\left[\begin{array}{c}[\begin{array}{c}C+(\sigma -u)\tan \phi \hfill \end{array}]y\prime \left(x\right)\hfill \\ +FS\left[\sigma +{\displaystyle \frac{q_y}{\cos {\beta }_{\mathrm{t}}}}-\gamma (1+K_{\mathrm{v}})(g(x)-y\left(x\right))\right]\hfill \end{array}\right]\mathrm{d}x=0,

\mathrm{\sum }M=0,

\left\{\begin{array}{c}{\begin{array}{c}\int \hfill \end{array}}_s\mathit{R}\times (\tau {\overrightarrow{t}}_{\mathrm{o}}+\sigma {\overrightarrow{n}}_{\mathrm{o}})\mathrm{d}s+\begin{array}{c}\int \hfill \end{array}\mathit{p}\times \mathit{q}\mathrm{d}l+\hfill \\ \int x\mathit{i}\times [-\gamma (1+K_{\mathrm{v}})(g(x)-y\left(x\right))\mathit{j}]\mathrm{d}s+\hfill \\ {\int }_{x_0}^{x_{\mathrm{m}}}\left\{\left[{\displaystyle \frac{1}{2}}(g(x)-y\left(x\right))\right]\mathit{j}\times [-\gamma K_{\mathrm{h}}(g(x)-y(x))]\mathit{j}\right\}\mathrm{d}x\hfill \end{array}\right\}=0,

(4){\int }_{x_0}^{x_{\mathrm{m}}}\left\{\begin{array}{c}[C+(\sigma -u)\tan \phi ][y\left(x\right)-xy\prime \left(x\right)]\hfill \\ -FS\left[\begin{array}{c}\sigma (x+y\left(x\right)y\prime (x\left)\right)+{\displaystyle \frac{1}{\cos {\beta }_{\mathrm{t}}}}(x{q}_y-q_{x}g\left(x\right))\hfill \\ -\gamma (1+K_{\mathrm{v}})(g(x)-y\left(x\right))+\gamma k_{\mathrm{h}}\left[{\displaystyle \frac{{g\left(x\right)}^2-{y\left(x\right)}^2}{2}}\right]\hfill \end{array}\right]\hfill \end{array}\right\}\mathrm{d}x=0,

where C is the cohesion of the material, u is the pore pressure, q_x and q_y are surcharges in x- and y-axis, tanf is the coefficient of internal friction, and y′(x) is the derivative of the failure surface with respect to the x-axis.

Subsequently, a functional G is required to group Eqs. (2)–(4), and the Newton-Cotes open compound integration method is applied (Midpoint composed rule). To apply this method, the Lagrange multipliers (λ₁ and λ₂) are used, that is, the functional is G=\mathrm{\sum }{F}_{\mathrm{h}}+{\lambda }_1\mathrm{\sum }{F}_{\mathrm{v}}+{\lambda }_2\mathrm{\sum }M, where the first, second and third terms represent the sum of horizontal forces, vertical forces, and moments, respectively. Simplifying Eq. (5), we obtain:

G=\left[\begin{array}{c}C+(\sigma -u)\tan \phi \hfill \\ -FS[\sigma y\prime \left(x\right)+k_{\mathrm{h}}\gamma (g(x)-y\left(x\right))-{\displaystyle \frac{q_x}{\cos {\beta }_{\mathrm{t}}}}]\hfill \end{array}\right]

+{\lambda }_1\left[\begin{array}{c}[\begin{array}{c}C+(\sigma -u)\tan \phi \hfill \end{array}]y\prime \left(x\right)\hfill \\ +FS\left[\sigma +{\displaystyle \frac{q_y}{\cos {\beta }_{\mathrm{t}}}}-\gamma (1+K_{\mathrm{v}})(g(x)-y\left(x\right))\right]\hfill \end{array}\right]

(5)+{\lambda }_2\left[\begin{array}{c}[\begin{array}{c}C+(\sigma -u)\tan \phi \hfill \end{array}][y\left(x\right)-xy\prime \left(x\right)]\hfill \\ -FS\left[\begin{array}{c}\sigma +(x+y\left(x\right)y\prime (x\left)\right)+{\displaystyle \frac{1}{\cos {\beta }_{\mathrm{t}}}}-(x{q}_y-q_{x}g\left(x\right))\hfill \\ -\gamma (1+k_{\mathrm{v}})(g\left(x\right)-y\left(x\right))+\gamma k_{\mathrm{h}}\left[{\displaystyle \frac{{g\left(x\right)}^2-{y\left(x\right)}^2}{2}}\right]\hfill \end{array}\right]\hfill \end{array}\right]\mathrm{.}

Besides, considering Euler’s equation, in this case, {\displaystyle \frac{\partial G}{\partial y}}={\displaystyle \frac{\mathrm{d}}{\mathrm{d}x}}\left({\displaystyle \frac{\partial G}{\partial y\text{'}}}\right) and calculating the derivatives of the functional G, we obtain:

(6)\begin{array}{c}[\tan \phi (x-x_{\mathrm{c}})+FS(y-y_{\mathrm{c}})]\left({\displaystyle \frac{\mathrm{d}\sigma }{\mathrm{d}x}}\right)\hfill \\ \\ \\ -(x-x_{\mathrm{c}})\left[\right({\displaystyle \frac{\mathrm{d}u}{\mathrm{d}x}})\tan \phi +\gamma (1+K_{\mathrm{v}})FS]+2C\hfill \\ \\ \\ +2(\sigma -u)\tan \phi +FS{K}_{\mathrm{h}}\gamma (y-y_{\mathrm{c}})=0\hfill \end{array},

where x_c and y_c are geometric centers dependent on the Lagrange multipliers (λ₁ and λ₂) and is related through, y_{\mathrm{c}}=-{\displaystyle \frac{\mathrm{1}}{{\lambda }_{\mathrm{2}}}}\mathrm{;}{x}_{\mathrm{c}}={\displaystyle \frac{{\lambda }_{\mathrm{1}}}{{\lambda }_{\mathrm{2}}}}.

Calculation of variations in mobile boundaries

In the problem of unfixed boundaries, we assume that one or both boundary points can move, allowing the admissible curves to expand. We then compare curves that have common points with the analyzed curve. Furthermore, curves with displaced border points can be considered. Note that in the planar failure mechanism, the failure passes through the base of the slope, point (0,0), and the equation of the straight line is represented by y(x) = mx, reaching one extreme. Therefore, the fundamental necessary extreme condition in the problem of unfixed boundaries or frontiers must be known, that is, the function y(x) must be a solution to Euler’s equation. Accordingly, the curve y = y(x), in which the extreme in the problem of unfixed frontiers is obtained, must be an extremal.

Additionally, for the surface or crest of the slope, the transversality condition shown in Ref. [²⁵] is used, which is represented by the following equation:

(7){[G+(g\prime \left(x\right)-y\prime \left(x\right))\frac{\partial G}{\partial y\prime }]}_{x=x_{\mathrm{m}}}=0,

where g′(x) is the derivative of the surface of the slope with respect to the axis of the abscissas and y′(x) is the derivative of the equation of the surface of failure. Substituting the functional G of Eq. (5) and (∂G/∂y′) for the mobile coordinate points (x_m, y_m), that is, the crest or surface of the slope, in Eq. (7), the boundary condition y (x_m) = g(x_m) can be determined. On simplifying this equation, we obtain the transversality equation:

\begin{array}{c}g\prime (x)\left[\begin{array}{c}-(C_{\mathrm{m}}+({\sigma }_{\mathrm{m}}-u_{\mathrm{m}})\tan {\phi }_{\mathrm{m}})(x_{\mathrm{m}}-x_{\mathrm{c}})\hfill \\ -FS{\sigma }_{\mathrm{m}}(y_{\mathrm{m}}-y_{\mathrm{c}})\hfill \end{array}\right]\hfill \end{array}

(8)\begin{array}{c}\begin{array}{c}\begin{array}{c}+(C_{\mathrm{m}}+({\sigma }_{\mathrm{m}}-u_{\mathrm{m}})\tan {\phi }_{\mathrm{m}})(y_{\mathrm{m}}-y_{\mathrm{c}})\hfill \end{array}\hfill \\ \begin{array}{c}\begin{array}{c}-FS{\sigma }_{\mathrm{m}}(x_{\mathrm{m}}-x_{\mathrm{c}})+\frac{FS{q}_x}{\cos {\beta }_{\mathrm{t}}}(g_{\mathrm{m}}-y_{\mathrm{c}})\hfill \end{array}\hfill \\ -\frac{FS{q}_y}{\cos {\beta }_{\mathrm{t}}}(x_{\mathrm{m}}-x_{\mathrm{c}})=0\mathrm{.}\hfill \end{array}\hfill \end{array}\hfill \end{array}

To determine the worst potential plane curve, we consider that a curve passes through the base of the slope and the coordinate on the crest is a moving boundary (see Fig. 2). Therefore, the transversality equation is used for the slope crest, Eq. (8). The stress and the corresponding coordinate point in the slope crest are determined using the equations, which matches additional terms in the expansion of the Taylor series, as observed in the following equation:

(9)f\text{'}\left(x_n\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\Delta }x}}\left(\mathrm{3}{f}_n-\mathrm{4}{f}_{n-\mathrm{1}}+f_{n-\mathrm{2}}\right),

(10)f\text{'}\left(x_n\right)={\displaystyle \frac{f_{n+\mathrm{1}}-f_{n-\mathrm{1}}}{{\displaystyle \frac{H}{\omega }}-\left(\left(r-\mathrm{2}\right)+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\right)\mathrm{\Delta }x}},

where H is the height of the slope, ω is the tanα slope of the failure surface, r is the number of slices in which the potentially sliding mass is divided, and Dx is the width of the slice. Equating Eqs. (9) and (10), we obtain:

f_{n+\mathrm{1}}=\left({\displaystyle \frac{{\displaystyle \frac{H}{\omega }}-\left(\left(r-\mathrm{2}\right)+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\right)\mathrm{\Delta }x}{\mathrm{2}\mathrm{\Delta }x}}\right)

(11)\left(\mathrm{3}{f}_n-\mathrm{4}{f}_{n-\mathrm{1}}+f_{n-\mathrm{2}}\right)+f_{n-\mathrm{1}},

where f_n+1 represents a function that allows to determine the stress on the crest of the slope.

Application of numerical methods and programming

Equilibrium equations (Eqs. (2)–(4)) are obtained by applying variational calculus to slope stability. These equations are solved through numerical methods using the numerical integration proposed by Newton-Cotes and the finite differences method. Nonlinear equations systems are derived and are used for programming and generating subroutines to feed the math software (Table 1). In this case, Maple is used.

Programming of the sum of horizontal forces Eq. (2)

for n from 3 by 1 to r do

f_{\mathrm{1}}:=f_{\mathrm{1}}

+hb\left(\begin{array}{l}C+\left({\sigma }_n-u_n\right)\hspace{0.17em}y-FS\hspace{0.17em}{\sigma }_n\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\Delta }x}}\right)\left(\mathrm{3}{y}_n-\mathrm{4}{y}_{n-\mathrm{1}}+y_{n-\mathrm{2}}\right)\hfill \\ -FS\hspace{0.17em}kh\hspace{0.17em}pu\left(g_n-y_n\right)\hfill \end{array}\right)

if n>d, then

f_{\mathrm{1}}:=f_{\mathrm{1}}

+hb\left(\begin{array}{l}C+\left({\sigma }_{n-d}-u_{n-d}\right)\hspace{0.17em}\psi \hfill \\ -FS\hspace{0.17em}{\sigma }_{n-d}\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\Delta }x}}\right)\left(-\mathrm{3}{y}_{n-d}+\mathrm{4}{y}_{n-d+\mathrm{1}}-y_{n-d+\mathrm{2}}\right)\hfill \\ -FS\hspace{0.17em}kh\hspace{0.17em}pu\left(g_{n-d}-y_{n-d}\right)\hfill \end{array}\right)

end if

end do

Programming of the sum of vertical forces Eq. (3)

for n from 3 by 1 to r do

f_{\mathrm{2}}:=f_{\mathrm{2}}

+hb\left(\begin{array}{l}\left(C+\left({\sigma }_n-u_n\right)\hspace{0.17em}\psi \right)\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\Delta }x}}\right)\left(\mathrm{3}{y}_n-\mathrm{4}{y}_{n-\mathrm{1}}+y_{n-\mathrm{2}}\right)\hfill \\ +FS\hspace{0.17em}{\sigma }_n-FS\hspace{0.17em}\hspace{0.17em}pu\left(\mathrm{1}+kv\right)\left(g_n-y_n\right)\hfill \end{array}\right)

if n>d, then

f_{\mathrm{2}}:=f_{\mathrm{2}}

+hb\left(\begin{array}{l}\left(C+\left({\sigma }_{n-d}-u_{n-d}\right)\hspace{0.17em}\psi \right)\cdot \left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\Delta }x}}\right)\hfill \\ \left(-\mathrm{3}{y}_{n-d}+\mathrm{4}{y}_{n-d+\mathrm{1}}-y_{n-d+\mathrm{2}}\right)\hfill \\ +FS\hspace{0.17em}{\sigma }_{n-d}-FS\hspace{0.17em}pu\left(\mathrm{1}+kv\right)\left(g_{n-d}-y_{n-d}\right)\hfill \end{array}\right)

end if

end do

Programming of the sum of moments Eq. (4)

for n from 3 by 1 to r do

f_{\mathrm{3}}:=f_{\mathrm{3}}+hb\left(\begin{array}{l}\left(C+\left({\sigma }_n-u_n\right)\hspace{0.17em}\psi \right)\hfill \\ \left(y_n-x_n\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\Delta }x}}\right)\left(\mathrm{3}{y}_n-\mathrm{4}{y}_{n-\mathrm{1}}+y_{n-\mathrm{2}}\right)\right)\hfill \\ -FS\hspace{0.17em}{\sigma }_n\left(x_n+y_n\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\Delta }x}}\right)\left(\mathrm{3}{y}_n-\mathrm{4}{y}_{n-\mathrm{1}}+y_{n-\mathrm{2}}\right)\right)\hfill \\ +FS\hspace{0.17em}\hspace{0.17em}{x}_n\hspace{0.17em}pu\left(\mathrm{1}+kv\right)\left(g_n-y_n\right)-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}FS\hspace{0.17em}pu\hspace{0.17em}kh\left(g_n^{\mathrm{2}}-y_n^{\mathrm{2}}\right)\hfill \end{array}\right)

if n>d, then

f_{\mathrm{3}}:=f_{\mathrm{3}}+hb\left(\begin{array}{l}\left(C+\left({\sigma }_{n-d}-u_{n-d}\right)\hspace{0.17em}\psi \right)\cdot \left(y_{n-d}-x_{n-d}\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\Delta }x}}\right)\left(-\mathrm{3}{y}_{n-d}+\mathrm{4}{y}_{n-d+\mathrm{1}}-y_{n-d+\mathrm{2}}\right)\right)\hfill \\ -FS\hspace{0.17em}{\sigma }_{n-d}\left(x_{n-d}+y_{n-d}\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\Delta }x}}\right)\left(-\mathrm{3}{y}_{n-d}+\mathrm{4}{y}_{n-d+\mathrm{1}}-y_{n-d+\mathrm{2}}\right)\right)\hfill \\ +FS\hspace{0.17em}\hspace{0.17em}{x}_{n-d}\hspace{0.17em}pu\left(\mathrm{1}+kv\right)\left(g_{n-d}-y_{n-d}\right)-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}FS\hspace{0.17em}pu\hspace{0.17em}kh\left(g_{n-d}^{\mathrm{2}}-y_{n-d}^{\mathrm{2}}\right)\hfill \end{array}\right)

end if

end do

Programming the condition of Euler Eq. (6)

for n from 3 by 1 to r do

f_{n+1}:=\left(\begin{array}{l}\left(\mathit{\Psi }\left(x_n-x_{\mathrm{c}}\right)\hspace{0.17em}+FS\left(y_n-y_{\mathrm{c}}\right)\right)\left(\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\Delta }x}}\right)\left(\mathrm{3}{\sigma }_n-\mathrm{4}{\sigma }_{n-1}+{\sigma }_{n-\mathrm{2}}\right)\right)\hfill \\ -\left(x_n-x_{\mathrm{c}}\right)\left(\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\Delta }x}}\right)\left(\mathrm{3}{u}_n-\mathrm{4}{u}_{n+1}+u_{n+\mathrm{2}}\right)\right)\psi +(1+kv)FSpu\hfill \\ +2C+2\left({\sigma }_n-{\sigma }_n\right)\psi +FSkhpu\left(y_n-y_{\mathrm{c}}\right)\hfill \end{array}\right)

end do

for n from 1 by 1 to 2 do

f_{n+1+r}:=\left(\begin{array}{l}\left(\mathit{\Psi }\left(x_n-x_{\mathrm{c}}\right)\hspace{0.17em}+FS\left(y_n-y_{\mathrm{c}}\right)\right)\left(\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\Delta }x}}\right)(-\mathrm{3}{\sigma }_n+\mathrm{4}{\sigma }_{n+1}-{\sigma }_{n+\mathrm{2}})\right)\hfill \\ -\left(x_n-x_{\mathrm{c}}\right)\left(\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{\Delta }x}}\right)\left(-\mathrm{3}{u}_n+\mathrm{4}{u}_{n+1}-u_{n+\mathrm{2}}\right)\right)\psi +(1+kv)FSpu\hfill \\ +2C+2\left({\sigma }_n-u_n\right)\psi +FSkhpu\left(y_n-y_{\mathrm{c}}\right)\hfill \end{array}\right)

end do

Condition of transversality Eq. (8)

\begin{array}{l}{f}_{r+4}:\psi p(-(C+({\sigma }_{r+1}-u_{r+1})\psi )(x_{r+1}-x_{\mathrm{c}})-{FS\sigma }_{r+1}(y_{r+1}-y_{\mathrm{c}}))\hfill \\ +(C+({\sigma }_{r+1}-u_{r+1})\psi )(y_{r+1}-y_{\mathrm{c}})-{FS\sigma }_{r+1}(x_{r+1}-x_{\mathrm{c}})\hfill \\ +{FSq}_x(g_{r+1}-y_{\mathrm{c}})-{FSq}_y(x_{r+1}-x_{\mathrm{c}})\hfill \end{array}

where \begin{array}{l}{x}_{r+\mathrm{1}}:=\left(x_r+h\right)\hfill \\ y_{r+\mathrm{1}}:=\left(\left({\displaystyle \frac{\mathrm{\Delta }x+h}{\mathrm{2}\mathrm{\Delta }x}}\right)\left(\mathrm{3}{y}_r-\mathrm{4}{y}_{r-\mathrm{1}}+y_{r-\mathrm{2}}\right)+y_{r-\mathrm{1}}\right)\hfill \\ {\sigma }_{r+\mathrm{1}}:=\left(\left({\displaystyle \frac{\mathrm{\Delta }x+h}{\mathrm{2}\mathrm{\Delta }x}}\right)\left(\mathrm{3}{\sigma }_r-\mathrm{4}{\sigma }_{r-\mathrm{1}}+{\sigma }_{r-\mathrm{2}}\right)+{\sigma }_{r-\mathrm{1}}\right)\hfill \\ u_{r+\mathrm{1}}:=\left(\left({\displaystyle \frac{\mathrm{\Delta }x+h}{\mathrm{2}\mathrm{\Delta }x}}\right)\left(\mathrm{3}{u}_r-\mathrm{4}{u}_{r-\mathrm{1}}+u_{r-\mathrm{2}}\right)+u_{r-\mathrm{1}}\right)\hfill \end{array}

Examples of normal stresses and stability analysis

The plane surface where failure occurs is known, the number of slices in which the sliding mass is divided is varied, and the stresses for each slice, SF, and the geometrical center are determined. Therefore, examples considering pore pressures, surcharge, and seism can be provided.

In this case, the slope of the curve of planar failure passing through the base of the slope is not known, and there is a mobile boundary on the surface of the slope. The most plane or unfavorable planar failure, minimum safety factor (SF_min), and slope of the planar failure (α critical) are determined. In addition, the normal stress for each slice and its geometric center are obtained.

Example of stress and stability analysis varying the number of slices

A simple slope is used as an example to calculate the stresses and stability for several slices. The dimensions and values of different parameters of the slope are as follows: height of the slope H = 6 m, angle that forms the free face of the slope with respect to the horizontal b = 90°, cohesion C = 20 kN/m², internal friction angle f = 22.6°, unit weight g = 20 kN/m³, failure angle α = 56.3°, horizontal seismic coefficient k_h = 0, vertical seismic coefficient k_v = 0, and pore pressure u = 0. The sliding mass is divided into four, six, eight, and ten slices, respectively.

For the sliding mass divided into four (Fig. 3(a)), six and eight (Fig. 3(b)), and ten slices, the stresses and stability are determined using the proposed programming for nonlinear equations obtained from variational calculus and the proposed subroutines of this research through Matlab, ESS, and Maple to compare results with the axis rotation method [²⁸] and the limit equilibrium method (Slide program).

Table 2 summarizes the normal stresses acting on the failure surface for each slice (s_n). SF is equal to 0.9998, and the geometric center values x_c and y_c, are different for the four slices. Table 3 lists the normal stresses for the four, six, eight, and ten slices obtained using only the Maple program. Table 4 presents the results obtained using the analytical method proposed by Ref. [²⁸].

Example of pore pressures, surcharge, seism, and stability analysis

This example is adopted from the book of the anchor manual [²⁶], in which it is desired to calculate the SF for rock excavation based on its geometric characteristics and resistant parameters. The dimensions and values of the different parameters of the slope are as follows: height of the slope H = 30 m, angle that forms the free face of the slope with respect to the horizontal b = 76°, cohesion C = 295 kN/m², internal friction angle f = 30°, unit weight g = 24 kN/m³, failure angle α = 45°, horizontal seismic coefficient k_h = 0.2, vertical seismic coefficient k_v = 0.1, water level height H₁ = 20 m, and surcharge q_y = 300 kN/m² .

The results obtained by applying Maple programming are listed in Table 5, and it passes through the slope toe whose equation is, y = xtan45° (Fig. 4). Moreover, Table 5 summarizes the distribution of normal stresses for the four slices; SF is equal to 1.156 and the geometric center values x_c and y_c, are different.

These results are practically the same as those in the anchors manual [²⁶].

Example of critical plane surface, analysis of stability and stresses

To determine which planar failure is the most unfavorable or critical, the minimum SF and normal stresses for each slice of the simple slope presented in [²¹] are revised. Here, the geometry and properties of the slope are as follows: height of the slope H= 11.54 m, angle that forms the free face of the slope with respect to the horizontal b = 30°, cohesion C=19.6 kN/m², internal friction angle f=10°, unit weight g=16.7 kN/m³, horizontal seismic coefficient k_h=0, vertical seismic coefficient k_v=0, and pore pressure u=0. The sliding mass is divided into four slices.

The Maple program is used to determine which planar failure surface is the most unfavorable or critical. The programming and subroutines of the Maple software are employed to investigate Eqs. (2)–(4), corresponding to the sum of the horizontal and vertical forces and the moments, and Eq. (6), corresponding to the Euler equation. To determine the slope of the most unfavorable planar failure surface, the transversality condition of the slope crest shown in Eq. (8) is used. The stress and the y coordinate are determined using the crest of the slope, which is obtained by applying Eq. (11) as follows:

(12){\sigma }_{n+\mathrm{1}}=\left({\displaystyle \frac{{\displaystyle \frac{H}{\omega }}-{\displaystyle \frac{\mathrm{5}}{\mathrm{2}}}\mathrm{\Delta }x}{\mathrm{2}\mathrm{\Delta }x}}\right)\left(\mathrm{3}{\sigma }_n-\mathrm{4}{\sigma }_{n-\mathrm{1}}+{\sigma }_{n-\mathrm{2}}\right)+{\sigma }_{n-\mathrm{1}},

(13)y_{n+\mathrm{1}}=\left({\displaystyle \frac{{\displaystyle \frac{H}{\omega }}-{\displaystyle \frac{\mathrm{5}}{\mathrm{2}}}\mathrm{\Delta }x}{\mathrm{2}\mathrm{\Delta }x}}\right)\left(\mathrm{3}{y}_n-\mathrm{4}{y}_{n-\mathrm{1}}+y_{n-\mathrm{2}}\right)+y_{n-\mathrm{1}},

where x_n+1 = H/ω is the ordinate value on the surface or crest of the slope, and ω = tanα is the slope of the plane surface of failure.

The solution is obtained by applying numerical methods for the system of nonlinear equations using the Maple program and iterating the x coordinate on the crest of the slope. From the aforementioned procedure, we can obtain the following.

1) The slope of the worst plane surface failure, i.e., ω = tanα = 0.285; thus, α = 15.9° (Fig. 5(a)).

2) The equation of the plane surface of failure that corresponds to the equation of a line passing through the foot of the slope (y = 0.285x).

3) The minimum SF, which is 1.95 (Fig. 5(b)).

4) The coordinate on the surface or slope top (x₅) (Fig. 5(c)), the width of the slice (Dx), and the normal stresses for each slice and the geometric center (x_c, y_c). The results are presented in Table 6.

Discussions

The results obtained using different methods for the normal stress and the SF are practically the same as those obtained using different mathematical programs (Table 2). However, when comparing the results from the calculus of variations and from the analytical method of the rotation of axis, it was found that the values for the normal stresses and the SF are practically equal (Tables 3 and 4). However, with respect to the geometric center values (x_c, y_c), a diverse set emerges when the system of nonlinear equations is solved using different programs and the number of slices or segments is varied (Tables 2 and 3). This discrepancy enables the relationship between x_c and y_c to be obtained. It should be noted that these values are related to a straight line with a constant slope (y_c/x_c), that is, the relation (−1/λ₁) is the slope of the line that results in a constant value. The line indicates all the probable values ​​that can reach x_c and y_c, which are related through the linear equation c = −2.583x_c + 2.134. Furthermore, the values of the geometric center of slope stability obtained in the present study can aid in the investigation of other types of failure, such as logarithmic spiral failure.

The critical planar failure surface determined in the present study was compared with the analytical method for SF calculation proposed in Ref. [²⁶]. It can be observed that in the present study, the SF of the planar failure is minimized, that is, the equation of the partial derivative of the SF is applied in relation to the α angle formed by the failure surface with respect to the horizontal, (∂SF/∂α) = 0. The observed values indicate that the minimum SF and the critical α are practically the same when the proposed calculus of variations technique is applied.

Conclusions

While applying the technique of calculus of variations in the stability of slopes, considering the mechanism of planar failure, some integrals are generated, which are solved using numerical methods. Simultaneously, there is a system of nonlinear equations that uses different programs and subroutines to identify and study the behavior of the aforementioned equations by varying the number of slices into which the sliding mass is divided. Slopes with simple geometry were revised and used with the proposed approach for comparison with traditional analytical methods, thereby allowing the proposed subroutines to be checked and providing the required validation.

With respect to the normal stresses (s_n) obtained using variational calculus for the different programs and the analytical method, we can conclude that the values were practically equal. These results encourage the practical application of the proposed programming approach.

The proposed programming and subroutines are innovative because they can be applied to slopes with more complicated geometries, such as parabolic, circular, and logarithmic spiral types, while considering the influence of surcharge, interstitial pressure, and pseudo-static load.

References

[1] Coggan J S, Stead D, Eyre J M. Evaluation of techniques for quarry slope stability assessment. Transactions of the Institution of Mining and Metallurgy. Section B. Applied Earth Science, 1998, 107: 139–147

[2] Shnorhokian S, MacNeil B, Mitri H. Volumetric analysis of rock mass instability around haulage drifts in underground mines. Journal of Rock Mechanics and Geotechnical Engineering, 2018, 10(1): 60–71

[3] Wyllie D C, Mah C. Rock Slope Engineering. 4th ed. New York: CRC Press, 2004

[4] Chen W F, Giger M W, Fang H Y. On the limit analysis of stability of slopes. Soil and Foundation, 1969, 9(4): 23–32

[5] Chen W F. Limit Analysis and Soil Plasticity. Amsterdam: Elsevier, 1975

[6] Chen J, Yang Z, Hu R, Zhang H. Study on the seismic active earth pressure by variational limit equilibrium method. Shock and Vibration, 2016, 2016: 1–10

[7] Terzaghi K. Mechanism of landslides. Application of geology to engineering. Geological Society of America, 1950, 83–123

[8] Baker R, Shukha R, Operstein V, Frydman S. Stability charts for pseudo-static slope stability analysis. Soil Dynamics and Earthquake Engineering, 2006, 26(9): 813–823

[9] Nouri H, Fakher A, Jones C. Evaluating the effects of the magnitude and amplification of pseudo-static acceleration on reinforced soil slopes and walls using the limit equilibrium horizontal slices method. Geotextiles and Geomembranes, 2008, 26(3): 263–278

[10] Zhang Y B, Chen G Q, Zen K, Kasama K, Dong S. Limit analysis of seismic slope stability based on tension-shear failure mechanism. In: The 89th International Symposium on Rock Slope Stability in Open Pit Mining and Civil Engineering. Vancouver, 2016, 77–90

[11] Liu G, Zhuang X, Cui Z. Three-dimensional slope stability analysis using independent cover based numerical manifold and vector method. Engineering Geology, 2017, 225: 83–95

[12] Rabczuk T, Zi G, Bordas S, Nguyen-Xuan H. A simple and robust three-dimensional cracking-particle method without enrichment. Computer Methods in Applied Mechanics and Engineering, 2010, 199(37–40): 2437–2455

[13] Ren H, Zhuang X, Cai Y, Rabczuk T. Dual-horizon peridynamics. International Journal for Numerical Methods in Engineering, 2016, 108(12): 1451–1476

[14] Ren H, Zhuang X, Rabczuk T. Dual-horizon peridynamics: A stable solution to varying horizons. Computer Methods in Applied Mechanics and Engineering, 2017, 318: 762–782

[15] Zhou S, Zhuang X, Zhu H, Rabczuk T. Phase field modelling of crack propagation, branching and coalescence in rocks. Theoretical and Applied Fracture Mechanics, 2018, 96: 174–192

[16] Zhou S, Zhuang X, Rabczuk T. Phase field modeling of brittle compressive-shear fractures in rock-like materials: A new driving force and a hybrid formulation. Computer Methods in Applied Mechanics and Engineering, 2019, 355: 729–752

[17] Zhou S, Zhuang X, Rabczuk T. Phase-field modeling of fluid-driven dynamic cracking in porous media. Computer Methods in Applied Mechanics and Engineering, 2019, 350: 169–198

[18] Kopácsy J. Distribution of stresses at failure shape of the sliding surface and theoretical slope height. In: Proceedings of the fifth International Conference on Soil Mechanics and Foundation Engineering. Paris: Dunod, 1961, 641–650

[19] Baker R, Garber M. Discussion of ‘On slip surface and slope stability analysis.’ by Chen W F, Snitbhan N. Soils and Foundations, 1977, 17(1): 65–68

[20] Baker R, Garber M. Theoretical analysis of the stability of slopes. Geotechnique, 1978, 28(4): 395–411

[21] Castillo E, Revilla J. An application of the calculus of variations to the slope stability. Soil Mechanics Transport Laboratory Newsletter José Luis Escario, 1976, 115: 3–23 (in Spanish)

[22] Luceño A, Castillo E. Critical analysis of variational methods applied to slope stability. Soil Mechanics Transport Laboratory Newsletter José Luis Escario, 1980, 45: 3–14 (in Spanish)

[23] Rojas S, Úcar R. Critical analysis of variational methods applied to slope stability. Science and Engineering, 2001, 22(1): 30–40

[24] Leshchinsky H, Huang C C. Generalized slope stability analysis: Interpretation, modification and comparation. Journal of Geotechnical Engineering, 1992, 118(10): 1559–1576

[25] Mac-Lennan J. Failure of a slope. Journal of Public Works: Organism Professional Engineers of Roads, Canals and Ports, 2004, 3446: 53–60

[26] Úcar R. Anchor manual in Civil Engineering. Madrid: Arias Montano Graphics’s, 2004 (in Spanish)

[27] Baker R. Variational slope stability analysis of materials with non-linear failure criterion. Electronic Journal of Geotechnical Engineering, 2005, 10: 1–22

[28] Belandría N, Úcar R, Bongiorno F. Determination of mathematical expressions for calculating the stresses applied to slope stability. Sciencie and Engineering, 2011, 32(3): 115–122 (in Spanish)

[29] Puig Adam P. Theoretical-practical Course of Differential Equations Applied to Physics and Technique. Madrid: Roberto Puig Alvarez, 1975 (in Spanish)

[30] Elsgoltz L. Differential Equations and Variational Calculus. Moscu: MIR, 1969 (in Spanish)

RIGHTS & PERMISSIONS

Higher Education Press