Rigorous Solution of Slopes ’ Stability considering Hydrostatic Pressure

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.


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][6][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][9][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][15][16][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, for c − φ 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.

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: 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: Because M is a point on the Mohr's circle, R can be rewritten as follows: According to (1.4) and (2): Solving the unary quadratic equation of p: When σ t > σ n , according to (1.3): 3) should be rewritten as follows:

5
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:  2 Complexity With θ on the noncharacteristic line stress boundary, (6) can be used to draw the stress line field near the stress boundary: 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: 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 γ ′ .
As shown in Figure 4(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.Such as x, y, σ e , and θ on the point 22 can be obtained by the points of 21 and 32, along line α:  4(b)).Setting an allowable value Δh, if Δh 1 ≤ Δh, the γ * of the points 21 and 22 all take γ ′ .If Δh 2 ≤ Δh, the γ * of the points 21 and 22 all take γ.If the Δh 1 > Δh and Δh 2 > Δ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) 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 x 22 and y 22 can be obtained.Then, substituting the first approximations of θ 22 and σ e22 into (10.1), the second approximations of x 22 and y 22 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.

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.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 (Figure 4).In the plastic regions ① and ②, stresses on the discontinuity line should obey the Mohr-Coulomb yield criterion.
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): The relationship of equivalent stresses on both sides of the discontinuity line is where σ + e and σ − e are equivalent stresses on both sides of the discontinuity line; η is given as follows: 3.1.2.Geometric Condition.As shown in Figure 6, the angle ∠BAD of the principal stress element is 2γ.The angle between the stress discontinuity line l s and the principal stress plane in the zone ① is ζ, and the angle in the zone ② Ground water table  4 Complexity is υ.From the geometric relationship shown in Figure 6, it can be obtained: According to the geometric relationship shown in Figure 2, (12) can be rewritten as follows: The relationship of θ on both sides of a stress discontinuity line can be obtained eliminating σ e1 and σ e2 from (19).
Substituting (18.1) and ( 18.2) into (20): The position of the stress discontinuity line can be determined according to (15) and ( 21) in the slip line field.

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]: where v x and v y are velocity components in the x and y directions.
Assuming that v α 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: Substituting (22.2) into (22.1): 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].

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 σ 0 ij has been got, the strain rate field ε * ij and the kinematically admissible velocity field v i 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 Figure 7, 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.
The calculation process (shown in Figure 8) is consisted of two parts: a statically admissible field and a kinematically admissible field.The specific calculation process is as follows.

Static Calculation
(1) The maximum principal stress on the noncharacteristic line stress boundary OA 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 point F 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 lines HF 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

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 point D 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.
(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 BC D 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 lines HF 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.

Example Verification
Geometric parameters and material properties are shown in Figure 10.The Fortran 95 compiler is utilized to compute the calculation program for the rigorous solution satisfying the static equilibrium and the kinematical requirement.7 Complexity When a slope reaches the limit state, the ultimate bearing capacity is composed of two parts (Figure 11).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.
The limit state can be seen as a status with safety factor Fs 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.
From Table 2, 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/m 3 , 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 Velocit unit (m/s) 8 Complexity 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.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.

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.

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  9 Complexity 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 (Figure 13).
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.

Figure 1 :
Figure 1: Basic coordinate system and noncharacteristic line stress boundary.

Figure 4 :
Figure 4: Cauchy boundary problem considering the impact of the ground water table.

Figure 6 :
Figure 6: Stress elements on both sides of stress discontinuity line.
plastic region III: transition plastic region V: top plastic region IV: rigid region AF: stress discontinuity line HF: velocity discontinuity line ABCD: velocity discontinuity line

Figure 8 :Figure 9 :
Figure 8: Flow chart of rigorous solution for slope stability.

Figure 10 :
Figure 10: Parameters of a calculation model.

Figure 11 :
Figure 11: Distribution of the sliding belt and ultimate loads under the limit state.

Figure 13 :
Figure 13: Comparison between this paper's method and the Spencer method.
Effective stress of Mohr's circle.
Substituting x 22 and y 22 in (

Table 1 :
Figure12: Velocity field of the slope under the limit state.Variation of the ultimate loads and the distribution of the sliding belt with the depth h.

Table 2 :
Variation of the ultimate loads and the distribution of the sliding belt with different soil parameters.