9512.net

# 1997-Ian B Donald-Slope stability analysis by the upper bound approach, fundamentals and methods

853

Slope stability analysis by the upper bound approach: fundamentals and methods
Ian B. Donald and Zuyu Chen

Abstract: A new method for stability analysis in soils and rocks is presented, based on the upper bound theorem of classical plasticity. The sliding mass is divided into a small number of discrete blocks, with linear interfaces between blocks and either linear or curved bases to individual blocks. By equating the work done by external loads and body forces to the energy dissipated in shearing, either a safety factor or a disturbance factor may be calculated. The rigorous theoretical background is established, from which it may be demonstrated that for several well-defined classical slope problems the equations for the multi-block solution reduce to the published closed-form solutions. Powerful optimization routines are provided in the computer program EMU to search for the critical failure mechanism giving the lowest factor of safety. Several examples are given to demonstrate that, for problems where the exact answers are known, the new method produces accurate values of safety factor and predictions of failure mechanism. Applications to practical problems have shown that the new method is as simple as the conventional limit equilibrium methods for practitioners. Key words: slope stability, upper bound theorem, energy method, factor of safety, methods of optimization. Résumé : L’on présente une nouvelle méthode d’analyse de stabilité dans les sols et les roches basée sur le théorème de la limite supérieure de la plasticité classique. La masse glissante est divisée en un petit nombre de blocs discrets avec des interfaces linéaires entre les blocs et des bases linéaires ou courbes pour les blocs individuels. En mettant en équation le travail fait par les charges extérieures et par les forces du poids propre avec l’énergie dissipée durant le cisaillement, un coefficient de sécurité ou un facteur de perturbation peuvent être calculés. L’on établit la base théorique rigoureuse est à partir de laquelle il peut être démontré que, pour plusieurs problèmes de pente classiques bien définis, les équations pour la solution de blocs multiples se réduit aux solutions exactes publiées. L’on fournit des routines efficaces d’optimisation dans le programme d’ordinateur, EMU, pour la recherche du mécanisme de rupture critique donnant le coefficient de sécurité le plus faible. L’on donne plusieurs exemples qui démontrent que, pour les problèmes pour lesquels les réponses exactes sont connues, la nouvelle méthode fournit des valeurs du coefficient de sécurité et des prédictions du mécanisme de rupture précises. Des applications à des problèmes pratiques ont montré que la nouvelle méthode est aussi simple pour les praticiens que les métohodes conventionnelles d’équilibre limite. Mots clés : stabilité des pentes, théorème de limite supérieure, méthode d’énergie, coefficient de sécurité, méthodes d’optimisation. [Traduit par la rédaction]

Introduction
Over several decades, the limit equilibrium method has almost dominated the profession for examining the stability of slopes, embankments, and other soil and rock structures. The method, originating from a basically empirical background (Fellenius 1936), has been greatly improved by Janbu (1973), Morgenstern and Price (1965), Spencer (1967), Sarma (1973), Fredlund and Krahn (1977), Chen and Morgenstern (1983), and other researchers to satisfy the complete requirements for force and moment equilibrium and to accommodate generalized slip surfaces. Recent updating of the method includes the automatic searching for the minimum factor of safety and its associated critical slip surface using the methods of optimization (Celestino and Duncan 1981; Nguyen 1985; Chen and Shao

Received August 20, 1996. Accepted July 5, 1997. I.B. Donald. Department of Civil Engineering, Monash University, Clayton VIC 3168, Australia. Z. Chen. China Institute of Water Resources and Hydropower Research, P.O. Box 366, Beijing 100044, China.
Can. Geotech. J. 34: 853–862 (1997).

1988; etc.). However, there is no theoretical reason which adequately explains the success of its extensive applications. On the other hand, the potential of extending the plasticity method to soil and rock stability analysis has been investigated by many researchers. Early efforts were made by Frontard (1922) and Jaky (1936), but neither of these methods was based on realistic failure surfaces or numerically tractable methods. Sokolovski (1954), Fang and Hirst (1970), Booker and Davis (1972), and Graham (1973) employed plasticity theory, mainly based on the slip-line fields method, and limited their studies to problems with simple geometries. Later work, mainly concerned with the upper bound analysis, was conducted by Finn (1967), Chen and Snitbahn (1975), Karal (1977a, 1977b), Chen and Chan (1984), and Izbicki (1981). Extensions of the upper bound solutions to nonlinear failure envelopes have been investigated by Baker and Frydman (1983), Zhang and Chen (1987), Drescher and Christopoulos (1988), and Collins et al. (1988). Gussmann (1982) presented the kinematical element method, which divides a soil mass into finite elements and calculates the velocity distribution and energy dissipations of these elements. Sloan (1988, 1989) used finite elements and linear programming to perform both lower

854 Fig. 1. Failure mechanism of the upper bound approach. The real failure mechanism is the shaded area.

Can. Geotech. J. Vol. 34, 1997

work has been further updated by the authors of this paper with a multiblock failure mechanism which provides a fully analytical formulation of upper bound solutions capable of producing a series of closed-form solutions originally offered by Sokolovski (1954). Extensions to bearing capacity and threedimensional slope stability analysis problems (Chen 1995) have shown the great potential of the new method in providing innovative tools to other areas where simple and complete numerical methods are still not available. Although the research findings have been presented in several papers (Giam and Donald 1991; Donald and Chen 1995; Chen 1995), this paper for the first time completely covers the theoretical background, numerical techniques, validations, and extensions of this new upper bound slope stability analysis method.

bound and upper bound analysis. His work represents an attempt to obtain “bound solutions” by numerical methods on a theoretically rigorous foundation. Similar approaches were presented by Chuang (1992a, 1992b), who employed a pair of primal–dual linear programs that encoded kinematic and static requirements, respectively, in a finite element discrete version. Rotational displacements were included in his work. Despite the great volume of research work, reports of its practical applications appear to be rare. It is not difficult to understand that the limit analysis method could not be used extensively in the profession unless the following problems have been properly solved. (1) Numerical tractability—Baker and Frydman (1983) and Zhang and Chen (1987), among others, employed calculus of variations to find the least upper bound and consequently limited themselves to only a few demonstrative examples that had simple slope geometries and material properties. As the majority of practical problems are non-analytic, some numerical techniques must be employed to accommodate the varied geometry and heterogeneity commonly encountered in geotechnical problems. It has been found, in the area of limit equilibrium methods, that the techniques of optimization have been successfully used to minimize the factor of safety. Their extensions to limit analysis methods certainly deserve investigation. (2) Rational failure mechanism—Partly because of the numerical difficulties mentioned previously, most researchers employed logarithm spirals as failure surfaces and assumed the sliding mass to be a rigid body within which the internal energy dissipations are totally ignored (Baker and Frydman 1983; Zhang and Chen 1987). Karal (1977a, 1977b) has correctly pointed out that the sliding mass is generally formed by a “composed mode” which includes both rigid and “soft” parts. Within the latter, energy dissipation cannot be ignored. It is suspected that these simplifications would be too great to produce sufficiently accurate solutions. (3) Demonstrations of feasibility and validity—The limit analysis method will not prove to be competitive with the limit equilibrium method until sufficient evidence has been gained to demonstrate its validity. Comparing the results with wellknown closed-form solutions is a means commonly employed (Chen 1975; Finn 1967). During the past decade, the first author and P. Giam worked intensively on a theoretically rigorous and numerically efficient upper bound limit analysis method (Giam and Donald 1991). A multiwedge failure mechanism was developed to include the energy dissipations both on the slip surface and within the sliding body. The early

Review of the numerical approach to the upper bound method
The upper bound theorem The extensions of the upper bound theorem in plasticity to solving geotechnical problems have been explored by Chen (1975). When applying the upper bound theory to slope stability problems, it is assumed that during failure, a slip surface Γ divides the slope into a plastic failure zone, in which the stress state at any point is either on or inside the yield surface, and an elastic zone, in which the displacement at any point is virtually negligible. Shear failure exclusively dominates along the slip surface and within the plastic zone (Fig. 1). The statement of the upper bound theorem, particularly concerned with slope stability analysis, can then be described as follows. Assign a kinematically admissible strain increment . in the plastic zone ??, and the velocity V* along the field ε? ij ? slip surface Γ , then the external surface load T* calculated by the equation .? dDs? = WV? + T?V? [1] ij εij d? + ∫ ∫ σ?
?? Γ?

will be either greater than or equal to the real surface load T, which is associated with a real failure mechanism represented by the plastic zone ? and the slip surface. Γ, where σ? ij is the stress in the plastic zone that produces ε? ij according to the normality flow rule, Ds is the energy dissipation on the slip surface, W is the weight of the sliding mass, and V is the velocity at the points at which W and T apply. The slope failure mechanism considered herewith generally offers a continuous stress as well as a velocity field except along the slip surface where the plastic velocity V changes abruptly to zero at the elastic zone through an infinitesimally thin shear band (Fig. 2). Based on the associated flow law and Mohr–Coulomb failure criterion, it can be shown (Chen 1975; Giam and Donald 1991) that the velocity at the shear band V inclines at an angle φ to the band and the energy dissipation along this slip surface per unit area is (Fig. 2) [2] dD = τVt + σnVn = (cos φ ? τ + sin φ ? σn)V = (c cos φ ? u sin φ)V where c and φ are shear strength parameters; Vt and Vn are components of V in tangential and normal directions, respectively;

Donald and Chen Fig. 2. Energy dissipation on a shear band. Fig. 3. A three-wedge failure mechanism.

855

u is the pore pressure applied on the shear band; σn and τ are the normal and shear stress on the shear band, respectively. Equation [2] indicates that the energy dissipation based on the Mohr–Coulomb criterion can be conveniently calculated without the knowledge of the stresses on the shear zone. This has contributed significantly to the success of the approach described herein. The procedure of solving structural problems includes assuming a series of compatible displacement patterns and calculating their respective load T* based on [1]. The one associated with the minimum T* will most likely be the true load that brings the structure into failure. The concept is referred to as the least upper bound approach. The approximated failure mechanisms: multiwedge systems For a slope concerned, the plastic zone is divided into a number of wedges which are created by inclined straight-line interfaces and bases. Figure 3 shows a three-wedge system. Each wedge moves with a velocity V that inclines at an angle φ to the slip surface, creating a relative velocity Vj to its neighbouring wedges along the interfaces. The wedge itself moves as a whole. No energy dissipation will develop within the wedge body. For a system containing n wedges and n – 1 interfaces, [1] is then replaced by its approximation
n?1

[4]

η = ηt =

T ? To To

This alternative, defined as alternative 1, can be conveniently used in bearing capacity problems. (2) The failure mechanism is brought about by a pseudohorizontal body force whose magnitude is determined by ηbW, where W is the self-weight of the plastic zone. The value of ηb that brings the failure is called the “coefficient of critical acceleration” (Sarma 1973). This alternative, defined as alternative 2, is preferred in slope stability analysis, since surface loads in many cases do not exist. Equation [3] can be replaced by
n?1

[5]

k=1

Djk

__ ? o ? ? + ∑ ?Ds k = WV + T V + ηbWV
k=1

n

__ where W is a force in the negative x axis direction with a magnitude of W. (3) The failure mechanism is created by decreasing the shear strength parameters by a coefficient called factor of safety, F, which provides a new cohesion ce and friction angle φe by the definitions [6] [7] ce = c F tan φ F

tan φe =

[3]

k=1

Djk

+∑
k=1

n

?Ds k

=

WV?

+

T?V?

Also based on [1], F can be obtained by solving the following equation using iterations:
n?1 n

where the first and second terms refer to the energy dissipations developed on the interfaces and slip surface, respectively, both calculated based on [2]. The disturbance that brings about the failure mechanism Most structural problems are concerned with the system that is not at failure and the question needing to be answered is how large an external disturbance need be to bring this system from a safe state to a state in which the failure mechanism appears. In slope stability analysis there are several alternatives for defining the external disturbance. (1) If there is a load To applying on the surface, then it can be increased to T at which failure happens. The disturbance factor η is then defined as

[8]

∑ Djek + ∑ ?Dsei = WV? + ToV?
k=1 i=1

The subscript e indicates that the corresponding terms are calculated based on ce and φe, which involve the unknown F. Since [8] is nonlinear, F will be obtained by trial and error or by the Newton–Raphson method. Sarma (1973) suggested an iteration method based on the critical acceleration concept. To facilitate presentation, the symbols ce and φe are invariably used for all three alternatives in the remaining part of the paper. This means that a factor of safety of unity is actually implied when alternatives 1 or 2 are concerned. Calculation of the compatible velocity field Let us observe the two adjoining wedges as shown in Fig. 4. The left and right wedges move with the absolute velocities Vl

856 Fig. 4. Velocity compatibility between adjacent blocks. The left wedge moves upward. (a) Velocities of the blocks. (b) Velocity hodograph.

Can. Geotech. J. Vol. 34, 1997 Fig. 5. Velocity compatibility between adjacent blocks. The left wedge moves downward. (a) Velocities of the blocks. (b) Velocity hodograph.

and Vr which incline at angles φel and φer to their bases, or θel and θer to the x axis, respectively. The relative velocity of the left wedge with respect to the right one along the interface is represented as Vj, which inclines at an angle φej or (π ? φej), depending on whether the left wedge moves upward or downward with respect to the right one, an important consideration that will be discussed subsequently. To allow the velocities assigned to the n-wedge failure mechanism to be kinematically compatible, the two adjoining wedges must not move to cause overlap or indentation. This implies that the velocity hodograph must be closed, i.e., [9] Vr + Vj = Vl sin(θl ? θj) sin(θr ? θj) sin(θr ? θl) sin(θr ? θj)

Solving [9], we obtain [10] Vr = Vl Vj = Vl

the first interface to the kth one that separates the kth and (k + 1)th wedges. It has been noted that the left wedge can move either upward or downward with respect to the right wedge. The former case (refer to Fig. 4), defined as case 1, is most usually encountered. However, if Vr lies over Vl, and consequently θr < θl as illustrated in Fig. 5, the left wedge will move downward instead of upward with respect to the right wedge. This situation, defined as case 2, occurs, for example, where the base of the left wedge is a weak band that offers a low friction angle and produces abrupt change in α along the slip surface. The relative velocity in this case would incline at an angle of (π ?φe) rather than φe to the interface. In general, we define [14] θl = π + αl ? φel [15] θr = π + αr ? φer and [16] θj = [17] θj = π ? δ + φej for case 1 2 3π ? δ ? φej for case 2 2

[11]

where θl, θr, and θj are measured from the positive x axis in a counterclockwise direction, i.e., 0 ≤ θ ≤ π, and δ is the inclination of the interface measured from the positive y axis to the positive x axis. Starting from the first interface, the velocity of any wedge at the right side of the kth interface, represented as V, can be expressed in terms of V1, the velocity of the first wedge by successive calculations based on the equation [12] V = κV1 where [13] κ=∏
i=1 k

It is worth noting that this argument applies to all methods that employ nonvertical slices, such as those proposed by Sarma (1979) or Kovari and Fritze (1984). If the direction of possible movement between two slices is not properly considered, these methods will also give incorrect answers on some occasions, particularly where the base failure surface contains sections of reverse curvature. Conditions of kinematic admissibility Determining the velocity field by [10] and [11] is subject to the conditions of kinematic admissibility which state as follows: [18] θl > θj > θr ? π for case 1 [19] θl < θj < θr + π for case 2

sin(αli ? φlei ? θji)
r j sin(αr i ? φei ? θi)

where α is the inclination of the base to the horizontal, and the superscripts l and r refer to the values at the left and right sides of the interfaces, respectively. The multiplication starts from

Donald and Chen Fig. 6. The multiblock failure mechanism.

857

[26] V = Vlk

sin(αlk ? φlek ? θj)
r sin(αr k ? φek ? θi)

dα ? ? x dξ × exp??∫ cot(α ? φe ? θj) dξ ? xk ? ? Since ? xk dα ? [27] Vlk = Vr cot(α ? φe ? θj) dξ k?1 exp ??∫ dξ ? ? xk?1 ? we have the following equation calculating the velocity field for the multiblock failure mechanism: Satisfaction of these conditions implies that the values of Vr and Vj calculated by [10] and [11] must be positive if Vr is positive. In summary, we present the following kinematic conditions for the multiwedge failure mechanism. The formal demonstration is given elsewhere (Donald and Chen 1992). Case 1: the left wedge moves upward, i.e., θr > θl, and it is required that θj take the definition of [16] and [20] θl ? θj > 0 Case 2: the left wedge moves downward, i.e., θr < θl, and it is required that θj take the definition of [17] and [21] and [22] θr ? θj > ? π The energy method formulated in terms of continuous media If a part of the slip surface has a smoothly curved shape with uniform shear strength parameters, i.e., φe is constant and α is continuous in an interval (xk, xk+1) (refer to Fig. 6), the velocity at any point can be calculated by integration rather than successive multiplications. Suppose the velocity at point A (refer to Fig. 6), with its abscissa value x, is V, which has an increment dV when A moves to B, which has the abscissa value (x + dx). Substituting V and V + dV for Vl and Vr, respectively, in [10], we have [23] V + dV = V sin(α ? φe ? θj) sin(α + dα ? φe ? θj) where dα ? ? x dξ [29] E(x) = κ exp ??∫ cot(α ? φe ? θj) dξ ? ? ? xo Equation [28] indicates that the velocity at any point of the slip surface, located right of the kth discontinuity of φe or α, can be calculated by direct integration within the interval (xo, x), based on Vo, the velocity at the left end of the slip surface. The effect of possible discontinuities in φe and α has been accounted for by the coefficient κ. Neglecting high-order terms, [11] can also be reduced to [30] Vj = ?V cosec(α ? φe ? θj)dα = ?Vo cosec(α ? φe ? θj)E(x)dα Equation [30] is applicable only when φe is constant and α is continuous. It indicates that within the interval (xk, xk+1), Vj is an infinitesimally small magnitude with the order of dα. If the slip surface is a straight line, i.e., dα/dx = 0, Vj will be zero. This means that the energy dissipation within a wedge body is zero if φe is uniform along the base of the wedge. With the approaches described herein, we are able to define a multiblock failure mechanism (refer to Fig. 6). The plastic zone of the slope is divided by a number of blocks that have curved bases with constant φe. Each block will be further divided into a number of slices, or wedges, whose interfaces incline at angles of δ from the vertical, which are determined by linear interpolation based on the corresponding values at the left and right interfaces of the block. The velocities along the slip surface and the interfaces can be readily determined

[28] V = Vr k?1

sin(αlk ? φlek ? θj)
r sin(αr k ? φek ? θj)

dα ? ? x dξ × exp??∫ cot(α ? φe ? θj) dξ ? xk?1 ? ? = Vo ∏
i=1 k

sin(αli ? φlei ? θji)
r j sin(αr i ? φei ? θi)

θl ? θj < 0

dα ? ? x dξ × exp??∫ cot(α ? φe ? θj) dξ ? xo ? ? x ? dα ? = κVo exp??∫ cot(α ? φe ? θj) dξ dξ ? ? xo ? = E(x)Vo

Neglecting magnitudes with the orders higher than dα, we have [24] ? dα dV = cot(α ? φe ? θj) dx V dx

Integration of [24] yields ? x dα ? dξ [25] V = Vr k exp ??∫ cot(α ? φe ? θj) dξ ? xk ? ? where Vr is the velocity at the right side of the kth interface, k and ξ is a dummy variable substituting for x. By virtue of [10], [25] can be written as

858

Can. Geotech. J. Vol. 34, 1997

by [28] and [30], except at the point of discontinuities where Vj is determined by [11]. Calculations of the work and energy dissipations With the velocity field obtained for a specified failure mechanism, the work done by the external loads and the internal energy dissipations of a block can be calculated by integration or summations. Considering a slice in a block with a base width dx, and taking Vo = 1, the various terms of work and energy dissipation for this slice can be determined as follows. Work done by the body forces The body forces applying on a slice include (i) self weight dW, (ii) pseudo-seismic force η′ dW, and (iii) pseudo-disturbance force ηb dW that brings about the failure mechanism. Work done by these body forces is [31] d(WV?) = dW[sin(α ? φe) + (η′ + ηb)cos(α ? φe)]E(x)

Internal energy dissipations The internal energy dissipations are contributed by two parts: (i) the energy dissipation along the base of a slice [33] dDs = (ce cos φe ? u sin φe)sec α E(x)dx and (ii) the energy dissipation along the interface between two adjoining slices at the point of discontinuity and within the blocks, which are calculated separately. If an interface is located at a point of the slip surface where α is continuous and φe is constant, the energy dissipation due to an increment in α between the two slices can be calculated based on [30] as [34] dDj = ?(cje cos φje ? uj sin φje)L cosec(α ? φe ? θj)E(x)dα where L is the length of the interface. If an interface is located at the point of the slip surface where α or φe changes abruptly, the energy dissipation at this particular kth interface is calculated based on [11] as [35] ?Djk = ?(cje cos φje ? uj sin φje)k Lk
l × cosec(αr ? φr e ? θj)k sin(?α ? ?φe)k E (xk)

Work done by the surface loads Surface loads should be transferred to the base of the slice where integration is made and are represented by dTx and dTy, respectively. For the case of disturbance alternative 1, another ηt dTx and ηt dTy are added. The work done by the surface loads and the pseudo-disturbance force is [32] d(TV?) = [(1 + ηt)dTy sin(α ? φe) + (1 + ηt)dTx cos(α ? φe)]E(x) The positive directions of the external forces that do work are defined to be opposite to the coordinate axis.

where ?α and ?φe are the increments in α and φ from the left side to the right side of the interface, respectively, and the superscript l represents the corresponding values at the left side of the interface. Formulations for the upper bound approach Substituting [31]–[35] into [3], [5], and [8], we obtain the formulas to calculate the disturbance factors for various alternatives with the following definitions:

[36] ?∫

xn ? dW dT y ? ? dW dTx ? ? ? G = ∫ ?(ce cos φe ? u sin φe)sec α ? ? + + sin(α ? φe)? ?η′ ? ? cos(α ? φe)?E(x)dx d x d x d x d x xo ? ? ? ? ? ? xn

xo

(cje

cos

φje

?

uj

sin

φje)

dα l L cosec(α ? φe ? θj) E(x)dx ? ∑(cje cos φje ? uj sin φje)k Lk cosec(αr ? φr e ? θj)k sin(?α ? ?φe)k E (xk) dx
k=1

n?1

[37]

Gb = ∫

xn xo

dW cos(α ? φe)E(x)dx dx

For alternative 3, the factor of safety is calculated by solving the equation [41] G = 0 in which F is involved in ce and φe. Iterations are necessary to find F.

xn dT dTx ? ? y sin(α ? φe) + cos(α ? φe)? E(x)dx [38] GT = ∫ ? dx xo ? dx ?

For alternative 1, the disturbance coefficient of the proportional pseudo-surface load is [39] G η = ηt = GT G Gb

Numerical techniques for finding the least upper bound
The method of optimization To determine the critical failure mode that gives the minimum factor of safety, various methods of optimization have been employed. A slip surface is divided by a number of nodal points A1, A2,..., Am (Fig. 7) with coordinates given by (i = 1, 2,..., m): ?x ? [42] Zi = ? i? y ? i?

For alternative 2, the coefficient for the critical acceleration is [40] η = ηb =

Note again that for alternatives 1 and 2 the symbols ce and φe involved in these equations represent the actual values c and φ, implying that F = 1.

Donald and Chen Fig. 7. Simulation of a generalized slip surface. Slip surfaces: 1, initial estimate; 2, critical. Fig. 8. Closed-form solution: weightless slope with a vertical surface load.

859

Each pair of contiguous nodal points is connected by a straight line or a smooth curve. The smooth curve is generated by spline functions and is generally preferred unless a weak band, such as the part A2A3 in Fig. 7, is simulated. At each nodal point, an interface is assigned with its inclination designated as δ. The value of factor of safety F, or the disturbance factor η, can then be determined by solving either [39], [40], or [41]. During the integrations, each block is further divided into a number of slices to accommodate possible changes in geometry and shear strength parameters. F and η are then expressed as functions with respect to x1, y1, x2, y2,..., xm, ym, δ1, δ2,..., δm for m nodal points: [43] F = F (Z, δ) = F (x1, y1, x2, y2,..., xm, ym, δ1, δ2,..., δm) [44] η = η (Z, δ) = η (x1, y1, x2, y2,..., xm, ym, δ1, δ2,..., δm) The task of evaluating the stability of a slope becomes a numerical problem of finding a set of variables Z and δ that gives the minimum F or η with the associated slip surface connected by the nodal points B1, B2,..., Bm. The technique of finding the minimum factor of safety is similar to that in conventional methods as discussed by a number of authors (Baker 1980; Celestino and Duncan 1981; Nguyen 1985; Li and White 1987; Sun 1984; Chen and Shao 1988). Chen (1992a, 1992b) has discussed the difficulties in finding the global minimum factor of safety. A method called random search has been proposed and proved to be helpful in approaching the critical failure mechanism. A computer program EMU (Energy Method Upper-bound) has been coded to perform various practical stability analysis problems, including the ship lock slope of the Three Gorges Project in China, with satisfactory results. Validations of the numerical method The validity and feasibility of the new method have been assessed by different approaches as described herewith. Comparison of the governing equation (eq. [39]) with the closed-form solution Figure 8 shows a slope composed of uniform material with cohesion c and friction angle φ. It is subjected to a uniform

vertical surface load q, but the unit weight of the material is zero. Based upon the slip-line fields method and the Sokolovski (1954) basic equations, it is possible to find the closedform solution to the limit load q: ? ? 1 + sin φ exp[(π ? 2χ) tan φ] ? 1? [45] q = c cot φ ? 1 ? sin φ ? ? where χ is the inclination of the slope from the horizontal. The critical slip surface consists of three lines. The straight lines AB and CD incline at an angle of ? to the slope surface and the vertical, respectively, where ? is defined as [46] ? = π φ ? 4 2

BC is a log-spiral line with its left and right interfaces BO and CO inclined at the angle ? to the slope surface and the vertical, respectively. It is possible to demonstrate that for this particular problem, with the known slip-lines, [45] is identical to [39]. For details, refer to Donald and Chen (1992). Comparisons of the numerical results with closed-form answers Figure 9 shows a test example which has closed-form solutions based on the slip-line field method (Sokolovski 1954). An inclined surface load is applied on a uniform, weightless slope with parameters c = 750 kPa, φ = 35°, χ = 35°, and δ′ = 24°. The closed-form solution gives q = 6228 kPa. In the numerical approach, the slip surface is simulated by splines connecting five nodal points. Using alternative 1, the value of η calculated by [39] for the initially estimated slip surface 1 was 0.07, and the value of Fo was found to be 1.034 if alternative 3 was subsequently employed. The results of 200 randomly generated slip surfaces suggested a slip surface 2 which gave η = 0.034. Starting with this slip surface, the optimization process gave the critical slip surface and interfaces shown as line 3 with ηm = 0.0033 and Fm = 1.007. The critical slip surface, interfaces, and minimum factor of safety obtained by the numerical upper bound approach are close to those anticipated by theory.

860 Fig. 9. Test problem for the closed-form solution shown in Fig. 8. Slip surfaces and interfaces: 1, initial estimate (F = 1.034); 2, random search result; 3, final solution (Fm = 1.007).

Can. Geotech. J. Vol. 34, 1997 Fig. 11. Force equilibrium of an inclined slice as shown in Fig. 6.

Fig. 10. Reevaluations of the Tianshenqiao landslide. Failure surfaces: I, actual; II, that obtained by Spencer’s method (Fm = 0.863); III, that obtained by EMU (Fm = 0.882).

Equivalence between the energy and force equilibrium methods
It has been argued that the new method employs the associated flow law, which might limit its applicability. This section demonstrates that the new method is equivalent to the conventional limit equilibrium method that employs nonvertical slices (Sarma 1979) and does not introduce the associated flow law. Figure 11 shows the forces applied on an inclined slice as shown in Fig. 6. Sarma (1979) assumes that the normal and shear forces on the base and the inclined interfaces all obey the Mohr–Coulomb criterion. In Fig. 11, P′ is the resultant of the normal effective force N′ on the slip surface and the frictional resistance generated by N′, which is tangent to the slip surface with a magnitude of N′ tan φe. P′ is therefore inclined at an angle of φe to the normal. A similar definition applies to G′, which refers to the interface. By projecting all the forces applied on the slice onto line AA′, which inclines at an angle of φe to the slip surface, the equilibrium equation for a slice can be formulated as [47] ?(?Ty + ?W)sin(α ? φe) ? (η′?W + ?Tx) cos(α ? φe) + (ce?x cos φe ? u?x sin φe)sec α ? ?[G′ cos(ψ ? + φje) ± cjeL sin ψ + ujL cos ψ] ? ?Q = 0 The symbols ± and ? + refer to the two possibilities of movement directions on the interface as discussed in the previous section. The upper symbol, for example + in ±, represents case 1, and the lower symbol, –, represents case 2, [48] ψ = α + δ ? φe
and ?Q refers to the disturbance factors and is determined by the following equations. For alternative 1,

Comparisons of the numerical results of EMU and other programs In 1989, a soil slope stability programs review project (Donald and Giam 1989) was organized by the Australian Association for Computer Aided Design (ACADS). A set of 10 standard test problems was issued to users and some authors of various slope stability programs, which included those developed by several invited specialists whose work has been well documented (Baker 1980; Fredlund 1977; Chen and Shao 1988). Donald and Giam (1989, 1992) gave a comprehensive review which confirmed that, for all 10 problems, EMU gave results equal to or very slightly higher than the “referee answers” which were based on conventional limit equilibrium methods and the answers given by the invited specialists. Figure 10 gives another example which evaluates the landslide of the Tianshenqiao Hydro-power Project documented by Chen and Shao (1988). The value of Fm obtained by EMU was 0.882 and that from the limit equilibrium method was 0.863. Although the critical slip surfaces differ in shape, they both suggest passing through the toe of the slope, a situation consistent with what had been observed in the field.

Donald and Chen

861

[49] [50]

?Q = ηt[?Ty sin(α ? φe) + ?Tx cos(α ? φe)] ?Q = ηb?W sin(α ? φe) dG′ dδ ? tan(ψ ? + φje) G′ = p(x)sec(ψ ? + φje) dx dx p(x) = p1(x) + p2(x) + p3(x)

and for alternative 2, Let ?x ? >0, and [47] becomes [51] where [52] ? dW dT x ? dW dTx ? [53] p1(x) = ?? ? dx + dx ?sin(α ? φe) ? ?η′ dx + dx ? ? ? ? ? × cos(α ? φe) + (ce cos φe ? u sin φe)sec α [54] p2(x) = ? + d j (c L sin ψ ± ujL cos ψ) dx e

For alternative 1, dTx ? ? dTy sin(α ? φe) + cos(α ? φe)? p3(x) = ?ηt? d x d x ? ? and for alternative 2, [55] [56] p3(x) = ?ηb dW sin(α ? φe) dx

least upper bound, which has enabled much wider applications to practical problems that might involve complicated slope contours and material heterogeneities. (3) The new method is supported by a sound theoretical background. A number of test examples show that it can give results as accurate as closed-form solutions. (4) The requirement for kinematic admissibility included in the new method allows more rational approaches for practical problems. The limitations of this method include a possible overestimate of factor of safety if the optimization routines fail to find the real or global minimum. The authors’ experience indicated that the random-search technique (Chen 1992b) was very useful for facilitating successful searches. Coupling the new approach with conventional methods such as the Morgenstern–Price method, which virtually implies a lower bound solution, will on most occasions give accurate solutions, since the gaps between the upper and lower bounds are generally very small, as shown in the examples in this paper. The potential for extending the method to other areas of geomechanics is highly promising.

Acknowledgements
The authors wish to acknowledge the early work in this field contributed by the late Dr. P. Giam. The financial support for this research came from the Australian Research Council and the China National Natural Science Foundation.

Solving the differential equation with relevant boundary conditions, solutions are obtained that are identical to [39] and [41]. The derivations, which are lengthy and complex, can be found in Chen and Donald (1993). The derivations demonstrate the following enlightening points: (1) Once a multiblock failure mode is established, the problem becomes statically determinate, and the answer is therefore unique, irrespective of the methods used to find it. (2) The energy approach employs the associated flow law as a working assumption to facilitate the solutions by [3]. The solutions have been mathematically demonstrated to be identical to those obtained by the force equilibrium approaches in which the associated law is not employed, but are obtained in a conceptually very simple way by establishing a work and energy balance equation. (3) Although both methods succeeded in finding the solution, the mathematics involved in the force equilibrium approach is complex and in fact would not have become available if the derivations were not directed by the known answers of [39] or [41]. This fact is even more significant if three-dimensional slope stability is concerned.

References
Baker, R. 1980. Determination of critical slip surface in slope stability computations. International Journal for Numerical and Analytical Methods in Geomechanics, 4: 333–359. Baker, R., and Frydman, S. 1983. Upper bound limit analysis of soil with non-linear failure criterion. Soils and Foundations, 23: 35–42. Booker, J.R., and Davis, E.H. 1972. A note on plasticity solutions for the stability of slopes in homogeneous clays. Géotechnique, 22: 509–513. Celestino, T.B., and Duncan, J.M. 1981. Simplified search for noncircular slip surface. In Proceedings, 10th International Conference on Soil Mechanics and Foundation Engineering, Stockholm, Vol. 3, pp. 391–394. Chen, W.F. 1975. Limit analysis and soil plasticity. Developments in Geotechnical Engineering No.7. Elsevier Scientific Publishing Co., New York. Chen, W.F., and Chan, S.W. 1984. Upper bound limit analysis of the stability of a seismic-infirmed earth slope. In Proceedings of an International Symposium on Geotechnical Aspects of Mass and Material Transportation, Bangkok, pp. 373–428. Chen, W.F., and Snitbahn, N.S. 1975. On slip surface and slope stability analysis. Soils and Foundations, 15: 41–49. Chen, Z.Y. 1992a. Experience with the search for minimum factor of safety of slopes. In Proceedings of the 6th Australian – New Zealand Conference on Geomechanics, Christchurch, New Zealand, pp. 426–431. Chen, Z.Y. 1992b. Random trials used in determining global minimum factors of safety of slopes. Canadian Geotechnical Journal, 29: 225–233. Chen, Z.Y. 1995a. Recent developments in the energy approach to slope stability analysis. In Modern developments in geomechanics, the Ian Boyd Donald Symposium. Edited by C.M. Haberfield. Monash University, Melbourne, Australia, June 7, 1995, pp. 51–63.

Conclusions
The new method proposed in this paper features the following: (1) It approaches the solution by establishing a compatible velocity field and using the upper bound theorem of plasticity. Compared to the conventional methods that employ force equilibrium conditions (Sarma 1979), this method is more easily formulated, since it includes only a scalar manipulation such as [1], [3], or [8] and presents more rational solutions by considering two possibilities of relative movement between slices. (2) It employs the method of optimization to approach the

862 Chen, Z.Y. 1995b. Recent developments in slope stability analysis. In Proceedings of the 8th International Congress on Rock Mechanics, September 25–30, Tokyo, Vol. 3, pp. 1041–1048. Chen, Z.Y., and Donald, I. 1993. On the equivalence between the energy and limit equilibrium method. Internal Report, Institute of Water Resources and Hydropower Research, Beijing, China. (In Chinese.) Chen, Z.Y., and Morgenstern, N.R. 1983. Extensions to the generalized method of slices for stability analysis. Canadian Geotechnical Journal, 20: 104–119. Chen Z.Y., and Shao, C. 1988. Evaluation of minimum factor of safety in slope stability analysis. Canadian Geotechnical Journal, 25: 735–748. Chuang, P.H. 1990a. Stability analysis in geomechanics by linear programming. I: formulation. ASCE Journal of Geotechnical Engineering, 118(GT11): 1696–1715. Chuang, P.H. 1990b. Stability analysis in geomechanics by linear programming. II: application. ASCE Journal of Geotechnical Engineering, 118(GT11): 1716–1726. Collins, I.F., Gunn, C.I.M., Pender, M.J., and Yan, W. 1988. Slope stability analysis for materials with a non-linear failure envelope. International Journal for Numerical and Analytical Methods in Geomechanics, 12: 533–550. Donald, I., and Chen, Z.Y. 1992. Slope stability analysis by an upper bound plasticity method. Internal Report, Monash University, Clayton, Australia. Donald, I., and Chen, Z.Y. 1995. Upper bound stability solutions in geomechanics. Computational plasticity, fundamentals and applications. In Proceedings of the 4th International Conference on Computational Plasticity, April, Barcelona, Pineridge Press, Vol. 2, pp. 1797–1808. Donald, I., and Giam, P. 1989. Soil slope stability programs review. ACADS – Association for Computer Aided Design, Australia, Report U255. Donald, I., and Giam, P. 1992. The ACADS slope stability programs review. In Proceedings of the 6th International Symposium on Landslides, Christchurch, New Zealand, February 10–14, Vol. 3, pp. 1665–1670. Drescher, A,. and Christopoulos, C. 1988. Limit analysis slope stability with non-linear yield condition. International Journal for Numerical and Analytical Methods in Geomechanics, 12: 341–345. Fang, H.Y., and Hirst, T.J. 1970. Application of plasticity theory to slope stability problems. Highway Research Record No. 233, pp. 26–38. Fellenius, W. 1936. Calculation of the stability of earth dams. In Proceedings of the 2nd International Congress on Large Dams, Washington, D.C., Vol. 4, p. 445. Finn, W. 1967. Application of limit plasticity in soil mechanics. ASCE Journal of the Soil Mechanics and Foundations Division, 89(SM5): 101–119. Fredlund, D.G., and Krahn, J. 1977. Comparison of slope stability methods of analysis. Canadian Geotechnical Journal, 14: 429–439. Frontard, M. 1922. Cycloides de glissement des terres. Compte rendus hebdomadaires de l’academie des sciences de Paris, 174: 526–528. Giam, P., and Donald, I. 1991. Upper bound solutions to soil stability

Can. Geotech. J. Vol. 34, 1997 problems via general wedge methods. In Proceedings of the 7th International Conference on Computer Methods and Advances in Geomechanics, Cairns, Australia, May 6–10, pp. 475–480. Graham, J. 1973. Plasticity solutions to stability problems in sand. Canadian Geotechnical Journal, 11: 238–247. Gussman, P. 1982. Kinematic elements for soils and rocks. In Proceedings of the 4th International Conference on Numerical Methods in Geomechanics, Edmonton, Alta., pp. 47–52. Izbicki, R.J. 1981. Limit plasticity approach to slope stability problems. ASCE Journal of Engineering Mechanics, 107(GT12): 228–233. Jaky, J. 1936. Stability of earth slopes. In Proceedings of the International Conference on Soil Mechanics, Cambridge, Mass., Vol. 2, pp. 125–129. Janbu, N. 1973. Slope stability computation. In Embankment dam engineering. Edited by R.C. Hirschfeld and S.J. Poulos. John Wiley and Sons, New York, pp. 47–86. Karal, K. 1977a. Application of energy method. ASCE Journal of the Engineering Mechanics Division, 103(GT5): 381–397. Karal, K. 1977b. Energy method for soil stability analysis. ASCE Journal of the Engineering Mechanics Division, 103(GT5): 431–455. Kovari, K., and Fritze, P. 1984. Recent developments in the analysis and monitoring of rock slopes. In Proceedings of the 4th International Symposium on Landslides, Vo1. 1, Toronto, University of Toronto Press, Toronto. Li, K.S., and White, W. 1987. Rapid evaluation of the critical slip surface in slope stability problems. International Journal for Numerical and Analytical Methods in Geomechanics, 11: 449–473. Morgenstern, N.R., and Price, V.E. 1965. The analysis of the stability of general slip surfaces. Géotechnique, 15: 79–93. Nguyen, V.U. 1985. Determination of critical slope failure surface. ASCE Journal of the Engineering Mechanics Division, 105: 238–250. Sarma, S.K. 1973. Stability analysis of embankments and slopes. Géotechnique, 23: 423–433. Sarma, S.K. 1979. Stability analysis of embankments and slopes. ASCE Journal of the Geotechnical Engineering Division, 105(GT12): 1511–1524. Sloan, S.W. 1988. Lower bound limit analysis using finite elements and linear programming. International Journal for Numerical and Analytical Methods in Geomechanics, 12: 61–77. Sloan, S.W. 1989. Upper bound limit analysis using finite elements and linear programming. International Journal for Numerical and Analytical Methods in Geomechanics, 13: 263–282. Sokolovski, V.V. 1954. Statics of soil media. Translated by Jones, D.H. and Schofield, A.N. 1960, Butterworth, London. Spencer, E. 1967. A method of analysis of the stability of embankments assuming parallel interslice forces. Géotechnique, 17: 11–26. Sun, J.S. 1984. The numerical analysis of strip partition method. Chinese Journal of Geotechnical Engineering, 6: 1–12. (In Chinese.) Zhang, X.J., and Chen, W.F. 1987. Stability analysis of slopes with general nonlinear failure criterion. International Journal for Numerical and Analytical Methods in Geomechanics, 11: 33–50.