Effect of undercut on the lower bound stability of vertical rock escarpment using finite element and power cone programming

Abstract

In the present study, the stability of a vertical rock escarpment is determined by considering the influence of undercut. Lower bound finite element limit analysis in association with Power Cone Programming (PCP) is applied to incorporate the failure of rock mass with the help of the Generalized Hoek-Brown yield criterion. The change in stability due to the presence of undercut is expressed in terms of a non-dimensional stability number (σ_ci/γH). The variations of the magnitude of σ_ci/γH are presented as design charts by considering the different magnitudes of undercut offset (H/v_u and w_u/v_u) from the vertical edge and different magnitudes of Hoek-Brown rock mass strength parameters (Geological Strength Index (GSI), rock parameter (m_i,), Disturbance factor (D)). The obtained results indicate that undercut can cause a severe stability problem in rock mass having poor strength. With the help of regression analysis of the computed results, a simplified design equation is proposed for obtaining σ_ci/γH. By performing sensitivity analysis for an undisturbed vertical rock escarpment, we have found that the undercut height ratio (H/v_u) is the most sensitive parameter followed by GSI, undercut shape ratio (w_u/v_u), and m_i. The developed design equation as well as design charts can be useful for practicing engineers to determine the stability of the vertical rock escarpment in the presence of undercut. Failure patterns are also presented to understand type of failure and extent of plastic state during collapse.

Graphical abstract

Keywords

undercut / vertical escarpment / stability / Hoek-Brown yield criterion / PCP

1 Introduction

Undercuts or notches are often observed at the toe of rock escarpments and cliffs in coastal regions. These features generally form due to erosion caused either by physical weathering or by chemical weathering from water and air (Fig.1(a)). Therefore, it becomes necessary to assess the stability of rock escarpments in the presence of undercut. In the past, many researchers [1–9] investigated the different factors corresponding to the formation and stability of undercut in steep slopes and cliffs.

In addition, several researchers [10–15] determined the case−specific stability of an overhanging rock escarpment and cliff in specific locations. Augustinus [10] investigated the effect of undercut in some glacial valley rock escarpments in New Zealand’s Southern Alps. Tsesarsky et al. [11] determined the stability of a 34 m height vertical rock slope in Haifa, Israel, having a notch offset of 11 m at toe, by using finite element analysis and discontinuous deformation analysis. Briaud [12] presented the different factors for rock cliff erosion and undercut formation in Normandy cliffs, France. Another case study in Haifa, Israel was presented by Tsesarsky and Hatzor [13], which reported the application of reinforcement of the vertical slope to avoid landslides due to the presence of undercut. Hayakawa and Matsukura [14] determined the stability of a waterfall cliff face in the presence of undercuts at Niagara Falls, USA, by using cantilever model analysis. Budetta [15] determined the stability of an undercut sea cliff in Campania, Southern Italy by using elastoplastic finite element analysis. Recently, for a vertical soil escarpment having undercut, Banerjee and Chakraborty [16] have presented the stability number in a series of designed charts. However, to the authors’ knowledge, no detailed study is available that deals with the stability of the vertical rock escarpment by considering the effect of undercut from a generalized point of view. Therefore, in the present study, a systematic analysis for the vertical rock escarpment is performed by considering the different extents of undercut (Fig.1(b)).

The failure of rock mass is modeled by considering the recently developed Generalized Hoek-Brown (GHB) yield criterion [17,18], which has been utilized by several researchers [19–25] to carry out different stability analyses. Amongst the available numerical methods, the lower bound finite element limit analysis (LBFELA) [26–28] does not assume any predefined failure surface, and most importantly, it always gives a safe stability number. Using LBFELA, some researchers [21,22,29,30] have shown that the conic programming technique is computationally very efficient in the implementation of the GHB yield criterion. Recently, Kumar and Rahaman [31] implemented Power Cone Programming (PCP) in a plane strain LBFELA framework for rock mass obeying the GHB yield criterion. Therefore, in this work, to determine the stability of a vertical rock escarpment, LBFELA in conjunction with PCP is employed. To express the variation of stability for a rock escarpment in the presence of undercut and rock mass strength properties, the obtained results are presented as design charts by using a non-dimensional stability number (σ_ci/γH). Furthermore, an empirical formulation is derived by utilizing the regression analysis on the computed results, so that the expression can be beneficial for engineering applications. In addition, a sensitivity analysis is also performed to check the significance effect of the input parameters on the obtained stability number.

The present study is outlined as follows. In Section 2, a brief review of the GHB yield criterion is discussed, and the variations of its different important parameters are presented. Section 3 contains the problem definition of the present study, with descriptions related to the selected finite element mesh as well as boundary conditions. In addition, the LBFELA framework in combination with the PCP is described. In Section 4, the developed code is employed in order to determine the variation in the magnitude of stability number and the failure patterns at collapse. The effects of different rock mass strength parameters on the stability number are discussed in Section 5. Next, the results of the present lower bound stability number are compared with the available literature in Section 6. A design equation is proposed, and a sensitivity analysis is performed for input parameters in Sections 7 and 8, respectively. Remarks and limitations related to the present study are described next in Section 9. Finally, conclusions are drawn in Section 10.

2 Generalized Hoek-Brown yield criterion

To model the rock mass at failure in the problem domain, the GHB yield criterion [17,18] is considered and is expressed as

(1)({\sigma }_1-{\sigma }_3)-{\left\{-{\sigma }_{\text{ci}}^{\left(\frac{1}{\alpha }-1\right)}{m}_{\text{b}}{\sigma }_1+s{\sigma }_{\text{ci}}\right\}}^{\alpha }⩽0,

where σ₁ and σ₃ are the major and minor principal stresses; σ_ci is the uniaxial compressive strength of intact rock mass; m_b, s, and α are the function of material constant (m_i), Geological Strength Index (GSI), and disturbance factor (D). The relationships between these parameters can be written as

\begin{array}{}\text{(2a)}& m_{\text{b}}=m_{\text{i}}\exp \left(\frac{GSI-100}{28-14D}\right),\end{array}

\begin{array}{}\text{(2b)}& s=\exp \left(\frac{GSI-100}{9-3D}\right),\end{array}

\begin{array}{}\text{(2c)}& \alpha =\frac{1}{2}+\frac{1}{6}\left({\text{e}}^{-\frac{GSI}{15}}-{\text{e}}^{-\frac{20}{3}}\right).\end{array}

Equation (1) is expressed by considering the normal tensile stress as positive. GSI values range between 10 and 100 to represent the rock mass behaviors between poor and intact; the magnitude of m_i varies between 1 and 35; whereas, D is selected between 0 and 1 based on disturbance of rock mass.

3 Problem definition, mesh, boundary details and methodology

An undercut having a horizontal width w_u, and vertical depth v_u is considered to be located in a rock escarpment of height (H), as shown in Fig.1(b). The rock mass is assumed to follow the GHB yield criterion. The aim of the present study is to evaluate the magnitude of stability number (σ_ci/γH) by considering the different influencing parameters, namely, the undercut height ratio (H/v_u), undercut shape ratio (w_u/v_u), and rock mass strength parameters (GSI, m_i, D). Here, γ is the unit weight of the rock mass. The stability number (σ_ci/γH) can be expressed as

(3)\frac{{\sigma }_{\text{ci}}}{\gamma H}=f\left(\frac{H}{v_{\text{u}}},\frac{w_{\text{u}}}{v_{\text{u}}},GSI,m_{\text{i}},D\right).

The chosen plane strain domain and associated boundary conditions are displayed in Fig.1(b). No surcharge load is considered in any of the free faces (OP, PQ, QR, RS, ST, and TM); therefore, the normal and shear stresses are assigned as zero (i.e., σ_n = 0 and τ_s = 0) along those boundaries. In addition, no boundary condition needs to be implemented along MN, NO, and OP as the shear strength of rock mass by default controls the stresses on these boundaries. The horizontal extent (L_h) and the vertical extent (T_r) of the domain are considered in a way that allows the boundary effect to be eliminated in the present study. The three-noded triangle elements (Fig.1(b)) are used to discretize the selected plane strain domain. A typical finite element mesh for GSI = 90, m_i = 15, D = 0, H/v_u = 3, w_u/v_u = 2, l_c/v_u = 1, d_c/v_u = 1 is demonstrated in Fig.1(c). The total numbers of elements, nodes, and discontinuities are described in Fig.1(c) by n_e, n_n, and n_d, respectively.

To showcase the mesh dependency on the stability number, five types of finite element meshes based on the number of elements are explored as shown in Tab.1. The different magnitudes of stability numbers are compared and it is found that the difference between the computed results for the very fine mesh and fine mesh is not significant. Moreover, the very fine mesh needs more computational time to obtain the stability number than fine mesh. Thus, the fine mesh is employed to perform the analysis.

By following the LBFELA formulation of Lysmer [26] and Sloan [27], the present study is performed to determine the stability of the vertical rock escarpment. Therefore, each node of the triangular elements has three unknown nodal stresses (σ_x, σ_y, and τ_xy). The magnitude of unit weight of rock mass (objective function) is maximized by applying the inequality and equality constraints. The inequality constraints are generated to satisfy the yield criterion at each node. As described earlier, the GHB yield criterion is considered in the problem domain to incorporate rock mass behavior at failure. In the present study, the Conic Programming technique PCP is utilized in the problem domain to model the GHB yield criterion [31] for rock mass. A detailed description regarding the application of the GHB yield criterion by using the PCP is presented. In addition, by following Sloan [27] as described below, the equality constraints are generated due to the requirement of fulfilling: 1) the equations of equilibrium throughout the domain; 2) the conditions of stress discontinuity along the edges of two adjacent elements; and 3) the conditions of stress boundary along the edges of boundary.

3.1 Equations of equilibrium

By considering the below equations, the stress equilibrium is satisfied in each element

\begin{array}{}\text{(4a)}& \frac{\mathrm{\partial }{\sigma }_x}{\mathrm{\partial }x}+\frac{\mathrm{\partial }{\tau }_{xy}}{\mathrm{\partial }y}=0,\end{array}

\begin{array}{}\text{(4b)}& \frac{\mathrm{\partial }{\tau }_{xy}}{\mathrm{\partial }x}+\frac{\mathrm{\partial }{\sigma }_y}{\mathrm{\partial }y}=\gamma ,\end{array}

where normal tensile stress is taken as positive. In LBFELA, the linear stress variation is anticipated within the each element. Therefore, by considering the concept of the linear shape function (N_i) in FEM, the unknown stresses within the each element can be written as

(5){\sigma }_x=\sum _{i=1}^3{N}_i{\sigma }_{x,i};{\sigma }_y=\sum _{i=1}^3{N}_i{\sigma }_{y,i};{\tau }_{xy}=\sum _{i=1}^3{N}_i{\tau }_{xy,i},

where i indicates the node number for an element. Now, by incorporating Eq. (5) into Eqs. (4a) and (4b), the below matrix form can be expressed

(6)[A_{\text{ce}}{]}_{2E\times 3N}\{{\sigma }^{\text{l}}{\}}_{3N\times 1}=\{b_{\text{ce}}{\}}_{2E\times 1},

where {\left\{{\sigma }^{\text{l}}\right\}}^{\mathrm{T}}=\{{\sigma }_{x,1}{\sigma }_{y,1}{\tau }_{xy,1}{\sigma }_{x,2}{\sigma }_{y,2}{\tau }_{xy,2}.........{\sigma }_{x,N}{\sigma }_{y,N}{\tau }_{xy,N}\}; vector {σ ^l} contains the nodal unknown stresses. The values of obtained matrix [A_ce] and vector {b_ce} are known.

3.2 Condition of discontinuity

In LBFELA, the continuity in shear and normal stresses is considered across the edges of two adjacent elements (j and k). Therefore, by utilizing the nodal pairs (1,3) and (2,4) (Fig.1(b)), four numbers of stress continuity equation can be written as

(7){\sigma }_{\text{n,1}}^j={\sigma }_{\text{n,3}}^k;{\sigma }_{\text{n,2}}^j={\sigma }_{\text{n,4}}^k;{\tau }_{\text{sh,1}}^j={\tau }_{\text{sh,3}}^k;{\tau }_{\text{sh,2}}^j={\tau }_{\text{sh,4}}^k.

On a plane having an angle ψ_h with the horizontal, the normal stress (σ_n) and shear stress (τ_sh) can be expressed as

\begin{array}{}\text{(8a)}& {\sigma }_{\text{n}}=({\sin }^2{\psi }_{\text{h}}){\sigma }_x+({\cos }^2{\psi }_{\text{h}}){\sigma }_y-(\sin 2{\psi }_{\text{h}}){\tau }_{xy},\end{array}

\begin{array}{}\text{(8b)}& {\tau }_{\text{sh}}=\left(-\frac{1}{2}\sin 2{\psi }_{\text{h}}\right){\sigma }_x+\left(\frac{1}{2}\sin 2{\psi }_{\text{h}}\right){\sigma }_y+(\cos 2{\psi }_{\text{h}}){\tau }_{xy}.\end{array}

By combining Eqs. (8) and (7), the discontinuity condition can be expressed for four nodes (1, 3 and 2, 4) in terms of twelve unknown stresses as

[A_{\text{dc}}^{\text{st}}{]}_{4\times 12}{\left\{{\sigma }_{\text{dc}}\right\}}_{12\times 1}={\left\{b_{\text{dc}}^{\text{st}}\right\}}_{4\times 1},

\begin{array}{c}[A_{\text{dc}}^{\text{st}}{]}_{4\times 12}=\left[\begin{array}{cccc}R& 0& -R& 0\\ 0& R& 0& -R\end{array}\right];\\ [R{]}_{2\times 3}=\left[\begin{array}{ccc}{\sin }^2{\psi }_{\text{h}}& {\cos }^2{\psi }_{\text{h}}& -\sin 2{\psi }_{\text{h}}\\ -{\displaystyle \frac{1}{2}}\sin 2{\psi }_{\text{h}}& {\displaystyle \frac{1}{2}}\cos 2{\psi }_{\text{h}}& \cos 2{\psi }_{\text{h}}\end{array}\right];\end{array}

(9)\begin{array}{rl}{\left\{{\sigma }_{\text{dc}}\right\}}_{12\times 1}=& \{{\sigma }_{x,1}^{\text{j}}{\sigma }_{y,1}^{\text{j}}{\tau }_{xy,1}^{\text{j}}{\sigma }_{x,2}^{\text{j}}{\sigma }_{y,2}^{\text{j}}{\tau }_{xy,2}^{\text{j}}\\ & {\sigma }_{x,3}^{\text{k}}{\sigma }_{y,3}^{\text{k}}{\tau }_{xy,3}^{\text{k}}{\sigma }_{x,4}^{\text{k}}{\sigma }_{y,4}^{\text{k}}{\tau }_{xy,4}^{\text{k}}\},\end{array}

where vector {σ_dc} comprises the nodal unknown stresses. The magnitudes of obtained matrix [A_{\text{dc}}^{\text{st}}] and vector \left\{b_{\text{dc}}^{\text{st}}\right\} are known.

3.3 Conditions for boundary

Along the boundary edge (m), the stress boundary conditions can be added as

(10){\sigma }_{\text{n,1}}^{\text{b}}=q_1^{\text{b}};{\sigma }_{\text{n,1}}^{\text{b}}=q_2^{\text{b}};{\tau }_{\text{sh,1}}^{\text{b}}=s_1^{\text{b}};{\tau }_{\text{sh,2}}^{\text{b}}=s_2^{\text{b}}.

where q_1^{\text{b}}, q_2^{\text{b}} and s_1^{\text{b}}, s_2^{\text{b}} are normal and shear stresses, respectively at the nodes 1 and 2 along the edge of boundary (m). If the boundary makes an angle ψ_b with horizontal, the boundary condition can be expressed in terms of matrix as

(11)[A_{\text{b}}^{\text{st}}{]}_{4\times 6}\{{\sigma }^{\text{b}}{\}}_{6\times 1}=\{b_{\text{b}}^{\text{st}}{\}}_{4\times 1},

where [A_{\text{b}}^{\text{st}}{]}_{4\times 6}=\left[\begin{array}{cc}{L}_{\text{b}}& 0\\ 0& L_{\text{b}}\end{array}\right]; [L_{\text{b}}{]}_{2\times 3}=\left[\begin{array}{ccc}{\sin }^2{\psi }_{\text{b}}& {\cos }^2{\psi }_{\text{b}}& -\sin 2{\psi }_{\text{b}}\\ -{\displaystyle \frac{1}{2}}\sin 2{\psi }_{\text{b}}& {\displaystyle \frac{1}{2}}\cos 2{\psi }_{\text{b}}& \cos 2{\psi }_{\text{b}}\end{array}\right];{\left\{{\sigma }^{\text{b}}\right\}}^{\text{T}}={\left\{\begin{array}{cccccc}{\sigma }_{x,1}^{\text{m}}& {\sigma }_{y,1}^{\text{m}}& {\tau }_{xy,1}^{\text{m}}& {\sigma }_{x,2}^{\text{m}}& {\sigma }_{y,1}^{\text{m}}& {\tau }_{xy,1}^{\text{m}}\end{array}\right\}}_{1\times 6}; \{b_{\text{b}}^{\text{st}}\}={\left\{\begin{array}{cccc}{q}_1^{\text{b}}& s_1^{\text{b}}& q_2^{\text{b}}& s_2^{\text{b}}\end{array}\right\}}_{1\times 4},

where vector {σ_b} comprises the nodal unknown stresses. The magnitudes of obtained matrix [A_{\text{b}}^{\text{st}}] and vector \left\{b_{\text{b}}^{\text{st}}\right\} are known.

3.4 Incorporation of yield condition

As described earlier, the GHB yield criterion [17,18] is considered to model the rock mass at failure and is written as

(12)({\sigma }_1-{\sigma }_3)-{\left\{\nu {\sigma }_1+s{\sigma }_{\text{ci}}\right\}}^{\alpha }⩽0\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}\nu =-{\sigma }_{\text{ci}}^{(\frac{1}{\alpha }-1)}{m}_{\text{b}}.

By following Kumar and Rahaman [31], a new variable (δ) is incorporated in Eq. (12) and presented as

(13)\delta ={\left\{\nu {\sigma }_1+s{\sigma }_{\text{ci}}\right\}}^{\alpha }\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}\delta ⩾0.

With the help of Eq. (13), Eq. (12) can be presented as

(14){\sigma }_1-{\sigma }_3⩽\delta .

For a plane strain problem, principal stresses can be described as

\begin{array}{}\text{(15a)}& {\sigma }_1=\frac{{\sigma }_x+{\sigma }_y}{2}+\frac{1}{2}\sqrt{{({\sigma }_x-{\sigma }_y)}^2+4{\tau }_{xy}^2},\end{array}

\begin{array}{}\text{(15b)}& {\sigma }_3=\frac{{\sigma }_x+{\sigma }_y}{2}-\frac{1}{2}\sqrt{{({\sigma }_x-{\sigma }_y)}^2+4{\tau }_{xy}^2}.\end{array}

From Eqs. (15a) and (15b), it can be expressed as

(16){\sigma }_1-{\sigma }_3=\sqrt{{({\sigma }_x-{\sigma }_y)}^2+4{\tau }_{xy}^2}.

From Eqs. (14) and (16), it can be written as

(17)\sqrt{{({\sigma }_x-{\sigma }_y)}^2+4{\tau }_{xy}^2}⩽\delta .

Equation (17) can be presented in the form of the standard second-order cone,

\begin{array}{}\text{(18a)}& \sqrt{{\kappa }_1^2+{\kappa }_2^2}⩽{\kappa }_3,\end{array}

where

\begin{array}{}\text{(18b)}& {\kappa }_1={\sigma }_x-{\sigma }_y;{\kappa }_2=2{\tau }_{xy};{\kappa }_3=\delta ,\end{array}

where κ₁, κ₂, and κ₃ are the auxiliary variables. By replacing the major principal stress (σ₁) as expressed in Eq. (15a) into Eq. (13), it can be written as

(19)\delta ={\left[\nu \left\{\frac{{\sigma }_x+{\sigma }_y}{2}+\frac{1}{2}\sqrt{{({\sigma }_x-{\sigma }_y)}^2+4{\tau }_{xy}^2}\right\}+s{\sigma }_{\text{ci}}\right]}^{\alpha }.

By expanding Eq. (17), the below form of Eq. (20) can be achieved.

(20)\begin{array}{rl}& {\left[\nu \left\{\frac{{\sigma }_x+{\sigma }_y}{2}+\frac{1}{2}\sqrt{{({\sigma }_x-{\sigma }_y)}^2+4{\tau }_{xy}^2}\right\}+s{\sigma }_{\text{ci}}\right]}^{\alpha }\\ ⩽& {\left[\nu \left\{\frac{{\sigma }_x+{\sigma }_y}{2}+\frac{1}{2}\delta \right\}+s{\sigma }_{\text{ci}}\right]}^{\alpha },\end{array}

where δ ≥ 0 and α > 0.

With the help of Eqs. (19) and (20), it can be written as

(21)\delta ⩽{\left[\nu \left\{\frac{{\sigma }_x+{\sigma }_y}{2}+\frac{1}{2}\delta \right\}+s{\sigma }_{\text{ci}}\right]}^{\alpha }.

Equation (21) can be presented in the form of the power cone equation, which is generally written as

\begin{array}{}\text{(22a)}& \sqrt{{\kappa }_6^2}⩽{\kappa }_4^{\alpha }{\kappa }_5^{(1-\alpha )},\end{array}

where

\begin{array}{}\text{(22b)}& {\kappa }_4=\nu \left\{\frac{{\sigma }_x+{\sigma }_y}{2}+\frac{1}{2}\delta \right\}+s{\sigma }_{\text{ci}};{\kappa }_5=1;{\kappa }_6=\delta ,\end{array}

where κ₄, κ₅, and κ₆ are the auxiliary variables. In the LBFELA framework, the obtained inequality Eqs. (18a) and (22a) are imposed at every nodes in terms of unknown stresses. The yield conditions can be expressed as

\begin{array}{}\text{(23a)}& [A_{\text{Scone}}^{\text{yield}}]\{{\sigma }^{\text{y}}\}⩽\{b_{\text{Scone}}^{\text{yield}}\},\end{array}

\begin{array}{}\text{(23b)}& [A_{\text{Pcone}}^{\text{yield}}]\{{\sigma }^{\text{y}}\}⩽\{b_{\text{Pcone}}^{\text{yield}}\},\end{array}

where the magnitude of matrixes [A_{\text{Scone}}^{\text{yield}}] and [A_{\text{Pcone}}^{\text{yield}}] and the magnitude of vectors \{b_{\text{Scone}}^{\text{yield}}\} and \{b_{\text{Pcone}}^{\text{yield}}\} are known. Note that vector {σ^y} contains the nodal unknown stresses.

Finally, to solve the problem by applying the PCP technique, all the constraint equations are assembled to form the global matrixes and vectors. The final form of the optimization problem can be expressed as

\begin{array}{}\text{(24a)}& \mathrm{M}\mathrm{a}\mathrm{x}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}\gamma ,\end{array}

\begin{array}{}\text{(24b)}& \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{e}\mathrm{d}\mathrm{t}\mathrm{o}[A_1]\left\{\overline{x}\right\}=\left\{b_1\right\},\end{array}

\begin{array}{}\text{(24c)}& [A_2]\left\{\overline{x}\right\}⩽\left\{b_2\right\},\end{array}

where {\overline{x}} is the global vector comprising the unknown stresses and auxiliary variables; [A₁] and {b₁} are the global matrix and vector of the equality constraints; [A₂] and {b₂} are the global matrix and vector of the inequality constraints. One can go through Kumar and Rahaman [31] and Chakraborty and Kumar [32] for the complete formulations of equilibrium, discontinuity, boundary, and yield constraints. An in-house computer code is written and implemented in MATLAB [33] to perform the analysis by applying LBFELA with the help of PCP. In the present study, the primal-dual interior point solver, MOSEK [34], is utilized to solve the optimization problem.

4 Results

4.1 Variation in the magnitude of σ_ci/γH

The change in stability due to the presence of the undercut is examined by varying: (1) GSI values (ranging from 10 to 90), (2) m_i values (ranging from 5 to 35), (3) w_u/v_u values (ranging from 1 to 5) (4) H/v_u values (ranging from 3 to 8), and (5) D values (ranging from 0 to 1). The results are presented by using a stability number σ_ci/γH in a series of design charts. The design charts are illustrated in Fig.2−Fig.4 for H/v_u = 3, 5, and 7. In addition, for H/v_u = 4, 6, and 8, the results are provided in Figs. S1−S3 (refer to the supplementary data). It needs to be described that a higher value of σ_ci/γH indicates a lower stability of the vertical rock escarpment.

In the present study, the magnitude of σ_ci/γH increases with the increase of w_u/v_u value. In contrast, the magnitude of σ_ci/γH decreases significantly with increase of the H/v_u value. Fig.2−Fig.4 and S1−S3 show that σ_ci/γH decreases with increase of the GSI value. However, the magnitude of σ_ci/γH increases considerably with the increase in m_i value. In addition, Fig.2−Fig.4 and S1−S3 show that the magnitude of σ_ci/γH increases with the D value.

4.2 Failure patterns

Fig.5 presents the failure patterns for the vertical rock escarpment in the presence of an undercut.The proximity of stress state at each node in rock mass with respect to the yield is expressed by a ratio a/d; where a = (σ₁ − σ₃), and d = {\sigma }_{\mathrm{c}\mathrm{i}}{\left(-m_{\mathrm{b}}{\displaystyle \frac{{\sigma }_1}{{\sigma }_{\mathrm{c}\mathrm{i}}}}+s\right)}^{\alpha }. The value of a/d becomes unity for the plastic state, and it becomes less than unity for the non-plastic state. Fig.5 illustrates the failure patterns for an undercut in the vertical rock escarpment with m_i = 5, D = 0, H/v_u = 5, w_u/v_u = 1, 3, and 5, and GSI = 50 and 90. It needs to be mentioned that the plastic zone determined by LBFELA only represents the initial failure having very minute deformations. It is observed that the plastic zone spreads away from the escarpment face with the increase of w_u/v_u value. In addition, it is found that for GSI = 50, the additional highly stressed zone (i.e., the dark zone) develops towards the escarpment face from near the upper edge of undercut. On the other hand, the highly stressed zone develops away from the escarpment face for GSI = 90. It is also noted that with the increase in w_u/v_u value, the plastic zone increases on the top surface and reduces towards the undercut base. For a vertical escarpment having H/v_u = 5, GSI = 50, m_i = 5, D = 0, the plastic zone increases on the top surface from nearly 0.82v_u to 1.51v_u with the increase in w_u/v_u value from 3 to 5 (Fig.5(c) and 5(e)). The effect of disturbance factor is shown in Fig.5(g) and 5(h). An increase of the plastic zone on the top surface is observed with the increase in D value.

5 Discussion

Due to the presence of the undercut, the variation of stability is examined in detail and presented in Fig.2−Fig.4 and S1−S3 by considering the effects of w_u/v_u, H/v_u, GSI, m_i, and D. For a clear understanding, Fig.6 presents the effects of parameters on stability number.

5.1 Effect of w_u/v_u on σ_ci/γH

For a vertical rock escarpment having H/v_u = 3, GSI = 50, m_i = 5, D = 0, the magnitude of σ_ci/γH increases from 72.33 to 1435.19 with the increse of w_u/v_u from 1 to 5. Therefore, the stability of the vertical rock escarpment decreases by a factor of approximately 20 when w_u/v_u increases from 1 to 5 for this specific case. While at H/v_u = 8, the magnitude of σ_ci/γH increases from 19.04 to 128.94 with the increase of w_u/v_u from 1 to 5. For this specific case, the stability of the vertical rock escarpment decreases by a factor of about 6.8 with the increase of w_u/v_u value from 1 to 5.

5.2 Effect of H/v_u on σ_ci/γH

For a vertical rock escarpment having w_u/v_u = 1, GSI = 50, m_i = 5, D = 0, the magnitude of σ_ci/γH decreases from 72.33 to 19.04 with the increase of H/v_u from 3 to 8. Therefore, the stability of the vertical rock escarpment increases by a factor of about 2.8 with the H/v_u from 3 to 8 for this specific case. While at w_u/v_u = 5, the magnitude of σ_ci/γH increases from 1425.19 to 128.94 with the increase of w_u/v_u from 1 to 5. For this specific case, the stability of the vertical rock escarpment increases by a factor of nearly 10 with the increase of w_u/v_u value from 1 to 5.

5.3 Effect of GSI on σ_ci/γH

For a vertical rock escarpment having w_u/v_u = 3, H/v_u = 5, m_i = 5, D = 0, the magnitude of σ_ci/γH decreases from 2910.25 to 8.33 with the increase of GSI from 10 to 90. Thus, stability of a vertical rock escarpment is increased by a factor of nearly 350 if GSI increases from 10 to 90. When m_i = 35, the magnitude of σ_ci/γH decreases from 18309.79 to 44.30 with the increase of GSI from 10 to 90. For this specific case, the stability of the vertical rock escarpment increases by a factor of nearly 412 with the increase of GSI value from 10 to 90.

5.4 Effect of m_i on σ_ci/γH

For a vertical rock escarpment having w_u/v_u = 3, H/v_u = 3, GSI = 50, D = 0, the magnitude of σ_ci/γH increases from 517.75 to 3547.62 with the increase of m_i from 5 to 35. Therefore, the stability of vertical rock escarpment decreases by a factor of around 5.8 when the m_i value changes from 5 to 35 for this specific case. While at H/v_u = 8, the magnitude of σ_ci/γH increases from 1204.86 to 6394.14 with the increase of m_i from 5 to 35. For this specific case, the stability of the vertical rock escarpment decreases by a factor of nearly 4.3 with increase of m_i from 5 to 35.

5.5 Effect of D on σ_ci/γH

For a vertical rock escarpment having w_u/v_u = 3, H/v_u = 4, GSI = 10, m_i = 5, the magnitude of σ_ci/γH increases from 4891.82 to 38908.53 with the increase of D from 0 to 1. For this specific case, the obtained result designates that the stability of the vertical rock escarpment drops by a factor of nearly 6.9 with the increase of D from 0 to 1. While at GSI = 90, the magnitude of σ_ci/γH increases from 13.36 to 16.67 with the increase of D from 0 to 1. Therefore, the stability of vertical rock escarpment decreases around 25% when the D value changes from 0 to 1 for this specific case.

6 Comparison

To the author’s knowledge, no stress-based study is available for the vertical rock escarpment based on considering the effect of undercut. Therefore, for the sake of comparison, a vertical rock escarpment having no undercut is considered in the present study. For the different magnitudes of GSI having m_i = 10, and D = 0, the calculated values of σ_ci/γH are compared (refer Fig.7) with the solutions of Michalowski and Park [35] for the upper bound limit analysis. It can be found that the magnitudes of σ_ci/γH provided by Michalowski and Park [35] are a little lower than the presently obtained lower bound solutions. However, the magnitudes of σ_ci/γH clearly show a similar trend to those of Michalowski and Park [35] for the different magnitudes of GSI.

7 Design equation for stability of vertical rock escarpment

In order to evaluate stability of vertical rock escarpment, a design equation is presented for undisturbed rock mass with the help of IBM SPSS STATISTICS 22. After several trial and error iterations, the below equation is proposed.

(25)\frac{{\sigma }_{\text{ci}}}{\gamma H}={\text{e}}^{F_1+F_2+F_3},

where F_1=\left(r_1{\left({\displaystyle \frac{w_{\text{u}}}{v_{\text{u}}}}\right)}^{r_2}{\left(GSI\right)}^{r_3}{\left(m_{\text{i}}\right)}^{r_4}{\left({\displaystyle \frac{H}{v_{\text{u}}}}-2\right)}^{r_5}\right); F_2= \left(r_6{\left({\displaystyle \frac{w_{\text{u}}}{v_{\text{u}}}}\right)}^{r_7}{\left(GSI\right)}^{r_8}{\left(m_{\text{i}}\right)}^{r_9}\right); F_3=\left(r_{10}{\left(GSI\right)}^{r_{11}}{\left(m_{\text{i}}\right)}^{r_{12}}{\left({\displaystyle \frac{H}{v_{\text{u}}}}-2\right)}^{r_{13}}\right).

where r_i (i varies from 1 to 13) are the constant coefficients. This equation is valid for undisturbed rock mass (D = 0) having H/v_u > 3. With the help of nonlinear least square regression analysis, the optimal values of the constant coefficients are obtained by minimizing the summation of the square of deviation between the results of the proposed equation and computed stability number. The optimum values of the constant coefficients are presented in Tab.2. To showcase the accuracy of the stability number obtained from the proposed equation, the coefficient of determination (R²) is calculated. The predicted magnitude of stability number is reasonably accurate as the obtained R² value is 99.6% (Fig.8).

8 Sensitivity analysis

In order to understand the sensitivity of the different parameters, What-If based sensitivity analysis is performed with the help of MS EXCEL by using the developed design equation. The input parameters are randomly selected and presented in Tab.3. Fig.9 presents the fluctuation of stability number as a Tornedo Chart. It is found that for undisturbed (D = 0) vertical rock escarpment, H/v_u is the most sensitive parameter followed, in order, by GSI, w_u/v_u, and m_i.

9 Remarks and limitations of the present study

In the present study, the effect of the undercut is numerically examined for the vertical rock escarpment by considering the structured mesh and the rectangular shape of undercut. Here, the original three dimensional problem is simplified to two dimensional plane strain problem. In addition, the outline of any undercut in reality can be irregular in shape. In that case, the obtained design charts can be useful by considering the equivalent undercut area (w_uv_u) for a particular magnitude of w_u/v_u because the presented design charts are developed by using w_u/v_u ratio. Here, w_u and v_u both are variable and are expressed as a dimensionless ratio.

It can be noted that for the case of w_u/v_u = 1, the stress concentration around the corner can be avoided by considering a arc-shaped undercut. Therefore, in future, the stability of vertical rock escarpment by considering the arc-shaped undercut can be performed. It should be cleared that the present study only determines the lower bound stability number; whereas, the exact stability number remains between the upper and lower bound solutions.

Sometimes, cracks can develop in a rock escarpment. In order to model the crack propagation, several studies were performed by considering cracking-particle method [36–38] as well as efficient and advanced remeshing techniques [39–41]. In addition, several advanced fracture modeling techniques like nonlocal operator methods [42] and explicit phase field methods [43] have also been incorporated to model the propagation of a crack. Therefore, the validity of the conclusions of this work can be inspected in future by incorporating analysis of crack propagation in the rock escarpment.

10 Conclusions

A general guideline for the stability of the vertical rock escarpment is carried out by implementing the numerical investigation in the presence of undercut. Therefore, the stability of the vertical rock escarpment is investigated in terms of a dimensionless stability number (σ_ci/γH) by considering different values of the undercut height ratio (H/v_u), undercut shape ratio (w_u/v_u), and rock mass parameters (GSI, m_i, D). The influences of the different parameters are illustrated as design charts by using σ_ci/γH. The LBFELA in conjunction with the PCP is employed to perform the numerical investigation. It is found that the stability of the vertical rock escarpment decreases with increase of w_u/v_u value. With increase of w_u/v_u from 1 to 5, the stability of the vertical escarpment decreases by a factor of nearly 11 for GSI = 50 having H/v_u = 5, m_i = 5, D = 0. The stability of the vertical rock escarpment increases significantly with the increase of H/v_u value. The stability of the vertical rock escarpment increases 54% and 88.9% with the increase of H/v_u value from 3 to 4, and 3 to 8, respectively, having GSI = 50, w_u/v_u = 3, m_i = 5, D = 0. The degradation of stability due to the weathering effect for a vertical rock escarpment can be calculated by comparing the obtained stability number in the design charts with real field measurements. It is expected that the developed design charts and design equation will be beneficial for practicing engineers to determine the stability of the vertical rock escarpment in the presence of undercut.

References

[1] Sunamura T. Geomorphology of rocky coasts. Coastal Morphology and Research, 1994, 3: 174–175

[2] Budetta P, Galietta G, Santo A. A methodology for the study of the relation between coastal cliff erosion and the mechanical strength of soils and rock masses. Engineering Geology, 2000, 56(3−4): 243–256

[3] Abderahman N. Evaluating the influence of rate of undercutting on the stability of slopes. Bulletin of Engineering Geology and the Environment, 2007, 66(3): 303–309

[4] Chu-Agor M L, Fox G A, Cancienne R M, Wilson G V. Seepage caused tension failures and erosion undercutting of hillslopes. Journal of Hydrology (Amsterdam), 2008, 359(3−4): 247–259

[5] Mirenkov V E. On probable failure of an undercut rock mass. Journal of Mining Science, 2009, 45(2): 105–111

[6] Young A P, Guza R T, O’Reilly W C, Hansen J E, Barnard P L. Short-term retreat statistics of a slowly eroding coastal cliff. Natural Hazards and Earth System Sciences, 2011, 11(1): 205–217

[7] Katz O, Mushkin A. Characteristics of sea-cliff erosion induced by a strong winter storm in the eastern Mediterranean. Quaternary Research, 2013, 80(1): 20–32

[8] Budetta P, Luca C D. Wedge failure hazard assessment by means of a probabilistic approach for an unstable sea-cliff. Natural Hazards, 2015, 76(2): 1219–1239

[9] Ukritchon B, Ouch R, Pipatpongsa T, Khosravi M H. Investigation of stability and failure mechanism of undercut slopes by three-dimensional finite element analysis. KSCE Journal of Civil Engineering, 2018, 22(5): 1730–1741

[10] Augustinus P C. Rock mass strength and the stability of some glacial valley slopes. Zeitschrift für Geomorphologie, 1995, 39(1): 55–68

[11] TsesarskyMHatzorY HLeviathanISaltzmanUSokolowksyM. Structural control on the stability of overhanging, discontinuous rock slopes. In: Alaska Rocks, 40th US Symposium on Rock Mechanics (USRMS). Anchorage, AK: American Rock Mechanics Association, 2005

[12] Briaud J L. Case histories in soil and rock erosion: Woodrow Wilson bridge, Brazos River Meander, Normandy Cliffs, and New Orleans Levees. Journal of Geotechnical and Geoenvironmental Engineering, 2008, 134(10): 1425–1447

[13] Tsesarsky M, Hatzor Y H. Kinematics of overhanging slopes in discontinuous rock. Journal of Geotechnical and Geoenvironmental Engineering, 2009, 135(8): 1122–1129

[14] Hayakawa Y S, Matsukura Y. Stability analysis of waterfall cliff face at Niagara Falls: An implication to erosional mechanism of waterfall. Engineering Geology, 2010, 116(1−2): 178–183

[15] Budetta P. Stability of an undercut sea-cliff along a cilento coastal stretch (Campania, Southern Italy). Natural Hazards, 2011, 56(1): 233–250

[16] Banerjee S K, Chakraborty D. Influence of undercut and surface crack on the stability of a vertical escarpment. Geomechanics and Engineering, 2017, 12(6): 965–981

[17] Hoek E, Carranza-Torres C, Corkum B. Hoek-Brown failure criterion—2002 edition. Proceedings of NARMS-Tac, 2002, 1(1): 267–273

[18] Hoek E, Brown E T. The Hoek–Brown failure criterion and GSI–2018 edition. Journal of Rock Mechanics and Geotechnical Engineering, 2019, 11(3): 445–463

[19] Merifield R S, Lyamin A V, Sloan S W. Limit analysis solutions for the bearing capacity of rock masses using the generalised Hoek–Brown criterion. International Journal of Rock Mechanics and Mining Sciences, 2006, 43(6): 920–937

[20] Li A J, Merifield R S, Lyamin A V. Effect of rock mass disturbance on the stability of rock slopes using the Hoek–Brown failure criterion. Computers and Geotechnics, 2011, 38(4): 546–558

[21] Kumar J, Mohapatra D. Lower-bound finite elements limit analysis for Hoek-Brown materials using semidefinite programming. Journal of Engineering Mechanics, 2017, 143(9): 04017077

[22] Kumar J, Rahaman O. Lower bound limit analysis of unsupported vertical circular excavations in rocks using Hoek−Brown failure criterion. International Journal for Numerical and Analytical Methods in Geomechanics, 2020, 44(7): 1093–1106

[23] Das S, Chakraborty D. Effect of interface adhesion factor on the bearing capacity of strip footing placed on cohesive soil overlying rock mass. Frontiers of Structural and Civil Engineering, 2021, 15(6): 1494–1503

[24] Das S, Chakraborty D. Effect of soil and rock interface friction on the bearing capacity of strip footing placed on soil overlying Hoek–Brown rock mass. International Journal of Geomechanics, 2022, 22(1): 04021257

[25] Das S, Halder K, Chakraborty D. Bearing capacity of interfering strip footings on rock mass. Geomechanics and Geoengineering, 2022, 17(3): 883–895

[26] Lysmer J. Limit analysis of plane problems in soil mechanics. Journal of the Soil Mechanics and Foundations Division, 1970, 96(4): 1311–1334

[27] Sloan S W. Lower bound limit analysis using finite elements and linear programming. International Journal for Numerical and Analytical Methods in Geomechanics, 1988, 12(1): 61–77

[28] Sloan S W. Geotechnical stability analysis. Geotechnique, 2013, 63(7): 531–571

[29] Makrodimopoulos A, Martin C. Lower bound limit analysis of cohesive-frictional materials using second-order cone programming. International Journal for Numerical Methods in Engineering, 2006, 66(4): 604–634

[30] Krabbenhøft K, Lyamin A V, Sloan S W. Three-dimensional Mohr−Coulomb limit analysis using semidefinite programming. Communications in Numerical Methods in Engineering, 2008, 24(11): 1107–1119

[31] Kumar J, Rahaman O. Lower bound limit analysis using power cone programming for solving stability problems in rock mechanics for generalized Hoek–Brown criterion. Rock Mechanics and Rock Engineering, 2020, 53(7): 3237–3252

[32] Chakraborty D, Kumar J. Stability of a long unsupported circular tunnel in soils with seismic forces. Natural Hazards, 2013, 68(2): 419–431

[33] MATLAB. Version 8.5. Natick, MA: MathWorks, 2015

[34] MOSEKApS. Version 9.0. Copenhagen: MOSEK, 2020

[35] Michalowski R L, Park D. Stability assessment of slopes in rock governed by the Hoek−Brown strength criterion. International Journal of Rock Mechanics and Mining Sciences, 2020, 127: 104217

[36] RabczukTZiGBordasSNguyen-XuanH. A simple and robust three-dimensional cracking-particle method without enrichment. Computer Methods in Applied Mechanics and Engineering, 2010, 199(37−40): 2437−2455

[37] Rabczuk T, Belytschko T. Cracking particles: A simplified meshfree method for arbitrary evolving cracks. International Journal for Numerical Methods in Engineering, 2004, 61(13): 2316–2343

[38] Rabczuk T, Belytschko T. A three-dimensional large deformation meshfree method for arbitrary evolving cracks. Computer Methods in Applied Mechanics and Engineering, 2007, 196(29−30): 2777–2799

[39] Areias P, Reinoso J, Camanho P P, de Sá J C, Rabczuk T. Effective 2D and 3D crack propagation with local mesh refinement and the screened Poisson equation. Engineering Fracture Mechanics, 2018, 189: 339–360

[40] Areias P, Rabczuk T. Steiner-point free edge cutting of tetrahedral meshes with applications in fracture. Finite Elements in Analysis and Design, 2017, 132: 27–41

[41] Areias P, Msekh M A, Rabczuk T. Damage and fracture algorithm using the screened Poisson equation and local remeshing. Engineering Fracture Mechanics, 2016, 158: 116–143

[42] Ren H, Zhuang X, Rabczuk T. A higher order nonlocal operator method for solving partial differential equations. Computer Methods in Applied Mechanics and Engineering, 2020, 367: 113132

[43] Zhuang X, Ren H, Rabczuk T. Nonlocal operator method for dynamic brittle fracture based on an explicit phase field model. European Journal of Mechanics. A, Solids, 2021, 90: 104380

RIGHTS & PERMISSIONS

Higher Education Press