Three Dimensional Slope Stability Analysis of Open Pit Mine

Abstract

The 3-dimensional slope stability analysis has been developing rapidly since the last decade, and currently a number of geomechanical researchers in the world have put forward ideas for optimization of slope design related to the economics and safety of mining operations. The 3-dimensional slope stability analysis methods has answered the assumption of spatial parameters in determining safety factors and the failure probability, thus the volume of failed material and the location of the most critical slopes can be determined. This chapter discusses two methods of 3-dimensional slope stability analysis, namely the limit equilibrium method (LEM) and finite element method (FEM). LEM 3D requires an assumption of failure type with the variable of analysis are the maximum number of columns, the amount of grid points, increment radius, and type of slip surface. On the other hand, FEM 3D requires an assumption of convergence type, absolute force and energy, with the variable of analysis are mesh type and maximum number of iterations. LEM 3D shows that the cuckoo algorithm is reliable in obtaining position and shape of slip surface. Meanwhile FEM 3D, the optimum iteration number needs to be considered to improve analysis efficiency and preserving accuracy.

1. Introduction

This chapter will mainly discuss about 3D slope stability analysis using Limit Equilibrium Method (LEM) and Finite Element Method (FEM). These two methods are widely developed by academics and have been applied by many mining geotechnical practitioners in slope stability analysis. In recent time most of the analysis are performed in 2D method because of its simplicity and lower operational cost [1]. However, 3D analysis is more justifiable to represent the actual geometry condition. Thus, to obtain 3D slope stability analysis has its own importance. In this regard, 2D and 3D analysis can both be performed for the same slope stability analysis problem to obtain a more convincing and realistic results. These two methods can be used to validate one to another [2]. The analysis of the 3D approach is carried out by [3, 4, 5] in the 1970’s. Ref. [5] modifying the slice model used in 2D analysis to become a column in 3D analysis. The study of 3D slope stability analysis model is further developed and evaluated by [2, 6, 7] states that 3D FoS value is always greater than 2D analysis, therefore 2D results are considered more conservative [8]. Most existing three-dimensional 3D slope stability analysis methods are based on simple extensions of corresponding two-dimensional 2D methods of analysis and a plane of symmetry or direction of slide is implicitly assumed. 3D asymmetric slope stability models based on extensions of Bishop’s simplified, Janbu’s simplified, and Morgenstern–Price’s methods are developed by [8].

2. Slope stability analysis

The stability of a slope can be determined by 2 criteria of considerations, which is the value of the safety factor (FoS) and also the value of the probability of failure (PoF), these two criteria are used to determine the optimal geometry of the pit opening. In order to obtain accurate analysis results, the information data regarding the geotechnical conditions of the model must be repetitive to the actual conditions. In general, the principal of the value of safety factor concept is the ratio between the shear strength along the slip surface required to maintain the slope at a stable condition, and the available shear strength [9]. The above definition can be mathematically described as:

$$\mathrm{FoS}=\frac{{\mathrm{S}}_{\mathrm{u}}}{{\mathrm{\tau }}_{\text{required}}}$$E1

The slope is assumed to be a model of an inclined plane [9]. By determining the resultant overall forces and moments acting in equilibrium state, the slope factor of safety value can be determined by comparing the amount of the resisting force to the driving force, or the resisting moment to overturning moment can see in Figure 1.

Figure 1.

Equilibrium force and moment in inclined plane [9].

$$\mathrm{FoS}=\frac{\text{Resisting force}}{\text{Mobilized force}}$$E2

$$\mathrm{FoS}=\frac{\text{Resisting moment}}{\text{Overturning moment}}$$E3

Based on the Mohr-Coulomb failure criteria, the definition of the value of the safety factor can be determined as follows:

$$\mathrm{FoS}=\frac{\mathrm{c}+\text{σtan}\mathrm{\varphi }}{{\mathrm{\tau }}_{\text{required}}}$$E4

Where c and ϕ are effective cohesion and internal friction angle.

The probabilistic failure is an approach that consider various input parameters that generates different Factor of Safety (FoS) values [10]. This information is based on the fact that every random input parameter has the same probability to yield a certain value of FoS. Regarding the difficulty and high expense of field and laboratory data collecting, this method is more attractive because of its representativeness. Figure 2 presents the concept of failure probability and the amount of uncertainty. Slope PoF is determined from the ratio between the area under the curve of the distribution of FoS <1 value to the distribution of FoS ≥1 value. The greater the range of distribution of FoS values, the higher the uncertainty of FoS values with the same PoF values.

Figure 2.

Concept of probabilistic failure.

By definition there is a linear relationship between the PoF value and the likelihood of failure, while this does not apply to the FoS relationship with the chance of failure [10]. A large FoS does not represent a more stable slope, because the implicit uncertainty is not captured by the FoS value. Slopes with FoS of 3 do not mean that they are 2 times more stable than FoS of 1.5, while slopes with a PoF value of 5% are 2 times more stable than slopes with a PoF value of 10%. Slope stability in general performed in two-dimensional analysis. But in modeling complex geometries, 2D analysis cannot simulate them properly. Therefore, the 3D analysis is considered to be able to describe the conditions in the field better than the 2D analysis. Analysis of slope stability with 3D limit equilibrium method starts by assuming the geometry of the sliding mass (Figure 3).

Figure 3.

Comparison between 3D and 2D single slope analysis [11].

The results of the calculation of slope stability can be expressed in safety factor. In this method, safety factor is not only influenced by the direction sliding, but also by the slip surface that safety factor is sensitive to critical slip surface locations. Therefore, the determination of a critical slip surface is very important. Safety factor can be obtained correctly if the determination of critical slip surface is accurate.

2.1 Slope design

Optimal slope geometry is obtained from a step-by-step assessment process [12] state there are 5 stages of the process, which is models, domains, design, analysis and implementation. The initial stage of the geotechnical model is determined by 4 parameters, namely the geological model, structural conditions, rock mass and hydrogeological model. At the domains stage, the failure modes are determined by 2 parameters, namely the strength of the material and the condition of the structure. A single slope design arrangement is determined by the regulations or standards used by the company and the capabilities of the equipment used. Determination of the haulage road width and also the overall slope angle is based on mine planning related to the economic aspects of the opening geometry made. Furthermore, the stability analysis of the slope geometry that has been designed refers to the parameters (structure, strength, groundwater, and in-situ stress). After the final design is obtained, a risk assessment is carried out to mitigate the potential for landslides that may occur. In the implementation stage, the functions of dewatering, blasting and monitoring of the progress of the design model and the movement of rock masses (Figure 4).

Figure 4.

Slope design process [12].

2.2 Limit equilibrium method 3D

These days the needs and pressure to analyze a slope 3 dimensionally is more sounds. This is because 2D analysis assumes that the width of slope is infinitely wide so then it neglects 3D effect [11]. In most of the cases the width to height ratio is not sufficiently long and varies perpendicular to the slide movement. Therefore, 3D analysis is considered important to be done to produce the representative FoS. Moreover, in 3D analysis the volume of failure can also be estimated while 2D analysis cannot. If the volume can be determinate, it can be useful as one of the considerations in giving failure prevention recommendation.

2.2.1 Safety factor for 3D slope stability

The 3-dimensional model is a refined version of the 2-dimensional by projecting the skid plane into a column and determining the resultant force, as well as the moment based on the x, y, and z directions. The equilibrium force and moments acting on the overall column mass are used to determine the following 3 possible direction of the slip plane:

1. The column moves in the same direction

2. The column moves towards one another

3. The column moves in the opposite direction

For the 3-dimensional analysis, the mass potential of the slip plane is divided into several columns. Ref. [8] give the equation of the Simplified Janbu method deduced from the Morgenstern-Price method to obtain a safety factor value of 3-dimensional analysis (Figure 5).

Figure 5.

3 dimensional column [8].

a_i is space angle for sliding direction with respect to the projected x – y plane, a_x, a_y are base inclination along x and y directions measure at the center of each column, E_xi, E_yi are inter-column normal forces in x and y directions, respectively, H_xi, H_yi are lateral inter-column shear forces in x and y directions, N’_i, U_i are effective normal and base pore watery force, P_vi, S_i is vertical external force, and base mobilized shear force, and X_xi, X_yi are vertical inter-column shear force in plane perpendicular to x and y directions.

With the mohr-coulomb collapse criteria, the safety factor is determined using the following equation:

$$\mathrm{F}=\frac{{\mathrm{S}}_{\mathrm{fi}}}{{\mathrm{S}}_{\mathrm{i}}}=\frac{{\mathrm{C}}_{\mathrm{i}}+\mathrm{N}{\prime }_{\mathrm{i}}tan\mathrm{\varphi }{\prime }_{\mathrm{i}}}{{\mathrm{S}}_{\mathrm{i}}}$$E5

where S_fi is ultimate resultant shear force available at the base of column i, N’_i is the effective base normal force, C_i is (c. Ai) and c and A_i are effective cohesive strength and the base area of the column. The base shear force S_i and normal base force N_i are expressed as the components of forces with respect to x, y, and z directions for column i.

$${\mathrm{S}}_{\mathrm{xi}}={\mathrm{f}}_1{\mathrm{S}}_{\mathrm{i}},{\mathrm{S}}_{\mathrm{yi}}={\mathrm{f}}_2{\mathrm{S}}_{\mathrm{i}},{\mathrm{S}}_{\mathrm{zi}}={\mathrm{f}}_3{\mathrm{S}}_{\mathrm{i}}$$E6

$${\mathrm{N}}_{\mathrm{xi}}={\mathrm{g}}_1{\mathrm{N}}_{\mathrm{i}},{\mathrm{N}}_{\mathrm{yi}}={\mathrm{g}}_2{\mathrm{N}}_{\mathrm{i}},{\mathrm{N}}_{\mathrm{zi}}={\mathrm{g}}_3{\mathrm{N}}_{\mathrm{i}}$$E7

where f₁, f₂, f₃ and g₁, g₂, g₃ = unit vectors in the direction of S_i and N_i. The projected shear angles a’ = same for all columns in the x – y plane in the present formulation, and by using this angle, the space shear angle ai found for each column.

$${\mathrm{a}}_{\mathrm{i}}={tan}^{-1}\left\{sin{\mathrm{\theta }}_{\mathrm{i}}/\left[cos{\mathrm{\theta }}_{\mathrm{i}}+\left(cos{\mathrm{a}}_{\mathrm{yi}}/tan\mathrm{a}\prime cos{\mathrm{a}}_{\mathrm{xi}}\right)\right]\right\}$$E8

$${\mathrm{\theta }}_{\mathrm{i}}={cos}^{-1}\left\{sin{\mathrm{a}}_{\mathrm{xi}}\bullet sin{\mathrm{a}}_{\mathrm{yi}}\right\}$$E9

An arbitrary intercolumn shear force function f (x, y) is assumed in the present analysis, and the relationships between the intercolumn shear and normal forces in the x- and y-directions are given as:

$${\mathrm{Xx}}_{\mathrm{i}}={\mathrm{Ex}}_{\mathrm{i}}\mathrm{f}\left(\mathrm{x},\mathrm{y}\right){\mathrm{\lambda }}_{\mathrm{x}}$$E10

$${\mathrm{Xy}}_{\mathrm{i}}={\mathrm{Ey}}_{\mathrm{i}}\mathrm{f}\left(\mathrm{x},\mathrm{y}\right){\mathrm{\lambda }}_{\mathrm{y}}$$E11

$${\mathrm{Hx}}_{\mathrm{i}}={\mathrm{Ey}}_{\mathrm{i}}\mathrm{f}\left(\mathrm{x},\mathrm{y}\right){\mathrm{\lambda }}_{\mathrm{xy}}$$E12

$${\mathrm{Hy}}_{\mathrm{i}}={\mathrm{Ex}}_{\mathrm{i}}\mathrm{f}\left(\mathrm{x},\mathrm{y}\right){\mathrm{\lambda }}_{\mathrm{yx}}$$E13

Where λ_x and λ_y are intercolumn shear force mobilization factors in x and y directions, respectively, and λ_xy and λ_yx are intercolumn shear force mobilization factors in xz and yz planes. Considering the vertical and horizontal force equilibrium for the ith column in the z, x, and y directions produces the following equations (Figure 6):

Figure 6.

Force equilibrium in columns [8].

$${\mathrm{F}}_{\mathrm{z}}=0={\mathrm{N}}_{\mathrm{i}}{\mathrm{g}}_{3\mathrm{i}}+{\mathrm{S}}_{\mathrm{i}}{\mathrm{f}}_{3\mathrm{i}}–\left({\mathrm{W}}_{\mathrm{i}}+{\mathrm{P}}_{\mathrm{vi}}\right)=\left({\mathrm{X}}_{\mathrm{xi}+1}–{\mathrm{X}}_{\mathrm{xi}}\right)+\left({\mathrm{X}}_{\mathrm{yi}+1}–{\mathrm{X}}_{\mathrm{yi}}\right)$$E14

$${\mathrm{F}}_{\mathrm{x}}=0={\mathrm{S}}_{\mathrm{i}}{\mathrm{f}}_{1\mathrm{i}}-{\mathrm{N}}_{\mathrm{i}}{\mathrm{g}}_{1\mathrm{i}}-{\mathrm{H}}_{\mathrm{xi}}+{\mathrm{H}}_{\mathrm{xi}+1}={\mathrm{E}}_{\mathrm{xi}+1}–{\mathrm{E}}_{\mathrm{xi}}$$E15

$${\mathrm{F}}_{\mathrm{y}}=0={\mathrm{S}}_{\mathrm{i}}{\mathrm{f}}_{2\mathrm{i}}-{\mathrm{N}}_{\mathrm{i}}{\mathrm{g}}_{2\mathrm{i}}-{\mathrm{H}}_{\mathrm{yi}}+{\mathrm{H}}_{\mathrm{yi}+1}={\mathrm{E}}_{\mathrm{yi}+1}–{\mathrm{E}}_{\mathrm{yi}}$$E16

Solving Equation, the base normal and shear forces can be expressed as

$${\mathrm{N}}_{\mathrm{i}}={\mathrm{A}}_{\mathrm{i}}+{\mathrm{B}}_{\mathrm{i}}{\mathrm{S}}_{\mathrm{i}}$$E17

$${\mathrm{S}}_{\mathrm{i}}=\frac{{\mathrm{C}}_{\mathrm{i}}+\left({\mathrm{A}}_{\mathrm{i}}-{\mathrm{U}}_{\mathrm{i}}\right)tan\mathrm{\varphi }{\prime }_{\mathrm{i}}}{\mathrm{F}\left[1-\left(\frac{{\mathrm{B}}_{\mathrm{i}}tan{\mathrm{\varphi }}_{\mathrm{i}}}{\mathrm{F}}\right)\right]}$$E18

$${\mathrm{A}}_{\mathrm{i}}=\frac{{\mathrm{W}}_{\mathrm{i}}+({\mathrm{P}}_{\mathrm{vi}}+\mathrm{\Delta }{\mathrm{Ex}}_{\mathrm{i}}{\mathrm{\lambda }}_{\mathrm{x}}+\mathrm{\Delta }{\mathrm{Ey}}_{\mathrm{i}}{\mathrm{\lambda }}_{\mathrm{y}})}{{\mathrm{g}}_{3\mathrm{i}}}$$E19

$${\mathrm{B}}_{\mathrm{i}}=-\frac{{\mathrm{f}}_{3\mathrm{i}}}{{\mathrm{g}}_{3\mathrm{i}}}$$E20

$$\mathrm{F}=\frac{{\mathrm{C}}_{\mathrm{i}}+\mathrm{N}{\prime }_{\mathrm{i}}tan\mathrm{\varphi }{\prime }_{\mathrm{i}}}{{\mathrm{S}}_{\mathrm{i}}}$$E21

Considering the overall force equilibrium in x-direction internal force E cancels out.

$$-\sum {\mathrm{H}}_{\mathrm{xi}}+\sum {\mathrm{N}}_{\mathrm{i}}{\mathrm{g}}_{1\mathrm{i}}–\sum {\mathrm{S}}_{\mathrm{i}}{\mathrm{f}}_{1\mathrm{i}}=0$$E22

Considering overall moment equilibrium in the x-direction

$$-\sum \left({\mathrm{W}}_{\mathrm{i}}+{\mathrm{P}}_{\mathrm{vi}}-{\mathrm{N}}_{\mathrm{i}}{\mathrm{g}}_{3\mathrm{i}}–{\mathrm{S}}_{\mathrm{i}}{\mathrm{f}}_{3\mathrm{i}}\right)\mathrm{RX}+\sum \left({\mathrm{N}}_{\mathrm{i}}{\mathrm{g}}_{1\mathrm{i}}-{\mathrm{S}}_{\mathrm{i}}{\mathrm{f}}_{1\mathrm{i}}\right)\mathrm{RZ}=0$$E23

Considering overall force equilibrium in the y-direction

$$-\sum {\mathrm{H}}_{\mathrm{yi}}+\sum {\mathrm{N}}_{\mathrm{i}}{\mathrm{g}}_{2\mathrm{i}}–\sum {\mathrm{S}}_{\mathrm{i}}{\mathrm{f}}_{2\mathrm{i}}=0$$E24

Considering overall moment equilibrium in the y-direction

$$-\sum \left({\mathrm{W}}_{\mathrm{i}}+{\mathrm{P}}_{\mathrm{vi}}-{\mathrm{N}}_{\mathrm{i}}{\mathrm{g}}_{3\mathrm{i}}–{\mathrm{S}}_{\mathrm{i}}{\mathrm{f}}_{3\mathrm{i}}\right)\mathrm{RY}+\sum \left({\mathrm{N}}_{\mathrm{i}}{\mathrm{g}}_{2\mathrm{i}}-{\mathrm{S}}_{\mathrm{i}}{\mathrm{f}}_{2\mathrm{i}}\right)\mathrm{RZ}=0$$E25

The directional safety factor F_x and F_y is determined as follows:

$${\mathrm{F}}_{\mathrm{x}}=\frac{\sum \left[{\mathrm{C}}_{\mathrm{i}}+\left({\mathrm{N}}_{\mathrm{i}}-{\mathrm{U}}_{\mathrm{i}}\right)tan{\mathrm{\varphi }}_{\mathrm{i}}\right]{\mathrm{f}}_{1\mathrm{i}}}{\sum {\mathrm{N}}_{\mathrm{i}}{\mathrm{g}}_{1\mathrm{i}}-\sum {\mathrm{H}}_{\mathrm{xi}}},\left(0<{\mathrm{F}}_{\mathrm{x}}<\infty \right)$$E26

$${\mathrm{F}}_{\mathrm{y}}=\frac{\sum \left[{\mathrm{C}}_{\mathrm{i}}+\left({\mathrm{N}}_{\mathrm{i}-}{\mathrm{U}}_{\mathrm{i}}\right)tan\mathrm{\varphi }{\prime }_{\mathrm{i}}\right]{\mathrm{f}}_{2\mathrm{i}}}{\sum {\mathrm{N}}_{\mathrm{i}}{\mathrm{g}}_{2\mathrm{i}}-\sum {\mathrm{H}}_{\mathrm{yi}}},\left(0<{\mathrm{F}}_{\mathrm{y}}<\infty \right)$$E27

Formulation 3D Bishop’s methods by considering the overall moment equilibrium equations in x or y direction and neglecting the inter-column vertical and horizontal shear forces.

$${\mathrm{F}}_{\mathrm{mx}}=\frac{\sum \left\{{\mathrm{K}}_{\mathrm{xi}}\left[{\mathrm{f}}_{2\mathrm{i}}{\mathrm{RZ}}_{\mathrm{i}}+{\mathrm{f}}_{3\mathrm{i}}{\mathrm{RY}}_{\mathrm{i}}\right]\right\}}{\sum \left({\mathrm{W}}_{\mathrm{i}}+{\mathrm{P}}_{\mathrm{vi}}\right){\mathrm{RY}}_{\mathrm{i}}+\sum {\mathrm{N}}_{\mathrm{i}}({\mathrm{g}}_{2\mathrm{i}}{\mathrm{RZ}}_{\mathrm{i}}-\left({\mathrm{g}}_{3\mathrm{i}}{\mathrm{RY}}_{\mathrm{i}}\right)}$$E28

$${\mathrm{F}}_{\mathrm{my}}=\frac{\sum \left\{{\mathrm{K}}_{\mathrm{yi}}\left[{\mathrm{f}}_{1\mathrm{i}}{\mathrm{RZ}}_{\mathrm{i}}+{\mathrm{f}}_{3\mathrm{i}}{\mathrm{RX}}_{\mathrm{i}}\right]\right\}}{\sum \left({\mathrm{W}}_{\mathrm{i}}+{\mathrm{P}}_{\mathrm{vi}}\right){\mathrm{RX}}_{\mathrm{i}}+\sum {\mathrm{N}}_{\mathrm{i}}({\mathrm{g}}_{1\mathrm{i}}{\mathrm{RZ}}_{\mathrm{i}}-\left({\mathrm{g}}_{3\mathrm{i}}{\mathrm{RX}}_{\mathrm{i}}\right)}$$E29

$${\mathrm{K}}_{\mathrm{xi}}=\frac{\left\{{\mathrm{C}}_{\mathrm{i}}+\left[({\mathrm{W}}_{\mathrm{i}}+{\mathrm{P}}_{\mathrm{vi}}\right)/\left({\mathrm{g}}_{3\mathrm{i}}-{\mathrm{U}}_{\mathrm{i}}\right)\right]tan{\mathrm{\varphi }}_{\mathrm{i}}\}}{1+\left({\mathrm{f}}_{3\mathrm{i}}tan{\mathrm{\varphi }}_{\mathrm{i}}/{\mathrm{g}}_{3\mathrm{i}}{\mathrm{F}}_{\mathrm{mx}}\right)}$$E30

$${\mathrm{K}}_{\mathrm{yi}}=\frac{\left\{{\mathrm{C}}_{\mathrm{i}}+\left[({\mathrm{W}}_{\mathrm{i}}+{\mathrm{P}}_{\mathrm{vi}}\right)/\left({\mathrm{g}}_{3\mathrm{i}}-{\mathrm{U}}_{\mathrm{i}}\right)\right]tan{\mathrm{\varphi }}_{\mathrm{i}}\}}{1+\left({\mathrm{f}}_{3\mathrm{i}}tan{\mathrm{\varphi }}_{\mathrm{i}}/{\mathrm{g}}_{3\mathrm{i}}{\mathrm{F}}_{\mathrm{my}}\right)}$$E31

Considering overall moment equilibrium about an axis passing through (x₀, y₀, z₀) and parallel to the z axis gives:

$$\sum \left(-{\mathrm{N}}_{\mathrm{i}}{\mathrm{g}}_{1\mathrm{i}}–{\mathrm{S}}_{\mathrm{i}}{\mathrm{f}}_{1\mathrm{i}}\right)\mathrm{RY}+\sum \left({\mathrm{N}}_{\mathrm{i}}{\mathrm{g}}_{2\mathrm{i}}-{\mathrm{S}}_{\mathrm{i}}{\mathrm{f}}_{2\mathrm{i}}\right)\mathrm{RX}=0$$E32

$${\mathrm{F}}_{\mathrm{mz}}=\frac{\sum \left[{\mathrm{K}}_{\mathrm{zi}}\left({\mathrm{f}}_{2\mathrm{i}}{\mathrm{RX}}_{\mathrm{i}}-{\mathrm{f}}_{3\mathrm{i}}{\mathrm{RY}}_{\mathrm{i}}\right)\right]}{\sum \mathrm{N}({\mathrm{g}}_{2\mathrm{i}}{\mathrm{RX}}_{\mathrm{i}}-\left({\mathrm{g}}_{1\mathrm{i}}{\mathrm{RY}}_{\mathrm{i}}\right)}$$E33

$${\mathrm{K}}_{\mathrm{zi}}=\frac{\left\{{\mathrm{C}}_{\mathrm{i}}+\left[({\mathrm{W}}_{\mathrm{i}}+{\mathrm{P}}_{\mathrm{vi}}\right)/\left({\mathrm{g}}_{3\mathrm{i}}-{\mathrm{U}}_{\mathrm{i}}\right)\right]tan{\mathrm{\varphi }}_{\mathrm{i}}\}}{1+\left({\mathrm{f}}_{3\mathrm{i}}tan{\mathrm{\varphi }}_{\mathrm{i}}/{\mathrm{g}}_{3\mathrm{i}}{\mathrm{F}}_{\mathrm{mz}}\right)}$$E34

For the 3D asymmetric Bishop’s method, at moment equilibrium point, the directional factors of safety, F_mx, F_my, and F_mz are equal to each other. Under this condition, the global factor of safety F_m based on moment can be determined as

$${\mathrm{F}}_{\mathrm{m}}={\mathrm{F}}_{\mathrm{mx}}={\mathrm{F}}_{\mathrm{my}}={\mathrm{F}}_{\mathrm{mz}}$$E35

Formulation 3D simplified Janbu’s methods by considering the overall force equilibrium equations and neglecting the inter-column vertical and horizontal shear forces.

$${\mathrm{A}}_{\mathrm{xi}}=\frac{\left\{{\mathrm{C}}_{\mathrm{i}}+\left[\left({\mathrm{W}}_{\mathrm{i}-}{\mathrm{P}}_{\mathrm{vi}}\right)/\left({\mathrm{g}}_{3\mathrm{i}-}{\mathrm{U}}_{\mathrm{i}}\right)\right]tan{\mathrm{\varphi }}_{\mathrm{i}}\right\}}{1+\left({\mathrm{f}}_{3\mathrm{i}}tan{\mathrm{\varphi }}_{\mathrm{i}}/{\mathrm{g}}_{3\mathrm{i}}{\mathrm{F}}_{\mathrm{sx}}\right)}$$E36

$${\mathrm{A}}_{\mathrm{yi}}=\frac{\left\{{\mathrm{C}}_{\mathrm{i}}+\left[\left({\mathrm{W}}_{\mathrm{i}-}{\mathrm{P}}_{\mathrm{vi}}\right)/\left({\mathrm{g}}_{3\mathrm{i}-}{\mathrm{U}}_{\mathrm{i}}\right)\right]tan{\mathrm{\varphi }}_{\mathrm{i}}\right\}}{1+\left({\mathrm{f}}_{3\mathrm{i}}tan{\mathrm{\varphi }}_{\mathrm{i}}/{\mathrm{g}}_{3\mathrm{i}}{\mathrm{F}}_{\mathrm{sy}}\right)}$$E37

$${\mathrm{F}}_{\mathrm{sx}}=\frac{\sum {\mathrm{A}}_{\mathrm{xi}}\left({\mathrm{f}}_{1\mathrm{i}}+{\mathrm{f}}_{3\mathrm{i}}{\mathrm{g}}_{1\mathrm{i}}/{\mathrm{g}}_{3\mathrm{i}}\right)}{\sum \left({\mathrm{g}}_{1\mathrm{i}}/{\mathrm{g}}_{3\mathrm{i}}\right)\left({\mathrm{W}}_{\mathrm{i}}+{\mathrm{P}}_{\mathrm{vi}}\right)}$$E38

$${\mathrm{F}}_{\mathrm{sy}}=\frac{\sum {\mathrm{A}}_{\mathrm{yi}}\left({\mathrm{f}}_{2\mathrm{i}}+{\mathrm{f}}_{3\mathrm{i}}{\mathrm{g}}_{2\mathrm{i}}/{\mathrm{g}}_{3\mathrm{i}}\right)}{\sum \left({\mathrm{g}}_{2\mathrm{i}}/{\mathrm{g}}_{3\mathrm{i}}\right)\left({\mathrm{W}}_{\mathrm{i}}+{\mathrm{P}}_{\mathrm{vi}}\right)}$$E39

For 3D asymmetric Janbu’s method, at the force equilibrium point, the directional factors of safety, F_sx, and F_sy are equal to each other. Under this condition, the global factor of safety F_f based on force is determined as follows:

$${\mathrm{F}}_{\mathrm{f}}={\mathrm{F}}_{\mathrm{sx}}={\mathrm{F}}_{\mathrm{sy}}$$E40

The safety factor is also used in vertical and 3D force equilibrium to achieve the simplified Janbu’s method.

2.2.2 Grid search in determination of slip surface

One of the methods that can be used to determine the critical slip surface is the grid search method [11]. In the grid search method, the first thing to do is to determine the size of the grid box with dimensions x, y, and z. After the grid box is available, the user determines the number of grid points that you want to use in the x, y, and z directions. This point serves as the center of rotation. Each center of rotation can produce a number of circles that are used as slip surfaces. The number of circles produced at each center of rotation from the minimum radius to the maximum radius is called the radius increment. Illustration of the number of grid points and radius increment can be seen in Figures 7 and 8.

Figure 7.

Illustration of grid point [11].

Figure 8.

Illustration of the radius increment in the grid search [11].

After assuming the field geomaterial failure, the next step is mass discretization of the sliding mass into a number of columns. Square nets are applied to the sliding mass so that the sliding mass is divided into columns. There are two kinds of columns; the active column where the column is inside the sliding mass boundary line, and the inactive column where these columns are outside the sliding mass boundary line. In the calculation, the inactive columns are ignored so that the discrete sliding mass is determined only from the sum of the active columns. Figure 9 shows the illustration of the discretization of the sliding mass using a square grid.

Figure 9.

Discretization of the sliding mass using a square grid [11].

After discretizing the sliding mass, internal and external forces in each column can be calculated based on moment equilibrium, force equilibrium, or both depending on what method of calculation is used (Figure 10).

Figure 10.

3D model for slope stability analysis [11].

The grid search method is used to find critical slip surface. The grid search method starts by specifying the grid box dimension. The location and dimensions of the grid box must cover the entire study area so that the search for critical slip surface can be performed optimally, for the influence of number of grid point and radius increment in determining safety factor result can be seen in Table 1. The result of 3D slope stability analysis using grid search can see in Figure 11.

Radius increment	Number of grid points	Factor of safety	Volume (m^3)	Location	Direction of sliding	Center of rotation
						X	Y	Z
	20 × 20 × 10	0.868	2,245,470	North-HW	246.7	371,070	9,586,190	470
10	30 × 30 × 15	0.874	1,998,950	North-HW	246.6	371,069	9,586,190	446
	40 × 40 × 20	0.873	2,009,660	North-HW	246.5	371,068	9,586,190	443
	20 × 20 × 10	0.868	2,245,470	North-HW	246.7	371,070	9,586,190	470
20	30 × 30 × 15	0.874	2,028,600	North-HW	246.5	371,069	9,586,200	429
	40 × 40 × 20	0.872	1,934,080	North-HW	246.4	371,069	9,586,200	428
	20 × 20 × 10	0.868	2,245,470	North-HW	246.7	371,070	9,586,190	470
30	30 × 30 × 15	0.906	1,628,780	North-HW	246.2	371,092	9,586,210	352
	40 × 40 × 20	0.867	1,937,390	North-HW	246.1	371,076	9,586,170	459
	20 × 20 × 10	0.868	2,245,470	North-HW	246.7	371,070	9,586,190	470
40	30 × 30 × 15	0.876	1,756,790	North-HW	246.2	371,076	9,586,200	404
	40 × 40 × 20	0.872	1,934,080	North-HW	246.4	371,069	9,586,200	428
	20 × 20 × 10	0.868	2,245,470	North-HW	246.7	371,070	9,586,190	470
50	30 × 30 × 15	0.876	2,068,950	North-HW	246.3	371,077	9,586,180	459
	40 × 40 × 20	0.879	1,968,990	North-HW	246.3	371,078	9,586,190	421

Table 1.

The influence of number of grid point and radius increment in determining safety factor.

Figure 11.

Grid search LEM 3D analysis result [11].

2.2.3 Cuckoo search in determination of slip surface

Grid Search is commonly used as slip surface searching method because the principle is simple and easy to understand [11]. However, this method can only calculate the circular slip surfaces, so it cannot represent the stability of slope in real condition. To make it more representative, non-circular slip surface option is also available in slope stability simulation software, and one of the searching methods in non-circular option is CS method.

CS is an algorithm which is used for solving the optimization problems. CS has been used in engineering field, such as welded beam and spring design optimization but there is still a few that use it for slope stability issue. Nowadays, the necessity to analyze 3D slope stability is more essential. The reason why slope stability problem should not be assumed 2 dimensionally that should be taken into account is the importance of determining volume of failure for risk management volume of failure is counted as one of the consequences, and it can be obtained by analyzing slope stability in 3D, indeed. Furthermore, there is a relationship between failure probability and volume of failure. Another issue exists when 3D analysis is performed using Grid Search on a vast area. The use of 3D analysis on a vast area can be complicated, and in real cases more advanced optimization methods are required. Therefore, CS is tried to be applied in order to determine the slip surface with minimum SF in 3D analysis. Thus, it can be suggested to be used as an alternative or other better option in 3D analysis.

CS means a metaheuristic optimization method that was developed by [13]. This method was inspired from cuckoo’s breeding behavior. In this research, CS is used as a slip surface search tool that has the lowest SF. CS is coupled with Lévy Flights random walk. There are few rules to use this algorithm as follows:

• Each cuckoo lay one egg at a time, and dumps it in a randomly chosen nest;

• The best nests with high quality of eggs (solutions) will carry over to the next generations;

• The number of available host nests is fixed, and a host can discover an alien egg with a probability pa ∈ [0, 1]. In this case, the host bird can either throw the egg away or abandon the nest so as to build a completely new nest in a new location.

The last rule can be approached using a fraction pa to determine the worst solutions of n nests that will be replaced with a new nest randomly. In order to solve the problem, it can be simply illustrated that every egg in a nest represents one new solution. The purpose is to use the new and potentially better solution to replace the current solution in the nest. In a certain condition, the nest may have 2 eggs (2 solutions) but this problem is simplified so one nest has only 1 solution.

The random walk as determinated by Lévy Flights can be described in the following formula:

$${\mathrm{x}}_{\mathrm{i}}^{\left(\mathrm{t}+1\right)}={\mathrm{x}}_{\mathrm{i}}^{\left(\mathrm{t}\right)}+\mathrm{\alpha }⨁\mathrm{L}\mathrm{é}\mathrm{vy}\left(\mathrm{\lambda }\right)$$E41

where α > 0 is the step size and Lévy(λ) is the position function from Lévy Flights (Figure 12).

Figure 12.

Cuckoo search LEM 3D analysis result [11].

$$\mathrm{L}\mathrm{é}\mathrm{vy}\sim \mathrm{u}={\mathrm{t}}^{-\mathrm{\lambda }},\left(1<\mathrm{\lambda }\le 3\right)$$E42

CS has been applied in many optimizations and computer intelligence with promising efficiency, this has been proved from design application in engineering field, scheduling problems, thermodynamic calculations, etc. Few examples of CS application in engineering field are designing spring, welded beam, and steel frame. The CS’s performance has also been compared with some metaheuristic algorithms such as PSO and GA, and the result shows that CS has higher success rate than. The result of the influence of max columns in x or y and max iteration in determining safety factor can be seen in Table 2.

Max iteration	Max columns in X or Y	FoS	Volume (m^3)
40	50	2.06	159,816
	60	2.02	179,523
	70	2.10	229,704
	80	2.10	303,186
	90	2.01	183,631
	100	2.01	187,898
80	50	2.03	179,192
	60	2.02	192,951
	70	2.03	167,961
	80	2.02	178,858
	90	2.02	171,821
	100	2.01	183,469
120	50	2.02	175,113
	60	2.02	180,573
	70	2.01	178,671
	80	2.01	193,606
	90	2.01	180,490
	100	2.02	190,336
160	50	2.00	194,125
	60	2.01	174,180
	70	2.03	170,430
	80	2.01	179,233
	90	1.00	19,008
	100	1.99	199,321
200	50	2.01	180,821
	60	2.01	179,535
	70	2.10	258,302
	80	2.19	270,897
	90	2.04	172,299
	100	1.99	192,171
400	50	2.00	195,539
	60	2.05	174,848
	70	2.02	199,719
	80	2.00	195,231
	90	2.00	194,572
	100	2.00	186,047
600	50	2.05	157,244
	60	2.04	171,805
	70	2.02	178,812
	80	2.01	184,185
	90	2.00	200,779
	100	2.00	195,602
800	50	0.82	38,239
	60	2.01	192,786
	70	2.01	204,299
	80	2.01	196,060
	90	2.01	187,613
	100	2.01	187,709
1000	50	2.02	185,006
	60	2.04	173,577
	70	2.03	175,154
	80	2.01	187,845
	90	2.01	189,787
	100	2.01	177,988

Table 2.

The influence of max columns in X or Y and max iteration in determining safety factor.

2.2.4 Analysis 3D limit equilibrium of open pit mine

3D limit equilibrium analysis method, data regarding 3D slope geometry model, material properties and 3D geological models are required. The 3D slope geometry model required for this method can be obtain from the reconstruction of pit surface model from either terrestrial or photogrammetric measurement methods. The next step is to create an external volume from the pit surface model which will then be analyzed to determine the position and shape of the slip surface, example of 3D slope geometry model can be seen in Figure 13.

Figure 13.

3D slope geometry model.

Limit equilibrium analysis method uses the properties values of materials obtained from laboratory tests results to calculate the value of the safety factor. An example of the input parameter data used in the analysis can be seen in Table 3.

Input parameter	Limonite	Saprolite	Bedrock	Unit
Unit Weight	17.19	15.77	29	γn (kN/m^3)
Cohesion	41.14	38.06	3270	C (KPa)
Internal Friction Angle	15.67	12.58	41.68	ϕ (°)
Modulus Young	20,000	20,000	150,000	ε (Kpa)
Poisson Ratio	0.30	0.40	0.20	ν

Table 3.

Material properties.

Geological modeling is the process of creating visual description geometry of rock lithology into software that represents the actual conditions. However, there are limitations in the modeling process, that related to the limited information on the data held and errors in data interpretation carried out. One of the methods that can be used to create the 3D geological model is interpolation the lithology information data from geotechnical drilling results. The 3D geological model will be used in the 3-dimensional slope stability analysis to determine the distribution characteristics of the rock lithology. In this analysis, the 3D geological floor model of limonite, saprolite and bedrock lithology is used and can be seen in Figure 14.

Figure 14.

3D geological floor model.

The results of the analysis using the Bishop Simplified method for the slip surface search method, both grid search and cuckoo search can be seen in Figures 15 and 16. The analysis results with the grid search show the safety factor value is 1.104, while the cuckoo search is 1.089. The position of the slip surface with the grid search and cuckoo search is the same position it indicates the accuracy of the cuckoo search, while the grid requires a grid box which must represent the actual slip surface position.

Figure 15.

Result 3D analysis using grid search.

Figure 16.

Result 3D analysis using cuckoo search.

2.3 Finite element method 3D

The finite element method has been widely applied by mining geotechnical practitioners in slope stability analysis, with the advantage that the stress–strain analysis in the material allows to determine the displacement and strain values acting on the model elements, but this method has weaknesses in the process. This analysis uses a large number of calculation matrices so that it requires a long computation time, especially if the analysis is carried out in 3D, the number of elements and nodes used will also be high so that the computation time can run for days depending on the number of nodes and types on mesh used. The computation time is closely related to the maximum number of iterations in the calculation of the force error or load imbalance (solid tolerance) in determining the convergence of model, the higher the maximum number of iterations, the calculation in determining the convergence level in the analysis result model will be more accurate, but unfortunately it will take a long time in the analysis.

Slope stability analysis using the finite element method takes into account the stress–strain analysis that works on each element of the model and focuses more on the analysis of the value of the deformation that occurs rather than the level of slope stability [2]. The advantages of analysis using the finite element method when compared to the limit equilibrium method are as follows:

• No need to assume the position of the slip surface, failure occur in zones where the shear strength of the material cannot maintain stability, due to the shear stress that works due to gravity.

• It does not require the concept of slice or columns, so it does not require an approach of forces acting on global equilibrium.

• Analyze the stress–strain so that it can see the deformation and the effective stress.

• The finite element method can monitor the movement of rock masses towards failure.

2.3.1 Strength reduction factor for 3D slope stability

The value of the safety factor (SF) in the finite element method is defined as the ratio of the actual material shear strength to the shear strength of the material when the model failure. This concept is similar to the limit equilibrium method, where the shear strength ratio of a material to its driving force [1]. The strength of the material at failure can be written in the following equation:

$${\mathrm{c}}_{\mathrm{f}}=\frac{\mathrm{c}}{\mathrm{SRF}}$$E43

$${\mathrm{\varphi }}_{\mathrm{f}}=arctan\frac{tan\mathrm{\varphi }}{\mathrm{SRF}}$$E44

With strength reduction factor (SRF) is the value of the factor for the decrease in the strength of the material, c_f is the cohesion value when the model failure and ϕ_f is the value of the internal friction angle when the model failure. In determining the strength of the material at failure, the technique of shear strength reduction (SSR) is used, where the actual material strength parameter is decreased step by step until the model become failure condition (non-convergent) [14]. The value of the material strength reduction factor for the Mohr-Coulomb criteria can be determined as follows:

$$\frac{\mathrm{\tau }}{\mathrm{SRF}}=\frac{\mathrm{c}}{\mathrm{SRF}}+\frac{tan\mathrm{\varphi }}{\mathrm{SRF}}$$E45

$$\frac{\mathrm{\tau }}{\mathrm{SRF}}={\mathrm{c}}_{\mathrm{f}}+tan{\mathrm{\varphi }}_{\mathrm{f}}$$E46

The systematic stage in the model analysis to determine the critical value of the material shear strength reduction factor (Critical SRF) is as follows:

• The first stage is to prepare the model for analysis using the finite element method, determine characteristics of strength and deformation properties of the material.

• The second stage is to increase the value of the material strength reduction factor (SRF) and calculate the material parameters using the Mohr-Coulomb criteria, then input the new material properties data into the model and recalculate and record the maximum total deformation value.

• Repeat the second stage using a systematic increase in the value of the material strength reduction factor, until the model becomes a non-convergent condition (failure). The critical SRF is determined at the highest srf value in the model to achieve the convergent condition.

3D model analysis is an extension of 2D analysis, but in modeling complex geometries 2D analysis cannot simulate them properly. The 2D analysis assumes that the slope width is infinitely wide [11]. However, in many cases the 2D analysis is considered more conservative because it results in a lower safety factor value compared to the 3-dimensional analysis. The weakness of the analysis which is carried out in 3 dimensions is that it takes a lot of time and money, which makes practitioners not want to switch. The results of the 3D analysis can be conservative if they are validated with a 2D analysis in the most critical areas [2]. For conservative results, the 2D and 3D analyzes should not be significantly different, however in many cases the 3D analysis provides a higher safety factor value. The advantage given if the analysis is carried out in 3D is that it can represent the actual slope conditions in real terms and can determine the critical position. The value of the safety factor in the finite element method analysis with the 3-dimensional model is determined by the use of the material shear strength reduction factor (SRF) technique. in order to obtain the true safety factor, the srf is gradually increased until the model become failure (non-convergent). When this critical value is found, the safety factor of the slope model is equal to the reduction in material strength (SF) ≈ (SRF).

The determination of the critical value of srf is determined by increasing and decreasing the SRF value step by step until the highest SRF value is obtained in the model to be able to achieve the convergence criteria, see Table 4, the slope model enters a non-convergent (failure) at srf 1.05, so that the critical value of srf is 1.04 and highlighted in bold color.

Step	SRF	Solid Tolerance	Convergence
1	1	0.0005	Yes
2	1.3	0.2813	No
3	1.14	0.0691	No
4	1.06	0.0152	No
5	1.02	0.0008	Yes
6	1.04	0.0009	Yes
7	1.05	0.0078	No

Table 4.

Critical SRF determination.

2.3.2 Mesh type

Finite element analysis method employs the utilization of elements and nodes to perform stress–strain analysis acting on each element. The process of establishing these elements is called meshing. Mesh type will affect the analysis result. Greater number of elements will lead to a more accurate analysis result. However, this will result a more time-consuming computation. There are various mesh types in 3-dimensional analysis, there are 4-nodes and 10-nodes with uniform and graded shape available to this method. The example of elements and nodes utilization in FEM are given in Figure 17.

Figure 17.

4-node graded and uniform.

To recapitulate the results of the analysis can be seen in Figure 18. From this figure provides information that there is an effect of using the mesh type on the srf value of the analysis results obtained, this is because the number of elements and nodes used is also different.

Figure 18.

SRF VS mesh type.

2.3.3 Maximum iteration optimization

In model analysis using the finite element method, the determination of the convergence criteria is limited by the analysis calculation (maximum iteration), the more iterations allowed, the more accurate the analysis results will be, but it also needs to be considered in terms of the slope model being analyzed, in the use of the maximum number. Optimal iteration related to the level of efficiency in the analysis so as not to waste too long.

At first, the model will be analyzed the stresses that work on each element due to the applied load, it will get the principal major stress and minor for each element, then enter the material properties for each material and start entering the material strength reduction factor (SRF) value for each element, at this stage the material properties value can be increased or decreased depending on the error value (solid tolerance) obtained. Furthermore, an Elasto-plastic analysis was performed using the Mohr-Coulomb failure criterion to obtain a new error value (solid tolerance). If the error value is still above the maximum value within the allowable iteration calculation limit, the SRF value will be lowered until the error value is below the maximum limit. The recapitulate the analysis results can be seen in Table 5.

Max iteration	SRF	Max displacement (m)	Computing time (s)
10,000	1.04	8.60	110,459
5000	1.04	8.74	57,990
2500	1.04	8.48	26,214
1000	1.04	24.29	14,622
500	1.04	14.34	6021
450	1.04	12.90	5422
400	1.04	12.72	4788
350	1.03	5.25	3650
300	1.03	5.42	3341
250	1.03	6.51	3159
200	1.01	1.70	2300
150	1.01	1.70	1969
100	1.00	1.67	1704
50	0.98	1.20	1541
25	0.68	1.09	1501

Table 5.

Influence of variable max iteration.

The analysis shows that the higher the maximum number of iterations, the SRF value also increases Table 4. As the SRF value increases, in the maximum number of iterations of 400 the SRF value constant at a value of 1.04, so this provides information that at a maximum of 400 iterations is the maximum optimal number of iterations. If seen from the effect of the number of iterations on the total displacement, the results will fluctuate, this is closely related to the srf value because the higher the srf value, the total displacement will also increase, but at 1.04 the srf value does not change/is consistent but the total displacement fluctuates because it is solid tolerance, the resulting solid tolerance also varies due to the use of different maximum iterations, this is because the slope conditions remain non-convergent (energy equilibrium is not achieved) above srf 1.04 in the number of iterations that have been set, however, solid tolerance will get closer to the maximum value, while for the computation time it is very clear that the higher the maximum number of iterations, the computation time will also increase, because it will take more time to search for convergence with the maximum iteration limit given, although the results will remain the same, that non-converging above 1.04 and the last converging value is 1.04.

2.3.4 Analysis 3D finite element method of open pit mine

Analysis using the finite element method uses the same data as the analysis carried out with boundary equilibrium, namely rock characteristics, 3D slope geometry and 3D geological models. Young and poisson ratio are input parameters in this analysis. If the boundary equilibrium analysis requires the position and shape of the slip plane, the analysis does not require these assumptions but uses the elements and nodes that are attached to the model. For an example of a 3D mesh model, see Figure 19. For this example, case a 4-node element type with a graded gradient is used.

Figure 19.

3D meshing model.

After the model is successfully meshed, it is necessary to determine the limit of the resistance that works in the field of analysis to determine the internal force that works on the restraint installation which plays a role in determining the deformation limit of the model. For an example of the restraint model can be seen in Figure 20.

Figure 20.

3D restraint model.

Analysis with the finite element method can see information about stress acting on the model and the total displacement (Figure 21).

Figure 21.

3D FEM analysis result.

3. Conclusions

Slope stability using LEM shows that the cuckoo algorithm is reliable in obtaining position and shape of slip surface. In finite element analysis method, the optimum iteration number needs to be considered to improve analysis efficiency and preserving accuracy.

Acknowledgments

The authors would like to thank PT X who has facilitated the collection of field data, as well as the Mining Computation Laboratory of Mining Engineering Department FTKE of University Trisakti which has facilitated the author in using the software used in modeling research data.