# Stability Analysis of Axisymmetric Concave Slopes Based on Two-Dimensional Limit Equilibrium Approach considering Additional Shear Resistance

Authors: Xiao-ming Liu; Rui Zhang; Jie Han; Sha Chen
Journal: Advances in Civil Engineering
 (2019)

---

## Abstract

Axisymmetric concave slopes, one special type of three-dimensional (3D) slopes, may be encountered in mining and civil engineering practice. Analysis of 3D slopes is generally complex and mostly relies on complicated numerical simulations. This paper proposes an elastoplastic solution for determining the additional shear resistances due to spatial effects of axisymmetric concave slopes. By incorporating the extra antislide forces, this paper proposes a simplified two-dimensional (2D) limit equilibrium procedure for the stability analysis of axisymmetric concave slopes. Combined with an iteration algorithm, the procedure can obtain the factors of safety for axisymmetric concave slopes in a simple and efficient way. Comparisons of the results from the proposed method and the numerical software FLAC3D are performed to demonstrate the validity of the proposed method for practical applications. Finally, the effects of several key parameters on the stability of axisymmetric concave slopes are investigated through a parametric study.


---

## Body

## 1. Introduction

Axisymmetric concave slopes, one special type of three-dimensional (3D) slopes, may be encountered in mining and civil engineering projects, such as open-pit excavations and earth-fill cofferdams. These slopes have a round or approximately round shape in a plan view but a 3D geometry in nature. Because of spatial 3D effects of stresses in the slope body, axisymmetric concave slopes cannot be directly treated as infinite long straight slopes with a uniform cross section; hence, typical two-dimensional (2D) limit equilibrium procedures are no longer applicable for axisymmetric concave slopes in a straightforward way. For this reason, further research is needed to investigate the stability of axisymmetric concave slopes.In the past, 2D limit equilibrium methods (LEM) have been commonly used for stability analyses of slopes in geotechnical engineering (e.g., [1–8]). In a typical 2D limit equilibrium method, a slip surface is often assumed a priori and the failure body, encompassed by the slip surface and the slope surface, is then divided into slices of soil mass. After that, force equilibrium and/or moment equilibrium conditions are established to solve the factor of safety (FoS), which is typically defined as the ratio of the shear strength (or resisting moment) to the shear stress (or driving moment).Considering that the actual failure surface of a slope is generally 3D in nature, a few researchers carried out stability analyses of slopes based on limit equilibrium approaches with 3D failure surfaces, e.g., Leshchinsky et al. [9], Lam and Fredlund [10], Chen et al. [11], and, more recently, Jiang and Zhou [12]. These studies assumed straight slopes with 3D failure surfaces (cylindrical, spherical, or others). However, as compared to straight slopes, axisymmetric concave slopes have a round or approximately round shape in the plan view, which has apparent spatial effects of stresses in the slope body. Ignoring these effects may lead to inaccurate or too conservative calculation of factors of safety for axisymmetric concave slopes.Numerical simulations are considered to be a powerful tool for analyzing the spatial (or 3D) effects developed in a concave slope. For example, Lorig [13] evaluated the stability of concave slopes using the numerical software FLAC3D, and Sun et al. [14] employed the displacement finite element software ABAQUS to develop stability charts for convex and concave slopes. In addition to numerical simulations, the method of characteristics has been used to analyze 3D stability of concave slopes, e.g., Jenike and Yen [15] and Jahanandish and Keshavarz [16]. However, these studies assumed homogeneous concave slopes.Although 3D limit equilibrium approaches or 3D numerical simulations may provide more reasonable results for stability analyses of slopes than 2D limit equilibrium approaches, building 3D models is generally time-consuming and usually requires special expertise of designers. In contrast, 2D limit equilibrium approaches are preferable for a preliminary analysis of slope stability in engineering practice because of their simplicity and long history of use. Generally, the FoS obtained by the 2D limit equilibrium approach is smaller than that based on 3D models. However, if additional shear resistance due to spatial effects of 3D slopes can be considered in the analysis, the results of the FoS based on the 2D limit equilibrium approach will be close to the actual value.Based on a 3D limit equilibrium approach, Zhang [17] proposed a practical method for 3D stability analysis of concave slopes in the plan view. In his method, the failure mass was divided into n vertical columns and the force and moment equilibrium conditions of the failure mass were established to calculate the FoS. The most important contribution of this work is that it considered the so-called “end force” P caused by the lateral pressure of soil acting on individual columns, i.e., accounted for the spatial (3D) effects of a concave slope. Because the end force P (normal to the radial line) is not perpendicular to the sliding direction, it provides additional shear resistance on the failure mass. Zhang [17] calculated the end force P empirically using the following equation:(1)P=KaγhVxz,where γ is the unit weight of soil, h is the height of a column, Ka is the active earth pressure coefficient, and Vxz is the projection of the base of each column V on the X-Z plane.In the authors’ opinion, the use of active earth pressure for the approximation of the end forceP may be too conservative in practical applications. This is because, when the concave slope fails, the sliding columns squeeze each other in the circumferential direction and the soil pressure in the circumferential direction will be much greater than active and static earth pressures. Therefore, a better method concerning the calculation of the soil pressure in the circumferential direction needs to be proposed for the stability analysis of concave slopes, which constitutes a major objective of this paper.Based on the above discussions, an elastoplastic solution is proposed in this study to give more reasonable estimations of the additional shear resistance due to the spatial effects of an axisymmetric concave slope. For this purpose, the concave slope is divided into a series of arch-shaped slices in the plan view (or thick cylinders in the space view). Then, the earth pressures acting on inner and outer faces of an individual cylinder are simplified as axisymmetric linearly distributed loads with depth, and the surcharge load on the slope surface is simplified as a uniformly distributed load acting on the upper face of the cylinder. Based on these simplified boundary conditions, an elastoplastic solution is deduced for the distributions of stressesσθ,σr,σz,andτrz in the cylinder. Finally, the solution of the circumferential stress σθ is used to calculate the additional shear resistance acting on individual column elements for a 2D limit equilibrium analysis for concave slopes. With the aforementioned elastoplastic solution, a simplified 2D limit equilibrium procedure for the stability analysis of concave slopes is proposed later in this paper. The proposed method is verified by comparing the calculated results with those from the numerical software FLAC3D.


## 2. Problem Descriptions

### 2.1. Geometry of a Concave Slope

Figure1 shows the geometry parameters of a typical axisymmetric concave slope with horizontal layers including (1) the slope height (H), (2) the slope angle (βc), and (3) the radius at the toe of the slope (Rc). In order to conduct a 3D limit equilibrium analysis, the concave slope can be divided into m × n individual vertical columns. This can be easily done by equally dividing the concave slope into m cylinder slices (in the plan view) along the radial direction and n wedge-shaped blocks along the circumferential direction. In the plan view, the overlapping area of the i-th cylinder slice and the j-th wedge-shaped block is designated as column (i, j), which is a fundamental element for the 3D limit equilibrium analysis of an axisymmetric concave slope.Figure 1
Typical geometry of an axisymmetric concave slope (1/4 model shown).


### 2.2. Definition of Additional Shear Resistance

Figure2(a) shows a typical column element (i, j) taken from an axisymmetric concave slope. Prior to defining the additional shear resistance, it is necessary to analyze all the internal and external forces acting on this element. The 3D limit equilibrium analysis of a general axisymmetric slope includes the following forces: (1) the normal force Ni and the shear force Si at the base of the column; (2) the normal force Pui due to pore-water pressure at the base of the column; (3) the external vertical forces Wi and Qi due to soil weight and surcharge load of the slope, respectively; (4) the intercolumn normal forces Eri−1 and Eri in the r-direction and Eθi in the θ-direction; and (5) the vertical intercolumn shear forces Xri−1 and Xri in the r-direction. Note that the subscript j for all these quantities has been omitted considering the inherent axisymmetry of the studied problem.Figure 2
Typical column element and its internal and external forces.
(a)
(b)The forcesNi, Si, Pui, Wi, Qi, Xri−1, Xri, and Eri play the same roles in a limit equilibrium analysis for both a straight slope and a concave slope. Therefore, the only difference between the 3D limit equilibrium analyses of a straight slope and a concave slope is that the force Eθ is not perpendicular to the sliding direction and can provide additional shear resistance on the column element (i, j). The additional shear resistances account for the spatial (3D) effects of the stresses in a concave slope and hence should be considered in practical applications.Figure2(b) shows the column element (i, j) and the corresponding intercolumn forces in the plan view. Figure 2(b) shows that the intercolumn forces Eθi can be decomposed into the component forces Eθri in the direction of sliding and the component forces Eθθi perpendicular to the direction of sliding. If the central angle dθ of the column element (i, j) is sufficiently small, the following relation can be easily obtained:(2)Eθri=Eθi⋅sindθ2≈Eθi⋅dθ2,and the resultant force Pri of the component forces Eθi acting on two sides of the column element (i, j) can be calculated as(3)Pri=2Eθri≈Eθidθ2.Based on the above discussions, the resultant forcePri acting in the opposite direction of movement is defined as the additional shear resistance for the specific column element (i, j). Because of the presence of the additional shear resistances, a concave slope is always more stable than a straight slope.


## 2.1. Geometry of a Concave Slope

Figure1 shows the geometry parameters of a typical axisymmetric concave slope with horizontal layers including (1) the slope height (H), (2) the slope angle (βc), and (3) the radius at the toe of the slope (Rc). In order to conduct a 3D limit equilibrium analysis, the concave slope can be divided into m × n individual vertical columns. This can be easily done by equally dividing the concave slope into m cylinder slices (in the plan view) along the radial direction and n wedge-shaped blocks along the circumferential direction. In the plan view, the overlapping area of the i-th cylinder slice and the j-th wedge-shaped block is designated as column (i, j), which is a fundamental element for the 3D limit equilibrium analysis of an axisymmetric concave slope.Figure 1
Typical geometry of an axisymmetric concave slope (1/4 model shown).


## 2.2. Definition of Additional Shear Resistance

Figure2(a) shows a typical column element (i, j) taken from an axisymmetric concave slope. Prior to defining the additional shear resistance, it is necessary to analyze all the internal and external forces acting on this element. The 3D limit equilibrium analysis of a general axisymmetric slope includes the following forces: (1) the normal force Ni and the shear force Si at the base of the column; (2) the normal force Pui due to pore-water pressure at the base of the column; (3) the external vertical forces Wi and Qi due to soil weight and surcharge load of the slope, respectively; (4) the intercolumn normal forces Eri−1 and Eri in the r-direction and Eθi in the θ-direction; and (5) the vertical intercolumn shear forces Xri−1 and Xri in the r-direction. Note that the subscript j for all these quantities has been omitted considering the inherent axisymmetry of the studied problem.Figure 2
Typical column element and its internal and external forces.
(a)
(b)The forcesNi, Si, Pui, Wi, Qi, Xri−1, Xri, and Eri play the same roles in a limit equilibrium analysis for both a straight slope and a concave slope. Therefore, the only difference between the 3D limit equilibrium analyses of a straight slope and a concave slope is that the force Eθ is not perpendicular to the sliding direction and can provide additional shear resistance on the column element (i, j). The additional shear resistances account for the spatial (3D) effects of the stresses in a concave slope and hence should be considered in practical applications.Figure2(b) shows the column element (i, j) and the corresponding intercolumn forces in the plan view. Figure 2(b) shows that the intercolumn forces Eθi can be decomposed into the component forces Eθri in the direction of sliding and the component forces Eθθi perpendicular to the direction of sliding. If the central angle dθ of the column element (i, j) is sufficiently small, the following relation can be easily obtained:(2)Eθri=Eθi⋅sindθ2≈Eθi⋅dθ2,and the resultant force Pri of the component forces Eθi acting on two sides of the column element (i, j) can be calculated as(3)Pri=2Eθri≈Eθidθ2.Based on the above discussions, the resultant forcePri acting in the opposite direction of movement is defined as the additional shear resistance for the specific column element (i, j). Because of the presence of the additional shear resistances, a concave slope is always more stable than a straight slope.


## 3. Determination of Additional Shear Resistances

### 3.1. Assumptions for Simplification

As discussed previously, a key issue for the stability analysis of an axisymmetric concave slope is the determination of the additional shear resistances. Since the shear resistances arein fact generated by the lateral pressures (stresses) in the circumferential directions of an axisymmetric concave slope, it is necessary to determine the inner circumferential stresses developing in the concave slope so as to more reasonably estimate the additional shear resistances. For this purpose, the axisymmetric concave slope is simplified as a set of thick cylinders (or cylinder slices in the plan view) of different inner and outer diameters. Based on the boundary conditions, the distributions of stresses in these cylinders can be derived via the elastoplastic theory.Figure3 shows the cross-sectional profile of the i-th thick cylinder element in an axisymmetric concave slope, which is also referred to as the i-th cylinder slice. To be clear, the following symbols are given: the height of the i-th cylinder is hi; the inner and outer diameters of the cylinder slice are ri−1 and ri, respectively; the surcharge load applied on the top of the i-th cylinder slice is qi; the average unit weight of the cylinder is γi; the soil pressures acting on the inner and outer faces of the cylinder are pi−1 and pi, respectively; and the vertical shear stresses acting on these two faces are τi−1 and τi, respectively. Before an elastoplastic solution is deduced for the distributions of stresses in the cylinder, some assumptions for simplification of the boundary conditions are proposed:(1)
The soil pressures (normal stresses),pi−1 and pi, acting on the inner and outer faces of the cylinder are axisymmetric and linearly distributed with depth.(2)
The vertical shear stresses,τi−1 and τi, acting on the inner and outer faces of the cylinder, are irrelevant to the calculation of the distribution of the circumferential stress σr, which will be demonstrated later in this paper.(3)
The thickness of the cylinder element is assumed to be sufficiently small; hence, the surcharge loadqi on the top of the cylinder can be treated as a uniform vertical load.Figure 3
Cross-sectional profile of thei-th cylinder element.


### 3.2. Elastoplastic Solution for Additional Shear Resistances

Based on the simplification assumptions, an elastoplastic solution is developed for the distributions of stressesσθ,σr,σz,andτrz in the cylinder element in this section. Then, the solution for the circumferential stress σr is used to calculate the additional shear resistances Pri for a specified vertical column element (i, j) in the limit equilibrium analysis.The stress equilibrium equation for thei-th cylinder can be written as(4)∂σr∂r+σr−σθr+∂τrz∂z=0,∂σz∂z+∂τrz∂r+τrzr−γi=0,and the simplified boundary conditions are expressed as(5)r=ri−1:σr=ki−1z,τrz=τi−1,r=ri:σr=kiz,τrz=τi,z=0:σz=qi,where ki−1 and ki are the gradients of the soil pressures acting on the inner and outer faces of the cylinder, respectively; τi−1 and τi are the vertical shear stresses acting on the inner and outer faces of the cylinder, respectively; and qi is the vertical surcharge load acting on top of the cylinder. For a specific column element (i, j), these quantities can be calculated using the following relations:(6)ki−1=2Eri−1hi2⋅ri−1dθ,(7)ki=2Erihi2⋅ridθ,(8)τi−1=Xri−1hi⋅ri−1dθ,(9)τi=Xri,jhi⋅ridθ,(10)qi=2Qibiri+ri−1dθ,where Eri−1 and Eri are the normal forces acting on the inner and outer faces of the element column (i, j); Xri−1 and Xri are the vertical shear forces acting on the inner and outer faces of the element column (i, j); Qi is the vertical surcharge load acting on the upper surface of the column element (i, j); bi is the thickness of the column element (i, j); and dθ is the central angle of the column element (i, j).According to the elastoplastic theory, the solution of equation (4) can be decomposed into a characteristic solution without a gravity force and a special solution with a gravity force. Therefore, a general solution is derived by the following equation:(11)∂σr∂r+σr−σθr+∂τrz∂z=0,∂σz∂z+∂τrz∂r+τrzr=0.Based on Love’s method [18], the stress components in equation (11) can be expressed as(12)σr=∂∂zμ∇2ϕ−∂2ϕ∂r2,(13)σθ=∂∂zμ∇2ϕ−1r∂ϕ∂r,(14)σz=∂∂z2−μ∇2ϕ−∂2ϕ∂z2,(15)τrz=∂∂r1−μ∇2ϕ−∂2ϕ∂z2,where σθ,σr,σz,andτrz are the stress components; ϕ is a potential function; and μ is Poisson’s ratio. The potential function ϕ must satisfy the following biharmonic condition:(16)∇2∇2ϕ=0,where ∇2 is the Laplacian operator.For a spatial axisymmetric problem, the potential function (ϕ) can be defined as(17)ϕ=A1z4+A2r4+A3z2r2+A4z3+A5z2lnr+A6zr2+A7r2lnr,where A1 to A7 are unknown constants. Combining equations (12)–(17) with the stress boundary condition (5), the constants A1 to A7 can be derived as(18)A1=2μ−112μ+1τi−1ri−1−τiriri−12−ri2+2−μ12μ+1kiri2−ki−1ri−12ri2−ri−12,(19)A2=2μ+1321+μτi−1ri−1−τiriri−12−ri2−μ321+μkiri2−ki−1ri−12ri2−ri−12,(20)A3=1−2μ61+μqi,(21)A4=−2μ4μ+1τi−1ri−1−τiriri−12−ri2−1−μ4μ+1kiri2−ki−1ri−12ri2−ri−12,(22)A5=ki−1−kiri−1,j2ri22ri2−ri−12,(23)A6=μ21+μqi,(24)A7=ki−1−kiri−12ri2μ+τi−1ri−τiri−1ri−1ri41−μri2−ri−12.By substituting the potential functionϕ in equation (17) and the constants in equations (18)–(24) into equations (12)–(15), the stress components can be further expressed as(25)σθ=kiri2−ki−1ri−12ri2−ri−12z−ki−1−kiri−12ri2zri2−ri−12r2,(26)σr=kiri2−ki−1ri−12ri2−ri−12z+ki−1−kiri−12ri2zri2−ri−12r2,(27)τrz=τi−iri−1−τiriri−12−ri2r+τi−1ri−τiri−1ri2−ri−12rri−1ri,(28)σz=qi−2τi−1ri−1−τiriri−12−ri2z.From equations (25) and (26), it can be easily seen that the shear stresses τi−1 and τi are irrelevant to the stresses σθ and σr, which is consistent with Assumption (3) in Section 3.1.Consider the field of gravity as a special solution of equation (4), the vertical component of which is(29)σz=γiz,and the general solution for the stress σz can be expressed as(30)σz=qi−2τi−1ri−1−τiriri−12−ri2z+γiz.Considering the soil satisfies the Mohr–Coulomb failure criterion, the failure condition can be expressed as(31)σθ−σr2=cz′tanφz′+σrsinφz′1−sinφz′,where cz′ and φz′ denote the effective cohesion and effective friction angle of the soil at depth z, respectively.From equations (25) and (26), it can be seen that a smaller value of r results in a lower value of σr and a higher value of σθ. Based on this fact and equation (31), the plastic zone first occurs at the radial distance of r = ri−1. Therefore, on the incipient failure of the concave slope, the plastic stresses σθ∗ and σr∗ can be deduced from equations (25) and (26) at the radial distance of r = ri−1 and be further simplified into(32)σθ∗=zri2−ri−122kiri2−ki−1ri−12−ki−1ri2,(33)σr∗=ki−1z.Substituting equations (32) and (33) into equation (31) leads to(34)ki−ki−1ri2ri2−ri−12z=cz′tanφz′+ki−1zsinφz′1−sinφz′.After some necessary manipulations, equation (34) can be transformed into the following relation:(35)kiz=ri2−ri−12ri2ki−1zri2ri2−ri−12+sinφz′1−sinφz′+cz′tanφz′sinφz′1−sinφz′.Finally, substituting equation (35) into equation (32) leads to(36)σθ∗=ki−1z1+sinφz′1−sinφz′+2ci′cosφz′1−sinφz′.Note that equation (36) is a linear expression with respect to depth z. Therefore, the lateral forces Eθi acting on two sides of the column element (i, j) can be calculated after combining equation (36) with equation (6):(37)Eθi=12⋅bi⋅hi⋅σθ∗z=hi=12⋅bi⋅hi⋅2Erihi⋅ridθ1+sinφz′1−sinφz′+2cz′cosφz′1−sinφz′,and the force Eθi acts on 1/3 the height of the column element (i, j).Finally, considering equation (3), the previously defined additional shear resistance Pri can be calculated as(38)Pri=12⋅bi⋅hi⋅2Erihi⋅ri1+sinφz′1−sinφz′+2ci′cosφz′1−sinφz′dθ.


## 3.1. Assumptions for Simplification

As discussed previously, a key issue for the stability analysis of an axisymmetric concave slope is the determination of the additional shear resistances. Since the shear resistances arein fact generated by the lateral pressures (stresses) in the circumferential directions of an axisymmetric concave slope, it is necessary to determine the inner circumferential stresses developing in the concave slope so as to more reasonably estimate the additional shear resistances. For this purpose, the axisymmetric concave slope is simplified as a set of thick cylinders (or cylinder slices in the plan view) of different inner and outer diameters. Based on the boundary conditions, the distributions of stresses in these cylinders can be derived via the elastoplastic theory.Figure3 shows the cross-sectional profile of the i-th thick cylinder element in an axisymmetric concave slope, which is also referred to as the i-th cylinder slice. To be clear, the following symbols are given: the height of the i-th cylinder is hi; the inner and outer diameters of the cylinder slice are ri−1 and ri, respectively; the surcharge load applied on the top of the i-th cylinder slice is qi; the average unit weight of the cylinder is γi; the soil pressures acting on the inner and outer faces of the cylinder are pi−1 and pi, respectively; and the vertical shear stresses acting on these two faces are τi−1 and τi, respectively. Before an elastoplastic solution is deduced for the distributions of stresses in the cylinder, some assumptions for simplification of the boundary conditions are proposed:(1)
The soil pressures (normal stresses),pi−1 and pi, acting on the inner and outer faces of the cylinder are axisymmetric and linearly distributed with depth.(2)
The vertical shear stresses,τi−1 and τi, acting on the inner and outer faces of the cylinder, are irrelevant to the calculation of the distribution of the circumferential stress σr, which will be demonstrated later in this paper.(3)
The thickness of the cylinder element is assumed to be sufficiently small; hence, the surcharge loadqi on the top of the cylinder can be treated as a uniform vertical load.Figure 3
Cross-sectional profile of thei-th cylinder element.


## 3.2. Elastoplastic Solution for Additional Shear Resistances

Based on the simplification assumptions, an elastoplastic solution is developed for the distributions of stressesσθ,σr,σz,andτrz in the cylinder element in this section. Then, the solution for the circumferential stress σr is used to calculate the additional shear resistances Pri for a specified vertical column element (i, j) in the limit equilibrium analysis.The stress equilibrium equation for thei-th cylinder can be written as(4)∂σr∂r+σr−σθr+∂τrz∂z=0,∂σz∂z+∂τrz∂r+τrzr−γi=0,and the simplified boundary conditions are expressed as(5)r=ri−1:σr=ki−1z,τrz=τi−1,r=ri:σr=kiz,τrz=τi,z=0:σz=qi,where ki−1 and ki are the gradients of the soil pressures acting on the inner and outer faces of the cylinder, respectively; τi−1 and τi are the vertical shear stresses acting on the inner and outer faces of the cylinder, respectively; and qi is the vertical surcharge load acting on top of the cylinder. For a specific column element (i, j), these quantities can be calculated using the following relations:(6)ki−1=2Eri−1hi2⋅ri−1dθ,(7)ki=2Erihi2⋅ridθ,(8)τi−1=Xri−1hi⋅ri−1dθ,(9)τi=Xri,jhi⋅ridθ,(10)qi=2Qibiri+ri−1dθ,where Eri−1 and Eri are the normal forces acting on the inner and outer faces of the element column (i, j); Xri−1 and Xri are the vertical shear forces acting on the inner and outer faces of the element column (i, j); Qi is the vertical surcharge load acting on the upper surface of the column element (i, j); bi is the thickness of the column element (i, j); and dθ is the central angle of the column element (i, j).According to the elastoplastic theory, the solution of equation (4) can be decomposed into a characteristic solution without a gravity force and a special solution with a gravity force. Therefore, a general solution is derived by the following equation:(11)∂σr∂r+σr−σθr+∂τrz∂z=0,∂σz∂z+∂τrz∂r+τrzr=0.Based on Love’s method [18], the stress components in equation (11) can be expressed as(12)σr=∂∂zμ∇2ϕ−∂2ϕ∂r2,(13)σθ=∂∂zμ∇2ϕ−1r∂ϕ∂r,(14)σz=∂∂z2−μ∇2ϕ−∂2ϕ∂z2,(15)τrz=∂∂r1−μ∇2ϕ−∂2ϕ∂z2,where σθ,σr,σz,andτrz are the stress components; ϕ is a potential function; and μ is Poisson’s ratio. The potential function ϕ must satisfy the following biharmonic condition:(16)∇2∇2ϕ=0,where ∇2 is the Laplacian operator.For a spatial axisymmetric problem, the potential function (ϕ) can be defined as(17)ϕ=A1z4+A2r4+A3z2r2+A4z3+A5z2lnr+A6zr2+A7r2lnr,where A1 to A7 are unknown constants. Combining equations (12)–(17) with the stress boundary condition (5), the constants A1 to A7 can be derived as(18)A1=2μ−112μ+1τi−1ri−1−τiriri−12−ri2+2−μ12μ+1kiri2−ki−1ri−12ri2−ri−12,(19)A2=2μ+1321+μτi−1ri−1−τiriri−12−ri2−μ321+μkiri2−ki−1ri−12ri2−ri−12,(20)A3=1−2μ61+μqi,(21)A4=−2μ4μ+1τi−1ri−1−τiriri−12−ri2−1−μ4μ+1kiri2−ki−1ri−12ri2−ri−12,(22)A5=ki−1−kiri−1,j2ri22ri2−ri−12,(23)A6=μ21+μqi,(24)A7=ki−1−kiri−12ri2μ+τi−1ri−τiri−1ri−1ri41−μri2−ri−12.By substituting the potential functionϕ in equation (17) and the constants in equations (18)–(24) into equations (12)–(15), the stress components can be further expressed as(25)σθ=kiri2−ki−1ri−12ri2−ri−12z−ki−1−kiri−12ri2zri2−ri−12r2,(26)σr=kiri2−ki−1ri−12ri2−ri−12z+ki−1−kiri−12ri2zri2−ri−12r2,(27)τrz=τi−iri−1−τiriri−12−ri2r+τi−1ri−τiri−1ri2−ri−12rri−1ri,(28)σz=qi−2τi−1ri−1−τiriri−12−ri2z.From equations (25) and (26), it can be easily seen that the shear stresses τi−1 and τi are irrelevant to the stresses σθ and σr, which is consistent with Assumption (3) in Section 3.1.Consider the field of gravity as a special solution of equation (4), the vertical component of which is(29)σz=γiz,and the general solution for the stress σz can be expressed as(30)σz=qi−2τi−1ri−1−τiriri−12−ri2z+γiz.Considering the soil satisfies the Mohr–Coulomb failure criterion, the failure condition can be expressed as(31)σθ−σr2=cz′tanφz′+σrsinφz′1−sinφz′,where cz′ and φz′ denote the effective cohesion and effective friction angle of the soil at depth z, respectively.From equations (25) and (26), it can be seen that a smaller value of r results in a lower value of σr and a higher value of σθ. Based on this fact and equation (31), the plastic zone first occurs at the radial distance of r = ri−1. Therefore, on the incipient failure of the concave slope, the plastic stresses σθ∗ and σr∗ can be deduced from equations (25) and (26) at the radial distance of r = ri−1 and be further simplified into(32)σθ∗=zri2−ri−122kiri2−ki−1ri−12−ki−1ri2,(33)σr∗=ki−1z.Substituting equations (32) and (33) into equation (31) leads to(34)ki−ki−1ri2ri2−ri−12z=cz′tanφz′+ki−1zsinφz′1−sinφz′.After some necessary manipulations, equation (34) can be transformed into the following relation:(35)kiz=ri2−ri−12ri2ki−1zri2ri2−ri−12+sinφz′1−sinφz′+cz′tanφz′sinφz′1−sinφz′.Finally, substituting equation (35) into equation (32) leads to(36)σθ∗=ki−1z1+sinφz′1−sinφz′+2ci′cosφz′1−sinφz′.Note that equation (36) is a linear expression with respect to depth z. Therefore, the lateral forces Eθi acting on two sides of the column element (i, j) can be calculated after combining equation (36) with equation (6):(37)Eθi=12⋅bi⋅hi⋅σθ∗z=hi=12⋅bi⋅hi⋅2Erihi⋅ridθ1+sinφz′1−sinφz′+2cz′cosφz′1−sinφz′,and the force Eθi acts on 1/3 the height of the column element (i, j).Finally, considering equation (3), the previously defined additional shear resistance Pri can be calculated as(38)Pri=12⋅bi⋅hi⋅2Erihi⋅ri1+sinφz′1−sinφz′+2ci′cosφz′1−sinφz′dθ.


## 4. Simplified 2D Limit Equilibrium Procedure for Stability Analysis of Concave Slopes considering Additional Shear Resistance

### 4.1. Equilibrium Conditions

A variety of 2D LEMs can be used to estimate the slope stability, such as Fellenius’ method, Bishop’s simplified method, Janbu’s simplified method, Spencer’s method, and Morgenstern–Price’s method. By incorporating additional shear resistance measured by equation (38) with any of these 2D LEMs, the stability of concave slopes can be assessed. In this section, Bishop’s method is taken as an example to illustrate the simplified 2D limit equilibrium procedure.Figure2 illustrates that a typical concave slope can be equally divided into n wedge-shaped blocks along the circumferential direction. Figure 4 shows the j-th wedge-shaped block extracted from the concave slope, which can be subsequently divided into m column elements.Figure 4
Thej-th wedge-shaped block.With this configuration, the equilibrium conditions of forces and moments for the wedge-shaped block are established as follows.After a failure surface is assumed, the weight of the column element (i, j) can be expressed as(39)Wi=12γihi−lisinαilicosαiri−1+ri2dθ,where αi is the inclination angle of the bottom face of the column element (i, j) to the horizontal plane and li is the length of the bottom face of the column element (i, j).The normal force caused by pore-water pressure at the base of the column element (i, j) is calculated as(40)Pui=Uiliri−1+ri2dθ,where Ui is the average pore water pressure.The shear resistanceSi due to soil cohesion and friction at the base of the column element (i, j) is calculated as(41)Si=PciF+Ni−Puitanφi′F,where F is the factor of safety, φi′ denotes the effective friction angle of the soil in the failure surface of the column element (i, j), and Pci is the shear resistance due to soil cohesion, which can be calculated as(42)Pci=ci′liri−1+ri2dθ,where ci′ denotes the effective cohesion of the soil in the failure surface of the column element (i, j). In equation (41), Ni is the normal force acting at the base of the column element (i, j). The quantity Ni can be determined by the force equilibrium in the z-direction as follows:(43)Wi−Sisinαi−Nicosαi+Xri−1−Xri=0.As in the Bishop algorithm, it is assumed that(44)Xri−1−Xri=0.Hence, combining equations (42)–(44) leads to the expression of Ni as follows:(45)Ni=1mαiWi−PcisinαiF+Puitanφi′sinαiF,where the parameter mαi can be calculated by(46)mαi=cosαi+tanφi′sinαiF.Substituting equation (45) into equation (41) leads to(47)Si=1F⋅mαiPcicosαi+tanφi′Wi−Puicosαi.Next, the force equilibrium in the radial (r) direction leads to(48)Eri=Eri−1+Sicosαi−Nisinαi+Pri,where Eri−1 and Eri are the lateral forces acting on the inner and outer faces of the column element (i, j) and Pri is the previously defined additional shear resistance, which acts at the height of hi/3 from the bottom face of the column element (i, j). Note the innermost lateral force Er0 at the toe of the concave slope is zero.Finally, the overall moment equilibrium of the failure body leads to(49)F=fF,Pri=∑R/mαiPci⋅cosαi+Wi−Pui⋅cosαitanφi′∑Wi⋅Rsinαi−∑Pri⋅Rri,where R is the distance from the bottom face of the column element (i, j) to the rotation centre O and Rri is the distance from the action line of the additional shear resistances Pri to the rotation centre O, which is calculated as(50)Rri=Rsinαi2+Rcosαi−hi32.


### 4.2. Iteration Algorithm for Solving Critical FoS

To solve the critical value of the FoS efficiently, an iteration algorithm is proposed in this section as follows:Step 1.
A potential failure surface is assumed and thej-th wedge-shaped block is subdivided into a series of column elements, as shown in Figure 4.Step 2.
The initial values of the additional shear resistances are set to zero, i.e.,Pri = 0, and an initial estimate of F is given, which is typically set to 1.Step 3.
As in Bishop’s algorithm, the value ofF is updated repeatedly using equations (39) to (40), (42), and (45)–(50) until the difference between the values of F obtained by the last two iterations is smaller than a certain tolerance (e.g., dF<0.01).Step 4.
The values ofEri,j are extracted from the analysis in Step 3, and then these values are substituted into equation (38) to update the values of the additional shear resistance Pri,j.Step 5.
Based on the updated values ofPri,j, repeat Steps 3 and 4 until the values of the additional shear resistances Pri,j converge to fixed values as the Euclid norm ||Pri,j|| is less than a specific tolerance (e.g., pri,j<0.001).Step 6.
The minimum value ofF is found by repeating Steps 1 to 5 for all possible failure surfaces.The minimum value ofF determined based on the above-mentioned steps is taken as the critical FoS for a concave slope.It is worth mentioning that the proposed method can be seen as a generalized traditional 2D limit equilibrium approach with the additional shear resistances considered in the analysis. The following section will demonstrate that this approach is simple and efficient for practical uses.


### 4.3. Verification of the Proposed Method

Based on the strength reduction technique incorporated in the numerical software FLAC3D, Zettler et al. [19] analyzed two 25-meter high slopes: (1) an axisymmetric concave slope with a 12-meter radius of curvature at the toe of the slope and (2) a straight slope. Table 1 lists the material and geometry parameters of the concave slope. For comparison, Table 2 lists the calculated results of the FoS obtained by Zettler et al. [19] and the proposed method, together with those from the Bishop algorithm for the straight slope.Table 1
Geometry and material parameters of the concave slope [19].
Material parameters
Geometry parameters
E (MPa)
υ
c (kPa)
φ (°)
25 kN/m3
tan(βc)
H (m)
R
c (m)
257
0.29
38
45
2500
2 : 1
25
12Table 2
Comparison of the FoS obtained by different methods.
Slope type
R
c/H
Methods
Proposed
Bishop
FLAC3D
Straight
∞
1.27
1.27
1.37
Concave
0.48
1.94
—
1.83Table2 shows that the results of the proposed method and Bishop’s method agree well for the straight slope. This is because no additional shear resistances are actually considered in the proposed method for a straight slope. In this case, the proposed method is basically the same as Bishop’s method. For the straight slope, the FoS calculated by FLAC3D using the strength reduction technique is 1.37. The difference in the results (FoS) between the proposed method and FLAC3D is less than 8%, which is small considering that these two methods are quite different. For the concave slope, the results of the FoS from the proposed method and FLAC3D are 1.94 and 1.83, respectively. Their difference is also less than 8%. The above discussions indicate that the proposed method is applicable for analyzing concave slopes in practical applications.


## 4.1. Equilibrium Conditions

A variety of 2D LEMs can be used to estimate the slope stability, such as Fellenius’ method, Bishop’s simplified method, Janbu’s simplified method, Spencer’s method, and Morgenstern–Price’s method. By incorporating additional shear resistance measured by equation (38) with any of these 2D LEMs, the stability of concave slopes can be assessed. In this section, Bishop’s method is taken as an example to illustrate the simplified 2D limit equilibrium procedure.Figure2 illustrates that a typical concave slope can be equally divided into n wedge-shaped blocks along the circumferential direction. Figure 4 shows the j-th wedge-shaped block extracted from the concave slope, which can be subsequently divided into m column elements.Figure 4
Thej-th wedge-shaped block.With this configuration, the equilibrium conditions of forces and moments for the wedge-shaped block are established as follows.After a failure surface is assumed, the weight of the column element (i, j) can be expressed as(39)Wi=12γihi−lisinαilicosαiri−1+ri2dθ,where αi is the inclination angle of the bottom face of the column element (i, j) to the horizontal plane and li is the length of the bottom face of the column element (i, j).The normal force caused by pore-water pressure at the base of the column element (i, j) is calculated as(40)Pui=Uiliri−1+ri2dθ,where Ui is the average pore water pressure.The shear resistanceSi due to soil cohesion and friction at the base of the column element (i, j) is calculated as(41)Si=PciF+Ni−Puitanφi′F,where F is the factor of safety, φi′ denotes the effective friction angle of the soil in the failure surface of the column element (i, j), and Pci is the shear resistance due to soil cohesion, which can be calculated as(42)Pci=ci′liri−1+ri2dθ,where ci′ denotes the effective cohesion of the soil in the failure surface of the column element (i, j). In equation (41), Ni is the normal force acting at the base of the column element (i, j). The quantity Ni can be determined by the force equilibrium in the z-direction as follows:(43)Wi−Sisinαi−Nicosαi+Xri−1−Xri=0.As in the Bishop algorithm, it is assumed that(44)Xri−1−Xri=0.Hence, combining equations (42)–(44) leads to the expression of Ni as follows:(45)Ni=1mαiWi−PcisinαiF+Puitanφi′sinαiF,where the parameter mαi can be calculated by(46)mαi=cosαi+tanφi′sinαiF.Substituting equation (45) into equation (41) leads to(47)Si=1F⋅mαiPcicosαi+tanφi′Wi−Puicosαi.Next, the force equilibrium in the radial (r) direction leads to(48)Eri=Eri−1+Sicosαi−Nisinαi+Pri,where Eri−1 and Eri are the lateral forces acting on the inner and outer faces of the column element (i, j) and Pri is the previously defined additional shear resistance, which acts at the height of hi/3 from the bottom face of the column element (i, j). Note the innermost lateral force Er0 at the toe of the concave slope is zero.Finally, the overall moment equilibrium of the failure body leads to(49)F=fF,Pri=∑R/mαiPci⋅cosαi+Wi−Pui⋅cosαitanφi′∑Wi⋅Rsinαi−∑Pri⋅Rri,where R is the distance from the bottom face of the column element (i, j) to the rotation centre O and Rri is the distance from the action line of the additional shear resistances Pri to the rotation centre O, which is calculated as(50)Rri=Rsinαi2+Rcosαi−hi32.


## 4.2. Iteration Algorithm for Solving Critical FoS

To solve the critical value of the FoS efficiently, an iteration algorithm is proposed in this section as follows:Step 1.
A potential failure surface is assumed and thej-th wedge-shaped block is subdivided into a series of column elements, as shown in Figure 4.Step 2.
The initial values of the additional shear resistances are set to zero, i.e.,Pri = 0, and an initial estimate of F is given, which is typically set to 1.Step 3.
As in Bishop’s algorithm, the value ofF is updated repeatedly using equations (39) to (40), (42), and (45)–(50) until the difference between the values of F obtained by the last two iterations is smaller than a certain tolerance (e.g., dF<0.01).Step 4.
The values ofEri,j are extracted from the analysis in Step 3, and then these values are substituted into equation (38) to update the values of the additional shear resistance Pri,j.Step 5.
Based on the updated values ofPri,j, repeat Steps 3 and 4 until the values of the additional shear resistances Pri,j converge to fixed values as the Euclid norm ||Pri,j|| is less than a specific tolerance (e.g., pri,j<0.001).Step 6.
The minimum value ofF is found by repeating Steps 1 to 5 for all possible failure surfaces.The minimum value ofF determined based on the above-mentioned steps is taken as the critical FoS for a concave slope.It is worth mentioning that the proposed method can be seen as a generalized traditional 2D limit equilibrium approach with the additional shear resistances considered in the analysis. The following section will demonstrate that this approach is simple and efficient for practical uses.


## 4.3. Verification of the Proposed Method

Based on the strength reduction technique incorporated in the numerical software FLAC3D, Zettler et al. [19] analyzed two 25-meter high slopes: (1) an axisymmetric concave slope with a 12-meter radius of curvature at the toe of the slope and (2) a straight slope. Table 1 lists the material and geometry parameters of the concave slope. For comparison, Table 2 lists the calculated results of the FoS obtained by Zettler et al. [19] and the proposed method, together with those from the Bishop algorithm for the straight slope.Table 1
Geometry and material parameters of the concave slope [19].
Material parameters
Geometry parameters
E (MPa)
υ
c (kPa)
φ (°)
25 kN/m3
tan(βc)
H (m)
R
c (m)
257
0.29
38
45
2500
2 : 1
25
12Table 2
Comparison of the FoS obtained by different methods.
Slope type
R
c/H
Methods
Proposed
Bishop
FLAC3D
Straight
∞
1.27
1.27
1.37
Concave
0.48
1.94
—
1.83Table2 shows that the results of the proposed method and Bishop’s method agree well for the straight slope. This is because no additional shear resistances are actually considered in the proposed method for a straight slope. In this case, the proposed method is basically the same as Bishop’s method. For the straight slope, the FoS calculated by FLAC3D using the strength reduction technique is 1.37. The difference in the results (FoS) between the proposed method and FLAC3D is less than 8%, which is small considering that these two methods are quite different. For the concave slope, the results of the FoS from the proposed method and FLAC3D are 1.94 and 1.83, respectively. Their difference is also less than 8%. The above discussions indicate that the proposed method is applicable for analyzing concave slopes in practical applications.


## 5. Parametric Study

This section presents a parametric study performed to investigate the influence of some important parameters on the stability of a concave slope, including the ratio of the radius of curvature at the toe of a concave slope to the slope heightRc/H, the slope angle βc, and the strength parameters soil cohesion c′ and soil friction angle φ′.Figure5 shows the variations of the calculated FoS with respect to the ratio Rc/H. For a comprehensive comparison, the results of the FoS obtained from FLAC3D are also shown in Figure 5. It needs to be mentioned that other necessary parameters used in FLAC3D are the same as those in Table 1. Figure 5 shows that the maximum error between the calculated FoS using the proposed method and that obtained by FLAC3D is 9.7%, and the error decreases with an increase of Rc/H. Considering the discrepancy between the two approaches, the difference in results is acceptable.Figure 5
Variations of FoS with respect to the ratioRc/H.Figure5 also shows that when Rc/H is greater than 10, the FoS calculated by the proposed method gradually decreases to the FoS given by Bishop’s method (1.274), indicating that a concave slope can be treated as an infinite long straight slope as long as the ratio Rc/H is sufficiently large.Figure6 shows the critical failure planes predicted by the proposed method with different Rc/H. As the Rc/H ratio decreases, the position of the critical failure planes moves upwards and their corresponding FoS increases. This phenomenon may be explained by the fact that the hoop stress effect is more apparent for a concave slope with a higher Rc/H, and hence, the stability of a concave slope is enhanced.Figure 6
Critical failure planes for differentRc/H ratios.To investigate the effects of the slope angleβc and the strength parameters c′ and φ′, a conversion factor Cf is defined as follows:(51)Cf=FoScFoSs,where FoSc and FoSs correspond to the factors of safety for a concave slope and a straight slope with the same characteristics, except for the analyzed parameter.Figure7 shows the variation of the conversion factor Cf with respect to Rc/H for different slope angles βc. It is clearly shown that a concave slope is more stable than a straight slope because Cf is always greater than 1.0. However, the contrast of stability between a concave slope and a straight slope becomes more apparent for steeper slopes (with higher values of tan βc), especially when Rc/H is less than 0.3. Also, it is interesting to note that the slope angle βc has an opposite impact on the factor Cf for the cases with Rc/H smaller or greater than 0.6, indicating that an obvious combined effect exists between the parameters βc and Rc/H. However, the reason for this phenomenon remains to be further discussed.Figure 7
C
f versus Rc/H for different slope angles.Figures8 and 9 show the variations of Cf with respect to Rc/H for different soil cohesions and friction angles, respectively. Figure 8 shows that the factor Cf increases as the soil cohesion increases, implying that the spatial effect in a concave slope is more apparent when the soil cohesion has a higher value. In contrast, Figure 9 shows that the factor Cf decreases as the soil friction angle increases, implying that the spatial effect in a concave slope weakens as the soil friction angle increases.Figure 8
C
f versus the Rc/H ratio for different soil cohesions.Figure 9
C
f versus Rc/H for different soil friction angles.


## 6. Conclusions

Because of its simplicity and efficiency, the 2D limit equilibrium method is still the most widely accepted approach for slope stability analysis in practice. However, the use of a 2D limit equilibrium method for stability analysis of concave slopes may lead to an inaccurate estimate of the FoS unless the additional shear resistances are considered in the analysis.This paper proposes an elastoplastic solution for calculating the additional shear resistances due to the spatial effects of stresses in an axisymmetric concave slope. Considering the additional shear resistances, a 2D limit equilibrium procedure combined with an efficient iteration algorithm is developed to calculate the FoS of concave slopes.The comparison of the calculated FoS from the proposed method and the numerical software FLAC3D demonstrated that the proposed method is accurate and efficient for stability analyses of concave slopes. Based on the proposed method, a parametric study was performed to study different effects of the parameters Rc,βc,c′,andφ′ on the stability of concave slopes, indicating that the 3D effects of stresses must be considered in the 2D analysis to give a reasonable estimate of the FoS of axisymmetric concave slopes.



---

*Source: 8491637-2019-07-21.xml*