Practical Slope Stability Analysis Using Finite Difference Codes

Practical Slope Stability Analysis Using Finite Difference Codes Loren Lorig, Itasca S.A. Pedro Varona, Itasca Consultores S.L. _.1 INTRODUCTION Numerical models are computer programs that attempt to represent the mechanical re- sponse of a rock mass subjected to a set of initial conditions (e.g., in-situ stresses, water levels), boundary conditions and induced changes (e.g., slope excavation). The result of a numerical model simulation is typically either equilibrium or collapse. If an equilibrium result is obtained, the resultant stresses and displacements at any point in the rock mass can be compared with measured values. If a collapse result is obtained, the predicted mode of failure is demonstrated. In either case, the factor of safety can be calculated. Numerical models divide the rock mass into elements. Each element is assigned a material model and properties. The material models are idealized stress/strain relations that describe how the material behaves. The simplest model is a linear elastic model, which uses the elastic properties (Young’s modulus and Poisson’s ratio) of the material. Elastic-plastic models use strength parameters to limit the shear stress that an element may sustain. The elements may be connected together (a continuum model) or separated by discontinuities (a discontinuum model). Discontinuum models allow slip and separation at explicitly located surfaces within the model. Numerical models tend to be general purpose in nature — that is, they are capable of solving a wide variety of problems. While it is often desirable to have a general-purpose tool 1 available, it requires that each problem be constructed individually. The elements must be arranged by the user to fit the limits of the geomechanical units and/or the slope geometry. Hence, numerical models often require more time to set up and run than special purpose tools (such as limit equilibrium methods.) There are several reasons why numerical models are used for slope stability studies. ƒ Empirical methods cannot confidently be extrapolated outside their databases. ƒ Other methods (e.g., analytic, physical, limit equilibrium) are not available or tend to oversimplify the problem, possibly leading to overly conservative solutions. ƒ Key geologic features, groundwater, etc. can be incorporated into numerical models, providing more realistic approximations of behavior of real slopes. ƒ Observed physical behavior can be explained. ƒ Multiple possibilities (e.g., hypotheses, design options) can be evaluated. _.2 EXPLICIT FINITE DIFFERENCE COMPUTER PROGRAMS (FLAC/FLAC3D, UDEC AND 3DEC) The specific numerical models discussed in this paper are called explicit finite-difference codes, or computer programs. Finite element programs are probably more familiar, but the finite difference method is perhaps the oldest numerical technique used to solve sets of differential equations. Both finite element and finite difference methods produce a set of algebraic equations to solve. While the methods used to derive the equations are different, the resulting equations are the same. The programs referred to in this paper use “explicit” time-marching schemes to solve the equations, whereas finite element methods usually solve systems of equations in matrix form. 2 Even though we usually are interested in a static solution to a problem, the dynamic equations of motion are included in the formulation of finite difference programs. One reason for doing this is to ensure that the numerical scheme is stable when the physical system being modeled is unstable. With non-linear materials, there is always the possibility of physical instability — for example, the failure of a slope. In real life, some of the strain energy is converted to kinetic energy. Explicit finite-difference programs model this process directly, because inertial terms are included. In contrast, programs that do not include inertial terms must use some numerical procedure to treat physical instabilities. Even if the procedure is successful at preventing numerical instability, the path taken may not be realistic. The consequence of including the full law of motion in finite difference programs is that the user must have some physical feel for what is going on. Explicit finite-difference programs are not black boxes that will “give the solution.” The behavior of the numerical system must be interpreted. FLAC (Fast Lagrangian Analysis of Continua) and UDEC (Universal Distinct Element Code) are two-dimensional finite-difference programs developed by Itasca Consulting Group (1998, 2000) specifically for geomechanical analysis. All codes discussed here can simulate varying loading and water conditions and have several pre-defined material models. They are unique in their ability to handle highly non-linear and unstable problems. FLAC is formulated to study continuum problems, although a limited number of discontinuities (in the form of interfaces) can be included. UDEC is formulated to study discontinuum problems involving large numbers (hundreds) of explicit discontinuities that divide the rock mass into blocks. The threedimensional equivalents of these codes are FLAC3D and 3DEC (Itasca 1997, 1999). 3 _.3 MODELING METHODOLOGY The primary objective of numerical modeling at any mine site is to understand the factors controlling present slope behavior and to make predictions about future slope behavior. There is usually concern about the ability to provide stable slopes at reasonable angles for increased mining depths. Numerical modeling attempts to address such concerns through simulation of the mining and mechanical response of the slopes. However, some uncertainty in numerical modeling arises from the generally poor knowledge of the material behavior of rock masses. Rock is a complex material that behaves in a non-linear and anisotropic manner. For problems such as the design of steel structures, where the material behavior is well known, it is usually assumed that a model can predict the behavior of the system accurately, without calibration. Owing to the complexity of a rock mass and the inherent uncertainty of its behavior, it is very important that any model to be used for design first be subjected to a detailed evaluation in which model predictions are compared with the observed and measured response of the slopes. Observed response includes appearance of tension cracks, toppling behavior, etc. Measured response includes instrumentation results from prisms, tape extensometers, inclinometers, etc. This process is known as calibration, or validation, of the modeling approach and is necessary before the model can be used with confidence in predictive studies. One problem associated with the calibration process is non-uniqueness — that is, more than one set of material parameters, for example, will satisfy the calibration conditions. A detailed discussion of modeling methodology is beyond the scope of this chapter. The interested reader is referred to a paper by Starfield and Cundall (1988), which includes three case studies that were examined using finite difference codes. 4 _.4 CREATING MODELS Modeling requires that the real problem be idealized (i.e., simplified) in order to fit the constraints imposed by available material models, computer capacity, etc. Analysis of rock mass response involves different scales. It is impossible — and undesirable — to include all features and details of rock-mass response mechanisms into one model. In addition, many of the details of rock mass behavior are unknown and unknowable; the approach to modeling is not as straightforward as it is, say, in other branches of mechanics. _.4.1 Two-Dimensional Analysis versus Three-Dimensional Analysis The first step in creating a model is to decide whether to perform two-dimensional or three-dimensional analyses. Until recently, three-dimensional analyses were relatively uncommon, but advances in personal computers have permitted three-dimensional analyses to be performed routinely. Strictly speaking, three-dimensional analyses are recommended/required if: (1) the direction of principal geologic structures does not strike within 20o to 30o of the strike of the slope; (2) the axis of material anisotropy does not strike within 20o to 30o of the slope; (3) the directions of principal stresses are thought to be not parallel or not perpendicular to the slope; (4) the distribution of geomechanical units varies along the strike of the slope; or (5) the slope geometry in plan cannot be represented by two-dimensional (i.e., axisymmetric or plain strain) analysis. 5 _.4.2 Continuum versus Discontinuum Model The next step is to decide whether to use a continuum code or a discontinuum code. This decision is seldom straightforward. There appear to be no ready-made rules for determining which type of analysis to perform. All slope stability problems involve discontinuities at one scale or another. However, useful analyses, particularly of global stability, have been made by assuming the rock mass can be represented as an equivalent continuum. Therefore, many analyses begin with continuum models. If the slope under consideration is unstable without structure, there is no point in continuing to discontinuum models. If, on the other hand, a continuum model appears to be reasonably stable, explicit incorporation of principal structures should give a more accurate estimation of slope behavior. Selection of joint geometry for input to a model is a crucial step in discontinuum analyses. Typically, only a very small percentage of joints can actually be included in a model in order to create models of reasonable size for practical analysis. Thus, the joint geometry data must be filtered to select only those joints that are most critical to the mechanical response. This is done by identifying those that are most susceptible to slip and/or separation for the prescribed loading condition. This may involve determining whether sufficient kinematic freedom is provided (especially in the case of toppling) and comparing observed behavior to model response (i.e., calibration). _.4.3 Selecting Appropriate Element Size The next step in the process is to select an appropriate element size. The finite difference zones in the programs listed previously assume that the stresses and strains within each element do not differ with position within the element — in other words, the elements are the lowest- 6 order elements possible. In order to accurately capture stress and strain gradients within the slope, it is necessary to use relatively fine discretizations. By experience, the authors have found that at least 20 (and preferably 30) elements are required over the slope height of interest. _.4.4 Initial Conditions Initial conditions are those conditions that existed prior to mining. The initial conditions of importance at mine sites are the in-situ stress field and the groundwater conditions. Since there are usually no measurements available regarding the initial or current stress field, different stress fields can be tried to determine their effect on model results. The effect of initial stresses is discussed separately in the last section of this chapter. _4.5 Boundary Conditions Boundaries are either real or artificial. Real boundaries in slope stability problems corre- spond to the natural or excavated ground surface and are usually stress-free. Artificial boundaries do not exist in reality. All problems in geomechanics (including slope stability problems) require that the infinite extent of a real problem domain be artificially truncated to include only the immediate area of interest. Figure _.1 shows typical recommendations for locations of the artificial far-field boundaries in slope stability problems. Artificial boundaries can be of two types: prescribed-displacement, or prescribed-stress. Prescribed displacement boundaries inhibit displacement in either the vertical director or horizontal direction, or both. Prescribed displacement boundaries are used to represent the condition at the base of the model and toe of the slope. Displacement at the base of the model is always fixed in both the vertical and horizontal directions to inhibit rotation of the model. Two assumptions can be made regarding the displacement 7 boundaries near the toe of any slope. One assumption is that the displacements near the toe are inhibited only in the horizontal direction. This is the mechanically correct condition for a problem that is perfectly symmetric with respect to the plane or axis representing the toe boundary. Strictly speaking, this condition only occurs in slopes of infinite length, which are modeled in two-dimensions assuming plane strain, or in slopes, which are axially symmetric (i.e., the pit is a perfect cone). In reality, these conditions are rarely satisfied. Therefore some models are extended laterally to avoid the need to specify any boundary condition at the toe of the slope. It is important to note that the difficulties with the boundary condition near the slope toe are usually a result of the two-dimensional assumption. In three-dimensional models, this difficulty usually does not exist. Insert Figure _-1 here. The far-field boundary location and condition must also be specified in any numerical model for slope stability analyses. The general notion is to select the far-field location so that it does not significantly influence the results. If this criterion is met, whether the stress is prescribed-displacement or prescribed-stress is not important. In most slope stability studies, a prescribed-displacement boundary is used. The authors have used a prescribed-stress boundary in a few cases and found no significant differences with respect to the results from a prescribeddisplacement boundary. The magnitude of the horizontal stress for the prescribed-stress boundary must match the assumptions regarding the initial stress state in order for the model to be in equilibrium. However, following any change in the model (i.e., an excavation increment), the prescribed-stress boundary causes the far-field boundary to displace toward the excavation while 8 maintaining its original stress value. For this reason, a prescribed-stress boundary is also referred to as a “following” stress, or constant stress, boundary, because the stress does not change and follows the displacement of the boundary. However, following stresses are most likely where slopes are cut into areas where the topography rises behind the slope. Even where slopes are excavated into an inclined topography, the stresses would flow around the excavation to some extent, depending on the effective width of the excavation perpendicular to the downhill topographic direction. In summary, ƒ A fixed boundary causes both stresses and displacements to be underestimated, whereas a stress boundary does the opposite. ƒ The two types of boundary condition “bracket” the true solution, so that it is possible to do tests with smaller models to get a reasonable estimate of the true solution by averaging the two results. A final point to be kept in mind is that all open-pit slope stability problems are threedimensional in reality. This means that the stresses acting in and around the pit are free to flow both beneath and around the sides of the pit. It is likely, therefore, that, unless there are very low strength faults parallel to the analysis plane, a constant stress or following stress boundary will overpredict the stresses acting horizontally. _.4.6 Rock Mass Material Models As noted previously, it is impossible to model all discontinuities in a large slope. This may be possible for a limited number of benches, but not large slopes. Therefore, much of the rock mass must be represented by an equivalent continuum in which the effect of the discontinui- 9 ties is to reduce the intact rock elastic properties and strength to those of the rock mass. This is true whether or not a discontinuum model is used. Thus, the discontinuities are not explicitly modeled; rather, they are assumed to be smeared throughout the rock mass. The process for initially estimating the rock mass properties is usually based on empirical relations as described, for example, by Hoek and Brown (1997). These initial properties are then modified, as necessary, through the calibration process. The rock mass material model used in most finite difference analyses is a linear elastic– perfectly plastic model. The shear strength is limited by the Mohr-Coulomb criterion. The tensile strength is limited by the specified tensile strength (taken to be 10% of the rock mass cohesion in most analyses). Using this model, the rock mass behaves in an isotropic manner. Strength anisotropy can be introduced through a ubiquitous joint model, which limits the shear strength according to a Mohr-Coulomb criterion in a specified direction. The direction often corresponds to a predominant jointing orientation. A more complete equivalent-continuum model that includes the effects of joint orientation and spacing is a micropolar (Cosserat) plasticity model. This model, as implemented in FLAC, is described in the context of slope stability by Dawson and Cundall (1996). The approach has the advantage of using a continuum model while still preserving the ability to explicitly consider realistic joint spacing. Unfortunately, the model has not yet been incorporated into any publicly available code. Real rock masses often appear to exhibit progressive failure — that is, the failure appears to progress over a period of time. Progressive failure is a complex process that is poorly understood and difficult to model. It may involve one or more of the following component mechanisms: 10 • gradual accumulation of strain or principal structures and/or within the rock mass; • increases in pore pressure with time; and • time-dependent deformation of material under constant load (i.e., creep). Each of these components is discussed briefly below in the context of slope behavior. Gradual accumulation of strain or principal structures within the rock mass usually results from excavation, and “time” is really related to the excavation sequence. In order to study the progressive failure effects due to excavation, one must either introduce characteristics of the post-peak or post–failure behavior of the rock mass into a strain–softening model or introduce similar characteristics into the explicit discontinuities. In practice, there are at least two difficulties associated with strain-softening rock mass models. The first is estimating the post-peak strength and the strain over which the strength reduces. There appear to be no empirical guidelines for estimating the required parameters. This means that the properties must be estimated through calibration. The second difficulty is that, for a simulation in which the response depends on shear localization and in which material softening is used, the results will depend on the element sizes. However, it is quite straightforward to compensate for this form of meshdependence. In order to do this, consider a displacement applied to the boundary of a body. If the strain localizes inside the body, the applied displacement appears as a jump across the localized band. The thickness of the band contracts until it is equal to the minimum allowed by the grid — i.e., a fixed number of element widths. Thus, the strain in the band is ε = u / n∆z where n is a fixed number, u is the displacement jump, and ∆z is the element width. 11 If the softening slope is linear (i.e., the change in a property value, ∆p , is proportional to strain), the change in property value with displacement is ∆p / ∆u = s / n ∆z where s is the input softening slope. In order to obtain mesh-independent results, we can input a scaled softening slope, such that: s = s′∆z where s′ is constant. In this case, ∆p / ∆u is independent of ∆z . If we define the softening slope by the critical s , then strain, ε crit s ∝ 1 / ∆z . ε crit For example, if we double the zone size, we must halve the critical strain for comparable results. Strain-softening models for discontinuities are much more common than similar relations for rock masses. Strain-softening relations for discontinuities are built into UDEC and 3DEC and can be incorporated into interfaces in FLAC and FLAC3D via a built-in programming language (i.e., FISH functions). Strain-softening models require special attention when computing safety factors. A technique to account for strain-softening behavior in safety factor calculations is discussed in the section on safety factors. Increases in pore pressure with time are not common in rock slopes for mines. More commonly, the pore pressures reduce due to deepening of the pit and/or drainage. However, there are cases in which the pore pressures do increase with time. In such cases, the slope may appear to progressively fail. 12 Time-dependent deformation of material under constant load (i.e., creep) is not commonly considered in the context of slope stability. It is much more common in underground excavations. Several material models are available to study creep behavior in rock slopes. These include classical viscoelastic models, power-law models, and the Burger-creep viscoplastic model. Work is under way to apply creep models to the study of slope behavior at Chuquicamata mine in Chile (Calderon, 2000). Of the models listed above, the Mohr-Coulomb plasticity model is used in approximately 90% of all slope analyses; the ubiquitous joint model is used approximately 10% of the time. Strain-softening models have been used, but rarely. _.4.7 Joint Material Models in UDEC and 3DEC The material models used most commonly to explicitly represent joints in UDEC and 3DEC is a linear elastic–plastic model. A peak and residual shear strength can be specified for the joints. The residual strength is used after the joint has failed in shear at the peak strength. The elastic behavior of the joints is specified by joint normal and shear stiffnesses. _.5 EFFECT OF FLOW ANALYSIS VERSUS WATER TABLE ON PREDICTED SLOPE STABILITY The effect of water pressure in reducing effective stresses and, hence, slope stability is well understood. However, the effect of various assumptions regarding specification of pore pressure distributions in slopes is not as well understood. Two methods are commonly used to specify pore pressure distributions within slopes. The most rigorous method is to perform a complete flow analysis and use the resultant pore pressures in the stability analyses. A less rig- 13 orous, but more common, method is to simply specify a water table. The resulting pore pressures are then assumed to be given by the product of the vertical depth below the water table, the water density and gravity. In this sense, the water table approach is equivalent to specifying a piezometric surface. Both methods use similar phreatic surfaces, but the second method underpredicts pore pressure concentrations that actually occur near the toe of a slope and slightly overpredicts the pore pressure behind the toe by ignoring the inclination of equipotential lines. Finally, seepage must be considered. The difference in water pressure that exists between two points at the same elevation (i.e., hydraulic gradient) results from seepage forces (or drag) as water moves through a porous medium. Flow analysis automatically accounts for seepage forces. To evaluate the error resulting from specifying a water table (without doing a flow analysis), two identical problems were run. In one case, a flow analysis was performed to determine the pore pressures. In the second case, the pressures were determined using only a piezometric surface that was assumed to be the phreatic surface taken from the flow analysis. The material properties and geometry for both cases are shown in Figure _.2. The right-hand boundary was extended to allow the far-field phreatic surface to coincide with the ground surface at a horizontal distance of two kilometers behind the toe. Permeability within the model was assumed to be homogeneous and isotropic. The error caused by specifying the water table can be seen in Figure _.3. The largest errors (up to 45%) are found just below the toe. Errors in pore pressure values behind the slope are generally less than 5%. The errors near the phreatic surface are not significant, as they result from the relatively small pore pressures just below the phreatic surface (i.e., small errors in small values result in large relative errors). 14 Insert Figure _.2 here. Insert Figure _.3 here. For a phreatic surface at the ground surface at a distance of 2 kilometers, a factor of safety of 1.1 is predicted using circular failure chart number 3 (Hoek and Bray 1981). The factor of safety determined by FLAC is approximately 1.15 for both cases. The FLAC analyses give similar safety factors because the distribution of pore pressures in the area where failure occurs (i.e., behind the slope) is very similar for the two cases. The conclusion drawn here is that there is no significant difference in predicted stability between a complete flow analysis and simply specifying a piezometric surface. However, it is not clear if this conclusion can be extrapolated to other cases involving, for example, anisotropic flow. _.6 EXCAVATION SEQUENCE AND INTERPRETATION OF RESULTS _.6.1 Excavation Sequence Simulating excavations in finite difference models poses no conceptual difficulties. However, the amount of effort required to construct a model depends directly on the number of excavation stages simulated. Therefore, most practical analyses seek to reduce the number of excavation stages. The most accurate solution is obtained using the largest number of excavation steps, since the real load path for any element in the slope will be followed closely. In theory, it is impossible to prove that the final solution is independent of the load-path followed. However, for many slopes, the stability seems to depend mostly on the condition of the slope (i.e., geome- 15 try, pore pressure distribution) at the time of analysis, and very little on the load path taken to get there. A reasonable approach regarding the number of excavation stages has evolved over the years. Using this approach, only a few (one, two or three) excavation stages are modeled. For each stage, two calculation steps are taken. In the first step, the model is run elastically to remove any inertial effects caused by sudden removal of a large amount of material. Next, the model is run allowing plastic behavior to develop. Following this approach, reasonable solutions to a large number of slope stability problems have been obtained. _.6.2 Interpretation of Results As noted in the introduction, finite difference programs are not black boxes that “give the solution.” The behavior of the numerical system must be interpreted, and results from finite difference models are interpreted in much the same way as prism data are interpreted. Finite difference programs record displacements and velocities at nominated points within the rock mass. During the analysis, the recorded values can be examined to see if they are increasing, remaining steady, or decreasing. Increasing displacements and velocities indicate an unstable situation; steady displacements and decreasing velocities indicate a stable situation. In addition, velocity and displacement vectors for every point in the model can be plotted. Fields of constant velocity and displacement indicate failure. The authors have found that velocities below 1e-6 indicate stability in FLAC and FLAC3D; conversely, velocities above 1e-5 indicate instability. Note that no units are given for velocities. This is because the velocities are not real, due to the damping and mass scaling used 16 to achieve static solutions. While the velocities are not real, the displacements are, although we can not say at what “time” the displacement occurs. The failure (plasticity) state (i.e., elastic, failed in tension, failed in shear) of points within the model can also be examined. Care must be used in examining the failure state indicators. For example, local overstressing (e.g., at the base of toppling columns) can appear to form a deep-seated slip surface when, in reality, it is just compressive failure of toppling columns. Therefore, the failure (plasticity) indicators must be reviewed in the context of overall behavior before any definitive conclusions can be drawn. _.6.3 Calculating Factor of Safety Evaluation of Safety Factor Using the Shear Strength Reduction Technique. For slopes, the factor of safety, F, is often defined as the ratio of the actual shear strength to the minimum shear strength required to prevent failure. A logical way to compute the factor of safety with a finite element or finite difference program is to reduce the shear strength until collapse occurs. The factor of safety is the ratio of the soil or rock’s actual strength to the reduced shear strength at failure. This shear-strength reduction technique was first used with finite elements by Zienkiewicz et al. (1975) to compute the safety factor of a slope composed of multiple materials. To perform slope stability analysis with the shear-strength reduction technique, simulations are run for a series of increasing trial factors of safety, f. Actual shear strength properties cohesion, c, and friction, φ, are reduced for each trial according to the equations: ctrial = (1/f)c (1) φtrial = arctan{(1/f)tanφ} (2) 17 If multiple materials and/or joints are present, the reduction is made simultaneously for all materials. The trial factor of safety is gradually increased until the slope fails. At failure, the safety factor equals the trial safety factor (i.e., f = F). Dawson et al. (1999) show that the shearstrength reduction factors-of-safety are generally within a few percent of limit equilibrium solutions when an associated flow rule is used. If a strain softening constitutive model is used, the softening logic should be turned off during the shear-strength reduction process or the factor of safety will be underestimated. When the slope is excavated, some zones will have exceeded their peak strength and some amount of softening will have taken place. During the strength reduction process, these zones should be considered as a new material with lower strength, but no further softening should be allowed due to the plastic strains associated to the gradual reduction of strength. The shear-strength reduction technique has two advantages over slope stability analyses with limit equilibrium. First, the critical failure surface is found automatically, and it is not necessary to specify the shape of the failure surface (e.g., circular, log spiral, piecewise linear, etc.) in advance. In general, the failure mode for slopes is more complex than simple circles or segmented surfaces. Second, numerical methods automatically satisfy translational and rotational equilibrium; not all limit equilibrium methods do. Consequently, the shear-strength reduction technique will usually determine a lower safety factor compared to other methods. For example, Zienkiewicz et al. (1975) show comparisons between factors of safety calculated by limitequilibrium and finite element methods for a slope composed of three materials. They report a factor of safety of 1.33 using Bishop’s method and 1.165 using the shear-strength reduction technique. Donald and Giam (1988) report results of factor of safety analyses for a homogeneous slope with a weak layer in the foundation. They report a factor of safety of 1.5 using 18 Bishop’s method with a simplex optimization and 1.34 using the shear-strength reduction technique. _.7 SOME USEFUL RESULTS _.7.1 Effect of Radius of Curvature on Slope Stability Most design analysis for slopes assumes a two-dimensional problem geometry (i.e., a unit slice through an infinitely long slope, under plane-strain conditions). In other words, the radii of both the toe and crest are assumed to be infinite. This is not the condition encountered in practice — particularly in open-pit mining, where the radii of curvature can have an important effect on safe slope angles. Concave slopes are inherently more stable than convex slopes due to the lateral restraint provided by material on either side of a potential failure in a concave slope. Despite its potential importance in slope stability, very little has been done to quantify this effect. Jenike and Yen (1961) presented the results of limit theory analysis of axisymmetric slopes in a rigid-perfectly plastic material, and determined S-shaped critical profiles that describe the theoretical failure shape for different radii of curvature. They showed that, as the radius of the slope increases, the profile of the stable slope in axial symmetry approaches the profile of the slope in plane strain. However, their method does not permit any evaluation of safety factors. In addition, their assumptions resulted in slopes that, in plane strain, took no account of cohesion. Hoek and Brown (1981) concluded that the analysis assumptions were not applicable to rock slope design. Piteau and Jennings (1970) studied the influence of plan curvature on the stability of slopes in four diamond mines in South Africa. As a result of caving from below the surface, slopes were all at incipient failure (i.e., a safety factor of one). The average slope height was 19 100 m. Piteau and Jennings found that the average slope angle for slopes with radius of curvature of 60 m was 39.5° as compared to 27.3° for slopes with a radius of curvature of 300 m. Hoek and Bray (1981) summarize their experience with the stabilizing effects of slope curvature as follows. When the radius of curvature of a concave slope is less than the height of the slope, the slope angle can be 10° steeper than the angle suggested by conventional stability analysis. As the radius of curvature increases to a value greater than the slope height, the correction should be decreased. For radii of curvature in excess of twice the slope height, the slope angle given by a conventional stability analysis should be used. To better quantify the effects of slope curvature, a series of generic analyses were performed. All analyses assumed a 500-m high dry slope consisting of an isotropic homogeneous material with 35o friction, 0.66 MPa cohesion and 2,600-kg/m3 density (i.e., similar to Figure _.2). Initial in-situ stresses are assumed to be lithostatic, and the excavation is made in 40-m decrements beginning from the ground surface. For these conditions, a factor of safety of 1.3 is predicted using circular failure chart number 1 (Hoek and Bray 1981). A factor of safety of 1.3 is a value that is frequently used in the design of slopes for open-pit mines. Two series of analyses were performed using FLAC. In the first series, the safety factor was calculated for axisymmetric conditions with various radii of curvature. The results are shown in Table _.1. Insert Table _.1 here. In the second series, the slope angle was increased until a safety factor of 1.3 was achieved. The results of this series of analyses are shown in Table _.2. 20 Insert Table _.2 here. The results shown for R = 500 m in Table _.2 confirm the notion that, when the radius of curvature of a concave slope is less than the height of the slope, the slope angle can be 10° steeper than the angle suggested by conventional (i.e., two-dimensional) stability analysis. One reason that designers are reluctant to take advantage of the beneficial effects of slope curvature is that the presence of discontinuities can often negate the effects. However, for massive rock slopes, or slopes with relatively short joint trace lengths, the beneficial effects of slope curvature should not be ignored — particularly in open-pit mines, where the economic benefits of steepening slopes can be significant. _.8 CASE STUDY: BOINAS EAST PIT Boinas East is one of the three open pits currently operated by Narcea Gold Mines in northwestern Spain. Clifford (2000) gives a general description of the mining operation. The east slope of the pit consists of a strong dolomite overlying a weak sandstone, weathered granite and the mineralized skarn, as shown in Figure _.4. The slope height is 280 m, and its radius of curvature at the base is 70 m. The slope angle is 65º in the dolomites and 45º in the sandstone and the granite. The pore pressure is negligible in the pit. The properties used for the analysis appear in Table _.3. Insert Figure _.4 here. Insert Table _.3 here. 21 A plane-strain FLAC analysis predicts a factor of safety of 1.05 with shear failure in the sandstone, and weathered granite and tensile failure in the dolomite (Figure _.5). In an axisymmetric analysis, the low radius of curvature relative to the height inhibits the tensile failure in the limestone, and the factor of safety increases to 1.70 (Figure _.6). Insert Figure _.5 here. Insert Figure _.6 here. However, the results of the axisymmetric analysis are too optimistic. While the geometry of the pit can be assumed to be conical, the sedimentary units dip east and do not follow an axisymmetric pattern. A more elaborate FLAC3D model (Figure _.7) shows a similar failure pattern (Figure _.8) in the east slope, but with a factor of safety of 1.35. Insert Figure _.7 here. Insert Figure _.8 here. _.9 IMPACT OF INITIAL STRESS STATE The role of stresses has been traditionally ignored in slope analyses. There are several possible reasons for this. ƒ Limit equilibrium analyses, which are widely used for stability analyses, can not include the effect of stresses in their analyses. Nevertheless, limit equilibrium analyses 22 are thought to provide reasonable estimates of stability in many cases, particularly where structure is absent (e.g., slopes in soil). ƒ Most stability analyses have traditionally been performed for soils, where the range of possible in-situ stresses is more limited than for rocks. Furthermore, many soil analyses have been performed for constructed embankments (e.g., dams), where in-situ stresses do not exist. ƒ Most slope failures are gravity-driven, and the effects of in-situ stress are thought to be minimal. ƒ In-situ stresses in rock masses are not routinely measured for slopes, and their effects are largely unknown. One particular advantage of stress analysis programs (i.e., numerical models) is their ability to include pre-mining initial stress states in stability analyses and to evaluate their importance. Five cases were run using a model similar to the model shown in Figure _.2 in order Tao evaluate the effects of in-situ stress state on stability. For each of the five cases, the slope angle was 60o degrees, the slope height was 400 m, the material density was 2,450 kg/m3, the friction angle was 32o, and the cohesion was 0.92 MPa. The results of FLAC analyses are shown in Table _4. Insert Table _.4 here. 23 Although, in general, it is impossible to say what effect the initial stress state will have on any particular problem, as behavior depends on several factors (i.e., orientation of major structures, rock mass strength, water conditions), some observations can be made. ƒ The larger the initial horizontal stresses are, the larger the horizontal elastic displacements will be. This is not much help, as displacements that are of elastic magnitude are not particularly important in slope studies. ƒ Initial horizontal stresses in the plane of analysis that are less than the vertical stresses tend to slightly decrease stability and reduce the depth of significant shearing with respect to a hydrostatic stress state. This observation may seem counterintuitive; we would expect smaller horizontal stresses to increase stability. The explanation lies in the fact that the lower horizontal stresses actually provide slightly decreased normal stress on potential shearing surfaces and/or joints within the slope. This observation was confirmed in UDEC analysis of a slope in Peru where in-situ horizontal stresses lower than the vertical stress led to deeper levels of joint shearing in toppling structures compared to cases involving horizontal stresses that were equal to or greater than the vertical stress. ƒ It is important to note that the regional topography may limit the possible stress states, particularly at elevations above regional valley floors. Three-dimensional models have been very useful in the past in addressing some regional stress issues. _.10 REFERENCES Calderon, A. 2000. CODELCO Chile, Division Chuquicamata. Personal communication. Clifford, D. 2000. Rio Narceas’s ongoing growth. Mining Magazine. 820-86. 24 Dawson, E. M., and P. A. Cundall. 1996. Slope stability using micropolar plasticity. Rock Mechanics Tools and Techniques (Proceedings of the 2nd North American Rock Mechanics Symposium, Montréal, June 1996), 1:551-558. M. Aubertin et al., Eds. Rotterdam: A. A. Balkema. Dawson, M. D., W. H. Roth, and A. Drescher. 1999. Slope stability analysis by strength reduction. Submitted for publication to Géotechnique. Donald, I. B., and S. K. Giam. 1988. Application of the nodal displacement method to slope stability analysis. Preprint Proceedings of the Fifth Australia-New Zealand Conference on Geomechanics — Prediction versus Performance (Sydney, 1988). 456-460. Hoek, E., and J. W. Bray. 1981. Rock Slope Engineering, 3rd Ed. London: The Institute of Mining and Metallurgy. Hoek, E., and E. T. Brown. 1997. Practical estimates of rock mass strength. Int. J. Rock Mech. Min.. Sci. & Geomech. Abstr. 34(8):1165-1186. Itasca Consulting Group, Inc. 1998. FLAC (Fast Lagrangian Analysis of Continua), Version 3.4. Minneapolis: ICG. Itasca Consulting Group, Inc. 1999. 3DEC (Three-Dimensional Distinct Element Code), Version 2.0. Minneapolis: ICG. Itasca Consulting Group, Inc. 2000. UDEC (Universal Distinct Element Code), Version 3.1. Minneapolis: ICG. Itasca Consulting Group, Inc. 1997. FLAC3D (Fast Lagrangian Analysis of Continua in 3 Dimensions), Version 2.0. Minneapolis: ICG. Jenike, A. W., and B. C. Yen. 1961. Slope stability in axial symmetry. Proceedings of the 5th Rock Mechanics Symposium, University of Minnesota, 689-711. London: Pergamon. 25 Piteau, D. R., and J. E. Jennings. 1970. The effects of plan geometry on the stability of natural slopes in rock in the Kimberley area of South Africa. Proceedings of the 2nd Congress of the International Society of Rock Mechanics (Belgrade), 3:Paper 7-4. Starfield, A. M., and P. A. Cundall. 1988. Towards a methodology for rock mechanics modelling. Int. J. Rock Mech. Min.. Sci. & Geomech. Abstr. 25(3):99-106. Zienkiewicz, O. C., C. Humpheson, and R. W. Lewis. 1975. Associated and non-associated visco-plasticity and plasticity in soil mechanics. Géotechnique. 25(4):671-689. 26 TABLE _.1 α = 45o Effect of radius of curvature on factor of safety for 45° slope R = 100 m F = 1.75 R = 250 m F = 1.65 R = 500 m F = 1.55 R=∞ F = 1.3 1 TABLE _.2 F = 1.3 Effect of radius of curvature on slope angle for safety factor = 1.3 R = 100 m α = 75o R = 250 m α = 60o R = 500 m α = 55o R=∞ α = 45o 2 TABLE _.3 Rock mass properties used in Boinas East Pit Study MATERIAL Strong dolomite Weak sandstone Weathered granite Marble Ore K(GPa) G(GPa) 7.92 1.31 1.17 2.45 1.11 5.09 0.75 0.66 1.41 0.62 φ (º) 43 21 20 33 19 COH (kPa) 1120 270 250 500 250 DEN(t/m3) 2.70 2.30 2.28 2.43 2.21 3 TABLE _.4 Effect of in-situ stress on slope stability In-plane Horizontal Stress σxx = σyy σxx = 2.0 * σyy σxx = 0.5 * σyy σxx = 2.0 *.σyy σxx = 0.5 * σyy Out-of-plane Horizontal Stress σzz = σyy σzz = 2.0 * σyy σzz = 0.5 * σyy σzz = 0.5 * σyy σzz = 2.0 * σyy Factor of Safety 1.30 1.30 1.28 1.30 1.28 4 FIGURE _.1 Minimum dimensions for slope analysis model JOB TITLE : Problem geometry and properties (*10^3) FLAC (Version 3.40) LEGEND 1.000 9-Mar- 0 9:08 step 0 -1.389E+02 <x< 2.639E+03 -1.264E+03 <y< 1.514E+03 Water table 45° 0.500 Grid plot 0 5E 2 Water Table 0.000 Density = 2600 kg/m3 Friction = 35 degrees Cohesion = 660 kPa -0.500 -1.000 Itasca Consulting Group, Inc. Minneapolis, Minnesota USA 0.250 0.750 1.250 (*10^3) 1.750 2.250 FIGURE_.2 Problem geometry and rock mass properties used to evaluate effect of pore pressure assumptions on slope stability FIGURE_.3 Error caused by specifying water table instead of performing a flow analysis JOB TITLE : Comparison of Pore Pressure Distributions (*10^3) FLAC (Version 3.40) LEGEND 8-Mar- 0 18:59 step 20000 Cons. Time 1.3380E+11 0.000E+00 <x< 2.100E+03 -1.000E+03 <y< 1.100E+03 0.800 Pore pressures underpredicted by water table in this area 0.400 Pore Pressure Error -4.50E-01 -3.50E-01 -2.50E-01 -1.50E-01 -5.00E-02 0.000 Contour interval= 5.00E-02 (zero contour omitted) Pore pressures overpredicted by water table in this area -0.400 -0.800 Itasca Consulting Group, Inc. Minneapolis, Minnesota USA 0.200 0.600 1.000 (*10^3) 1.400 1.800 JOB TITLE : (*10^2) FLAC (Version 3.40) 4.500 LEGEND Strong Dolomite 3.500 3-Jan- 0 13:43 step 14810 -3.544E+01 <x< 6.733E+02 -2.182E+02 <y< 4.905E+02 friction 1.900E+01 2.000E+01 2.100E+01 3.300E+01 4.300E+01 Grid plot 0 2.500 Weak Sandstone Weathered Granite 1.500 0.500 2E 2 -0.500 Ore -1.500 Marble Itasca Consulting Group, Inc. Minneapolis, Minnesota USA 0.500 1.500 2.500 3.500 (*10^2) FIGURE_.4 Geometry of Boinas East Pit with axial symmetry 4.500 5.500 6.500 JOB TITLE : (*10^2) FLAC (Version 3.40) 4.500 LEGEND 3.500 3-Jan- 0 16:20 step 131299 -3.544E+01 <x< 6.733E+02 -2.182E+02 <y< 4.905E+02 2.500 Boundary plot 0 2E 2 Plasticity Indicator * at yield in shear or vol. X elastic, at yield in past o at yield in tension Max. shear strain-rate 1.500 0.500 Contour interval= 5.00E-06 Minimum: 0.00E+00 Maximum: 9.50E-05 -0.500 -1.500 Itasca Consulting Group, Inc. Minneapolis, Minnesota USA 0.500 1.500 2.500 3.500 (*10^2) FIGURE _.5 Failure mechanism of Boinas East Pit with plane-strain assumption 4.500 5.500 6.500 JOB TITLE : (*10^2) FLAC (Version 3.40) 4.500 LEGEND 3.500 3-Jan- 0 21:05 step 89907 -3.544E+01 <x< 6.733E+02 -2.182E+02 <y< 4.905E+02 2.500 Boundary plot 0 2E 2 Plasticity Indicator * at yield in shear or vol. X elastic, at yield in past o at yield in tension Max. shear strain-rate 1.500 0.500 Contour interval= 5.00E-07 Minimum: 0.00E+00 Maximum: 2.60E-05 -0.500 -1.500 Itasca Consulting Group, Inc. Minneapolis, Minnesota USA 0.500 1.500 2.500 3.500 (*10^2) FIGURE _.6 Failure mechanism of Boinas East Pit with axial symmetry 4.500 5.500 6.500 FLAC3D 2.00 Step 31966 Model Projection 09:40:39 Fri Mar 10 2000 Center: X: 1.169e+003 Y: 1.026e+004 Z: 4.630e+002 Dist: 3.906e+003 Rotation: X: 39.196 Y: 359.485 Z: 23.312 Size: 1.300e+003 Block Group OD H T GR MIN SK DOL AGR USS LSS Alfredo Varela Rio Narcea Gold Mines SA FIGURE _.7 Boinas East Pit FLAC3D model FLAC3D 2.00 Step 31966 Model Projection 10:05:21 Fri Mar 10 2000 Center: X: 1.210e+003 Y: 1.026e+004 Z: 3.733e+002 Dist: 3.906e+003 Rotation: X: 0.000 Y: 0.000 Z: 0.000 Size: 1.146e+003 Plane Origin: X: 6.610e+002 Y: 1.028e+004 Z: 3.040e-013 Plane Orientation: Dip: 90.000 DD: 0.000 Contour of Velocity Mag. Plane: on 0.0000e+000 to 2.0000e-006 4.0000e-006 to 6.0000e-006 8.0000e-006 to 1.0000e-005 1.2000e-005 to 1.4000e-005 1.6000e-005 to 1.8000e-005 2.0000e-005 to 2.2000e-005 2.4000e-005 to 2.6000e-005 2.8000e-005 to 3.0000e-005 3.2000e-005 to 3.4000e-005 3.6000e-005 to 3.8000e-005 4.0000e-005 to 4.2000e-005 4.4000e-005 to 4.6000e-005 4.8000e-005 to 5.0000e-005 5.2000e-005 to 5.4000e-005 5.6000e-005 to 5.8000e-005 5.8000e-005 to 6.0000e-005 Interval = 2.0e-006 Alfredo Varela Rio Narcea Gold Mines SA FIGURE _.8 Cross-section through Boinas East Pit FLAC3D model showing failure mode 4