- Email: [email protected]
Safety analysis using finite elements
Computers and Georechnics, Vol. 22, No. 2, pp. 165-181, 1998
PII:SO266-352X(98)-00005-6
0 1998 Elsevier Science Ltd. All rights reserved Printed in Great Britain 0266-352X/98/$19.00+0.00
ELSEVIER
Safety Analysis Using Finite Elements M. M. Fariaf*
& D. J. Nay108
“Departamento de Engenharia Civil, Universidade de Brasilia, Brasilia, DF, CEP 70.910-900, Brazil bDepartment of Civil Engineering, University of Wales Swansea, Singleton Park, Swansea, U.K., SA2 8PP
(Received
20 September
1996; revised version received accepted 16 February 1998)
12 February
1998;
ABSTRACT A finite element based method for determining safety factors which has certain advantages over conventional limit equilibrium methods is described. The stress field produced by a jinite element analysis is used in conjunction with a conventional limit method to determine the safety factor. It is shown how the stress values at theJinite element integrating points may be interpolated to points on a potential slip surface selected by trial as in the conventional method. This involves three basic problems: (a) ident@cation of the element which contains a chosen point on the slip surface; (b) determination of the local co-ordinates of this point; and (c) determination of the element nodal stress values. Procedures to solve these problems are described. The problem of a vertical cut in clay is analysed and the results are compared with analytical solutions. 0 1998 Elsevier Science Ltd. All rights reserved
INTRODUCTION The occurrence of failures in embankments and excavated slopes is not uncommon. Collapses in these cases may, and often do, bring serious consequences to both life and property. The risks involved justify the need to adequately assess and assure safety. The choice of suitable methods of safety analysis is an important step. Limit equilibrium methods provide a powerful means of evaluating safety. They are based on the conventional sliding wedge or curved slip surface *To whom correspondence
should be addressed. 165
166
M. M. Farias and D. J. Naylor
methods. There are many different methods with different degrees of acceptability. One of the earliest is that of Fellenius [l] often referred to as the Swedish slip circle method. Bishop [2] made an important contribution to the study of safety of circular surfaces. Several methods were later developed for general slip surface analysis. The most widely used of these are Janbu’s method [3], Morgenstern and Price’s method [4] and Sarma’s method [5]. All methods of safety analysis have some qualifications associated with their use. Yet some limit equilibrium methods are so widely used that these qualifications are in danger of being forgotten. There are situations in which they can be misleading and should be treated with caution. This is because they generally provide either a rigorous or an approximate upper bound, so that the actual safety factor may be less than that computed. The analyses are therefore unconservative. Davis and Booker [6] have shown that the slip circle method may overestimate safety by a factor of 4 times when applied to a footing on level ground with cohesion that increases with depth. A further shortcoming of limit methods is that they provide no information about the failure mechanism. These methods do not take account of the stress-strain behaviour of the soil. Assumptions have to be made about the distribution of stresses and the possible shape of potential slip surfaces. There is no indication where yield starts or how it develops. The strength of finite element based methods lies in its potential to answer these questions. It is possible to use the finite element method to obtain both the safety factor and information about the collapse mechanism. It is, however, not easy to obtain a safety factor accurate to within the confidence limits achievable by limit equilibrium methods. A fine mesh is required, code capable of giving reliable results with the Mohr-Coulomb elasto-plastic model for loading states close to failure is needed, and finally it is usually necessary to carry out a set of analyses with c and tan4 progressively reduced by a factor which will become the safety factor when failure is eventually reached. These analyses become progressively more expensive as this factor is increased. Naylor [7] referred to this method as the “direct” finite element procedure. Naylor [7] then described what he called the “enhanced limit” method which provided an alternative way for obtaining the safety factor. Its essential feature is that the stress field produced by a finite element analysis is used in conjunction with a conventional limit method. It offers two major advantages over the direct procedure. First, only one finite element analyses is required. Second, for safety factors significantly greater than one the analysis will usually be relatively cheap and free of convergence problems. This paper looks into the detail of the application of the enhanced limit method. It shows how it can provide an economical means of calculating the
Safety analysis usingfinite elements
167
safety factor and at the same time provide information about the failure mechanism. It even shows that meaningful information can be obtained from a linear elastic analysis, although it should be emphasised at the outset that an elastic analysis should only be used to give a preliminary indication and should not replace a non-linear analysis.
ENHANCED
LIMIT METHOD
Figure 1 illustrates the method. Consider the trial slip surface AB in diagram (a). The dashed line in diagram (b) shows the variation of shear strength, r-, along AB. The actual shear stress distribution is indicated by the full line for non-linear behaviour. A properly constituted non-linear analysis should produce such a line lying entirely under or on the strength line. However, for linear analyses, the shear stresses represented by the chain-dotted line may exceed the strength line in regions of potential yield. The safety factor for this trial surface is defined as the ratio of the area under the strength curve to the area under the computed shear stress curve, i.e.:
(1)
This definition is equivalent to that used in conventional limit analyses where, at least for circular slips, it is usually expressed as the ratio of righting to disturbing moment. The stress (effective in the present context) at a point along the potential slip surface AB is defined by the stress components (a,, aY, rXY).These must LCCSI
F-i T
7,
@@mm
k-4 .-.-._
1”
(nm-!inear)
‘.
‘.._
.-.-.- ..__._.
_.
.
. . . . _. .._._._... -.--
7” (linear elastic)
A
(a) Fig. 1. Shear stress and strength
B
(W distribution
along a trial slip surface.
I
168
M. M. Farias and D. J. Naylor
be transformed by the following equations in order to give the normal and shear stress components along the surface: a, = o;, sin2 0 + aYcos2 8 - 2t,, sin Ocos 8
(2)
rn = (OX- aY) sin 8 cos 19+ rXY(cos2 8 - sin2 0)
(3)
The shear strength at the point may be computed using Mohr-Coulomb: rf= c’+cr’tanf$
(4)
The method is conceptually simple. The main difficulty lies in the determination of the stress components for a finite number of points along the trial surface. Naylor [7] used contour plots of the three stress components, over which the trial slip surface was superimposed. Then the stress components were interpolated manually. This process involved considerable work. To make the method practical, automation was clearly necessary. This was done by Farias [8]. The procedure is as follows.
DETERMINATION
OF STRESS COMPONENTS
The stresses at the points defining a potential slip surface could be computed directly from the displacement formulation during the finite element analysis in the same way the Gauss point stress values are computed. In this case the slip surface has to be known a priori which in general it will not be. Indeed, the main advantage of the finite element analysis is the identification of the failure mechanism as well as the location of the most feasible zone inside which the potential slip surface lies. Besides, the displacement formulation can give misleading stress results for points very close to the element borders. The solution is to adopt some kind of “post-processing” technique to compute the stresses along the surface. The stress at a given point P(x, y) can be computed by interpolating the nodal stresses in the finite element to which the point belongs. Let the stress to be computed at point P be represented by a row matrix with m components, e.g. cp = [ax, aY, rX,,]. Let c?” be an [n x m] matrix of nodal stresses, where IZ is the number of nodes of the element. The following interpolating can be used:
where N=[Ni,Nz,... , N,] is a row matrix of interpolating (6, r) are the local co-ordinates of point P.
functions and
Safety analysis using finite elements
In order to perform the interpolation solved: l
l
l
169
in Eqn (5), three problems have to be
given a general point expressed by its Cartesian co-ordinates, the element to which the point belongs must first be identified; once the element is identified, the local co-ordinates of the point must be computed in order to evaluate the interpolation functions; finally the element nodal values of the stress components must be found, since these are generally computed at the Gauss integration points.
The solution to the problems is explained separately in the following sections. Element identification
Given a finite element mesh and a general point P(x, y), it is required to identify the element inside which the point lies. A detailed discussion of various approaches to the problem of identification of a point with respect to a polygon has been given by Nordbeck and Rydsteadt [9]. Sloan [lo] describes the implementation of an improved version of the Nordbeck and Rydsteadt algorithm for linear-sided polygons of arbitrary shape. Sloan’s method is based on the so called “enlarged orientation theorem” which states that: 1. A point is inside a polygon if it is closer to a side than to a vertex, and the triangle generated by the vertices of this side and the point taken in an anti-clockwise order has a positive area. 2. A point is inside a polygon if it is closer to a vertex than to a side and this vertex is concave. All elements used in finite element analysis are convex polygons. If only straight-sided finite elements were used, then assertion (1) is sufficient and a simpler implementation could be used. Figure 2 illustrates the process for a quadrilateral element. For each side of a given element, a triangle connecting the corner nodes, Ci(xi, yi) and Ci+i(xi+i, yi+i) and point P(x, _y) can be defined. Two possibilities are depicted: (a) If the point lies inside the element the connectivity of all triangles previously defined will be listed in an anti-clockwise order [Fig. 2(a)]; (b) If the point lies outside the element the connectivity of at least one of the triangles so defined will be listed in a clockwise order [Fig. 2(b)]. The connectivity order is given by the sign of the triangle area, A=-
1 Xi-X 2 X-Y
Xi+1 -X
Yi+1-Y
(6)
M. 44. Farias and D. J. Naylor
170
(a)
(b)
Fig. 2. (a) Point inside element;
(b) point outside element.
If the area is positive, the triangle is listed in an anti-clockwise order. The opposite is true for a negative area. If the area is zero, the point lies on one of the element’s sides. For domains with curvilinear boundaries only assertion (2) of Sloan’s theorem applies and the algorithm may break down for points close to the element sides. Angelov and Manoach [l l] have extended Sloan’s method to account for this case for application to finite element mesh generation problems. Local co-ordinates The transformation from local (et q) to global co-ordinates (x, Y) is easily achieved by interpolation of the element’s nodal co-ordinates (xi, . . . , x,), (Yl,. . . ,un):
X =
2 j=l
Nj(t,
T)Xj
(7)
Y =
2
Nj(67 17)Yj
(8)
j=l
where II is the number of element nodes and Nj(& q) is the shape function for node j evaluated at (6, n). However, the inverse transformation from global since the interpola(x, y) to local (6, q ) co - or d’ma t es is not straightforward, tion functions are not expressed in terms of the global co-ordinates (x, y). A change in global position (Ax, Ay) may be approximated to the first order by:
Safety analysis usingJinite elements
171
(11) The matrix in the relation above represents the transpose of the Jacobian matrix, JT. Inversion gives explicit expressions for the change in local co-ordinates:
An iterative procedure can be used to determine the local co-ordinates (c, II) corresponding to the global co-ordinates (x, u). This is outlined in Table 1. The proposed algorithm was tested for the curved eight noded isoparametric elements shown in Fig. 3. The sides of the elements are described by parabolas according to the eight noded quadrilateral element shape functions. The co-ordinates of 36 points in the upper right quarter of the element were first computed using Eqn (7) and Eqn (8), in which the values of 6 and n were varied to from 0.0 to 1.0 in intervals of 0.2. The local co-ordinates were than computed back using the proposed algorithm and the
Iterative Step Step Step Step Step
1: 2: 3: 4: 5:
TABLE 1 scheme for the determination
of local co-ordinates
i=O Initialise. Choose initial guess. (4i, Vi) = (O,O) Evaluate shape functions at ({i, vi). Compute the global co-ordinates (xi, yi), using Eqns (7) and (8). Compute the necessary change in global position in order to approximate
the point
P(X, _V).: (AX, Ay) = (X - XL,Y-9 J’i) Step 6: Step 7: Step 8:
Step 9:
EXIT IF(Ax2 + AJ’~)~‘~< tol; Check convergence. Compute the corresponding change in local co-ordinates (A.$, An), using Eqn (12). Update the local co-ordinates: {i+i = 6i + At Q+I = ni + An Make i = i+ 1 and go to Step 3.
172
M. M. Farias and D. J. Naylor
pmposad algorithm Iterations
1 x 0
0 1
5:
A5
Fig. 3. Test of the proposed
algorithm.
number of iterations was recorded. A tolerance of lop6 was used for the computation of the global co-ordinates. The number of iterations necessary for convergence is also shown in Fig. 3. For points along a straight line 0, = constant for 6 = 0 or x = constant for r] = 0, in this example) only one iteration is necessary, since the relation between global and local co-ordinates becomes linear. Therefore no iterations will be necessary for linear sided elements. However, the number of iterations necessary for convergence increases as the line joining points corresponding to constant values of < or q gets more curved. Stress interpolation Because of the nature of the displacement variation assumed in the finite element method, the stresses computed are generally discontinuous across element boundaries. They should of course not be in zones of constant or continuously varying material properties, but could be when there are discontinuities in the properties across the boundaries. For numerically integrated elements, the best stress values are generally those computed at the integration points [12]. It is convenient, however, to have the stresses at the nodes as the element shape functions can then be used to calculate the stresses at any required (e, n) location (on the trial surface in this instance) within the element. These nodal stresses are best not calculated directly as the nodes are known to be poor stress sampling points. Instead a smoothing scheme is preferred as will now be outlined. To avoid discontinuities at element boundaries analysts have often resorted to the averaging of the nodal values computed directly from the finite element displacement formulation. Although economical and simple this solution can be highly unsatisfactory in same cases. Hinton and
Safety analysis usingjnite
elements
173
Campbell [13] have reported errors of orders of magnitude for the analysis of an infinitely long thick incompressible plate. In many cases, it is possible to get a better stress picture using a least square stress smoothing of the stresses calculated at the integrating points. This has been adopted by Hinton [14] for linear, parabolic and cubic isoparametric two and three dimensional elements. A similar approach has been applied to the smoothing of experimental data by Hinton and Irons [15]. Oden [16] and Oden and Brauchli [17] have applied the method to the calculation of the so called “consistent conjugate stresses”. The process is described as a mixed u - 0 formulation by Zienkiewicz and Taylor [18], who also call it “projection” or “variational recovery”. This smoothing technique is presented next. Let the smoothed stresses, D*, at any point within an element be expressed as CT*= N,cT
(13)
where 3 is a [n x m] matrix with m smoothed stress components for n element nodes, and N, = [Ni, N2, . . . , NJ ’ is a vector of smoothing shape frictions evaluated at the point. The smoothed nodal stresses 6 are the unknowns of the problem. The smoothing functions N,, may be of different order from the shape functions N used in the finite element analysis. The smoothed stresses (T*should approximate in a weak sense the stresses, 6, computed from the displacement formulation, using the shape functions N. This approximation may be written as
J- N;(c*
&)dS2 = 0
n
(14)
in which the integral is taken over the volume of the element. Substituting Eqn (13) into Eqn (14) gives
(W or MC=P
(16)
The “smoothing matrix” on the left-hand side of Eqn (16) is often found in dynamics problems as the “mass matrix”. Both the mass matrix and the
M. M. Farias and D. J. Naylor
174
“force vector” on the right-hand side can easily be evaluated using numerical integration. They can be assembled into a “global mass matrix” and a “global force vector” and solved for the whole domain. This process produces unique nodal stress values and is generally referred to as “global smoothing”. An iterative process to solve Eqn (16) using a “diagonalized” or lumped mass matrix is described in [ 181. If the smoothing is performed over individual elements it is referred to as “local smoothing”. This is equivalent to an exact least square fit to the unsmoothed integration point values [ 131. For a 2 x 2 Gauss integration rule the method can be interpreted as a bilinear fit through the Gauss point stresses. The locally smoothed nodal values are not unique and average values are often computed. Several error measures may be defined to give an indication of the order of the difference between computed and smoothed stress. These are often used for mesh refinements in adaptive procedures [19]. A simple and direct measure is the so called L2 norm. For a stress error (e, = D* - 6) the L2 error norm is defined as:
(17)
SEARCH
FOR
THE
CRITICAL
The mobilisation of shear strength stress ratio” (OSR) defined as:
SLIP
SURFACE
at a point can be indicated
OSR=s
4
by the “over-
(18)
where 0d = oi - ai and 4 is its failure value for the same elective mean stresses ai = (a; + 4)/2. The over-stress ratio is generally a little higher than the alternative “stress level” measure defined in the same way but with o-$ taken at the corresponding ai rather than ai. For a properly constituted non-linear program the OSR should not exceed unity. For linear elastic analysis the over-stress ratio can be greater than one. Even this case important information about the development of failure can be obtained. A potential slip surface would be expected to pass through the zones of highest OSR. Contour plots of over-stress ratio can therefore be
Safety analysis usingjinite elements
175
used to confine the search of a critical surface to a well defined region of the domain. Further important information on the development of failure surfaces can be obtained from potential slip directions plots. The slip directions are inclined at f(45” - 9/2) to the major principal stress, dt, where \I, is the dilatancy angle. For safety factors close to 1.0 the failure surface should be tangential to the slip directions. However, if the safety factor is high these will not generally coincide with the direction of the eventual slip surface. This is because the principal stress directions change as the soil strength is mobilised to cause failure. The use of OSR contours and slip direction plots provide a powerful means of detecting potential slip surfaces. An interactive graphics program greatly reduces the effort involved in these operations.
DETECTION
OF PHYSICAL
FAILURE
As failure approaches the yield zones spread and eventually a collapse mechanism is developed. Large displacements are expected for points in the slip area. In finite element analysis the proximity of failure can be associated with nonconvergence of iterative processes. The number of iterations tends to increase and the size of the load increments may be drastically reduced in the case of automatic step algorithms. Some limit must be imposed on the number of iterations and step size. For tangential formulations in which no iterations are used failure is generally indicated by excessive residual values. It is obviously important to identify when a numerical breakdown or program stop corresponds to a real failure of the structure. The following procedures can help to identify a physical failure: (a) over-stress contours; (b) deviatoric strain contours; (c) displacement increment (velocity) plots; (d) eigenvalue analysis. Over-stress contours can identify failure by indicating yield zones which extend to the borders of the domain and sometimes isolate two portions of the body. Deviatoric strain plots can reinforce this picture, as can velocity plots which typically show zones with big displacements indicating a rigid body movement of part of the mesh. The analysis of the eigenvalues and eigenvectors of the stiffness matrix prior to breakdown can give further indications of physical failure. If a mechanism is developing the lowest eigenvalue of the stiffness matrix should be small. The eigenvectors of the stiffness matrix represent possible displacement modes. The displacement modes associated with the lowest eigenvalues may disclose the nature of the failure mechanism.
M. M. Farias and D. J. Naylor
176
A TEST CASE A vertical cut in clay was chosen to demonstrate the points previously discussed. The problem is illustrated in Fig. 4. The height of the cut is H = 10 m. The clay was assumed to be a Tresca material with a cohesion c = 50 kPa, a friction angle 4 = 0” and a bulk unit weight y = 20 kN rnp3. The domain was discretized in 63 eight-noded elements numbered from the bottom left to the top right corner. A rather wobbly mesh was chosen to fully test the proposed procedures in the context of a real problem. A geostatic initial stress state was assumed with K, = 1.0. The excavation was simulated by the application of a surface traction representing the initial horizontal stresses along the vertical face. Theoretical safety factor
A slip plane at 45” was defined by a sequence of points spaced 1.Om apart, starting at the bottom left point A (Fig. 4). The elements inside which the points lie, as well as the local co-ordinates were successfully computed using the proposed procedures. An analysis of the static equilibrium, of the 45” wedge gives an upper bound solution to the safety factor: (19)
where H,, is the critical height, Q = 4.0 and y is the total unit weight of the soil. The value of c/yH was chosen as 0.25, so as to have an incipient failure SCALE(m) 012
3
Smooth
Rough
Fig. 4. Vertical cut in clay.
Safety analysis using finite elements
177
for the upper bound value of CLThe exact value of (Yis not yet known but is known to lie in the range 3.64 -Ca! < 3.83 [20,21] for different shapes of slip surfaces. The upper bound of 3.83 is associated with a curved slip surface as compared with the value of 4.0 associated with the 45” plane. Linear elastic analysis
A preliminary linear elastic analysis was carried out. The elastic parameters were taken as E = 10000 kPa and u = 0.30. Figure 5 shows the distribution of shear stress and shear strength along the 45” slip surface. The stress components were computed using the three schemes previously discussed, i.e.: (a) local smoothing; (b) local smoothing plus nodal averaging; and (c) global smoothing. The over-stress ratio exceeds unity at the toe of the cut thus indicating potential local yielding. However, an interesting picture arises from the safety factor values: (a) F, = 1.088 with local smoothing; (b) F, = 1.085 for local smoothing plus nodal averaging; and (c) F, = 1.098 for global smoothing. This illustrates the insensitivity of the result to the type of smoothing used, but more significantly suggests that even an elastic linear analysis can provide meaningful information about the safety factor and failure mechanism. The change in safety factor due to different smoothing techniques proved to be irrelevant in this example. Perhaps in an example with different material properties, these differences could be more significant. The global 140 120 -Average
100 f s B
-A-
Fs=l.O85
Global Fs=l.o98
80
E 111
60--
i f
>c 40 --
x
20 -06 0
2
4
6
10
8
12
14
Distance (m) Fig. 5. Distribution
of shear stress and shear strength
along the slip plane. Linear elastic.
M. M. Farias and D. J. Naylor
178
smoothing solution showed some slight differences near the toe when compared to the other schemes. This is due to the stress concentration at the bottom of the cut as a consequence of the fixity conditions at point A. This concentration is also indicated by the least square error as measured by Eqn (17). The error for element 1 was the largest. Elasto-plastic
analysis
A Tresca elasto-plastic model was used. This was applied using a version of the critical state model (c.s.m.) described by Farias [S]. Advantage was taken of the versatility of this model to adapt it to incorporate the Tresca yield surface with c = 50 kPa and @= 0”. The load was applied incrementally using an automatic step size algorithm. The algorithm failed to converge when 99% of the total load had been applied. At this stage the load increments were in the order of lO-5 of the total load. Figure 6(a) shows a plot of OSR contours at the last converged load step prior to failure. The yield zone extends along a shear band inclined roughly at 45”. In this rather wide zone the over-stress ratio is within the range 0.95 to 1.O. A narrower band can be obtained from the deviatoric strain contours in Fig. 6(b). It would be possible to choose from a number of surfaces within this zone, but the velocity vectors prior to failure shown in Fig. 7, while generally aligned along the 45” slip plane, show a preference for a slightly curved surface. This is particularly evident at the toe of the cut. Figure 7 shows the displacement increments at the last converged time step. These concentrate almost exclusively in the upper wedge. The displacements are roughly uniform except along the slip surface where the strain localisation is evident. There are clear indications that a mechanism is developing.
e
(‘4
(4 Fig. 6. (a) Over-stress
ratio and (b) deviatoric
deformation
at end of analysis.
Safety analysis using finite elements
Fig. 7. Displacement
increments
179
prior to failure. Elasto-plastic
analysis.
The distribution of shear stress and shear strength along the 45” slip plane is shown in Fig. 8. These are plotted for three different smoothing techniques: (a) local smoothing; (b) local smoothing plus nodal averaging; and (c) global smoothing. Although some local variations in the plotted curves can be observed the safety factor hard1 changes when these are inte rated. The computed safety factors were $L 1.081, I$” = 1.083 and 2 ‘) = 1.089 respectively. It can be noted that the shear strength was neirly totally mobilised along most of the surface, except for a small region close to the 60
t
0,30
-e
110
-c-Average
b
Local Fs=l.O61 Fs=l.O63
ii
P) 20 c 10
+
Gbbal Fs=l.O69
+
Strength
04 0
2
4
6
6
10
12
14
Distance (m) Fig. 8. Distribution
of shear stress and shear strength
along the slip plane. Elasto-plastic.
180
M. M. Farias and D. J. Naylor
toe. In some points the shear stress curve slightly exceeds that of the shear strength. This is due to the smoothing techniques, rather than to any violation of the yield criteria at Gauss points where the stresses are actually computed. Clearly the actual safety factor at the end of the analysis was F,= 1.0since failure was well defined. The value of F,= 1.08 is attributed to the assumption of the 45” plane. It is interesting that it is within 1% of the value obtained from the linear elastic analysis.
CONCLUSIONS The iterative approach proposed for computing the local co-ordinates of a point on the slip surface was found to be extremely efficient. No significant difference in the safety factor resulted from the three different approaches based on a least squares fit for the extrapolation of the Gauss point stresses to the nodes. These approaches were: (a) local smoothing; (b) local smoothing plus nodal averaging; and (c) global smoothing. The local smoothing technique gives the fastest results. The commonly used procedure of taking nodal averages to avoid discontinuities across element borders proved to be more expensive than the use of a global smoothing technique. An approximate identification of the failure mechanism can be obtained from over-stress ratio contours. These indicate yielded zones which at failure extend to the borders of the domain and sometimes isolate two portions of the body. Velocity plots can show zones associated with relatively big displacements and thereby identify a rigid body movement. Deviatoric strain plots show the localization of shear bands as failure approaches. The test case of a vertical cut in clay showed that elasto-plastic analyses can accurately reproduce failure. This study also shows that a remarkably good estimate of the safety factor can in this case be obtained from a linear elastic analysis, which suggests that such an analysis has a useful preliminary role in assessing safety. The example cited above shows the great potential of using results of finite element analyses to assess safety. It allows not only the safety factor to be quantified accurately, but more importantly it provides vital information about the failure mechanisms. REFERENCES 1. Fellenius, W., Calculation of the stability of earth dams. Proc. 2nd Congress on Large
Dams,
1936, 4, 445.
Safety analysis using finite eiements
181
2. Bishop, A. W., The use of the slip circle in the stability analysis of slopes. Geotechnique, 1955, 5(l), 7-17. 3. Janbu, N., Application of composite slip surfaces for stability analysis. Proc. European Conference on Stability of Earth Dams, 1954, 3, 43349. 4. Morgenstern, N. R. and Price, V. E., The analysis of the stability of general slip surfaces. Geotechnique, 1965, 15(l), 79-93. 5. Sarma, S. K., Stability analysis of embankments slopes. Journal of the Geotechnical Engineering Division-ASCE, 1979, 105(GT12), 15 1l-l 524. 6. Davis, E. H. and Booker, J. R., The effect of increasing strength with depth on the bearing capacity of clay. Geotechnique, 1989, 23(4), 551-563. 7. Naylor, D. L, Finite elements and slope stability. In Numerical Methods in Geomechanics, D. Reidel Publishing Company, J. B. Martins (eds), 1982, 229244. 8. Farias, M. M., Numerical analysis of clay core dams. Ph.D thesis, University of Wales, University College of Swansea, 1993. 9. Nordbeck, S. and Rydsteadt, B., Computer cartography point-in-polygon programs. BIT, 1967, 7, 39. 10. Sloan, S. W., A point-in-polygon program. Advances in Engineering Software, 1985, 7, 45-47. 11. Angelov, T. A. and Manoach, E. S., A point-in-domain identification program. Advances in Engineering Software, 1989, 11(2), 99-105. 12. Too, J. M., Two dimensional plate and finite prism isoparametric elements and their applications. Ph.D thesis, University of Wales, University College of Swansea, 197 1. 13. Hinton, E. and Campbell, J. S., Local and global smoothing of discontinuous finite element functions using a least square method. International Journal for Numerical Methods in Engineering, 1974, 8, 461480. 14. Hinton, E., Least squares analysis using finite elements. M.Sc. thesis, University of Wales, University College of Swansea, 1968. 15. Hinton, E. and Irons, B. M., Least squares smoothing of experimental data using finite elements. Journal of the British Society for Strain Measurement, 1968, 14,24-27. 16. Oden, J. T., Finite Elements of Non-Linear Continua. McGraw-Hill, New York, 1971. 17. Oden, J. T. and Brauchli, J., On the calculation of consistent stress distributions in finite element applications. International Journal for Numerical Methods in Engineering, 197 1, 3, 3 17-325. 18. Zienkiewicz, 0. C. and Taylor, R. L., The Finite Element Method, 4th edn. McGraw-Hill, New York, 1990. 19. Zienkiewiez, 0. C. and Zhu, J. Z., A simple error estimator and adaptive procedure for practical engineering analysis. International Journal for Numerical Methods in Engineering, 1987, 24, 337-357. 20. Josselin de Jong, G., Application of the calculus of variations to the vertical cut off in cohesive frictionless soil. Geotechnique, 1980, 30(l), l-16. 21. Pastor, J., Analyse limite: determination de solutions statiques completes. Revice de Mechanique Appliquee, 1978, 2(2), 167-197.
PII:SO266-352X(98)-00005-6
0 1998 Elsevier Science Ltd. All rights reserved Printed in Great Britain 0266-352X/98/$19.00+0.00
ELSEVIER
Safety Analysis Using Finite Elements M. M. Fariaf*
& D. J. Nay108
“Departamento de Engenharia Civil, Universidade de Brasilia, Brasilia, DF, CEP 70.910-900, Brazil bDepartment of Civil Engineering, University of Wales Swansea, Singleton Park, Swansea, U.K., SA2 8PP
(Received
20 September
1996; revised version received accepted 16 February 1998)
12 February
1998;
ABSTRACT A finite element based method for determining safety factors which has certain advantages over conventional limit equilibrium methods is described. The stress field produced by a jinite element analysis is used in conjunction with a conventional limit method to determine the safety factor. It is shown how the stress values at theJinite element integrating points may be interpolated to points on a potential slip surface selected by trial as in the conventional method. This involves three basic problems: (a) ident@cation of the element which contains a chosen point on the slip surface; (b) determination of the local co-ordinates of this point; and (c) determination of the element nodal stress values. Procedures to solve these problems are described. The problem of a vertical cut in clay is analysed and the results are compared with analytical solutions. 0 1998 Elsevier Science Ltd. All rights reserved
INTRODUCTION The occurrence of failures in embankments and excavated slopes is not uncommon. Collapses in these cases may, and often do, bring serious consequences to both life and property. The risks involved justify the need to adequately assess and assure safety. The choice of suitable methods of safety analysis is an important step. Limit equilibrium methods provide a powerful means of evaluating safety. They are based on the conventional sliding wedge or curved slip surface *To whom correspondence
should be addressed. 165
166
M. M. Farias and D. J. Naylor
methods. There are many different methods with different degrees of acceptability. One of the earliest is that of Fellenius [l] often referred to as the Swedish slip circle method. Bishop [2] made an important contribution to the study of safety of circular surfaces. Several methods were later developed for general slip surface analysis. The most widely used of these are Janbu’s method [3], Morgenstern and Price’s method [4] and Sarma’s method [5]. All methods of safety analysis have some qualifications associated with their use. Yet some limit equilibrium methods are so widely used that these qualifications are in danger of being forgotten. There are situations in which they can be misleading and should be treated with caution. This is because they generally provide either a rigorous or an approximate upper bound, so that the actual safety factor may be less than that computed. The analyses are therefore unconservative. Davis and Booker [6] have shown that the slip circle method may overestimate safety by a factor of 4 times when applied to a footing on level ground with cohesion that increases with depth. A further shortcoming of limit methods is that they provide no information about the failure mechanism. These methods do not take account of the stress-strain behaviour of the soil. Assumptions have to be made about the distribution of stresses and the possible shape of potential slip surfaces. There is no indication where yield starts or how it develops. The strength of finite element based methods lies in its potential to answer these questions. It is possible to use the finite element method to obtain both the safety factor and information about the collapse mechanism. It is, however, not easy to obtain a safety factor accurate to within the confidence limits achievable by limit equilibrium methods. A fine mesh is required, code capable of giving reliable results with the Mohr-Coulomb elasto-plastic model for loading states close to failure is needed, and finally it is usually necessary to carry out a set of analyses with c and tan4 progressively reduced by a factor which will become the safety factor when failure is eventually reached. These analyses become progressively more expensive as this factor is increased. Naylor [7] referred to this method as the “direct” finite element procedure. Naylor [7] then described what he called the “enhanced limit” method which provided an alternative way for obtaining the safety factor. Its essential feature is that the stress field produced by a finite element analysis is used in conjunction with a conventional limit method. It offers two major advantages over the direct procedure. First, only one finite element analyses is required. Second, for safety factors significantly greater than one the analysis will usually be relatively cheap and free of convergence problems. This paper looks into the detail of the application of the enhanced limit method. It shows how it can provide an economical means of calculating the
Safety analysis usingfinite elements
167
safety factor and at the same time provide information about the failure mechanism. It even shows that meaningful information can be obtained from a linear elastic analysis, although it should be emphasised at the outset that an elastic analysis should only be used to give a preliminary indication and should not replace a non-linear analysis.
ENHANCED
LIMIT METHOD
Figure 1 illustrates the method. Consider the trial slip surface AB in diagram (a). The dashed line in diagram (b) shows the variation of shear strength, r-, along AB. The actual shear stress distribution is indicated by the full line for non-linear behaviour. A properly constituted non-linear analysis should produce such a line lying entirely under or on the strength line. However, for linear analyses, the shear stresses represented by the chain-dotted line may exceed the strength line in regions of potential yield. The safety factor for this trial surface is defined as the ratio of the area under the strength curve to the area under the computed shear stress curve, i.e.:
(1)
This definition is equivalent to that used in conventional limit analyses where, at least for circular slips, it is usually expressed as the ratio of righting to disturbing moment. The stress (effective in the present context) at a point along the potential slip surface AB is defined by the stress components (a,, aY, rXY).These must LCCSI
F-i T
7,
@@mm
k-4 .-.-._
1”
(nm-!inear)
‘.
‘.._
.-.-.- ..__._.
_.
.
. . . . _. .._._._... -.--
7” (linear elastic)
A
(a) Fig. 1. Shear stress and strength
B
(W distribution
along a trial slip surface.
I
168
M. M. Farias and D. J. Naylor
be transformed by the following equations in order to give the normal and shear stress components along the surface: a, = o;, sin2 0 + aYcos2 8 - 2t,, sin Ocos 8
(2)
rn = (OX- aY) sin 8 cos 19+ rXY(cos2 8 - sin2 0)
(3)
The shear strength at the point may be computed using Mohr-Coulomb: rf= c’+cr’tanf$
(4)
The method is conceptually simple. The main difficulty lies in the determination of the stress components for a finite number of points along the trial surface. Naylor [7] used contour plots of the three stress components, over which the trial slip surface was superimposed. Then the stress components were interpolated manually. This process involved considerable work. To make the method practical, automation was clearly necessary. This was done by Farias [8]. The procedure is as follows.
DETERMINATION
OF STRESS COMPONENTS
The stresses at the points defining a potential slip surface could be computed directly from the displacement formulation during the finite element analysis in the same way the Gauss point stress values are computed. In this case the slip surface has to be known a priori which in general it will not be. Indeed, the main advantage of the finite element analysis is the identification of the failure mechanism as well as the location of the most feasible zone inside which the potential slip surface lies. Besides, the displacement formulation can give misleading stress results for points very close to the element borders. The solution is to adopt some kind of “post-processing” technique to compute the stresses along the surface. The stress at a given point P(x, y) can be computed by interpolating the nodal stresses in the finite element to which the point belongs. Let the stress to be computed at point P be represented by a row matrix with m components, e.g. cp = [ax, aY, rX,,]. Let c?” be an [n x m] matrix of nodal stresses, where IZ is the number of nodes of the element. The following interpolating can be used:
where N=[Ni,Nz,... , N,] is a row matrix of interpolating (6, r) are the local co-ordinates of point P.
functions and
Safety analysis using finite elements
In order to perform the interpolation solved: l
l
l
169
in Eqn (5), three problems have to be
given a general point expressed by its Cartesian co-ordinates, the element to which the point belongs must first be identified; once the element is identified, the local co-ordinates of the point must be computed in order to evaluate the interpolation functions; finally the element nodal values of the stress components must be found, since these are generally computed at the Gauss integration points.
The solution to the problems is explained separately in the following sections. Element identification
Given a finite element mesh and a general point P(x, y), it is required to identify the element inside which the point lies. A detailed discussion of various approaches to the problem of identification of a point with respect to a polygon has been given by Nordbeck and Rydsteadt [9]. Sloan [lo] describes the implementation of an improved version of the Nordbeck and Rydsteadt algorithm for linear-sided polygons of arbitrary shape. Sloan’s method is based on the so called “enlarged orientation theorem” which states that: 1. A point is inside a polygon if it is closer to a side than to a vertex, and the triangle generated by the vertices of this side and the point taken in an anti-clockwise order has a positive area. 2. A point is inside a polygon if it is closer to a vertex than to a side and this vertex is concave. All elements used in finite element analysis are convex polygons. If only straight-sided finite elements were used, then assertion (1) is sufficient and a simpler implementation could be used. Figure 2 illustrates the process for a quadrilateral element. For each side of a given element, a triangle connecting the corner nodes, Ci(xi, yi) and Ci+i(xi+i, yi+i) and point P(x, _y) can be defined. Two possibilities are depicted: (a) If the point lies inside the element the connectivity of all triangles previously defined will be listed in an anti-clockwise order [Fig. 2(a)]; (b) If the point lies outside the element the connectivity of at least one of the triangles so defined will be listed in a clockwise order [Fig. 2(b)]. The connectivity order is given by the sign of the triangle area, A=-
1 Xi-X 2 X-Y
Xi+1 -X
Yi+1-Y
(6)
M. 44. Farias and D. J. Naylor
170
(a)
(b)
Fig. 2. (a) Point inside element;
(b) point outside element.
If the area is positive, the triangle is listed in an anti-clockwise order. The opposite is true for a negative area. If the area is zero, the point lies on one of the element’s sides. For domains with curvilinear boundaries only assertion (2) of Sloan’s theorem applies and the algorithm may break down for points close to the element sides. Angelov and Manoach [l l] have extended Sloan’s method to account for this case for application to finite element mesh generation problems. Local co-ordinates The transformation from local (et q) to global co-ordinates (x, Y) is easily achieved by interpolation of the element’s nodal co-ordinates (xi, . . . , x,), (Yl,. . . ,un):
X =
2 j=l
Nj(t,
T)Xj
(7)
Y =
2
Nj(67 17)Yj
(8)
j=l
where II is the number of element nodes and Nj(& q) is the shape function for node j evaluated at (6, n). However, the inverse transformation from global since the interpola(x, y) to local (6, q ) co - or d’ma t es is not straightforward, tion functions are not expressed in terms of the global co-ordinates (x, y). A change in global position (Ax, Ay) may be approximated to the first order by:
Safety analysis usingJinite elements
171
(11) The matrix in the relation above represents the transpose of the Jacobian matrix, JT. Inversion gives explicit expressions for the change in local co-ordinates:
An iterative procedure can be used to determine the local co-ordinates (c, II) corresponding to the global co-ordinates (x, u). This is outlined in Table 1. The proposed algorithm was tested for the curved eight noded isoparametric elements shown in Fig. 3. The sides of the elements are described by parabolas according to the eight noded quadrilateral element shape functions. The co-ordinates of 36 points in the upper right quarter of the element were first computed using Eqn (7) and Eqn (8), in which the values of 6 and n were varied to from 0.0 to 1.0 in intervals of 0.2. The local co-ordinates were than computed back using the proposed algorithm and the
Iterative Step Step Step Step Step
1: 2: 3: 4: 5:
TABLE 1 scheme for the determination
of local co-ordinates
i=O Initialise. Choose initial guess. (4i, Vi) = (O,O) Evaluate shape functions at ({i, vi). Compute the global co-ordinates (xi, yi), using Eqns (7) and (8). Compute the necessary change in global position in order to approximate
the point
P(X, _V).: (AX, Ay) = (X - XL,Y-9 J’i) Step 6: Step 7: Step 8:
Step 9:
EXIT IF(Ax2 + AJ’~)~‘~< tol; Check convergence. Compute the corresponding change in local co-ordinates (A.$, An), using Eqn (12). Update the local co-ordinates: {i+i = 6i + At Q+I = ni + An Make i = i+ 1 and go to Step 3.
172
M. M. Farias and D. J. Naylor
pmposad algorithm Iterations
1 x 0
0 1
5:
A5
Fig. 3. Test of the proposed
algorithm.
number of iterations was recorded. A tolerance of lop6 was used for the computation of the global co-ordinates. The number of iterations necessary for convergence is also shown in Fig. 3. For points along a straight line 0, = constant for 6 = 0 or x = constant for r] = 0, in this example) only one iteration is necessary, since the relation between global and local co-ordinates becomes linear. Therefore no iterations will be necessary for linear sided elements. However, the number of iterations necessary for convergence increases as the line joining points corresponding to constant values of < or q gets more curved. Stress interpolation Because of the nature of the displacement variation assumed in the finite element method, the stresses computed are generally discontinuous across element boundaries. They should of course not be in zones of constant or continuously varying material properties, but could be when there are discontinuities in the properties across the boundaries. For numerically integrated elements, the best stress values are generally those computed at the integration points [12]. It is convenient, however, to have the stresses at the nodes as the element shape functions can then be used to calculate the stresses at any required (e, n) location (on the trial surface in this instance) within the element. These nodal stresses are best not calculated directly as the nodes are known to be poor stress sampling points. Instead a smoothing scheme is preferred as will now be outlined. To avoid discontinuities at element boundaries analysts have often resorted to the averaging of the nodal values computed directly from the finite element displacement formulation. Although economical and simple this solution can be highly unsatisfactory in same cases. Hinton and
Safety analysis usingjnite
elements
173
Campbell [13] have reported errors of orders of magnitude for the analysis of an infinitely long thick incompressible plate. In many cases, it is possible to get a better stress picture using a least square stress smoothing of the stresses calculated at the integrating points. This has been adopted by Hinton [14] for linear, parabolic and cubic isoparametric two and three dimensional elements. A similar approach has been applied to the smoothing of experimental data by Hinton and Irons [15]. Oden [16] and Oden and Brauchli [17] have applied the method to the calculation of the so called “consistent conjugate stresses”. The process is described as a mixed u - 0 formulation by Zienkiewicz and Taylor [18], who also call it “projection” or “variational recovery”. This smoothing technique is presented next. Let the smoothed stresses, D*, at any point within an element be expressed as CT*= N,cT
(13)
where 3 is a [n x m] matrix with m smoothed stress components for n element nodes, and N, = [Ni, N2, . . . , NJ ’ is a vector of smoothing shape frictions evaluated at the point. The smoothed nodal stresses 6 are the unknowns of the problem. The smoothing functions N,, may be of different order from the shape functions N used in the finite element analysis. The smoothed stresses (T*should approximate in a weak sense the stresses, 6, computed from the displacement formulation, using the shape functions N. This approximation may be written as
J- N;(c*
&)dS2 = 0
n
(14)
in which the integral is taken over the volume of the element. Substituting Eqn (13) into Eqn (14) gives
(W or MC=P
(16)
The “smoothing matrix” on the left-hand side of Eqn (16) is often found in dynamics problems as the “mass matrix”. Both the mass matrix and the
M. M. Farias and D. J. Naylor
174
“force vector” on the right-hand side can easily be evaluated using numerical integration. They can be assembled into a “global mass matrix” and a “global force vector” and solved for the whole domain. This process produces unique nodal stress values and is generally referred to as “global smoothing”. An iterative process to solve Eqn (16) using a “diagonalized” or lumped mass matrix is described in [ 181. If the smoothing is performed over individual elements it is referred to as “local smoothing”. This is equivalent to an exact least square fit to the unsmoothed integration point values [ 131. For a 2 x 2 Gauss integration rule the method can be interpreted as a bilinear fit through the Gauss point stresses. The locally smoothed nodal values are not unique and average values are often computed. Several error measures may be defined to give an indication of the order of the difference between computed and smoothed stress. These are often used for mesh refinements in adaptive procedures [19]. A simple and direct measure is the so called L2 norm. For a stress error (e, = D* - 6) the L2 error norm is defined as:
(17)
SEARCH
FOR
THE
CRITICAL
The mobilisation of shear strength stress ratio” (OSR) defined as:
SLIP
SURFACE
at a point can be indicated
OSR=s
4
by the “over-
(18)
where 0d = oi - ai and 4 is its failure value for the same elective mean stresses ai = (a; + 4)/2. The over-stress ratio is generally a little higher than the alternative “stress level” measure defined in the same way but with o-$ taken at the corresponding ai rather than ai. For a properly constituted non-linear program the OSR should not exceed unity. For linear elastic analysis the over-stress ratio can be greater than one. Even this case important information about the development of failure can be obtained. A potential slip surface would be expected to pass through the zones of highest OSR. Contour plots of over-stress ratio can therefore be
Safety analysis usingjinite elements
175
used to confine the search of a critical surface to a well defined region of the domain. Further important information on the development of failure surfaces can be obtained from potential slip directions plots. The slip directions are inclined at f(45” - 9/2) to the major principal stress, dt, where \I, is the dilatancy angle. For safety factors close to 1.0 the failure surface should be tangential to the slip directions. However, if the safety factor is high these will not generally coincide with the direction of the eventual slip surface. This is because the principal stress directions change as the soil strength is mobilised to cause failure. The use of OSR contours and slip direction plots provide a powerful means of detecting potential slip surfaces. An interactive graphics program greatly reduces the effort involved in these operations.
DETECTION
OF PHYSICAL
FAILURE
As failure approaches the yield zones spread and eventually a collapse mechanism is developed. Large displacements are expected for points in the slip area. In finite element analysis the proximity of failure can be associated with nonconvergence of iterative processes. The number of iterations tends to increase and the size of the load increments may be drastically reduced in the case of automatic step algorithms. Some limit must be imposed on the number of iterations and step size. For tangential formulations in which no iterations are used failure is generally indicated by excessive residual values. It is obviously important to identify when a numerical breakdown or program stop corresponds to a real failure of the structure. The following procedures can help to identify a physical failure: (a) over-stress contours; (b) deviatoric strain contours; (c) displacement increment (velocity) plots; (d) eigenvalue analysis. Over-stress contours can identify failure by indicating yield zones which extend to the borders of the domain and sometimes isolate two portions of the body. Deviatoric strain plots can reinforce this picture, as can velocity plots which typically show zones with big displacements indicating a rigid body movement of part of the mesh. The analysis of the eigenvalues and eigenvectors of the stiffness matrix prior to breakdown can give further indications of physical failure. If a mechanism is developing the lowest eigenvalue of the stiffness matrix should be small. The eigenvectors of the stiffness matrix represent possible displacement modes. The displacement modes associated with the lowest eigenvalues may disclose the nature of the failure mechanism.
M. M. Farias and D. J. Naylor
176
A TEST CASE A vertical cut in clay was chosen to demonstrate the points previously discussed. The problem is illustrated in Fig. 4. The height of the cut is H = 10 m. The clay was assumed to be a Tresca material with a cohesion c = 50 kPa, a friction angle 4 = 0” and a bulk unit weight y = 20 kN rnp3. The domain was discretized in 63 eight-noded elements numbered from the bottom left to the top right corner. A rather wobbly mesh was chosen to fully test the proposed procedures in the context of a real problem. A geostatic initial stress state was assumed with K, = 1.0. The excavation was simulated by the application of a surface traction representing the initial horizontal stresses along the vertical face. Theoretical safety factor
A slip plane at 45” was defined by a sequence of points spaced 1.Om apart, starting at the bottom left point A (Fig. 4). The elements inside which the points lie, as well as the local co-ordinates were successfully computed using the proposed procedures. An analysis of the static equilibrium, of the 45” wedge gives an upper bound solution to the safety factor: (19)
where H,, is the critical height, Q = 4.0 and y is the total unit weight of the soil. The value of c/yH was chosen as 0.25, so as to have an incipient failure SCALE(m) 012
3
Smooth
Rough
Fig. 4. Vertical cut in clay.
Safety analysis using finite elements
177
for the upper bound value of CLThe exact value of (Yis not yet known but is known to lie in the range 3.64 -Ca! < 3.83 [20,21] for different shapes of slip surfaces. The upper bound of 3.83 is associated with a curved slip surface as compared with the value of 4.0 associated with the 45” plane. Linear elastic analysis
A preliminary linear elastic analysis was carried out. The elastic parameters were taken as E = 10000 kPa and u = 0.30. Figure 5 shows the distribution of shear stress and shear strength along the 45” slip surface. The stress components were computed using the three schemes previously discussed, i.e.: (a) local smoothing; (b) local smoothing plus nodal averaging; and (c) global smoothing. The over-stress ratio exceeds unity at the toe of the cut thus indicating potential local yielding. However, an interesting picture arises from the safety factor values: (a) F, = 1.088 with local smoothing; (b) F, = 1.085 for local smoothing plus nodal averaging; and (c) F, = 1.098 for global smoothing. This illustrates the insensitivity of the result to the type of smoothing used, but more significantly suggests that even an elastic linear analysis can provide meaningful information about the safety factor and failure mechanism. The change in safety factor due to different smoothing techniques proved to be irrelevant in this example. Perhaps in an example with different material properties, these differences could be more significant. The global 140 120 -Average
100 f s B
-A-
Fs=l.O85
Global Fs=l.o98
80
E 111
60--
i f
>c 40 --
x
20 -06 0
2
4
6
10
8
12
14
Distance (m) Fig. 5. Distribution
of shear stress and shear strength
along the slip plane. Linear elastic.
M. M. Farias and D. J. Naylor
178
smoothing solution showed some slight differences near the toe when compared to the other schemes. This is due to the stress concentration at the bottom of the cut as a consequence of the fixity conditions at point A. This concentration is also indicated by the least square error as measured by Eqn (17). The error for element 1 was the largest. Elasto-plastic
analysis
A Tresca elasto-plastic model was used. This was applied using a version of the critical state model (c.s.m.) described by Farias [S]. Advantage was taken of the versatility of this model to adapt it to incorporate the Tresca yield surface with c = 50 kPa and @= 0”. The load was applied incrementally using an automatic step size algorithm. The algorithm failed to converge when 99% of the total load had been applied. At this stage the load increments were in the order of lO-5 of the total load. Figure 6(a) shows a plot of OSR contours at the last converged load step prior to failure. The yield zone extends along a shear band inclined roughly at 45”. In this rather wide zone the over-stress ratio is within the range 0.95 to 1.O. A narrower band can be obtained from the deviatoric strain contours in Fig. 6(b). It would be possible to choose from a number of surfaces within this zone, but the velocity vectors prior to failure shown in Fig. 7, while generally aligned along the 45” slip plane, show a preference for a slightly curved surface. This is particularly evident at the toe of the cut. Figure 7 shows the displacement increments at the last converged time step. These concentrate almost exclusively in the upper wedge. The displacements are roughly uniform except along the slip surface where the strain localisation is evident. There are clear indications that a mechanism is developing.
e
(‘4
(4 Fig. 6. (a) Over-stress
ratio and (b) deviatoric
deformation
at end of analysis.
Safety analysis using finite elements
Fig. 7. Displacement
increments
179
prior to failure. Elasto-plastic
analysis.
The distribution of shear stress and shear strength along the 45” slip plane is shown in Fig. 8. These are plotted for three different smoothing techniques: (a) local smoothing; (b) local smoothing plus nodal averaging; and (c) global smoothing. Although some local variations in the plotted curves can be observed the safety factor hard1 changes when these are inte rated. The computed safety factors were $L 1.081, I$” = 1.083 and 2 ‘) = 1.089 respectively. It can be noted that the shear strength was neirly totally mobilised along most of the surface, except for a small region close to the 60
t
0,30
-e
110
-c-Average
b
Local Fs=l.O61 Fs=l.O63
ii
P) 20 c 10
+
Gbbal Fs=l.O69
+
Strength
04 0
2
4
6
6
10
12
14
Distance (m) Fig. 8. Distribution
of shear stress and shear strength
along the slip plane. Elasto-plastic.
180
M. M. Farias and D. J. Naylor
toe. In some points the shear stress curve slightly exceeds that of the shear strength. This is due to the smoothing techniques, rather than to any violation of the yield criteria at Gauss points where the stresses are actually computed. Clearly the actual safety factor at the end of the analysis was F,= 1.0since failure was well defined. The value of F,= 1.08 is attributed to the assumption of the 45” plane. It is interesting that it is within 1% of the value obtained from the linear elastic analysis.
CONCLUSIONS The iterative approach proposed for computing the local co-ordinates of a point on the slip surface was found to be extremely efficient. No significant difference in the safety factor resulted from the three different approaches based on a least squares fit for the extrapolation of the Gauss point stresses to the nodes. These approaches were: (a) local smoothing; (b) local smoothing plus nodal averaging; and (c) global smoothing. The local smoothing technique gives the fastest results. The commonly used procedure of taking nodal averages to avoid discontinuities across element borders proved to be more expensive than the use of a global smoothing technique. An approximate identification of the failure mechanism can be obtained from over-stress ratio contours. These indicate yielded zones which at failure extend to the borders of the domain and sometimes isolate two portions of the body. Velocity plots can show zones associated with relatively big displacements and thereby identify a rigid body movement. Deviatoric strain plots show the localization of shear bands as failure approaches. The test case of a vertical cut in clay showed that elasto-plastic analyses can accurately reproduce failure. This study also shows that a remarkably good estimate of the safety factor can in this case be obtained from a linear elastic analysis, which suggests that such an analysis has a useful preliminary role in assessing safety. The example cited above shows the great potential of using results of finite element analyses to assess safety. It allows not only the safety factor to be quantified accurately, but more importantly it provides vital information about the failure mechanisms. REFERENCES 1. Fellenius, W., Calculation of the stability of earth dams. Proc. 2nd Congress on Large
Dams,
1936, 4, 445.
Safety analysis using finite eiements
181
2. Bishop, A. W., The use of the slip circle in the stability analysis of slopes. Geotechnique, 1955, 5(l), 7-17. 3. Janbu, N., Application of composite slip surfaces for stability analysis. Proc. European Conference on Stability of Earth Dams, 1954, 3, 43349. 4. Morgenstern, N. R. and Price, V. E., The analysis of the stability of general slip surfaces. Geotechnique, 1965, 15(l), 79-93. 5. Sarma, S. K., Stability analysis of embankments slopes. Journal of the Geotechnical Engineering Division-ASCE, 1979, 105(GT12), 15 1l-l 524. 6. Davis, E. H. and Booker, J. R., The effect of increasing strength with depth on the bearing capacity of clay. Geotechnique, 1989, 23(4), 551-563. 7. Naylor, D. L, Finite elements and slope stability. In Numerical Methods in Geomechanics, D. Reidel Publishing Company, J. B. Martins (eds), 1982, 229244. 8. Farias, M. M., Numerical analysis of clay core dams. Ph.D thesis, University of Wales, University College of Swansea, 1993. 9. Nordbeck, S. and Rydsteadt, B., Computer cartography point-in-polygon programs. BIT, 1967, 7, 39. 10. Sloan, S. W., A point-in-polygon program. Advances in Engineering Software, 1985, 7, 45-47. 11. Angelov, T. A. and Manoach, E. S., A point-in-domain identification program. Advances in Engineering Software, 1989, 11(2), 99-105. 12. Too, J. M., Two dimensional plate and finite prism isoparametric elements and their applications. Ph.D thesis, University of Wales, University College of Swansea, 197 1. 13. Hinton, E. and Campbell, J. S., Local and global smoothing of discontinuous finite element functions using a least square method. International Journal for Numerical Methods in Engineering, 1974, 8, 461480. 14. Hinton, E., Least squares analysis using finite elements. M.Sc. thesis, University of Wales, University College of Swansea, 1968. 15. Hinton, E. and Irons, B. M., Least squares smoothing of experimental data using finite elements. Journal of the British Society for Strain Measurement, 1968, 14,24-27. 16. Oden, J. T., Finite Elements of Non-Linear Continua. McGraw-Hill, New York, 1971. 17. Oden, J. T. and Brauchli, J., On the calculation of consistent stress distributions in finite element applications. International Journal for Numerical Methods in Engineering, 197 1, 3, 3 17-325. 18. Zienkiewicz, 0. C. and Taylor, R. L., The Finite Element Method, 4th edn. McGraw-Hill, New York, 1990. 19. Zienkiewiez, 0. C. and Zhu, J. Z., A simple error estimator and adaptive procedure for practical engineering analysis. International Journal for Numerical Methods in Engineering, 1987, 24, 337-357. 20. Josselin de Jong, G., Application of the calculus of variations to the vertical cut off in cohesive frictionless soil. Geotechnique, 1980, 30(l), l-16. 21. Pastor, J., Analyse limite: determination de solutions statiques completes. Revice de Mechanique Appliquee, 1978, 2(2), 167-197.