# Rigorous Solution of Slopes’ Stability considering Hydrostatic Pressure

Authors: Chengchao Li; Pengming Jiang; Aizhao Zhou
Journal: Complexity
 (2018)

---

## Abstract

According to characteristics of soils in failure, a sliding mechanism of slopes in limit state is divided into five parts, for building a slip line field satisfying all possible boundary conditions. An algorithm is built to obtain the rigorous solution approaching upper and lower bound values simultaneously, which satisfies the static boundary and the kinematical boundary based on the slip line field, while stress discontinuity line and velocity discontinuity line are key points. This algorithm is copared with the Spencer method to prove its feasibility with a special example. The variation of rigorous solution, including an ultimate load and a sliding belt the rigid body sliding along rather than a single slip surface for friction-type soils, is achieved considering hydrostatic pressure with soil parameters changing.


---

## Body

## 1. Introduction

The stability of slopes has been regarded as a classic and difficult problem for engineers because of less boundary constraints, compared with the earth pressure of retaining wall and the bearing capacity of foundation. Over the past few years, many investigators have evaluated slope stability, thereby developing many methods for meeting engineering requirements, such as the limit equilibrium method (LEM) [1, 2], the finite element method (FEM) [3, 4], and the limit analysis method (LAM) [5–7]. The LEM captures the static equilibrium of rigid blocks on a particular slip surface, while not considering the plastic deformation of soils. The equilibrium equation accounts for the whole slice but not guarantees each point in soils. The strength reduction method (SRM) [8–10] is the main finite element slope stability method currently employed, by which stress field and displacement of soils in slopes can be calculated with an elastic-plastic constitutive model to get the safety factor. Although the displacement mutation, the numerical calculation of nonconvergence, and other criteria can determine the slope instability, the slope displacement calculated by the SRM is not the actual displacement of soils and the safety factor required is an approximation. For the slope stability, it is sometimes not necessary to get the variation of stress field and strain field, only to get the ultimate load. Based on the extremum principal [11], the lower bound (LB) [12, 13] solution can be got by static analysis for limit equilibrium problems and the upper bound (UB) [14–17] solution can be got by dynamic analysis. If the lower bound solution satisfies all the kinematic conditions or the upper bound solution satisfies all the static conditions, the solution will be the rigorous one. Compared with obtaining the safety factor, the solution calculated by the upper and lower bound theorems is closer to the real condition because the ultimate load approaches the upper bound solution and the lower bound solution at the same time, satisfying all possible boundary conditions.When a slope is on the verge of collapse, forc−φ type, a sliding belt is emerged within the slope due to the friction between soils. And the sliding belt is not a single slip surface, but a thin shear zone. This paper builds the slip line field satisfying the static boundary condition and the velocity boundary condition according to the characteristics of the stress discontinuity line and the velocity discontinuity line and compiles an algorithm to gain the distribution of the sliding belt and the ultimate capacity of the slope considering the hydrostatic pressure. It also indicates the variation of the ultimate bearing capacity and the sliding belt with different parameters. The application of this algorithm is proved by comparing with the Spencer’s method.


## 2. Slip Line Field

To analyze the slope stability based on the LAM, the slip line field is constructed according to the stress boundary conditions; one of which is a noncharacteristic line stress boundary with normal stress and shear normal stress, and the other one is the interface of a rigid region and a plastic region or a stress discontinuity surface [18]. Under the limit state, only the noncharacteristic line boundary can be used to construct the slip line field, while the other is unknown. Basic coordinate and noncharacteristic line stress boundary is shown in Figure 1. θ is the angle between the direction of the maximum principal stress σ1, and the y-axis and Γ is the noncharacteristic line stress boundary. φ1−π<φ1≤π is the angle between the tangential direction t of the boundary Γ and the y-axis, so the angle between t and σ1 is θ−φ1. According to σn and τn, two Mohr’s circles tangent to the Coulomb failure line can be drawn (Figure 2). When soils are in the extreme state, σt, σn, and τn on the boundary Γ are expressed, respectively, as follows:
(1.1)σt=p±Rcos2θ−φ1,(1.2)σn=p∓Rcos2θ−φ1,(1.3)τn=Rsin2θ−φ1,(1.4)R=p+σcsinφ,where σn, τn, and σt are normal stress, shear stress, and tangent normal stress on the boundary Γ, respectively. p is average stress; R is radius of Mohr’s circle; c is soil cohesion; and φ is internal friction angle. σc is cohesive internal stress, which is given by: σc=ccotφ.Figure 1
Basic coordinate system and noncharacteristic line stress boundary.Figure 2
Mohr’s circles on noncharacteristic line stress boundary.BecauseM is a point on the Mohr’s circle, R can be rewritten as follows:
(2)R=σn−p2+τn2.According to (1.4) and (2):
(3.1)p2cos2φ−2σn+csinφcosφp+σn2+τn2−c2cos2φ=0.Solving the unary quadratic equation ofp:
(3.2)p1,2=sec2φσn+csinφcosφ±σnsinφ−ccosφ2−τn2cos2φ.Whenσt>σn, according to (1.3):
(4)θ+=φ1+12arcsinτnR=φ1+12arcsinτnp1+σcsinφ.Whenσt<σn, (1.3) should be rewritten as follows:
(5)τn=Rsinπ−2θ−φ1=Rsin2π2−θ−φ1,θ−=φ1−12arcsinτnR+π2=φ1−12arcsinτnp2+σcsinφ+π2,where θ+ and θ− are the maximum principal direction angles.According to the above analysis, the analytic expressions of noncharacteristic line stress boundary conditions can be uniformly written as follows:(6)p1,2=sec2φσn+csinφcosφ±σnsinφ+ccosφ2−τn2cos2φ,θ+,−=φ1±12arcsinτnp1,2+σcsinφ∓π4+π4,σt1,2=σn±2p1,2+σc2sin2φ−τn2.Withθ on the noncharacteristic line stress boundary, (6) can be used to draw the stress line field near the stress boundary:
(7)α=θ−μalonglineα,β=θ+μalonglineβ,where α and β are the angles between stress characteristic lines through the point M and the y-axis; μ is the angle between the stress characteristic line and the direction of major principal stress and given as follows:
(8)μ=π4−φ2.In this case, stress boundary conditions are achieved. In order to facilitate the following calculations, new variablesσe and θ are used instead of σx, σy, and τxy (Figure 3). σe is the effective stress, and θ is the direction angle of the maximum principal stress, written as the (9). In addition, σe and θ of each point in slip line field can be calculated by solving the basic boundary value problems. In the theory of hyperbolic partial differential equations, there are three basic boundary value problems [19, 20], of which Cauchy problem is the key issue the others can be solved by. Here, only the Cauchy problem is introduced and the rest of the problems (Riemann problem and mixed boundary value problem) can be referred to the literature [19]. Taking into account of the hydrostatic pressure, the soils above the ground water take unit natural weight γ and the soils below ground water take unit floating weight γ′.
(9)σe=12σx+σy+σc,σxσy=σe1±sinφcos2θ−σc,τxy=σesinφsin2θ.Figure 3
Effective stress of Mohr’s circle.As shown in Figure4(a), OA is a smooth and continuous stress boundary line, a nonslip line, in which the coordinates x,y of each point, effective stress σe, and the direction angle θ are known. By solving the Cauchy problem, the distribution of slip line field in the triangular OAB surrounded by the lines α and β through the points of O and A, respectively, could be obtained.Figure 4
Cauchy boundary problem considering the impact of the ground water table.
(a)
(b)Such asx, y, σe, and θ on the point 22 can be obtained by the points of 21 and 32, along line α:
(10.1)σe22−σe21−σe22+σe21+2σctanφθ22−θ21=−γ∗sinθ22+θ21/2+μy22−y21cosφcosθ22+θ21/2−μy22−y21=x22−x21cotθ22+θ212−μ,along line β:
(10.2)σe22−σe32−σe22+σe32+2σctanφθ22−θ32=−γ∗sinθ22+θ32/2+μy22−y32cosφcosθ22+θ32/2+μy22−y32=x22−x32cotθ22+θ322+μ,when soils are above ground water, γ∗=γ. When soils are below ground water, γ∗=γ′.If the ground water table is between two points, it should be simplified to calculate (Figure4(b)). Setting an allowable value Δh, if Δh1≤Δh, the γ∗ of the points 21 and 22 all take γ′. If Δh2≤Δh, the γ∗ of the points 21 and 22 all take γ. If the Δh1>Δh and Δh2>Δh, γ∗ of the points 21 and 22 take γ and γ′, respectively. In fact, the ground water table is not a straight line. The smaller the grid in the slip line field is, the more similar to a straight line the water line is.The above two sets of nonlinear equations ((10.1) and (10.2)) need to be solved by an iterative method. It can be assumed that the initial values of point 22 is
(11)x22=x21+x322,y22=y21+y322,orθ22=θ21+θ322,σe22=σe21+σe322.Substitutingx22 and y22 in (11) into (10.2), the approximations of θ22 and σe22 can be obtained. Similarly, substituting θ22 and σe22 in (11) into (10.1), the approximations of x22 and y22 can be obtained. Then, substituting the first approximations of θ22 and σe22 into (10.1), the second approximations of x22 and y22 can be gained, so repeatedly until the accuracy of the calculation to meet the precision. And x, y, σe, and θ of each point can be achieved in the entire OAB. Finally, according to the (9), σx, σy, and τxy in the slip line field could be gained.


## 3. Stress Discontinuity Line and Velocity Discontinuity Line

When soils reach the limit state, the stress discontinuity field and the velocity discontinuity field will be generated under complex boundaries [21]. With different boundary conditions, the stress discontinuity line and the velocity discontinuity line will also change. Therefore, the distribution and the characteristic of the two discontinuity lines are significant for the slope stability.


### 3.1. Characteristics of Stress Discontinuity Line

#### 3.1.1. Stress Equation

The stress discontinuity line is actually a thin transition zone, where the stress changes rapidly. It divides the element of the discontinuity line crossing through two plastic regions ① and ②, while only the tangential stress can be interrupted and the normal stress and the shear stress stay the same on both sides (Figure4). In the plastic regions ① and ②, stresses on the discontinuity line should obey the Mohr-Coulomb yield criterion.
(12)σt1≠σt2,σn1=σn2,τn1=τn2,(13)ω=θ++δ+=θ−+δ−.where θ+ and θ− are the maximum principal direction angles of the plastic regions ① and ②; w is the tangential direction angle of the discontinuity line; and δ+ and δ− are the angles between the maximum principal stress direction and tangential direction of stress discontinuity line and are given as follows (Figure 5):
(14)δ+=arctan−1±1+cos2φtan2θ+−θ−tanθ+−θ−1+sinφ,δ−=arctan1±1+cos2φtan2θ+−θ−tanθ+−θ−1+sinφ.Figure 5
Stress discontinuity line.The relationship of equivalent stresses on both sides of the discontinuity line is(15)σe+=ησe−,where σe+ and σe− are equivalent stresses on both sides of the discontinuity line; η is given as follows:
(16)η=1+cos2φtan2θ+−θ−+sinφ2sin2θ+−θ−cos2φ.


#### 3.1.2. Geometric Condition

As shown in Figure6, the angle ∠BAD of the principal stress element is 2γ. The angle between the stress discontinuity line ls and the principal stress plane in the zone ① is ζ, and the angle in the zone ② is υ. From the geometric relationship shown in Figure 6, it can be obtained:
(17)2γ=ζ+υ,(18.1)θ+=π2−ζ+ω,(18.2)θ−=υ+ω.Figure 6
Stress elements on both sides of stress discontinuity line.According to the geometric relationship shown in Figure2, (12) can be rewritten as follows:
(19)σe11−sinφcos2θ+−w=σe21−sinφcos2θ−−w,σe1sin2θ+−w=σe2sin2θ−−w.The relationship ofθ on both sides of a stress discontinuity line can be obtained eliminating σe1 and σe2 from (19).
(20)cos2ω−θ+−θ−+sinφcosθ−−θ+=0.Substituting (18.1) and (18.2) into (20):
(21)sinζ−υ−sinφsinζ+υ=0.The position of the stress discontinuity line can be determined according to (15) and (21) in the slip line field.


#### 3.1.3. Kinematic Equation

Because the stress discontinuity line is emerged by the rigid region reducing to the limit state, the tangential velocity on the discontinuity line remains unchanged [21]:
(22.1)dvxcosw+dvysinw=0,where vx and vy are velocity components in the x and y directions.Assuming thatvα and vβ are the velocity components in the direction of the characteristic line, a rigid kinematic equation of the stress discontinuity line will be obtained:
(22.2)vx=vαsinθ+μ−vβsinθ−μsecφ,vy=−vαsinθ+μ+vβsinθ−μsecφ.Substituting (22.2) into (22.1):
(23)sinθ+μ−ωdvα−sinθ−μ−ωdvβ=−vαcosθ+μ−ω+vβcosθ−μ−ωdθ.


### 3.2. Characteristics of Velocity Discontinuity Line

Obeying the associated flow rule, the velocity discontinuity line is a slip line or an envelope of slip line. The boundary between a rigid region and a plastic region is a velocity discontinuity line. So the velocity discontinuity line is the slip line at both ends of the stress discontinuity line, dividing the slope into rigid zone and plastic zone. For the Mohr-Coulomb material, the angle between the velocity direction and the slip line isφ [22, 23].


## 3.1. Characteristics of Stress Discontinuity Line

### 3.1.1. Stress Equation

The stress discontinuity line is actually a thin transition zone, where the stress changes rapidly. It divides the element of the discontinuity line crossing through two plastic regions ① and ②, while only the tangential stress can be interrupted and the normal stress and the shear stress stay the same on both sides (Figure4). In the plastic regions ① and ②, stresses on the discontinuity line should obey the Mohr-Coulomb yield criterion.
(12)σt1≠σt2,σn1=σn2,τn1=τn2,(13)ω=θ++δ+=θ−+δ−.where θ+ and θ− are the maximum principal direction angles of the plastic regions ① and ②; w is the tangential direction angle of the discontinuity line; and δ+ and δ− are the angles between the maximum principal stress direction and tangential direction of stress discontinuity line and are given as follows (Figure 5):
(14)δ+=arctan−1±1+cos2φtan2θ+−θ−tanθ+−θ−1+sinφ,δ−=arctan1±1+cos2φtan2θ+−θ−tanθ+−θ−1+sinφ.Figure 5
Stress discontinuity line.The relationship of equivalent stresses on both sides of the discontinuity line is(15)σe+=ησe−,where σe+ and σe− are equivalent stresses on both sides of the discontinuity line; η is given as follows:
(16)η=1+cos2φtan2θ+−θ−+sinφ2sin2θ+−θ−cos2φ.


### 3.1.2. Geometric Condition

As shown in Figure6, the angle ∠BAD of the principal stress element is 2γ. The angle between the stress discontinuity line ls and the principal stress plane in the zone ① is ζ, and the angle in the zone ② is υ. From the geometric relationship shown in Figure 6, it can be obtained:
(17)2γ=ζ+υ,(18.1)θ+=π2−ζ+ω,(18.2)θ−=υ+ω.Figure 6
Stress elements on both sides of stress discontinuity line.According to the geometric relationship shown in Figure2, (12) can be rewritten as follows:
(19)σe11−sinφcos2θ+−w=σe21−sinφcos2θ−−w,σe1sin2θ+−w=σe2sin2θ−−w.The relationship ofθ on both sides of a stress discontinuity line can be obtained eliminating σe1 and σe2 from (19).
(20)cos2ω−θ+−θ−+sinφcosθ−−θ+=0.Substituting (18.1) and (18.2) into (20):
(21)sinζ−υ−sinφsinζ+υ=0.The position of the stress discontinuity line can be determined according to (15) and (21) in the slip line field.


### 3.1.3. Kinematic Equation

Because the stress discontinuity line is emerged by the rigid region reducing to the limit state, the tangential velocity on the discontinuity line remains unchanged [21]:
(22.1)dvxcosw+dvysinw=0,where vx and vy are velocity components in the x and y directions.Assuming thatvα and vβ are the velocity components in the direction of the characteristic line, a rigid kinematic equation of the stress discontinuity line will be obtained:
(22.2)vx=vαsinθ+μ−vβsinθ−μsecφ,vy=−vαsinθ+μ+vβsinθ−μsecφ.Substituting (22.2) into (22.1):
(23)sinθ+μ−ωdvα−sinθ−μ−ωdvβ=−vαcosθ+μ−ω+vβcosθ−μ−ωdθ.


## 3.1.1. Stress Equation

The stress discontinuity line is actually a thin transition zone, where the stress changes rapidly. It divides the element of the discontinuity line crossing through two plastic regions ① and ②, while only the tangential stress can be interrupted and the normal stress and the shear stress stay the same on both sides (Figure4). In the plastic regions ① and ②, stresses on the discontinuity line should obey the Mohr-Coulomb yield criterion.
(12)σt1≠σt2,σn1=σn2,τn1=τn2,(13)ω=θ++δ+=θ−+δ−.where θ+ and θ− are the maximum principal direction angles of the plastic regions ① and ②; w is the tangential direction angle of the discontinuity line; and δ+ and δ− are the angles between the maximum principal stress direction and tangential direction of stress discontinuity line and are given as follows (Figure 5):
(14)δ+=arctan−1±1+cos2φtan2θ+−θ−tanθ+−θ−1+sinφ,δ−=arctan1±1+cos2φtan2θ+−θ−tanθ+−θ−1+sinφ.Figure 5
Stress discontinuity line.The relationship of equivalent stresses on both sides of the discontinuity line is(15)σe+=ησe−,where σe+ and σe− are equivalent stresses on both sides of the discontinuity line; η is given as follows:
(16)η=1+cos2φtan2θ+−θ−+sinφ2sin2θ+−θ−cos2φ.


## 3.1.2. Geometric Condition

As shown in Figure6, the angle ∠BAD of the principal stress element is 2γ. The angle between the stress discontinuity line ls and the principal stress plane in the zone ① is ζ, and the angle in the zone ② is υ. From the geometric relationship shown in Figure 6, it can be obtained:
(17)2γ=ζ+υ,(18.1)θ+=π2−ζ+ω,(18.2)θ−=υ+ω.Figure 6
Stress elements on both sides of stress discontinuity line.According to the geometric relationship shown in Figure2, (12) can be rewritten as follows:
(19)σe11−sinφcos2θ+−w=σe21−sinφcos2θ−−w,σe1sin2θ+−w=σe2sin2θ−−w.The relationship ofθ on both sides of a stress discontinuity line can be obtained eliminating σe1 and σe2 from (19).
(20)cos2ω−θ+−θ−+sinφcosθ−−θ+=0.Substituting (18.1) and (18.2) into (20):
(21)sinζ−υ−sinφsinζ+υ=0.The position of the stress discontinuity line can be determined according to (15) and (21) in the slip line field.


## 3.1.3. Kinematic Equation

Because the stress discontinuity line is emerged by the rigid region reducing to the limit state, the tangential velocity on the discontinuity line remains unchanged [21]:
(22.1)dvxcosw+dvysinw=0,where vx and vy are velocity components in the x and y directions.Assuming thatvα and vβ are the velocity components in the direction of the characteristic line, a rigid kinematic equation of the stress discontinuity line will be obtained:
(22.2)vx=vαsinθ+μ−vβsinθ−μsecφ,vy=−vαsinθ+μ+vβsinθ−μsecφ.Substituting (22.2) into (22.1):
(23)sinθ+μ−ωdvα−sinθ−μ−ωdvβ=−vαcosθ+μ−ω+vβcosθ−μ−ωdθ.


## 3.2. Characteristics of Velocity Discontinuity Line

Obeying the associated flow rule, the velocity discontinuity line is a slip line or an envelope of slip line. The boundary between a rigid region and a plastic region is a velocity discontinuity line. So the velocity discontinuity line is the slip line at both ends of the stress discontinuity line, dividing the slope into rigid zone and plastic zone. For the Mohr-Coulomb material, the angle between the velocity direction and the slip line isφ [22, 23].


## 4. Algorithm

For slope stability problems, the equilibrium equation and the plastic flow equation should be settled simultaneously to obtain the rigorous solution, which is difficult to achieve in mathematics and only depends on numerical methods. If a statically admissible stress fieldσij0 has been got, the strain rate field ε̇ij∗ and the kinematically admissible velocity field vi will be obtained by the stress field based on the associated flow rule. If the strain rate field and the velocity field are nothing less than the kinematical field, in this case, the plastic region corresponding to such a statically admissible stress field and a kinematically admissible velocity field must be the sliding belt in the extreme state and the external load must be the ultimate bearing capacity. Based on the upper and lower bound theorems, the numerical algorithm is established to get the distribution of the sliding belt and the ultimate bearing capacity considering the hydrostatic pressure.When a slope is on the verge of collapse (shear failure), not all the points in sliding mass reach to the yield state, but the mixture of the plastic region and the rigid region is emerged. As shown in Figure7, the entire slope is divided into five areas. The stress discontinuity line divides the slope foot into the strong plastic region and the weak plastic region, passive, and active, respectively. The boundary between two regions is the velocity discontinuity line, which is the slip line at both ends of the stress discontinuity line. So the stress discontinuity line and the velocity discontinuity line determine the distribution of the plastic region and the rigid region. The stress and deformation in each region are different as well as kinematic features. Only to meet all the possible kinematically admissible conditions, the solution will be exact.Figure 7
Calculation model.The calculation process (shown in Figure8) is consisted of two parts: a statically admissible field and a kinematically admissible field. The specific calculation process is as follows.Figure 8
Flow chart of rigorous solution for slope stability.


### 4.1. Static Calculation

(1)
The maximum principal stress on the noncharacteristic line stress boundaryOA is 0, and its direction angle is θ1, which is perpendicular to the slope surface. From the boundary OA, the slip line of the active region I is calculated in a very dense grid by solving the Cauchy problem.
(2)
To get the stress discontinuity line In the active region, the specific approach is to firstly select the pointF in the grid. The stress discontinuity line is drawn by (15) and (21). A group of the maixmum principal stress direction angles of the stress discontinuity line in the active region is assumed as θ+. Then, the direction angle θ−, the tangential direction angle of the discontinuity line w, and the coordinate of each point in the line AF in the passive region can be obtained by the (13) and (15). And the slip line field of the region II can be obtained by solving the Cauchy problem.
(3)
According to the direction angle of the maximum principal stressθ artificially assumed and the coordinate of each point in the line FE, the effective stress σe is obtained. From the lines FB and FE, the slip line field of region III can be got by solving the Riemann problem. On the basis of maximum principal stress direction angle θ2=90°, the slip line field of region V can be got by solving the mixed boundary value problem.
(4)
The numerical integration is carried out along the slip linesHF and FE to get the forces in the x and y directions and the moments of each force to the point O. Then the vertical stress on the boundary ED is calculated.
(5)
At this point, the static calculation is over. Then, check whether the force and the moment of the rigid block satisfy the equilibrium condition and whether the points within the rigid block satisfy the yield condition. If not satisfied, the positon of pointF is modified to recalculate. The whole calculation process is shown in Figure 9.Figure 9
Flow chart of static calculation.


### 4.2. Velocity Field Analysis

According to the associated flow rule, the dynamic analysis is carried out on the basis of the slip line field got by the static calculation.(1)
The absolute velocity of pointD on slope top is assumed as 1 m/s. The velocity component vα and vβ in the slip line ABCD is obtained according to the velocity characteristic equation.
(24)vα=vxcosθ−μ+vysinθ−μ,vβ=vxcosθ+μ+vysinθ+μ.
(2)
The velocity field of region II can be obtained by solving the mixed boundary value problem for the condition that the tangential velocity on the discontinuity line remains unchanged (23). From the slip lines BCD and BF, the velocity field in regions III and V and the velocity of each point in the line FE can be got by solving the Riemann problem. Then, the velocity vα− and vβ− of each point in the line AF in the passive region II is translated into the velocity vα+ and vβ+ in the active region I. Finally, the velocity field in the region I and the velocity in line FH are obtained by solving the Cauchy problem from the line AF.
(3)
At this point, the velocity analysis is over. Then, check whether the velocity in the linesHF and FE satisfy the kinematically rigid condition of the sliding block IV. If not satisfied, the input values of θ on the slip line FE and the position of the point F are modified to recalculate all.


## 4.1. Static Calculation

(1)
The maximum principal stress on the noncharacteristic line stress boundaryOA is 0, and its direction angle is θ1, which is perpendicular to the slope surface. From the boundary OA, the slip line of the active region I is calculated in a very dense grid by solving the Cauchy problem.
(2)
To get the stress discontinuity line In the active region, the specific approach is to firstly select the pointF in the grid. The stress discontinuity line is drawn by (15) and (21). A group of the maixmum principal stress direction angles of the stress discontinuity line in the active region is assumed as θ+. Then, the direction angle θ−, the tangential direction angle of the discontinuity line w, and the coordinate of each point in the line AF in the passive region can be obtained by the (13) and (15). And the slip line field of the region II can be obtained by solving the Cauchy problem.
(3)
According to the direction angle of the maximum principal stressθ artificially assumed and the coordinate of each point in the line FE, the effective stress σe is obtained. From the lines FB and FE, the slip line field of region III can be got by solving the Riemann problem. On the basis of maximum principal stress direction angle θ2=90°, the slip line field of region V can be got by solving the mixed boundary value problem.
(4)
The numerical integration is carried out along the slip linesHF and FE to get the forces in the x and y directions and the moments of each force to the point O. Then the vertical stress on the boundary ED is calculated.
(5)
At this point, the static calculation is over. Then, check whether the force and the moment of the rigid block satisfy the equilibrium condition and whether the points within the rigid block satisfy the yield condition. If not satisfied, the positon of pointF is modified to recalculate. The whole calculation process is shown in Figure 9.Figure 9
Flow chart of static calculation.


## 4.2. Velocity Field Analysis

According to the associated flow rule, the dynamic analysis is carried out on the basis of the slip line field got by the static calculation.(1)
The absolute velocity of pointD on slope top is assumed as 1 m/s. The velocity component vα and vβ in the slip line ABCD is obtained according to the velocity characteristic equation.
(24)vα=vxcosθ−μ+vysinθ−μ,vβ=vxcosθ+μ+vysinθ+μ.
(2)
The velocity field of region II can be obtained by solving the mixed boundary value problem for the condition that the tangential velocity on the discontinuity line remains unchanged (23). From the slip lines BCD and BF, the velocity field in regions III and V and the velocity of each point in the line FE can be got by solving the Riemann problem. Then, the velocity vα− and vβ− of each point in the line AF in the passive region II is translated into the velocity vα+ and vβ+ in the active region I. Finally, the velocity field in the region I and the velocity in line FH are obtained by solving the Cauchy problem from the line AF.
(3)
At this point, the velocity analysis is over. Then, check whether the velocity in the linesHF and FE satisfy the kinematically rigid condition of the sliding block IV. If not satisfied, the input values of θ on the slip line FE and the position of the point F are modified to recalculate all.


## 5. Example Verification

Geometric parameters and material properties are shown in Figure10. The Fortran 95 compiler is utilized to compute the calculation program for the rigorous solution satisfying the static equilibrium and the kinematical requirement.Figure 10
Parameters of a calculation model.When a slope reaches the limit state, the ultimate bearing capacity is composed of two parts (Figure11). While the load on the rigid region reaches 77.54 kPa, the other on the plastic region changes from 102.05 kPa to 112.45 kPa, considering the hydrostatic pressure. The rigid body rotates around point E, slipping along the sliding belt. The velocity of each point in the velocity discontinuity line ABCD increases gradually from top to foot (Figure 12). The angle between a velocity discontinuity line and a slip line is φ. The position of points E and D calculated determines the general distribution of the sliding belt. In order to facilitate the appearance, the figures (Figures 11 and 12) show only part of the slip lines and the actual existence is hundreds of groups.Figure 11
Distribution of the sliding belt and ultimate loads under the limit state.Figure 12
Velocity field of the slope under the limit state.The limit state can be seen as a status with safety factorFs equal to 1.0, and the external load at this time is the ultimate bearing capacity. When considering the hydrostatic pressure, the soils below ground water take floating weight resulting in the slide force falling and the slide resistance remains unchanged. From the Table 1, the ultimate load has become lager compared with the load without considering the ground water. As the depth increases causing the slide force to rise, the ultimate loads in both regions decrease, of which the one on the rigid region is influenced greater than the other by the ground water. While changing the depth, the distribution of sliding belt hardly alters. As the water level rises, the width of the sliding belt increases just a little.Table 1
Variation of the ultimate loads and the distribution of the sliding belt with the depthh.
Ground water depthh (m)
Ultimate load (kPa)
OE (m)
ED (m)
Rigid region
Plastic region (average)
0
77.80
107.69
2.32
0.87
0.5
77.54
107.25
2.32
0.87
1.0
76.69
106.74
2.31
0.87
1.5
73.34
105.39
2.31
0.86
2.0
70.78
104.38
2.30
0.86
2.5
64.09
101.67
2.29
0.86
3.0
59.83
99.93
2.28
0.85
4.0
43.89
93.51
2.27
0.85
No ground water
33.85
86.98
2.25
0.85From Table2, when weight of soils increases gradually, the ultimate load on the rigid region decreases corresponding to the other on plastic region increasing. The distribution of sliding belt undergoes great changes due to the reducing width of ED. For γ=0 kN/m3, the rigid region in the slope disappears under the limit state and the whole slip field becomes plastic. Changes of cohesive only bring about the increase of ultimate loads on both regions, and the sliding belt has a little change. Internal friction angle has a great impact on the load on the plastic region and the sliding belt. With friction angle rising, the load on plastic region increases sharply and the width of sling belt fills out gradually. When the angle φ reduces to 0 resulting the plastic region disappeared, the sliding belt transforms into a signal slip surface and the bearing capacity of a slope depends on the cohesive. At this time, the ultimate load calculated by the algorithm is considered as the external load on the slope top and this status can be calcualted by the Spencer method to get the slip surface and the safety factor. As shown in Figure 13, the critical slip surface obtained by this paper’s method is basically consistent with the Spencer method. The safety factor Fs calculated by the Spencer method is 0.99, when the ultimate load is 65.53 kPa. Generally speaking, the comparison proves the feasibility of the algorithm proposed by this paper.Table 2
Variation of the ultimate loads and the distribution of the sliding belt with different soil parameters.
c (kPa), φ (°), γ (kN/m3), h (m)
Ultimate load (kPa)
OE (m)
ED (m)
Rigid region
Plastic region
c=15φ=20h=0.5
γ=10
83.23
84.86
2.25
1.14
γ=15
79.99
93.97
2.29
0.99
γ=20
77.54
107.25
2.32
0.87
γ=25
71.91
125.62
2.34
0.77
γ=20φ=20h=0.5
c=10
53.52
78.61
2.33
0.85
c=15
77.54
107.25
2.32
0.87
c=20
110.57
191.98
2.31
0.89
c=25
198.88
331.25
2.30
0.92
γ=20c=15h=0.5
φ=15
54.22
75.89
2.39
0.79
φ=20
77.54
107.25
2.32
0.87
φ=25
108.52
198.55
2.28
1.03
φ=30
163.35
387.97
2.26
1.12Figure 13
Comparison between this paper’s method and the Spencer method.With the hydrostatic pressure, the difference between the two loads on both regions becomes smaller compared with the state of no ground water. The vairable amplitude of the ultimate loads and the sliding belt with the different parameters is also diminished.


## 6. Discussion

During the process of searching the rigorous solution of slopes stability, it is necessary to find a slip line field satisfying all possible static and kinematic conditions. Although the calculation is so difficult that the iterative computation takes thousands of times, the author used the Fortran 95 compiler to compile the calculation program, which only takes less than one minute to obtain the rigorous solution and greatly simplified the calculation. Since the grid does not meet the requirements of infinite subdivision, it is not each point on the boundary but the intersection of the grid and the boundary that satisfies the stress boundary condition and the velocity boundary condition.With the calculation model, the program is compiled based on the algorithm developed to obtain the ultimate bearing capacity and the distribution of the sliding belt. The influence of the ground water table on the ultimate bearing capacity and the sliding belt of a slope is discussed by an example. The results show that the ultimate load increases considering the hydrostatic pressure, while decreasing with the descending water level. However, the sliding belt remains constant. The soil weight and the internal friction angle have a tremendous impact on the belt, especially for the later. The ultimate load is influenced by the three parameters, of which the friction angle mainly controlled the load on the plastic region. When the internal friction angle reduces to zero, the sliding belt will translate into a traditional slip surface. Considering the ultimate load in this status as the external load on the slope top, the critical slip surface obtained by the Spencer method is basically consistent with the one got by this paper and the safety factor obtained by the Spencer method is 0.99 very close to 1.0. So the feasibility of this algorithm is verified by the specific example. Due to the hydrostatic pressure, the difference of loads on the two regions and the variation of sliding belt become smaller compared with not considering the ground water.On the basis of this study, the next research focus is to add the seepage field into the slip field, which are not completely coincide. The authors hope that others can be motivated to consider adding the penetration force into the equilibrium equation of slip line to gain the numerical solution and the influence of seepage on the sliding belt and the ultimate load.


## 7. Conclusion

Based on the upper and lower bound theorems, the rigorous solution satisfying all static boundaries and kinematic boundaries is obtained, which approaches lower and upper bound solutions at the same time and is considered as the actual solution. An algorithm for calculating the ultimate bearing capacity and the distribution of sliding belt has been developed, and the computer program has been accomplished. The application of this algorithm is verified by a specific example (Figure13).Considering the hydrostatic pressure, the variation of the ultimate load and the sliding belt under different parameters has been gained. The ultimate load decreases with the rising water level and become the minimum with no ground water. The sliding belt is controlled by the internal friction angle and the soil weight, while cohesive has a obvious impact on the ultimate loads. The difference of ultimate loads on the two regions and the variation of sliding belt become smaller because of the soils taking the floating weight.



---

*Source: 2829873-2018-06-12.xml*