TOTAL STRESS DYNAMIC ANALYSIS OF COQUITLAM DAM by CHARISSA W. DHARMASETIA B.A.Sc, The University of British Columbia, 1995 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in THE FACULTY OF GRADUATE STUDIES Department of Civil Engineering We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA July, 2000 © Charissa Dharmasetia, 2000 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of <?*<f-'A/££Z/Ai& The University of British Columbia Vancouver, Canada Date 3/ AOfST, JooO DE-6 (2/88) ABSTRACT Earthquake-induced liquefaction of loose saturated sands can cause large deformations to occur resulting in flow slides. The near catastrophic failure of the Lower San Fernando Dam due to the 1971 earthquake is one of the most well known examples of a flow slide. However, the concern is not limited to flow slide situations in which the driving stresses are greater than the residual strength after liquefaction. Significant deformations can occur for a non-flow slide condition. These are caused by a reduction in stiffness of liquefied materials as well as from inertia forces caused by the earthquake motion. Limit equilibrium type analyses have been used reliably in the prediction of the occurrence of a flow slide. However, difficulty has been encountered in the prediction of earthquake-induced deformation in the event that a flow slide is not predicted to occur. The problems have arisen from the modeling of the behaviour of liquefied soils. Soil exhibits stiff behaviour under cyclic loading and the strains associated up to the point of triggering of liquefaction are small. Upon liquefaction triggering, soil behaves like a liquid and initially strains under very small shear stresses. However, upon further deformation, the soil dilates and regains stiffness and strength. Most deformation analysis procedures overly simplify the effects of earthquake inertia forces and do not adequately model the stress-strain behaviour of liquefied soil. The proposed total stress procedure attempts to take account of both of the above effects. The procedure is separated into three main phases of cyclic-induced liquefaction behaviour of sands: the pre-triggering, triggering, and post-triggering response of soils. The triggering of liquefaction in each element of a soil structure is predicted by weighting the cyclic shear stresses induced by a prescribed base motion. Upon triggering of liquefaction, liquefied stress-strain parameters are assigned to zones predicted to liquefy as they occur. The pre-triggering and triggering phases of the procedure were verified using S H A K E analyses. Similarly, the post-triggering phase of the procedure was compared with results obtained from Bartlett and Youd's empirical equation. In both cases, reasonable agreement was found. The method was finally applied to the Coquitlam Dam and the results compared with two of the more commonly used deformation analyses such as variations of the Modified Modulus method and Jitno and Byrne's extended Newmark method. The predicted results from all three methods are in reasonable agreement. TABLE OF CONTENTS ABSTRACT » TABLE OF CONTENTS iv LIST OF TABLES ix LIST OF FIGURES x LIST OF SYMBOLS xiv ACKNOWLEDGEMENTS xvii CHAPTER 1 1 INTRODUCTION 1 CHAPTER 2 4 BEHAVIOUR OF SATURATED UNDRAINED SOILS 4 2.1 Introduction 4 2.2 Monotonic Loading Behaviour 4 2.3 Cyclic Loading Behaviour 7 2.3.1 Cyclic Resistance 12 2.4 Post-Liquefaction Cyclic Loading Behaviour 15 2.5 Post-Liquefaction Monotonic Loading Behaviour 19 2.5.1 Residual Strength 24 2.5.2 Limiting Shear Strain 26 2.6 Post-Liquefaction Volumetric Deformation 27 2.7 Summary 30 CHAPTER 3 31 iv EVALUATION OF EARTHQUAKE INDUCED DEFORMATIONS OF EARTH STRUCTURES ..31 3.1 Introduction 31 3.2 Empirical Methods 32 3.2.1 Hamadaetal. 32 3.2.2 Bartlett and Youd 33 3.3 Pseudo Static Methods 34 3.3.1 Newmark's Method 35 3.3.2 Byrne's Extended Newmark Approach 38 3.4 Two Dimensional Methods 41 3.4.1 Modified Modulus Approach 42 3.4.2 Strain Potential / Dynamic Stress Path Approach 42 3.4.3 Jitno and Byrne's Approach 43 3.5 Non-Linear Effective Stress Analysis 44 3.6 Summary 45 CHAPTER 4 46 PROPOSED TOTAL STRESS METHOD 46 4.1 Introduction 46 4.2 Triggering of Liquefaction 47 4.2.1 Cyclic Stress Ratio 47 4.2.2 Cyclic Resistance Ratio 49 4.3 Proposed Method 49 4.3.7 Pre-Liquefaction Trigger 51 4.3.2 Liquefaction Trigger 52 4.3.3 Post-Liquefaction Trigger 54 4.3.4 Mohr-Coulomb Model 59 4.3.4.1 Introduction 59 4.3.4.2 Yield Sudace 60 v 4.3.4.3 Non-associative Flow Rule 60 4.4 Required Parameters 61 4.4.1 Pre-liquefaction Parameters 61 4.4.2 Post-Liquefaction Parameters 62 4.4.2.1 Residual Strength 62 4.4.2.2 Limiting Shear Strain 62 4.5 Summary 63 CHAPTER 5 64 VERIFICATION OF TOTAL STRESS PROCEDURE 64 5.1 Introduction 64 5.2 Triggering Verification Using SHAKE 65 5.2.1 Description of FLAC 65 5.2.2 Case Analyzed. 68 5.2.3 Results 69 5.3 Post-Triggering Verification Using Bartlett and Youd Method 70 5.3.1 Case Analyzed. 71 5.3.2 Results 71 5.3.2.1 Effect of G/Gmax ratio 74 5.3.2.2 Effect of Critical Damping Ratio 74 5.3.2.3 Effect of water table 75 5.3.2.4 Effect of Earthquake Time Histories and Magnitude 75 5.4 Summary 82 CHAPTER 6 83 PAST ANALYSES OF COQUITLAM DAM 83 6.1 Introduction 83 6.2 Field and Laboratory Testing 85 vi 6.2.1 Field Exploration 85 6.2.1.1 Mud Rotary Drilling 86 6.2.1.2 In-Situ Testing 88 6.2.1.2.1 Dynamic Cone Penetration Testing 88 6.2.1.2.2 Pressuremeter Testing 88 6.2.1.2.3 Electric Cone Penetration Testing 89 6.2.2 Laboratory Testing 91 6.2.2.1 Index Tests 91 6.2.2.2 Resonant Column Tests 92 6.2.2.3 Consolidated Undrained Triaxial Test 92 6.2.2.4 Cyclic Triaxial Tests 94 6.2.2.5 Post Cyclic Monotonic Tests 94 6.3 Seismic Stability Analyses 95 6.3.1 1979 Analyses 95 6.3.2 1984 Analyses 98 6.4 Other Studies 101 6.5 Summary 101 CHAPTER 7 102 TOTAL STRESS DYNAMIC ANALYSES OF COQUITLAM DAM 102 7.1 Introduction 102 7.2 Pre-Earthquake Analysis 102 7.2.1 Construction Analysis 103 7.2.1.1 Model Parameters 103 7.2.1.2 Results 106 7.2.2 End of Reservoir Filling Analysis 107 7.2.2.1 Flow Model Parameters 109 7.2.2.2 Results 110 vii 7.3 Total Stress Dynamic Analysis 113 7.3.1 Introduction 113 7.3.2 Liquefaction Triggering Analysis 114 7.3.3 Flow Slide Analysis 116 7.3.4 Deformation Analysis 117 7.3.4.1 Model Parameters 117 7.3.4.2 Results 119 7.4 Proposed Total Stress Procedure 121 7.4.1 Model Parameters 121 7.4.2 Results 123 7.5 Deformation Due to Consolidation 127 7.6 Summary 129 CHAPTER 8 131 CONCLUSIONS AND RECOMMENDATIONS 131 BIBLIOGRAPHY 134 APPENDIX 1 140 viii LIST OF TABLES Table 3. 1 Post-Liquefaction Stress-Strain Parameters 41 Table 4. 1 Equivalent Number of Cycles 48 Table 5. 1 Summary of Earthquake Histories used in Bartlett-Youd Verification 75 Table 6. 1 Summary of Mud Rotary Drill Holes 87 Table 6. 2 Summary of Laboratory Tests 91 Table 7. 1. Summary of Dam Material Parameters used in Static Analysis 106 Table 7. 2. Summary of Dam Material Properties for Seepage Analysis 110 Table 7. 3 Corrected Cyclic Resistance Ratio for Coquitlam Dam Materials 114 Table 7. 4 S H A K E input parameters 115 Table 7. 5 Results of Liquefaction Induced Displacement of Coquitlam Dam 119 Table 7. 6 Summary of Parameters Used in Proposed Model 122 ix LIST OF FIGURES Figure 2.1 Characteristic Undrained Monotonic Response of Saturated Soils (Chern, 1985) 5 Figure 2. 2 Characteristic Undrained Response of Saturated Sands Under Cyclic Loading (Vaid and Chern, 1985) 9 Figure 2. 3 Stress Ratio versus Number of Cycle to Liquefaction (De Alba, 1976) 11 Figure 2. 4 Ko versus Effective Overburden Pressure (Pillai and Byrne, 1994) 12 Figure 2. 5 Correction Factor for Static Shear ( N C E E R , 1997) 13 Figure 2. 6 Post-Liquefaction Stress-Strain Curves (Kuerbis, 1989) 17 Figure 2. 7 Surface Acceleration versus Relative Displacement (Byrne and Mclntyre., 1994) 18 Figure 2. 8 Approximate Stress-Strain Relationship during Shaking Table Tests (Sasaki et al., 1992) 18 Figure 2. 9 Post-Liquefaction Monotonic Behaviour of Saturated Sands (Seed, 1979) a. Sacramento River Sand, D R = 40%; b. Mine Tailings, D R = 95% 20 Figure 2. 10 Effect of Excess Pore Pressure on Post-Liquefaction Behaviour of Saturated Sands (Yasuda, 1994) 22 Figure 2.11 Characterization of Post-Liquefaction Curve (Thomas, 1992) 23 Figure 2. 12 Shear strength of strain-softening Syncrude Sand (Vaid et al. 1998) 25 Figure 2. 13 Relationship Between Corrected "Clean Sand" Blowcount (N^o-os and Undrained Strength (Sr) (Seed and Harder, 1990) 26 Figure 2. 14 Volumetric Strain versus Maximum Amplitude Shear Strain for Different Relative Densities (Ishihara and Yoshimine, 1992) 28 Figure 2 .15 Volumetric Strain versus Factor of Safety Against Liquefaction (Ishihara and Yoshimine, 1992) 29 Figure 3. 1 Idealized Potential Sliding Slope (Newmark, 1965) 35 Figure 3. 2 Work Energy Approach to Newmark Method (Byrne, 1990) 37 Figure 3. 3 Idealized Pre- and Post- Liquefaction Behaviour of Sand (Haile et al., 1996) 38 x Figure 3. 4 Work Energy Approach to Extended Newmark (Jitno, 1995) 38 Figure 3. 5 Linear and Non-Linear Stress-Strain Curves (Byrne, 1990) 39 Figure 4.1 Cyclic Resistance Ratio Based on S P T Blow Count, M 7.5 (NCEER, 1997) 50 Figure 4. 2 Cyclic Induced Liquefaction Behaviour of Sands 51 Figure 4. 3 Cyclic-Induced Shear Stress History 52 Figure 4. 4 Relationship Between Normalized Shear Stress and Number of Cycles to Liquefaction (Byrne, 1990) 53 Figure 4. 5 Modeled Cyclic Induced Behaviour of the Soil Element 55 Figure 4. 6 Stress-Strain Response of Soil 55 Figure 4. 7 Actual and Idealized Post-Liquefaction Behaviour of Saturated Sands Under Cyclic Loading 56 Figure 4. 8 Definition of Loading and Unloading in Proposed Model 57 Figure 4. 9 Comparison Between Actual and Modelled Hysteretic Loops 58 Figure 4. 10 Linear Elastic-Plastic Behaviour of Soils 59 Figure 5. 1 F L A C Calculation Cycle (FLAC 3.3 manual, 1995) 65 Figure 5. 2 Variation of Normalized Critical Damping Ratio with Frequency (FLAC 3.3 manual, 1995) 67 Figure 5. 3 Shear Moduli and Damping Attenuation Curves 68 Figure 5. 4 Comparison of Dynamic Response Between F L A C and S H A K E 70 Figure 5. 5 Comparison of Displacements Predicted by Proposed Model and Bartlett-Youd 73 Figure 5. 6 Comparison of Response to Different Earthquake Records Scaled to P G A 0.15g ....78 Figure 5. 7 Acceleration Time Histories of Wildlife and Caltechb Earthquakes 79 Figure 5. 8 Acceleration Time Histories of Gilroy and San Fernando (a) Earthquakes 80 Figure 5. 9 Comparison of Response to Different Earthquake Records Scaled to 0.27g 81 Figure 6. 1 Plan and Location of Coquitlam Dam 84 xi Figure 6. 2 Typical Cross Section of Coquitlam Dam (not to scale) 85 Figure 6. 3 Typical Cone Penetration Resistance Profile of Core Material 90 Figure 6. 4 C U Triaxial Test Results of Core Material 93 Figure 6. 5 Cyclic Resistance of Core Material of Coquitlam Dam 94 Figure 6. 6 C U Triaxial Test Results of Liquefied Core Material 96 Figure 6. 7 1979 Computed Cyclic Stress Ratio Contours Based on Hand Calculations 97 Figure 6. 8 1979 Computed Cyclic Stress Ratio Contours Based on Dynamic Analyses 98 Figure 6. 9 1984 Computed Cyclic Stress Contours Based on Hand Calculations 99 Figure 6. 10 S O I L S T R E S S Deformation Analysis Based on Loss in Stiffness 100 Figure 7. 1 Grid Used For Analysis of Coquitlam Dam 105 Figure 7. 2 Material Zones For Analysis of Coquitlam Dam 105 Figure 7. 3 Mohr-Coulomb Parameter Fit to Hyperbolic Parameter 106 Figure 7. 4 End of Construction Displacement Vectors 108 Figure 7. 5 Vertical Displacements Along Centreline of Dam 108 Figure 7. 6 End of Construction Vertical Stress Contours 108 Figure 7. 7 End of Reservoir Filling Flow Vectors 111 Figure 7. 8 End of Reservoir Filling Pore Pressure Contours 111 Figure 7. 9 End of Reservoir Filling Total Head Contours 111 Figure 7. 10 End of Reservoir Filling Displacement Vectors 112 Figure 7. 11 End of Reservoir Filling Vertical Stress Contours 112 Figure 7. 12 Liquefaction Triggering of Coquitlam Dam 118 Figure 7. 13 Limit Equilibrium Flow Slide Analyses 118 Figure 7. 14 Typical Post-Liquefaction Displacement Vector Pattern 120 Figure 7 .15 Comparison of Velocity History of Crest of Dam 121 Figure 7. 16 Zones Predicted to Liquefy by Proposed Method 123 Figure 7. 17 Post-Liquefaction Displacement Pattern Predicted by Proposed Procedure 125 xii Figure 7.18 Post-Liquefaction Displacement Pattern Assuming Liquefiable Zones Triggered at the Same Time 125 Figure 7. 19 Cyclic Shear-Strain Plot Predicted in Dam Core Material 126 Figure 7. 20 Displacement Pattern Due to Volumetric Deformations Only 128 Figure 7. 21 Distorted Grid Due to Post-Liquefaction Deformation and Settlement 128 xiii LIST OF SYMBOLS (N^ S P T blow counts normalized to an overburden stress of 1 tsf (98 kPa) (Ni)eo (Ni) normalized to 60 percent of hammer energy (Ni)6o-cs Equivalent clean sand blow count (tdy)max Maximum dynamic shear stress a Ratio of shear stress on horizontal plane to effective stress, direction of loading with repect to bedding plane B Bulk modulus c Cohesion C R R Cyclic resistance ratio C S R Critical stress ratio (Chapter 2) or Cyclic stress ratio (subsequent Chapters) D 1 0 Effective grain size D 5 0 Mean grain size D m a X Maximum damping ratio A N e q Increment in equivalent number of cycles D r Relative density Au Excess pore water pressure increment e a Axial strain E V Volumetric strain F Force <|) Friction angle of soil FL or Fug Factor of safety against liquefaction fs Shear yield function g Acceleration due to gravity G Shear modulus y Shear Strain YL Limiting shear strain Gliq Liquefied shear modulus Ymax Maximum amplitude of shear strain xiv G m a x Maximum shear modulus g s Shear potential function Ysat Saturated unit weight of soil k Horizontal or vertical seismic coefficient k 2 Shear modulus number (dynamic analysis) k2,max Maximum shear modulus number (dynamic analysis) K a C R R correction factor due to effect of static shear bias kb Bulk modulus number kg Shear modulus number (Mohr-Coulomb, static analysis) KCT C R R correction factor for confining stress X plastic modulus number m Bulk modulus exponent N Number of cycles n Shear modulus exponent P a Atmospheric pressure P G A Peak ground acceleration PT Phase transformation state rd Stress reduction factor due to soil depth Pdry Bulk dry density of soil a ' 3 Minor effective principal stress a ' c Effective confining stress a ' 0 Initial effective confining stress a 'vo Effective overburden stress G'I Major effective principal stress CTd Deviator Stress S P T Standard penetration test S r Residual strength S S Steady state a v 0 Total overburden stress CTxx Normal stress in x direction (% Normal stress in y direction xv o-zz Normal stress in z direction T Shear Stress TIS Cyclic shear stress causing liquefaction in 15 cycles T D Y Dynamic or cyclic shear stress i s t Static shear stress txy Shear stress on horizontal plane v Dilation angle co. Angular frequency com i n Minimum angular frequency E,j Critical damping ratio c^ min Minimum critical damping ratio xvi ACKNOWLEDGEMENTS I would like to express my sincere gratitude to my supervisor, Dr. P.M. Byrne for his support and guidance during the course of my studies and research. I would also like to thank my colleagues in the Geotechnical Engineering Research Group at the University of British Columbia, particularly, Dr. Michael Beaty for reviewing this thesis and Dr. Humberto Puebla for making suggestions for solving problems encountered in F L A C . Financial support provided by BC Hydro and the valuable work experience gained through the Professional Partnership Program was sincerely appreciated. Thank you to Dr. Michael Lee of BC Hydro, who also reviewed this thesis and provided helpful suggestions and edits. Finally, I thank God and Our Lady for the strength needed to help me complete this and also to my family for their loving support, understanding, and most importantly patience in dealing with my banterings. I would also like to thank my office-mates for their friendship and support. Last, but not least, I am very grateful to Karen Savage and Saman Vazinkhoo for their cooperation, encouragement and motivation. xvii Chapter 1. Introduction CHAPTER 1 INTRODUCTION Dynamic loading of saturated granular soils causes an increase in porewater pressure. Depending on the duration and intensity of earthquake motion and the permeability of the soil, the pore pressures may rise to equal the initial overburden stress in which case liquefaction is said to occur. This rise in porewater pressure causes a reduction in the strength and stiffness of the soil. The reduced strength is generally referred to as the residual strength. Liquefaction may cause large displacements depending on the earthquake intensity and duration and on the initial static shear stress acting on the structure. There are many case histories of catastrophic failures of earth structures which have occurred as a result of earthquake-induced liquefaction. One of the more well known examples of such behaviour is the 1978 failure of two tailings dams associated with the Mochikoshi gold mine in Japan in which a large volume of tailings materials was released. Near catastrophic failures have also occurred as a result of liquefaction due to earthquake loading. In 1971, an earthquake hit the San Fernando Valley in California. The earthquake caused major damage to the area, in particular to the Upper and Lower San Fernando Dams. An upstream flow slide which extended to the upper portion of the downstream slope left the Lower San Fernando Dam with only approximately 1.5 m of freeboard. Although the damage to the Upper San Fernando Dam was less severe, the crest of the dam moved 1.5 m downstream and settled about 0.8 m. The movements of the Lower San Fernando Dam and Mochikoshi Tailings Dams occurred because the driving forces within the structures exceeded the residual strengths of the soils. The movements come to a stop when large changes in geometry have caused a reduction in the difference between the driving stresses and residual strengths. If the changes in geometry 1 Chapter 1. Introduction cannot cause the stresses and strengths to be balanced, unlimited deformation, such as those observed in the Mochikoshi tailings dams, may occur. Although flow slides may not occur, situations in which the driving stresses are less than the residual strength can cause be a major concern. Large deformations can be induced from the reduction in stiffness of liquefied materials as well as from inertia forces caused by the earthquake base motions. The damage to the Upper San Fernando is an example of such behaviour. Limit equilibrium analyses have long been used to assess the stability of earth structures under static and dynamic conditions. These types of analyses have been shown to be reliable in predicting the flow slide potential of earth structures. The flow slide potential is assessed by a factor of safety in which the driving stresses are compared to material strengths. If a factor of safety of less than unity is computed then a flow slide may occur. If the factor of safety is greater than unity then a flow slide is not likely to occur. However, the deformations associated with liquefaction may still be significant. Newmark was the first to advance the concept of assessing seismic response by computing displacements rather than a factor of safety. Newmark proposed a single degree of freedom model in which soil is assumed to exhibit rigid plastic behaviour. The method works well in situations where deformations occur on a distinct surface. However, the method does not consider the displacements due to the reduced strength and stiffness of liquefied soils. More sophisticated two-dimensional procedures have been developed as a result of the increase in understanding of the behaviour of liquefied materials. The two more commonly-used methods are Lee's Modified Modulus Method and Jitno and Byrne's extended Newmark Method. Both procedures use finite elements to compute deformations. The Modified Modulus approach assumes that displacements are due to the reduction in strength and stiffness of the liquefied soil and does not consider inertia forces induced by earthquake loading. Jitno and Byrne extended Newmark's approach to a multi degree of freedom system in which the loss in strength and stiffness of liquefied soils is considered. The method has proven to be reliable, however, because the displacements are computed based on energy balance of the system, the computed 2 Chapter 1. Introduction displacements are sensitive to the relative stiffness of adjacent layers as well as to the grid geometry. Rigorous non-linear effective stress methods are more fundamental than the simple total stress procedures. The rise in porewater pressure during earthquake motion can be modeled and liquefaction is triggered in different zones at different times. However, these analyses are complicated and have not been fully validated. A proposed total stress procedure is presented in this thesis. The procedure considers the two main effects of liquefaction induced deformations: the reduction in strength and stiffness and the inertia forces due to earthquake loading. The method combines triggering and deformation analyses together into one procedure. Triggering of liquefaction in each element can be predicted by weighting the cyclic shear stress amplitude during a prescribed base motion. Post-liquefaction parameters are assigned when liquefaction is considered to have been triggered. The triggering of liquefaction was approximately verified against the results using S H A K E . Furthermore, the resulting displacements were validated against the predictions using Bartlett and Youd's empirical approach. The proposed method was then applied to Coquitlam Dam. In this thesis, the results using the proposed procedure on Coquitlam Dam were compared to the results using the more conventionally used methods to predict liquefaction induced displacements. 3 Chapter 2 Behaviour of Saturated Undrained Soils CHAPTER 2 BEHAVIOUR OF SATURATED UNDRAINED SOILS 2.1 Introduction Cohesionless soils can be susceptible to excessive deformations under undrained loading conditions. This phenomenon, liquefaction, is due to porewater pressure build up and decrease in effective stresses and can be initiated by static or dynamic loading or by changes in groundwater conditions. Undrained response due to static or monotonic loading as opposed to cyclic loading is usually considered separately. Interest in monotonic loading is related to failure associated with flow slide. The characteristic behaviour is large deformations due to low shear resistances. Interest in cyclic loading is related to the accumulation of deformations under repeated loading such as in earthquake shaking and/or a flow slide. 2.2 Monotonic Loading Behaviour Monotonic loading behaviour of saturated sands has been investigated by a number of researchers [eg. Castro(1969), Lee and Seed(1970), Chern(1985), Thomas(1992), and Sivathayalan(1994)]. Typical triaxial compression undrained response can be characterized by the 3 curves shown in Fig 2.1. Different responses are associated with the relative density and initial confining stress prior to shearing. At a given initial effective confining stress, the behaviour changes from type 1 to 3 with increasing relative density. Type 1 response represents a strain softening or contractive behaviour wherein the sand exhibits a large reduction in shear resistance under large strains. This type of response is termed liquefaction (Castro, 1969; Seed, 1979; Vaid and Chem,1985). Initial development of strains and pore pressure is slow until an effective stress ratio corresponding to a peak deviator stress is 4 Chapter 2 Behaviour of Saturated Undrained Soils Normal Effective Stress Figure 2.1 Characteristic Undrained Monotonic Response of Saturated Soils (Chern, 1985) 5 Chapter 2 Behaviour of Saturated Undrained Soils attained. This effective stress ratio, termed the critical stress ratio (CSR) (Vaid and Chern,1983,1985), has been found to be unique for a given sand in undrained compression loading (Vaid and Chern, 1985; Vaid et al., 1989). In contrast, the C S R in extension or in simple shear is not unique and is lower than in compression. The C S R increases with an increase in void ratio but is always smaller than that in compression mode. Once the C S R has been attained, the soil experiences a sudden decrease in resistance along with a sudden increase in pore pressure with strain. The soil resistance eventually reaches a minimum constant value, termed the steady state (SS) strength (Castro, 1969). The S S strength is related to initial void ratio (Castro, 1969). However, at a given void ratio, the S S strength in compression is greater than that in extension (Vaid et al., 1989). In type 2 response, strain softening behaviour occurs over a limited strain range then transforms into strain hardening behaviour after a minimum strength is reached. This type of behaviour has been termed limited liquefaction (Castro, 1969). Like type 1 response, strain softening is initiated at the C S R and continues until a minimum undrained strength is reached. However, upon further deformation, the soil begins to regain strength as the pore pressure decrease occur. The state where the soil changes from contractive to dilative behaviour is termed phase transformation (PT) (Ishihara et al., 1975). In the effective stress diagram, this state is represented by the point of sharp reversal of the effective stress path. The friction angle mobilized in the phase transformation phase is identical to that mobilized in steady state in type 1 response (Vaid and Chern, 1985). Upon further straining past the PT state, the effective stress path follows the ultimate failure line. This line is unique for a given water deposited sand (Vaid and Chern, 1985; Vaid and Thomas, 1994). Type 3 response represents strain hardening behaviour with no loss in shear resistance. The shear resistance increases with increasing strains. In contrast to type 2 response, the reversal in effective stress path is gradual once it has crossed the P T / S S line and approaches the ultimate failure line. 6 Chapter 2 Behaviour of Saturated Undrained Soils 2.3 Cyclic Loading Behaviour Undrained cyclic loading causes a cumulative increase in pore pressure which may lead to the development of large shear strains. The soil would then be considered as having liquefied. This strain development can be due to either liquefaction, limited liquefaction, or cyclic mobility, depending on the initial state of the sand. A close link exists between the response of sands under monotonic and cyclic loading. It has been shown that liquefaction and limited liquefaction occur in the same manner in cyclic as in monotonic loading (Castro, 1969; Vaid and Chern, 1985). Liquefaction-type behaviour in cyclic loading is shown in Fig. 2.2a and b. As in monotonic loading, pore pressure and strain development are slow until the effective stress state crosses the C S R line. Upon intersecting the C S R line, the effective stress path reaches the PT or S S line, at which point the soil resistance decreases to a minimum value upon further straining. If the initial static shear stress level is greater than the steady state strength, unlimited deformation may occur. The C S R and P T / S S lines are the same in cyclic as in monotonic loading, implying that they are unique (Vaid and Chern, 1985). Figures 2.2c and d show a limited liquefaction-type behaviour under cyclic loading. Strain softening behaviour is initiated when the effective stress path reaches the C S R line. Large strains occur as the stress path moves from the C S R to PT lines. As in monotonic behaviour, a sharp turn around is observed when the stress path intersects the PT line. This results in the development of large strains (Fig 2.2c). Subsequent unloading causes a large increase in pore pressure leading to a transient state of zero effective stress. Strain recovery in this process is small (Fig 2.2c). Further unloading in the extension region of cyclic loading causes large strains to develop as the sand deforms under virtually zero stiffness. Upon subsequent re-loading in the compression region, the sand exhibits very soft behaviour until it reaches the strain obtained prior to the first zero effective stress condition, after which it regains much of its stiffness. Cyclic mobility due to cyclic loading is shown in Figures 2.2e and f. Strain and pore pressure accumulate slowly with increasing number of cycles. Deformation is not due to marked 7 Chapter 2 Behaviour of Saturated Undrained Soils strain softening at any loading stage. Upon reaching the PT line, the effective stress path sharply turns around, accompanied by the development of significant strains. Further unloading causes large pore pressures to develop, bringing the effective stress close to a state of zero effective stress. Subsequent re-loading finally causes a transient zero effective stress condition and the sample undergoes large deformation in the process. Similar to the limited liquefaction response, the concern with cyclic mobility is with the accumulation of limited deformation with continued cyclic loading. The type of response under cyclic undrained loading depends on factors such as deposition, grain angularity, initial void ratio, effective confining stress, initial shear stress level, induced cyclic stress and number of loading cycles. Because liquefaction can be loosely defined as the development of large strains, it can be caused by any of the above three mechanisms. Cyclic resistance to liquefaction is defined as the cyclic stress required to cause a specified strain in a specific number of cycles. In other words, if the cyclic resistance is exceeded, the response of the soil will change from small strain to large strain. 8 Chapter 2 Behaviour of Saturated Undrained Soils Figure 2. 2 Characteristic Undrained Response of Saturated Sands Under Cyclic Loading (Vaid and Chern, 1985) 9 Chapter 2 Behaviour of Saturated Undrained Soils 2.3.1 Cyclic Resistance to Triggering of Liquefaction Cyclic triaxial tests, cyclic simple shear tests, and cyclic torsional tests have been used to evaluate the liquefaction resistance of saturated sands by many researchers. Liquefaction resistance has been found to be dependent on initial confining stress, shaking intensity, number of loading cycles, and relative density. It is generally considered that the ratio of cyclic shear stress to initial confining stress, referred to as the cyclic resistance ratio (Seed, 1984) is an important parameter in liquefaction analysis. For the triaxial loading condition, this ratio is defined as GD/ZG'O, where CTd is a single amplitude of cyclic axial stress and a ' 0 is the initial effective confining stress. The results of cyclic triaxial tests are usually plotted in terms of cyclic stress ratio (CSR) versus number of cycles to liquefaction. Liquefaction in laboratory tests is achieved when a specifed double amplitude axial strain of usually 5 percent is obtained. Typical results are shown in Fig 2.3. The cyclic resistance ratio is usually multiplied by 0.65 in order to represent simple shear loading conditions. This correction is considered conservative since it can range from 0.66 to about 0.78 depending on relative density and confining stress (Sivathayalan, 1994). 10 Chapter 2 Behaviour of Saturated Undrained Soils N u m b e r of C y c l e s , N c O00 Figure 2. 3 Stress Ratio versus Number of Cycle to Liquefaction (De Alba, 1976) 2.3.1.1 Effect of Confining Stress The effects of confining stress on the cyclic resistance of sands have been investigated by means of a cyclic triaxial device on reconstituted Fraser River sands (Thomas, 1992) and on undisturbed samples from frozen Duncan Dam foundation material under high confining stresses (Pillai and Byrne, 1994). The results were plotted in terms of a correction factor, K o , defined as the cyclic resistance ratio of soil at a ' 0 divided by the ratio at a ' 0 = 100 kPa, versus effective confining stress, a ' 0 . As shown in Fig. 2.4 (Pillai and Byrne, 1994), the values for K C T decreases with increasing effective confining stress, implying that the cyclic resistance ratio decreases with increase in effective confining stress. However, the increase in confining stress is associated with an increase in relative density. As a result, the K„ values simulate an increase in confining stress as well as in relative density. The results from both tests plot above the values proposed by Seed and Harder(1990). This indicates that the values previously used may have been conservative. 11 Chapter 2 Behaviour of Saturated Undrained Soils 1.2 1 0.8 0.6 0.4 0.2 0 A — o o © : • X N ^ \ : : .8 \. ~=^y^m^ A _ . ' —- --ai-r-j-rt \ \ -=ft ,, „ "v. \ Duncan Dam X ' — — lings o Dr=60% -\ ' " \ •—-s . Seed and Harder (1990) ^ ^ a 1- Dr=65% ~Dr=60% o - — -Df=66%-O Dr=70% _ « 4 ~Dr=20% \ ~° i . i ' i Dr=40% Dr=60% 500 1000 1500 Effective Confining Stress (kPa) Figure 2. 4 Ko versus Effective Overburden Pressure (Pillai and Byrne, 1994) 2.3.1.2 Effect of Relative Density The effect of relative density on the cyclic resistance of sands has been well documented. Liquefaction resistance increases with increase in relative density at all levels of confining stress although it is more pronounced at lower confining stresses. 2.3.1.3 Effect of Static Shear Bias The effects of static shear on the liquefaction potential of sands have been investigated by a number of researchers (Vaid and Finn, 1979; Vaid and Chern, 1983; Vaid and Chern, 1985; Pillai and Stewart, 1992). It has been found that the effects of static shear depends on the relative density of sands as well as on the confining stress. As a result, a correction factor, K H , defined as the ratio of the C R R at any initial static shear and the C R R at zero static shear has been 12 Chapter 2 Behaviour of Saturated Undrained Soils developed, is usually plotted versus a, the ratio of shear stress on a horizontal plane to the effective normal stress. As shown, in Fig. 2.5, for dense sands, the cyclic resistance increases with increase in initial static shear. In contrast, the cyclic resistance decreases with increase in static shear for fairly loose sands. 2.0 1.5 Ka i.o 0.5 h T D r « 55-70% (N,)6o= 14-22 D r « 3 5 % (N,)6o«4-6 A. T avo' < 3 tsf D r « 45-50% (N,)«o«8-12 0.1 0.2 0.3 a=(Ts/aVf5) 0.4 Figure 2. 5 Correction Factor for Static Shear (NCEER, 1997) 2.3.1.4 Effect of Fines Content Based on field observations of liquefaction and non-liquefaction in sands, it has been considered conservative to ignore the presence of fines contents greater than 5% in silty sands (Seed et al., 1985). As a result, research has concentrated on liquefaction in clean sands. However, recent studies in sand and silt mixures by means of cyclic triaxial device (Kuerbis et al., 13 Chapter 2 Behaviour of Saturated Undrained Soils 1988; Ishihara and Kosecki, 1989; Koester, 1992) have shown that the presence of non-cohesive or non-plastic fines may show liquefaction potential as great or greater than that in clean sands. Liquefaction in silty sands and silts is usually defined by the development of 5% or 10% double amplitude strains rather than the attainment of 100% pore pressure in a cyclic triaxial device. This is because it has been found that significant deformations develop prior to significant development of pore pressure. This is especially true in undisturbed silts, in which, with the exception of loose silts, the development of pore pressure is similar to that in medium dense sands (Zhou et al., 1995). However, in general, over 80% pore pressure is developed in reconstituted silty specimens when liquefaction is considered to have been triggered. The fines content and the plasticity of fines have been found to most strongly affect the cyclic resistance of silty sands. For a given relative density, cyclic strength decreases with increase in the content of low plasticity fines of up to 20 to 30 percent (Kuerbis et al., 1988; Koester, 1992; Singh, 1994; Erten and Maher, 1995). However, if the relative density is expressed in terms of sand skeleton, the cyclic strength increases only slightly with increasing fines content to 20 percent. This suggests that the silts in silty sands only occupy the void spaces in the sand skeleton and the behaviour is essentially controlled by the sand skeleton void ratio (Kuerbis et al., 1988). In contrast, the cyclic strength increases with increase in silt content greater than about 30 percent although never exceeds that of the clean parent sand (Koester, 1992; Singh, 1994). Most researchers have found that the cyclic strength increases with increase in plasticity (Ishihara and Koseki, 1989; Prakash and Sandoval, 1992). However, Zhou et al. (1995) have found that the cyclic resistance decreases with increase in plasticity in reconstituted specimen although the reverse is true in undisturbed specimens. This may be due to the increase in interparticle cementation between fines particles with time, resulting in more difficult particle separation during cyclic loading. To a lesser extent, the cyclic resistance of silty sands is influenced by the overconsolidation ratio and time. Ishihara et al. (1989) compiled a data base of case histories of 14 Chapter 2 Behaviour of Saturated Undrained Soils liquefaction of silty sands. They found that cyclic resistance increases with time and with increase in overconsolidation ratio. 2.4 Post-Liquefaction Cyclic Loading Behaviour The cyclic loading behaviour past the liquefaction state has been investigated in the laboratory by a number of researchers by means of a triaxial device (Seed and Lee, 1966; Vaid and Chern, 1985; Kuerbis, 1989) and a triaxial torsion shear device (Towhata and Ishihara, 1985). Typical stress-strain response under isotropic consolidated triaxial loading conditions are shown in Fig. 2.6. Limited liquefaction behaviour is shown in Fig 2.6a. Axial strains prior to liquefaction are small compared to those upon liquefaction. Once liquefaction has occurred in the extension mode, large strains develop. Loading in the compression region causes the sample to deform under low stiffness over a limited amount of strain after which strain hardening behaviour is exhibited. Subsequent loading in the extension region causes the sample to deform under almost zero stiffness for a greater amount of strain. Cyclic mobility behaviour is shown in Fig 2.6b. Similar response as that in limited liquefaction is observed, although stiffness is regained over a smaller strain range. The post-liquefaction strain increment initially increases with each cycle until a certain number of cycles is reached after which the strain increment decrease slightly. There seems to be a maximum cyclic strain which additional load cycles do not significantly alter. This is consistent with the findings of De Alba et al. (1976) in the large simple shear device. Laboratory test results show that, after initial deformation under zero stiffness, stiffness of liquefied soils increases with increase in strains under cyclic loading. In fact, even loose soils under low static shear can exhibit dilative behaviour. This implies that liquefied soils can transmit earthquake induced shear waves to overlying layers. Field evidence in terms of recorded acceleration data from the Wildlife site in the 1987 Superstition Hills earthquake (Byrne and Mclntyre, 1994) and results from shaking table tests (Sasaki et al., 1992) confirm this. 15 Chapter 2 Behaviour of Saturated Undrained Soils The Wildlife site, California, is a gently sloping area located near the Alamo river consisting of loose saturated silty sands which has liquefied a number of times in past earthquakes. Prior to the 1987 Superstition Hills earthquake, it had been heavily instrumented to monitor the possible re-occurrence of liquefaction. Holzer et al. (1989) presented the field behaviour. Byrne and Mclntyre (1994) analyzed the recorded earthquake accelerations at ground surface and at the top of the unliquefied layer. The accelerations at both locations were integrated to obtain displacements. The surface acceleration was then plotted versus the relative displacements between the two layers. The resulting plot, shown in Fig. 2.7, is similar to a stress versus strain plot since shear stress and strain are proportional to acceleration and relative displacement respectively. The plot shows a substantial decrease in stiffness at some point during earthquake loading indicating the triggering of liquefaction in some layers between the surface and nonliquefied base. Similar to the plot shown in Fig 2.6, significant deformation is developed upon liquefaction although stiffness increases upon further deformation. A series of shaking table tests were conducted by Sasaki et al. ( 1992) in order to investigate the deformation mechanism of liquefied sands. Acceleration, excess pore pressure, and lateral displacement were recorded and plotted during cyclic loading. In addition, similar to Byrne and Mclntyre, acceleration versus lateral displacement at the surface relative to the shaking table was plotted as shown in Fig. 2.8. 16 Chapter 2 Behaviour of Saturated Undrained Soils AXIAL STRAIN (PERCENT) Figure 2. 6 Post-Liquefaction Stress-Strain Curves (Kuerbis, 1989) Chapter 2 Behaviour of Saturated Undrained Soils 200 Relative Displacement (cm) Figure 2. 7 Surface Acceleration versus Relative Displacement (Byrne and Mclntyre., 1994) Figure 2. 8 Approximate Stress-Strain Relationship during Shaking Table Tests (Sasaki et al., 1992) 18 Chapter 2 Behaviour of Saturated Undrained Soils The plot shows that after a number of cycles, the soil stiffness decreases noticeably causing significant deformations to develop in the sand. However, the soil regains stiffness upon further straining. In addition, consistent with the observations of De Alba et al. (1976), there do not seem to be any significant difference between the strain increments in the 20th and 40th cycles. 2.5 Post-Liquefaction Monotonic Loading Behaviour Failures of the Lower San Fernando Dam and Dam No. 2 of the Mochikoshi tailings dam indicate that large deformations can occur after cessation of earthquake shaking. Monotonic loading tests performed on cyclically liquefied material are used to investigate the behaviour of soils under such conditions. A comparison between the pre- and post-liquefaction monotonic stress-strain behaviour of two types of sand is shown in Fig. 2.9 (Seed, 1979). In both sands, a pore pressure ratio of 100% is developed during cyclic loading prior to monotonic loading. As shown in Fig. 2.9a., the pre-liquefied sand gains strength steadily to attain a peak deviator stress of about 7 kg/cm 2 at approximately 20% axial strain. In contrast, the liquefied sand loses most of its stiffness over a large strain before the sample dilates at about 20% axial strain and exhibits strain hardening behaviour. The sample regains most of its strength when an axial strain of about 40% has been developed. 19 Chapter 2 Behaviour of Saturated Undrained Soils Figure 2. 9 Post-Liquefaction Monotonic Behaviour of Saturated Sands (Seed, 1979) a. Sacramento River Sand, DR = 40%; b. Mine Tailings, DR = 95% As shown in Fig. 2.9b, similar stress-strain behaviour is observed in dense mine tailings sand. In contrast to the looser sample, the sample deforms under zero stiffness over a smaller strain range of nearly 10% and regains most of its strength when an axial strain of 30% has been developed. This implies that the rate of stiffness build-up increases with increasing relative density which is consistent with the findings of Thomas (1992). Simple shear post-liquefaction tests on undisturbed samples of Duncan Dam foundation material have shown that samples with initial static shear stress are stiffer than those without static shear stress (Salgado and Pillai, 1993). Post-liquefaction stiffnesses were reduced by about 2 to 50 times the pre-liquefaction stiffness. It is important to note that samples were cyclically loaded to a maximum shear strain of 4% and not necessarily until a pore pressure of 100% has been developed. 20 Chapter 2 Behaviour of Saturated Undrained Soils The effects of pore pressure ratio developed during cyclic loading on post-liquefaction monotonic behaviour have been investigated by means of triaxial and torsional shear devices (Yasuda et al., 1991,1994; Thomas, 1992). In addition, the effects of additional cyclic loading past the liquefaction stage on the post-liquefaction response and expressed as a factor of safety has been studied (Yasuda et al., 1994). As shown in Fig. 2.10, the liquefied stiffness of soil decreases with increasing excess pore pressure ratio. Stiffnesses were reduced by as much as 500 to 1000 times the original stiffnesses when 100% pore pressure is developed. Stiffnesses are decreased further when samples are cycled past the initial liquefaction stage. Recent studies on the undrained response of loose Syncrude sand under various loading conditions by Vaid et al. (1998) have shown that the liquefied strength and stiffness of the soil is dependent on the direction of loading relative to the direction of deposition. The liquefied sand is strongest when the load is applied in the same direction as the direction of deposition and weakest when applied perpendicular to the direction of deposition. Strengths were reduced by as much as a factor of 5. Thomas (1992) showed that the post-liquefaction undrained stress-strain curves on Fraser River Sand can be characterized into three distinct regions as shown in Fig. 2.11. Region 1 is the region with near zero stiffness. Region 2 is a region of increasing stiffness and can be approximated by a parabolic curve. The stress-strain curve in region 3 is nearly linear. The stress-strain response has been found to become stiffer with increasing density. The length of region 1 increases with decrease in relative density and the slope of region 3 increases with increase in relative density. For loose sands, a similar trend is observed as the stiffness of the response increases with increase in confining stress. However, the effect of confining stress is not as apparent at higher relative densities. 21 Chapter 2 Behaviour of Saturated Undrained Soils (a) 20 30 40 50 60 Shear strain , r (96) 80 (b) Dr=46.4~49.5% (kgf/cm1) o 5 o 4 CO* to <J CD •> l_ 55 2 t_ o o 1 > o ° 0 I 1 ' [ — — i 1 Triaxial compression test —i 1 1 . Toyoura Sand ' •du/oc'=0.38 ^ u / O e = 0 . 2 6 ^ / / / / o '=0.5kgf/m2 . 7/* / / ^u/o ev0.94 / / . 4 u / o c ' e 1 . 0 l i i 4 6 Axial Strain, c (%) 8 10 Figure 2.10 Effect of Excess Pore Pressure on Post-Liquefaction Behaviour of Saturated Sands(Yasuda, 1994) 22 Chapter 2 Behaviour of Saturated Undrained Soils Figure 2.11 Characterization of Post-Liquefaction Curve (Thomas, 1992) 23 Chapter 2 Behaviour of Saturated Undrained Soils 2.5.1 Residual Strength The residual strength and limiting shear strain are the most sought after parameters when testing soils under post-liquefaction monotonic loading. The residual strength value for soils exhibiting contractive behaviour is the steady state strength. However, for soils exhibiting limited liquefaction behaviour, the post liquefaction strength depends on the strain level. It has been suggested that the residual strength is the strength achieved after the dilation process is complete. This residual strength value may be very high if the soil exhibits strong dilative response. However, it is commonly assumed that the residual strength is the strength at the phase transformation state. This value is considered to be conservative. Based on laboratory testing, the residual strength has been found to be influenced by relative density, fines content, and direction of loading. Because the residual strength is also dependent on vertical consolidation pressure, a'vo, it is often expressed as a fraction of a'vo- In general, the residual strength decreases with decrease in relative density and increase in fines content. Triaxial compression and extension, simple shear, and hollow cylinder tests on Syncrude sand (Vaid et al., 1998) have shown that the residual strength is anisotropic and that the residual strength ratio, s r / a ' v 0 , decreases as much as by a factor of 5 as loading changes from compression to extension. Compression loading occurs when the major principal stress is perpendicular to the plane of deposition and is defined as a = 0°. Simple shear loading corresponds to a = 45° and extension loading to a = 90°. As shown in Fig. 2.12, for a void ratio of about 0.70, the values range from s/a'vo ~ 0.3 to 0.06. The residual strength can be determined directly from laboratory testing or indirectly by correlating equivalent clean sand S P T (N^o values to a residual strength using the empirical chart developed by Seed and Harder (1990) and shown in Fig. 2.13. The residual strength can also be obtained by using empirically derived equations such as the ones developed by Stark and Mesri (1992) and Byrne (1990) as shown below: a) Byrne (1990) 24 Chapter 2 Behaviour of Saturated Undrained Soils Sr =0.0284 Pa-e(°-173<JV')»-) > 0.087 o r , , / Eq. 2. 1a b) Stark and Mesri (1992) = 0.0055 - ( J V , ) 1 ^ 60-e.v Eq. 2. 2b 400 a ' ^ 4 0 0 kPa, b=0.5, e.=0.790 amt b and a r held constant during shear a0 = 0 30 — 4 5 -90 8 10 1^ - £3 ( 55) Figure 2.12 Shear strength of strain-softening Syncrude Sand (Vaid et al. 1998) However, the residual strength ratio of Duncan Dam foundation material has been found to be 0.21, which is almost 3 times greater than that predicted by Stark and Mesri's empirical relationship (Salgado and Pillai, 1993). Lab test results on Tia Juana silty sand and Lagunillas sandy silt yielded residual strength ratios varying from about 0.08 to 0.18 (Ishihara, 1993). Back-calculated residual strength ratios from case histories of silty sands which have liquefied yielded values in the same range. 25 Chapter 2 Behaviour of Saturated Undrained Soils 2 0 0 0 1600 o z ui cr cc < X io Q Ul z < cr Q z 3 _l < g ui CE 1 2 0 0 8 0 0 4 0 0 • E A R T H Q U A K E - IN0UCE0 LtOUEFACTION ANO SLIOINC CASE HISTORIES WHERE S P T OATA AND R E S I D U A L STRENGTH P A R A M E T E R S HAVE S E E N M E A S U R E D . o • E A R T H Q U A K E - INOUCEO LIQUEFACTION ANO SLIOING CASE HISTORIES WHERE S P T OATA A N D R E S I D U A L STRENGTH P A R A M E T E R S HAVE B E E N ESTIMATED. CONSTRUCTION - INOUCEO L I Q U E F A C T I O N AND SLIOING C A S E HISTORIES. LOWER SAN F E R N A N D O DAM O 4 8 12 16 2 0 24 E Q U I V A L E N T C L E A N S A N D S P T B L O W C O U N T . ( N , ) 6 0 . c s Figure 2.13 Relationship Between Corrected "Clean Sand" Blowcount (N^eo-cs and Undrained Strength (Sr) (Seed and Harder, 1990). 2.5.2 Limiting Shear Strain The limiting shear strain is the limited amount of shear strain that could be developed during cyclic loading regardless of the cyclic stress ratio and number of cycles. It is often assumed to be the strain corresponding to the point at which the residual strength has been developed. The limiting shear strain can be determined from laboratory testing or by correlating (Ni)6o-cs values with limiting strain values as proposed by Seed et al. (1985). Alternatively, the limiting strains can be obtained from the empirical relationship developed by Byrne (1990): ^ = 1 0 ( 2 2 - 0 . 0 5 . < / V , W E q 2 2 The limiting shear strain depends on the relative density of the soil. Loose soils exhibiting contractive behaviour may have limiting shear strain values ranging from 10 to possibly over 30%. 26 Chapter 2 Behaviour of Saturated Undrained Soils On the other hand, the limiting shear strain of medium dense to dense soils may only range from 2 to 5%. 2.6 Post-Liquefaction Volumetric Deformation Excess pore pressures are generated when undrained saturated soils are subjected to cyclic loading caused by earthquake shaking. Pore pressures dissipate during and at the cessation of earthquake loading causing volumetric strains. These strains manifest on the surface as ground settlement. Although post-liquefaction settlement may not be as damaging as lateral ground displacements, field evidence have shown that settlement should not be ignored. Lee and Albaisa (1974) investigated the cyclic induced settlement behaviour of saturated sands by means of cyclic triaxial tests. They found that the reconsolidation volumetric strains due to dissipation of excess pore pressure increases with increasing grain size of the soil, decreasing relative density, and increasing excess pore pressure generated during cyclic loading. In addition, they found that the volumetric change increases with additional dynamic loads past the condition of development of 100% pore pressure ratio. Tatsuoka et al. (1984) studied reconsolidation volumetric strains after 100% excess pore pressure has been generated by cyclic undrained simple shear loading and found that the amount of settlement strongly depends on relative density and maximum shear strain developed in the soil. This agrees well with the findings of Lee and Albaisa. They also found that settlement is insensitive to effective pressure prior to cyclic loading. Tokimatsu and Seed (1987) plotted available data in terms of volumetric strain versus relative density. Their plot shows that volumetric strains are strongly dependent on relative density and cyclic shear strain. However, their data show that there is a maximum cyclic shear strain above which cyclic shear strain does not affect the reconsolidation volumetric strains. Ishihara and Yoshimine (1992) plotted available data in a series of curves in terms of volumetric strain versus maximum shear strain for several values of relative density. The plot, shown in Fig. 2.14, supports the idea that a maximum volumetric strain due to consolidation 27 Chapter 2 Behaviour of Saturated Undrained Soils exists. Interestingly, the value is attained when the maximum amplitude of shear strain is at least 8% for all relative densities. The data can also be plotted in terms of factor of safety against liquefaction. As a result, Ishihara and Yoshimine prepared a chart to approximate the volumetric strains if both the factor of safety against liquefaction and the relative density are known. This chart is shown in Fig. 2.15. Figure 2.14 Volumetric Strain versus Maximum Amplitude Shear Strain for Different Relative Densities (Ishihara and Yoshimine, 1992) 28 Chapter 2 Behaviour of Saturated Undrained Soils Figure 2.15 Volumetric Strain versus Factor of Safety Against Liquefaction (Ishihara and Yoshimine, 1992) 29 Chapter 2 Behaviour of Saturated Undrained Soils 2.7 Summary The understanding gained from laboratory tests have allowed engineers to develop methods to model liquefaction behaviour and predict response. This chapter reviewed the behaviour of soils under saturated undrained monotonic and cyclic loading based on extensive laboratory testing by numerous researchers. Deformations associated with liquefaction can be the result of a loss in strength and stiffness upon liquefaction or the dissipation of the excess porewater pressures generated during liquefaction. Deformations obtained prior to triggering of liquefaction are small compared to those that occur after liquefaction. Upon liquefaction, the soil temporarily behaves like a liquid and may deform significantly under a driving stress. The amount of strain the soil must undergo before recovering strength and stiffness depends on its relative density. Despite the liquid-like behaviour over a relatively large strain range, it has been shown that liquefied soil can transmit earthquake induced shear stresses. 30 Chapter 3 Evaluation of Earthquake Induced Deformations of Earth Structures CHAPTER 3 EVALUATION OF EARTHQUAKE INDUCED DEFORMATIONS OF EARTH STRUCTURES 3.1 Introduction There are three basic concerns when dealing with seismic analyses of earth structures. The first is whether or not cyclic loading due to earthquake shaking will induce liquefaction. This procedure is well documented and usually involves a triggering analysis wherein the the factor of safety against triggering of liquefaction is assessed by comparing the cyclic resistance ratio with the cyclic stress ratio caused by the design earthquake. A factor of safety of less than unity means that liquefaction is likely to occur. The second concern is whether the residual strength of the liquefied materials is adequate to prevent the occurrence of a flow slide. This involves a limit equilibrium analysis where the stability of a potential sliding mass is assessed by a factor of safety. The factor of safety is the ratio of the shear strength of the soil to the driving shear stress. A residual strength is assigned to those zones predicted to have liquefied. A factor of safety of less than unity implies that the earth structure is not stable and that large displacements may occur. If the factor of safety is greater than unity, flow failure may not occur, although the displacements may not be acceptable. The magnitude of these displacements is the third concern in liquefaction analysis. Prediction of earthquake-induced displacements of earth structures is a more difficult problem. A reliable estimate should take account of displacements due to inertia forces as well as those due to the reduced liquefied soil stiffness under gravity loads. In the past 30 years, various methods have been developed to estimate seismic deformations of earth dams, embankments and natural slopes. Methods vary from mechanics-based simple one-dimensional methods or rigorous effective stress methods to empirical based formulas. 31 Chapter 3 Evaluation of Earthquake Induced Deformations of Earth Structures Sophisticated effective stress methods are more fundamental and complicated than the more simple one-dimensional methods. However, in some cases, displacements predicted from simpler methods are more appropriate. The various methods developed are summarized along with their respective advantages and disadvantages in the following paragraphs from the most simple to the more complicated. 3.2 Empirical Methods There are mainly two empirical equations used to predict potential seismic-induced deformations. Both equations were developed using field data of liquefaction-induced displacements from past earthquakes and have considered soil conditions and topography. 3.2.1 Hamadaetal. Hamada et al. (1987) compiled data of ground displacements observed from sites liquefied during the 1964 Niigata, 1971 San Fernando, and 1983 Nihonkai-Chubu earthquakes, which have earthquake magnitudes of greater than M7.0. The database involved mainly of gently sloping sites consisting of loose, clean, uniform sand of medium grain size. The data showed that the liquefaction-induced displacements are most strongly influenced by the ground slope and thickness of liquefied layer. As a result, Hamada et al. proposed the following empirical relationship to estimate displacements: i \_ D = 0 . 7 5 / / 2 # 3 E q . 3 . 1 where D = displacement, meters H = thickness of liquefied layer, meters 9 = the larger of the ground slope or slope of base of liquefied layer, percent 32 Chapter 3 Evaluation of Earthquake Induced Deformations of Earth Structures A comparison between displacements estimated from the empirical equation and displacements observed from the field data on which the equation was based found that most of the displacements ranged from about half to twice the predicted value. This large scatter is probably because many key factors such as soil type, soil density, and level of shaking intensity were not accounted for in the equation. The equation implies that soils with S P T (N^o values of 4 and 20 would give the same displacement assuming they have both liquefied. 3.2.2 Bartlett and Youd Bartlett and Youd (1995) used a database of liquefaction induced lateral spread case histories from eight different major earthquakes in Japan and the United States to develop their more comprehensive empirical method. The earthquakes ranged in magnitude from M 6.4 to M 9.2. They used regression analysis in order to determine which factors such as earthquake magnitude, soil conditions, and topography, most strongly influence lateral spread displacement. Bartlett and Youd found that two models with different parameters were required to predict lateral spread displacement for free face and ground slope conditions. Their models are as follows: a. Free Face Model log(DH+0.01) = -16.366 + 1.178M - 0.927 log R - 0.013R + 0.657 logW + 0.348 logT 1 5 + 4.527 log(100-Fi5) - 0.922 D50 1 5 Eq. 3. 2a b. Ground Slope Model log(DH+0.01) = -15.787 + 1.178M - 0.927 log R - 0.013R + 0.429 log S + 0.348 logT 1 5 +4.527 log(100-F 1 5) - 0.922 D50 1 5 Eq. 3. 2b where D H = horizontal displacement, meters M = moment magnitude of earthquake R = horizontal distance to nearest seismic source, kilometers W = 100 * Height of free face/L, distance from channel S = ground slope, percent 33 Chapter 3 Evaluation of Earthquake Induced Deformations of Earth Structures T 1 5 = cummulative thickness of liquefied layer with (N^o < 15, meters F15 = average fines content in layer T 1 5 , percent D50 1 5 = average mean grain size in layer T15, millimeters Although Bartlett and Youd's model is superior to Hamada's empirical equation, the results are very sensitive to earthquake magnitude, distance to seismic source, and to the average fines content. Because the model was developed from primarily western U.S. and Japanese data, it is most applicable to sites with similar characteristics. For example, the model yields reasonable predictions for regions with earthquake magnitude ranging from 6.0 to 8.0, underlain by shallow layers of sandy material (F 1 5 < 50%) with (N^o values less than 15. In addition, the depth to the liquefied layer should be less than 15 meters. Bartlett and Youd compared the displacements measured from their data set and compared them with those predicted by their equation. About ninety percent of their estimated displacements are within half to twice the observed displacements. 3.3 Pseudo-Static Methods The seismic stability of earth structures were initially assessed by a pseudo-static limit equilibrium analysis. This type of analysis is similar to that used to evaluate flow slide potential. However, in contrast, the effects of earthquake shaking are accounted for and represented by constant horizontal and/or vertical accelerations. The pseudostatic accelerations produce inertial forces which act through the centroid of the failure mass. A factor of safety of less than unity may be acceptable because the actual vibrational earthquake loads act only for a finite time. This type of analysis gives an index of stability but does not give any estimation of displacements. Newmark (1965) was the first to develop a procedure wherein displacements can be predicted. 34 Chapter 3 Evaluation of Earthquake Induced Deformations of Earth Structures 3.3.1 Newmark's Method Newmark proposed a simple method to predict earthquake-induced displacements of earth structures. A potential sliding soil wedge is modelled as a rigid block resting on an inclined plane as shown in Fig. 3.1. water table m.g.A(t) r7~77 m.g.N earthquake motion a. Potential sliding surfaces of an earth dam. b. Idealized potential sliding mass. Figure 3.1 Idealized Potential Sliding Slope (Newmark, 1965) The soil is assumed to exhibit rigid plastic behaviour where deformations occur when earthquake driving forces exceed the yield resistance of the block. The yield resistance is the force which causes the factor of safety against sliding to equal unity and has a value of mgN where m is the mass of the block, g is the gravity acceleration, and N is the yield acceleration in terms of fraction of gravity. Displacements of the block for any time history of base motion can be computed by double integrating the difference between the base acceleration and yield acceleration at times where the base acceleration is greater than the yield acceleration to obtain relative velocity. The relative velocity is then integrated to obtain displacements. Newmark applied his procedure to estimate the displacements of a sliding block due to four western US earthquakes which were normalized to the same peak acceleration and velocity. He considered cases where the soil resistance is similar uphill as downhill (symmetrical) and 35 Chapter 3 Evaluation of Earthquake Induced Deformations of Earth Structures where the yield acceleration in one direction is significantly higher that the other (unsymmetrical). The following equations were proposed as a result of his analyses: a. Symmetrical Resistance V2 N D = - — • ( ! - - ) E q . 3 . 2 b. Unsymmetrical Resistance 2gN A V2 A \ A N/A > 0.5 D = (1 ) • — Eq. 3. 3a. 2gN A N V2 A 0.13 < N/A < 0.5 D = Eq. 3.3b. 2gN N N/A < 0.13 D = ^— Eq. 3.3c. 2gN where V = velocity of the mass, taken as peak ground velocity A = maximum earthquake acceleration in terms of fraction of g. Newmark's equation can also be developed from energy principles where the work done by the external forces (W e x t) minus the work done by the stress field (W i n )) is equal to the change in kinetic energy. This principle is expressed as : We,-Wint=±M-(Vf2-V2) = -±M-V2 Eq. 3. 4 where the final resting velocity, V f, is equal to zero, and V is the specified initial velocity of the block which comes from earthquake forces greater than the yield resistance of the potential sliding block. 36 Chapter 3 Evaluation of Earthquake Induced Deformations of Earth Structures The external force on the block, as shown in Fig. 3.2. is the gravity driving force which is constant with displacement. The work done by the external force is the area beneath the driving force line. The work done by the internal force is the area beneath the soil resistance line. If only one earthquake pulse is considered, Eq. 3.5. can manipulated and reduced to V2 D = Eq. 3. 5 2gN If six velocity pulses were considered this equation would be identical to Newmark's formula in Eq. 3.4c for unsymmetrical case N/A < 0.13. V = Velocity M = Mass of the block D = Seismic displacement ^ V Force (b) oil Resistance Driving force, Mg sin Ot Displacement, D Figure 3. 2 Work Energy Approach to Newmark Method (Byrne, 1990) This method has proven to be very useful for predicting earthquake-induced displacements. However, its assumption that soils behave in a rigid-perfectly plastic manner is too simplistic and may cause the displacements to be underestimated. Also, because yield acceleration depends on effective stress, hence on pore pressure, it would be too difficult to obtain a value for yield acceleration because pore pressures change considerably during earthquake loading (Seed, 1966). Finally, the method only gives one displacement value along a potential sliding surface rather than an overall displacement pattern. 37 Chapter 3 Evaluation of Earthquake Induced Deformations of Earth Structures 3.3.2 Byrne's Extended Newmark Approach In contrast to the assumption made in Newmark's approach, soil does not exhibit rigid perfectly plastic behaviour when liquefied but rather a large reduction in stiffness is observed when the pore pressure rises to cause the effective stresses to drop down to zero. The stress-strain curves more ideally resemble the curves shown in Fig. 3.3. As a result, Byrne (1990) extended Newmark's approach to incorporate the essentials of the stress-strain characteristics of loose saturated sands into a work-energy approach, as shown in Fig. 3.4. b H O DZ m ID I •p Pre-liquefaction S„/a'„ Shear Strain, y Figure 3. 3 Idealized Pre- and Post- Liquefaction Behaviour of Sand (Haile et al., 1996) Strain Figure 3. 4 Work Energy Approach to Extended Newmark (Jitno, 1995) 38 Chapter 3 Evaluation of Earthquake Induced Deformations of Earth Structures Point P represents the pre-earthquake stress state of a soil element in an earth structure. Upon liquefaction, the stress state drops from point P to point Q, which occurs at a low strain. The resistance increases with strain until it reaches a residual strength value of S r . Because the driving stresses from the slope are constant, the system accelerates as it deforms. A s a result, the system has a velocity when it reaches point R, where the resistance is equal to the driving stress. The system continues to deform until the external work done by the driving force is equal to the internal work done by the system. If, in addition, the system has a velocity due to earthquake shaking it will continue to deform to point T. As shown in Fig. 3.4. Newmark's method does not consider the displacements from point P to S. Byrne developed equations to estimate liquefaction induced displacements by modelling a slope as crust lying on a layer of liquefied soil. The liquefied layer is represented as a block resting on an inclined plane. The liquefied layer can be assumed to exhibit linear elastic plastic or non-linear elastic plastic behaviour as shown in Fig 3.5. 4-* Linear Soil Resistance Non-linear x, •a Shear Strain, y Figure 3. 5 Linear and Non-Linear Stress-Strain Curves (Byrne, 1990) The developed equations are the following: 39 Chapter 3 Evaluation of Earthquake Induced Deformations of Earth Structures a. Linear Elastic Plastic D<D L D=DB+L* + *?L Eq.3.6a 1 K D > D L D = D„ + - ( £ > , - £ > „ ) + ^ - Eq. 3. 7b b. Non-Linear Elastic Plastic D < D U D = — ( ^ - ^ - - M V 2 ) E q . 3 . 7 a (\MV2-\KLD2 +SRDL) D > D L D = ^ Eq. 3. 8b where D = liquefaction induced displacement D s t = static displacement corresponding to point R in Fig. 3.5 D L = displacement corresponding to mobilisation of residual strength, S r K L = stiffness of liquefied soil, equal to T s t / y L M = soil mass V = velocity of soil mass Tst = initial static shear stress YL = limiting shear strain The above equations were developed considering one velocity pulse. This was justified because the displacements due to pulses prior to liquefaction are small compared to those that occur upon liquefaction. The displacements due to velocity pulses after liquefaction are accounted for by the introduction of a factor of safety against liquefaction (Jitno, 1995). Table 3.1 was proposed by Byrne (1996) in which the limiting shear strain and residual strength ratio can be determined if the S P T (N^eo value and the factor of safety against liquefaction are known. The 40 Chapter 3 Evaluation of Earthquake Induced Deformations of Earth Structures effect of the accumulation of strain with velocity pulses is accounted for in the residual strain. For example, a strain increment of loose sand may be 2 or 3% per cycle. Ten cycles would cause 20 to 30% shear strain. Table 3.1 Post-Liquefaction Stress-Strain Parameters (NOeo s,/a'o 7L F T R I G « 1 . 0 FTRIG ® 0.5 0-4 0.05-0.10 2 5 - 5 0 > 100 4- 10 0.10-0.20 1 0 - 2 5 3 0 - > 100 10- 15 0.15-0.40 8 -15 2 0 - 3 5 15-20 0.30 - 0.50 5 -10 15-25 >20 >0.50 <5 < 15 The displacements predicted from this method were compared to displacements observed from North American case histories (Jitno, 1995). The predicted displacements using a linear stress-strain assumption were close to the observed displacements. In contrast, most of the predicted displacements using a non-linear stress-strain assumption were within half to twice the measured values. This suggests that a linear stress-strain assumption may be more appropriate for North American case histories. Although the methods described above consider inertia forces and the loss in stiffness due to liquefaction, they predict displacements on only a distinct failure surface. As a result, they are not able to predict overall deformation of earth structures. 3.4 Two Dimensional Methods One-dimensional methods are only able to predict liquefaction-induced displacements on a potential failure surface. In the field, however, deformations may result from the accumulation of strain increments throughout the earth structure. As a result, more realistic displacements can be predicted by analyses performed on a more general field rather than on a specific surface. 41 Chapter 3 Evaluation of Earthquake Induced Deformations of Earth Structures Simple two-dimensional methods using the finite element approach have been proposed by a number of researchers, a few of which will be summarized below. 3.4.1 Modified Modulus Approach Lee (1974) and Yasuda et al. (1992) among others have proposed a simple procedure for predicting liquefaction-induced displacements of earth structures. Permanent displacements occur under the pre-earthquake static stresses in the structure as a result of a loss of strength and decrease in stiffness of the liquefied soil. A hyperbolic stress-strain relationship is assumed for the soil. The procedure requires that a simple gravity turn on finite element analysis be performed twice: once using the pre-earthquake stiffnesses to establish reference deformations and compute pre-earthquake stresses; the second time using the reduced soil stiffness under the stresses computed previously kept constant. The difference in deformation between both analyses is computed and assumed to represent the deformation due to earthquake. Although this method is able to give an overall deformation pattern of an earth structure and accounts for reduction in stiffness due to liquefaction, it may underestimate the earthquake-induced displacements. This is because it does not consider the effect of inertia forces or deformation due to dissipation of excess pore pressures. However, there are field examples such as the failure of the Lower San Fernando Dam and Dam No. 2 of the Mochikoshi tailings dam, where deformations occur after cessation of earthquake shaking and inertia forces due to earthquake loading need not be considered. 3.4.2 Strain Potential / Dynamic Stress Path Approach To account for the effects of inertia forces, Seed and his coworkers proposed a procedure for predicting earthquake-induced displacements using a strain potential concept. The procedure involves a series of steps which must be more or less rigorously observed and approached with caution in order to obtain reasonable results (Seed, 1979). The steps can be summarized as follows: 42 Chapter 3 Evaluation of Earthquake Induced Deformations of Earth Structures a) Compute pre-earthquake stresses in embankment using finite element method. b) Compute the stresses and dynamic response of the structure induced by a selected design earthquake using a dynamic analysis procedure. c) Subject representative samples of the embankment materials to the combined effects of the initial static stresses and the superimposed dynamic stresses and determine their effects in terms of the generation of pore water pressures and the development of strains. Perform a sufficient number of these tests to permit similar evaluations to be made, by interpolation, for all elements comprising the embankment. d) Use strains obtained from laboratory induced by the combined effects of dynamic and static stresses to estimate the overall deformation of the dam. In this method, effects of the earthquake are represented by equivalent static stresses which can be converted to nodal point forces in a finite element mesh. The equivalent static stresses are read off from static stress-strain curves of soils comprising each element of the embankment and represent the stresses required to cause the strain potentials as determined in step b. Although this method accounts for inertia forces in the embankment, it has been found to underestimate displacements in Upper San Fernando Dam. This may be due to the fact that it assumes that the post-earthquake curves can be represented by the pre-cyclic stress-strain curve. This procedure may be applicable to elements which are not predicted to liquefy but is not valid for those elements which have. As shown in Chapter 2, the initial shear modulus of post-liquefaction stress-strain curves are significantly less than that of the pre-cyclic stress-strain curve. As a result, this method only partly accounts for the effect of a reduction in stiffness. 3.4.3 Jitno and Byrne's Approach In order to obtain an overall deformation pattern, Jitno(1995) extended Byrne's one-dimensional method as described above to a two-dimensional system where a pseudo-dynamic finite element approach is used. The procedure is similar to the modified modulus approach in 43 Chapter 3 Evaluation of Earthquake Induced Deformations of Earth Structures that two finite element analyses are performed: the first using pre-earthquake parameters to establish reference deformations; the second using post-earthquake stress-strain parameters. However, in this procedure, the displacements in the second analysis must satisfy the energy principles as described in section 3.3.1. The displacements can be computed from the solution of the following equation: where [K] = global stiffness matrix of system {A} = nodal displacement {F} = static load vector {AF} = additional load vector to produce energy balance of Eq. 3.5 The internal work done in a multi-degree-of-freedom system is the work done by the element stresses and strains and the external work corresponds to the work done by a static load vector, {F}-{A}T. The additional force, AF, is equal to the product of a horizontal or vertical seismic coefficient, k, and the weight of a soil element, w. However, k does not correspond to the peak ground acceleration but is rather the result of an iterative procedure to balance the energy equation. 3.5 Non-Linear Effective Stress Analysis The best approach to accurately predict earthquake-induced displacement is a non-linear effective stress analysis wherein the soil response during an earthquake is simulated by means of a constitutive model developed from experimental observation. Many constitutive models have been developed which can more or less accurately predict the cyclic behaviour of soils although no one model can predict the behaviour of all types of soils. Most of the proposed models are based on complicated plasticity theories which incorporate different hardening properties. A non-associative flow rule can be used which fully couples shear stress and volumetric strain. If the volume is constrained, shear-induced Eq. 3. 8 44 Chapter 3 Evaluation of Earthquake Induced Deformations of Earth Structures contractive volumetric strains cause an increase in pore water pressures thus a decrease in effective stresses. Displacements can be computed by solving the equations of motion using a step-by-step integration procedure. Changes in effective stress can be considered at every step. Unlike the simpler methods, in which liquefiable zones are assumed to trigger at the same time, both the triggering of liquefaction in different elements and the post-liquefaction displacement can be predicted. Although the procedure is more fundamental than other methods, it is complex, time consuming, and as a result costly. In some cases the simpler methods can equally predict the earthquake induced displacements. An effective stress analysis is preferable to the simpler methods particularly when drainage effects are important (Byrne, 1996). 3.6 Summary The various available methods for predicting earthquake-induced displacements have been summarized. Most of the methods either use simplifying assumptions which neglect key parameters to a reliable prediction of liquefaction-induced displacements or are so complex that they become time-consuming, expensive, and inefficient for most applications. There is room for improvement in the procedures for predicting earthquake-induced displacements. The following chapter presents a proposed procedure which attempts to capture the essential features of soil behaviour under cyclic loading in a total stress approach. 45 Chapter 4 Proposed Two-Dimensional Method CHAPTER 4 PROPOSED TOTAL STRESS METHOD 4.1 Introduct ion As discussed in the previous chapter, a reliable simple procedure to predict liquefaction-induced displacements should identify liquefied zones and account for displacements due to inertia forces and those due to the reduced stiffness of the liquefied soil. The present available simple methods use assumptions which overly simplify the effects of earthquake inertia forces and fail to adequately model the stress-strain behaviour of liquefied soil. As a result, a total stress dynamic analysis procedure is proposed. In the proposed procedure, stiff pre-liquefaction stress-strain parameters are used initially. The time history of shear stress for a prescribed base input motion is computed for each element. Each pulse is weighted and when or if sufficient pulses occur, liquefaction is triggered for that element. Soft post-liquefaction stress-strain parameters are assigned to the element and the analysis is continued. In this manner, the various elements of the soil structure liquefy and soften at different times and allow both the amount and pattern of displacement to be predicted as shaking proceeds. In this chapter, the proposed two-dimensional method is presented. The procedure by which liquefaction is triggered in each element as well as the constitutive model used to represent the stress-strain behaviour of the liquefied soil are discussed. Finally, the key parameters to be used in the analysis are presented along with recommended values. First, however, the traditional triggering analysis procedure is summarized and presented. 46 Chapter 4 Proposed Two-Dimensional Method 4.2 Triggering of Liquefaction The determination of liquefied zones has traditionally been assessed by using the method proposed by Seed and his coworkers (1990). A factor of safety against liquefaction, F L , which compares the cyclic resistance ratio, C R R , with the cyclic stress ratio, C S R , is computed: F L = C R R / C S R Eq. 4. 1 A factor of safety of less than unity implies that liquefaction would likely occur. 4.2.1 Cyclic Stress Ratio The C S R is usually evaluated using a total stress equivalent linear dynamic analysis program such as S H A K E (Schnabel et al., 1972) or F L U S H (Lysmer et al., 1974). Both programs use a frequency domain approach in which a base motion is represented by its harmonics. The response in terms of shear stresses or accelerations is the sum of the responses from each harmonic. The response is based on the solution of the wave equation. The nonlinearity of soil behaviour in terms of shear modulus and damping is accounted for by the use of equivalent linear soil properties. An iterative procedure to obtain shear modulus and damping values compatible with the effective strains in each layer is used. S H A K E performs a one-dimensional dynamic response analysis based on a continuous solution of the wave equation. In contrast, F L U S H performs a two or an approximate three-dimensional analysis by assuming a lumped mass system. Two-dimensional analyses can generally provide sufficient accuracy except in cases where the structure is subject to significant three-dimensional effects such as high dams in narrow canyons wherein 3-D analyses may be preferable (Seed and Harder, 1990). One-dimensional dynamic analyses are generally not recommended for dams although they can provide reasonably accurate estimates of cyclic stresses within embankments, especially ones with flat slopes and high crest length to dam height ratio. However, one- and two-dimensional analyses of the Lower San Fernando Dam have shown that one-dimensional 47 Chapter 4 Proposed Two-Dimensional Method analyses underestimate both accelerations and cyclic shear stresses near the crest and upper faces of the embankment (Seed and Harder, 1990). The C S R is evaluated from the maximum induced shear stress computed by the dynamic analysis from the following relationship: C S * = f t 6 5 ' ( r * ) n " E q . 4 . 2 where ( T d y ) m a x = maximum shear stress computed from the dynamic analysis, CT'VO = the initial effective overburden stress The factor, 0.65, is an equivalent loading factor which converts the random earthquake induced shear stresses into an equivalent uniform cyclic shear stress acting over a number of cycles. The equivalent number of cycles depend on the earthquake magnitude as depicted in Table 4.1 (Seed et al., 1975), based on a study by Seed and Idriss (1967). For example, an equivalent uniform shear stress equal to 0 .65 (T d y ) m a x would act over 15 cycles for an earthquake magnitude of 7.5. Table 4.1 Equivalent Number of Cycles Magnitude Number of Earthquake Cycles 5'U 2-3 6 5 6 3 / 4 10 A 15 8 1 / 2 26 Alternatively, the C S R can be evaluated from an empirical relationship developed in conjunction with the simplified liquefaction evaluation procedure proposed by Seed and Idriss (1971) based on the results of S H A K E analyses, as follows: CSR = 0 . 6 5 - - ^ r d Eq. 4. 3 where a = peak ground acceleration 48 Chapter 4 Proposed Two-Dimensional Method g = acceleration due to gravity a'vo, a'vo = total and effective overburden stress respectively rd= stress reduction factor due to soil depth 4.2.2 Cyclic Resistance Ratio One of the best methods of evaluating C R R is by testing undisturbed samples of representative materials. The number of cycles to cause liquefaction is determined for different cyclic stress ratios and the results are plotted on a semi-logarithmic scale as shown earlier in Fig. 2.3. The C R R is usually taken as the cyclic shear stress required to cause liquefaction in 15 cycles. However this method is limited to relatively large and important projects due to the difficulty and high cost associated with obtaining high quality undisturbed samples. The C R R is generally determined by using S P T testing and Seed's liquefaction assessment chart, shown in Fig. 4.1, which empirically correlates the C R R with ( N ^ o values based on past liquefaction experience. The C R R from Seed's chart is for an earthquake magnitude of 7.5 corresponding to 15 cycles. As discussed in Chapter 2, the C R R must be corrected for overburden stress, static shear bias, and earthquake magnitude. Alternatively, the C R R can also be evaluated from empirical charts correlating C R R and cone penetration (CPT) tip resistance. Stark and Olsen (1995) have compiled cone penetration liquefaction assessment charts for different types of soils ranging from clean sands, silts to gravelly soils. 4.3 Proposed Method As shown in Chapter 2, there are three main phases in the cyclic-induced liquefaction behaviour of sands which are depicted in Fig. 4.2. In the pre-liquefaction triggering phase the soil exhibits stiff behaviour as it is subjected to cyclic loading. Because the strains developed at this phase are relatively small, the soil shear moduli are usually determined from small strain tests such as resonant column tests. Upon the triggering of liquefaction, the second phase, the soil 49 Chapter 4 Proposed Two-Dimensional Method behaves like a liquid and the shear stresses drop down to zero as the shear stresses try to reverse directions. Further straining in the post-liquefaction phase, the third phase, causes the soil to deform significantly under nearly zero shear stiffness before finally regaining strength upon even further deformation. The shear stiffnesses in the post-liquefaction phase are reduced by factors of 25 up to 500 of the pre-liquefaction stiffness (Haile et al., 1996). 0.6 0.5 o > ca GO o U 0.4 0.3 0.2 0.1 Percent Fines = 35 15 < 5 A|0 I -CRR curves for 5,15, and 35 percent fines, respectively I FINES CONTENT ^ 5% ijf0 Modified Chinese Code Proposal (clay content = 5%) ® 3^0 Adjustment Recommended By Workshop M a r g i n a l N o Liquefac t ion L i q u e f a c t i o n L ique fac t ion Pan - American data Japanese data Chinese data • o A 10 20 30 40 Corrected Blow Count, (Ni)^ 50 Figure 4.1 Cyclic Resistance Ratio Based on SPT Blow Count, M 7.5 (NCEER, 1997) 50 Chapter 4 Proposed Two-Dimensional Method x, Liquefaction Trigger f y Post-Liquefaction Trigger Figure 4. 2 Cyclic Induced Liquefaction Behaviour of Sands The above three phases are captured simply in the proposed procedure. The manner in which the three phases are modeled are presented and discussed. 4.3.1 Pre-Liquefaction Trigger The pre-liquefaction analysis of the proposed method is similar to that of the current methods in that the soil is assumed to exhibit linear behaviour and the cyclic induced shear stress is compared to the cyclic resistance. Unlike the current methods, however, the proposed method does not consider the nonlinear soil behaviour with respect to shear modulus and damping through an iterative procedure. Parameters correspond to the uniform Rayleigh stiffness-proportional viscous damping and constant shear moduli which are estimated as a fraction of the maximum shear modulus, G m a x , parameters are selected prior to analysis. These pre-liquefaction parameters are maintained in each element during the dynamic shear stress history until the element is predicted to liquefy. 51 Chapter 4 Proposed Two-Dimensional Method 4.3.2 Liquefaction Trigger The triggering of liquefaction is evaluated by tracking the dynamic shear stress history for each element. An example of a dynamic shear stress history for a prescribed base motion is shown in Fig. 4.3. The cyclic shear stress, tdy, is computed at each time step and is defined as the difference between the current shear stress, xxy, and the static bias, tst. The static bias, T s„ is equal to the shear stress calculated from static analysis. A maximum value of T D Y defines a half-cycle or a cyclic pulse, so that the dynamic shear stress history for an element can be represented as a series of half-cycles or pulses which cumulatively contribute to the triggering of liquefaction. Liquefaction is predicted to occur when a number of pulses of sufficient amplitude, has been accumulated. A / <^Jy ' } ^ / ^ ^ W r r r r i ! ? \ /* : : Half eyck • t ime j Figure 4 . 3 Cyclic-Induced Shear Stress History As shown in Fig. 4.3, the dynamic shear stress history for an element is non-uniform and must be converted to an equivalent uniform cyclic history in order for liquefaction to be adequately predicted. The amplitude of the uniform cyctic history is taken to be T 1 5 , the cyclic shear stress required to cause liquefaction in 15 cycles. Each pulse is weighted to an equivalent number of cycles, AN e q , at T 1 5 . Liquefaction is predicted when the summation of A N e q is equal to 15. Weighting of each cyclic pulse is carried out using the relationship between x d y and number of cycles required to cause liquefaction. The cyclic shear stress is normalized with 52 Chapter 4 Proposed Two-Dimensional Method respect to T 1 5 . This relationship, as shown in Fig. 4.4, is based on laboratory test data. Application of the relationship is as follows. For T D Y equal to x 1 5 , the corresponding half cycle is weighted such that A N e q is equal to 0.5 and liquefaction is triggered in 30 half-cycles or 15 full cycles in that element. On the otherhand, if r d y is equal to 1.5T 1 5, then A N e q is approximately equal to 2.5 and liquefaction is triggered in 6 half-cycles or 3 full cycles of x d y = 1.5T 1 5, as schematically shown in Fig. 4.4. Similarly, if i d y is equal to about 0.75x1 5, then A N e q is equal to 0.075 liquefaction is triggered in 100 full cycles of the above amplitude in that element. Because the dynamic shear stress histories in all elements are not equal, liquefaction in each element is triggered at different times. 2.0 1.8 1.6 1.4 1.2 in ~> 1.0 0.8 0.6 0.4 0.2 0.0 I I I I I I I I I i l l ! ! J I \ •\ \ Approximate F i t Curve •A I I I I I I T 1 I I I I I I 10 100 Number of Cycles to Liquefaction I I I I II 1000 Figure 4. 4 Relationship Between Normalized Shear Stress and Number of Cycles to Liquefaction (Byrne, 1990) The relationship shown in Fig. 4.5 is mathematically expressed as follows: T D Y / T 1 5 > 1 0 (—1 . 0 ) (0.5-1*) A A ^ = 0 . 5 - 1 0 r ' 5 Eq. 4. 4a 53 Chapter 4 Proposed Two-Dimensional Method (I*-1.0) (0.40-^2-) T d y / x 1 5 < 1.0 ANeq = 0 . 5 - 1 0 r , s Eq. 4.4b where T D Y = half-cycle stress amplitude T N 5 = cyclic shear stress required to cause liquefaction in 15 cycles A N e q = equivalent number of cycles at T 1 5 N 1 5 = 15 cycles 0.5 = half cycle. It has been stated that in order for the equivalent number of cycle weighting process to be consistent, Fig. 4.4 should be represented by a single linear relationship on a log-log plot (Idriss, 1998). However, a study, as shown in the next chapter, showed that this procedure is consistent. Upon liquefaction, the soil temporarily behaves as a liquid and has very low resistance to shear stresses. In this procedure, this transition phase into a liquid state was conservatively ignored and a condition of 100% pore pressure was immediately imposed upon the element as soon as liquefaction is triggered, as shown in Fig. 4.5. This effect is simulated in the numerical model by setting the horizontal stresses equal to the vertical stress (CTXX = a z z = a y y) and the shear stresses equal to zero (T^ = 0) after the summation of A N e q exceeds 15. In addition, post-liquefaction stress-strain parameters are assigned to those zones which are estimated to have liquefied. 4.3.3 Post-Liquefaction Trigger As discussed earlier, loading after liquefaction has been triggered causes the soil to deform considerably under very low shear stresses. However, after a significant amount of strain has been developed, the liquefied soil eventually regains strength and stiffness. The stiffness of the liquefied soil has been found to be reduced by a factor of 25 to 500 times that of pre-liquefied soil as shown in Fig. 4.6. 54 Chapter 4 Proposed Two-Dimensional Method Figure 4. 5 Modeled Cyclic Induced Behaviour of the Soil Element b H CJ V> IS) e in Pre-liquefaction S,/o'„ Shear Strain, y Figure 4. 6 Stress-Strain Response of Soil In the proposed procedure, the liquefied soil is assumed to exhibit bi-linear elastic plastic behaviour. The behaviour observed in the laboratory as well as that idealized by the model are shown in Fig. 4.8. Although the non-linear behaviour of the liquefied soil in loading mode appears 55 Chapter 4 Proposed Two-Dimensional Method to be overly simplified by using a linear stress-strain model, analyses on case histories using both linear and non-linear assumptions have shown that a linear assumption predicts response as well as or better than a non-linear assumption for monotonic post-liquefaction loading (Jitno, 1995). A observed soil ber model approximation I • / / .—1— / j aviour ! • j ^ Y Ylim Y = 1/10 y,im Figure 4. 7 Actual and Idealized Post-Liquefaction Behaviour of Saturated Sands Under Cyclic Loading The plot shown in Fig. 4.2 as well as the plots shown in Chapter 2 show that liquefied soil exhibits stiffer behaviour in unloading than in loading mode. As a result, different loading and unloading bi-linear shear moduli are required to capture the response of the idealized soil behaviour under cyclic loading. An examination of the laboratory data indicates that the ratio of the unloading and loading modulus can vary from about 7 to over 12. Hence, an average constant value of 10 is assumed. This is also shown in Fig. 4.8. Loading and unloading modes need to be defined in order to differentiate the shear moduli assigned to a zone. This is done by storing, updating, and comparing the maximum earthquake-induced shear stress, ( T d y ) m a x , in each zone to the 'current' zone shear stress, T d y , at 56 Chapter 4 Proposed Two-Dimensional Method each time step. If t d y is greater than (T d y) m ax, then the loading shear modulus is assigned to the zone and the ( T d y ) m a x is then updated to equal T d y. In contrast, unloading mode is defined by the case where i d y is less than 0.99 ( T d y ) m a x , and is depicted in Fig. 4.8. If x d y stays the same sign as ( T d y ) m a x , but starts to increase in value, then the unloading modulus is still used in the zone. When t d y changes sign relative to ( T d y ) m a x , then the loading modulus is used and the maximum shear stress is reset to the current shear stress value and updated at each time step. t T k i ^^\\dea\\T.ed pre-liquefaction curve ~ r * ^ / . 0 . 9 9 { t 0 „ (^)_ ( t j „ w 0.99(T„)4 j 1 L Idealized post- liquefaction curve Figure 4. 8 Definition of Loading and Unloading in Proposed Model The use of loading and unloading shear moduli creates shear stress versus strain loops when modeling soil under dynamic loading conditions. These loops, known as hysteresis loops, dissipate energy and act to damp the soil system. Although damping is somewhat accounted for in the system, it is not accurate because the modeled loops only approximate the hysteresis loops observed from laboratory tests and case histories. A comparison between hysteresis loops observed from the model and those observed from laboratory data is shown in Fig. 4.9 and indicates that the damping applied in the model may be greater than that experienced under field and laboratory conditions, especially when hysteretic damping is combined with viscous damping. 57 Chapter 4 Proposed Two-Dimensional Method Different loading and unloading shear moduli were used only after liquefaction has been triggered. This is because equivalent linear properties, which have been obtained in S H A K E and F L U S H , are used in the pre-liquefaction analysis. The loading liquefied shear modulus used is computed from the ratio of the residual shear strength to limiting shear strain. This is depicted in Fig. 4.10. 40 20 CO CO co O o H -20 -40 -15 Actual Hysteresis Loop Modeled HysteresisLoop -10 0 S h e a r S t ra in , y 10 15 Figure 4. 9 Comparison Between Actual and Modeled Hysteretic Loops 58 Chapter 4 Proposed Two-Dimensional Method Figure 4.10 Linear Elastic-Plastic Behaviour of Soils As shown in Chapter 2, soil residual strength is anisotropic and can display up to five times as much strength in compression loading than in extension loading. This was not considered in this version of the proposed model because it does not significanlty affect the model's application presented in this thesis. However, it is recommended that strength anisotropy be incorporated in a more refined version of this proposed model. 4.3.4 Mohr-Coulomb Model 4.3.4.1 Introduction The total stress analysis described here cannot be performed using programs based on a frequency domain approach. Because the cyclic shear stresses must be assessed for at least each half cycle, a computer program based on a time domain approach is needed. The computer code F L A C (Fast Lagrangian Analysis of Continua) which is a two-dimensional explicit finite difference program was selected for the analyses. A subroutine which performs the proposed procedure was written for use in F L A C . A more detailed description of F L A C is presented in Chapter 5. F L A C has several built-in constitutive models, one of which is the simple and verified Mohr-Coulomb model. The procedure uses a modified Mohr-Coulomb model to simulate the soil behaviour under static and dynamic loading conditions. 59 Chapter 4 Proposed Two-Dimensional Method The Mohr-Coulomb model assumes soil exhibits simple linear-elastic plastic behaviour, as shown in Fig. 4.8. In terms of plasticity, the model features a yield surface, a flow rule, and has been modified to include a definition for loading and unloading for dynamic problems. 4.3.4.2 Yield Sudace The yield surface is defined by the Mohr-Coulomb shear yield function: f,=&-o'3N4-r2cJW;<0 Eq . 4 .5 where fs = shear yield function <j> = friction angle CT'I = major effective principal stress c = cohesion 1 + sin d> a ' 3 = minor effective principal stress = 1 - sin (j> To determine loading states, trial incremental principal effective stresses are assumed to be completely due to elastic principal strain increments and related by Hooke's law. Principal effective stresses are updated at each step by summing the incremental principal effective stresses and the previous stresses. If the new stresses exceed the failure criterion, they are modified using plasticity theory. Plastic strain increments are related to the principal stresses by a flow rule. 4.3.4.3 Non-associative Flow Rule The flow rule has the following form: As? = E q 4 6 where Aej P = plastic principal strain increment 5gs= shear potential function corresponding to a non-associative flow rule 5o'i = principal effective stress increment The shear potential function has the form: 60 Chapter 4 Proposed Two-Dimensional Method 8, = ° " i - o - 3 W v Eq. 4. 7 where Af,, = 1 + sin v 1 - sin v , v = dilation angle 4.4 Required Parameters The key parameters in the total stress procedure are the common Mohr-Coulomb variables such as cohesion, c, friction angle, shear modulus, G , and bulk modulus, B. For pre-liquefaction dynamic analyses, the key parameters are the shear modulus and damping. For post-liquefaction analyses, the key parameters are the residual strength, expressed as a cohesion value, and the loading and unloading shear moduli. 4.4.1 Pre-liquefaction Parameters The main parameter in the pre-liquefaction portion of the analysis are the equivalent elastic shear modulus, G , which is taken to be a fraction of G m a x , and damping. The shear modulus at low shear strains, G m a x , can be determined from laboratory tests such as resonant column tests, from empirical correlations with cone penetration data, in-situ shear wave velocities, or from a data base of values. The fraction of shear modulus value to use is more difficult to determine and depends mainly on the input dynamic motion and the nonlinear stress-strain properties of the soil. The value was taken based on the average of the values output in S H A K E or F L U S H analyses. Similarly, the damping value is very difficult to determine. F L A C has the option of using Rayleigh- or "local"-type damping. Unfortunately, there is no database of values for Rayleigh-type damping parameters. As shown in Chapter 5, the value was taken as the average of the values output in S H A K E or F L U S H analyses. 61 Chapter 4 Proposed Two-Dimensional Method 4.4.2 Post-Liquefaction Parameters 4.4.2.1 Residual Strength As presented in Chapter 2, the residual strength can be determined directly from laboratory testing or indirectly by correlating equivalent clean sand S P T (N^o values to a residual strength using the empirical chart developed by Seed and Harder (1990). The residual strength can also be obtained by using empirically derived equations such as the ones developed by Stark and Mesri (1992) and Byrne (1990). It was shown in Chapter 2 that the residual strength depends significantly on the pore pressure developed during liquefaction, relative density, fines content, and direction of loading, i.e. direction of major principal stress. For example, based on Vaid et al.(1998), a change in the direction of loading can change the residual strength by a factor of up to 5. Because the direction of major principal stress can range in the field condition from extension to simple shear to compression (Byrne et al., 1998), it is recommended that the residual strength be varied with direction of loading or that the lowest, most conservative, residual strength be used. 4.4.2.2 Limiting Shear Strain The "limiting" shear strain, as defined in this procedure, is the strain required to mobilize the residual strength. This limiting shear strain is mainly required to define the liquefied shear modulus, G i iq , along with the residual strength. The limiting shear strain value can be determined from laboratory testing and can vary significantly depending on the value of the relative density of the soil as well as on the liquefaction behaviour. Most empirical shear strain values, such as those proposed by Seed et al. (1985), refer to a maximum shear strain value which can not be exceeded rather than as the strain required to mobilize the residual strength although they have been used in that context. In addition, these empirical values were based on pseudo-dynamic post-liquefaction procedures. Because of this, the range of values to be used in the proposed procedure should be determined by comparing the results of analyses using the proposed procedure with field experience. 62 Chapter 4 Proposed Two-Dimensional Method 4.5 Summary A relatively simple total stress procedure for predicting earthquake-induced behaviour has been presented. In this procedure, liquefaction is triggered in each zone or element individually and post-liquefaction parameters are assigned to liquefied zones as they occur. The triggering of liquefaction is predicted by weighting the cyclic stress pulses as they occur to obtain an equivalent number of cycles and accumulating cycles for a prescribed base motion. Unlike other total stress methods, earthquake time histories are directly incorporated in the displacement analysis. Displacements due to the effects of earthquake inertia forces and reduction in soil stiffness are accounted for. The analysis can be performed using a computer program based on a time domain approach such as F L A C . A subroutine incorporating the procedure has been written for use in F L A C . 63 Chapter 5. Verification of Total Stress Procedure CHAPTER 5 VERIFICATION OF TOTAL STRESS PROCEDURE 5.1 Introduction As in any other new method for predicting liquefaction induced displacements, the procedure must be validated against case histories in order to be credible. The proposed procedure was verified against Bartlett and Youd's relationship for ground slope conditions. As summarized in Chapter 3, Bartlett and Youd's equation is based on observations of lateral displacement case histories. The procedure is also compared to two state-of-practice methods in predicting liquefaction induced displacements. The three different methods were applied to analyses of Coquitlam Dam. The triggering process is comparable to the equivalent linear method. The equivalent linear method, which is embodied in the program S H A K E , has been shown to give reasonably accurate results. The method essentially captures the variation of shear modulus and damping ratio with strain. The method approximates the nonlinear behaviour of soil by using an equivalent elastic shear modulus, G , that is compatible with the average strain level for the particular input motion. S H A K E uses a different ratio of G / G m a x for each layer. A similar procedure can be used in F L A C , but for simplicity a constant G / G m a x value is assigned to all zones. In addition, the type of damping available in the computer code F L A C is different from that used in S H A K E . This chapter presents and compares the results of dynamic analyses performed using F L A C and S H A K E to validate the pre-triggering procedure. The post-triggering procedure is verified against Bartlett and Youd's empirical equation. These results are also presented and discussed. The comparison between the proposed procedure and other methods applied to Coquitlam Dam will be presented in Chapter 7. Available data and past analyses of Coquitlam Dam will be summarized in Chapter 6. 64 Chapter 5. Verification of Total Stress Procedure 5.2 Triggering Verification Against SHAKE 5.2.1 Description of FLAC F L A C (Fast Lagrangian Analysis of Continua) is a two-dimensional explicit finite difference program initially developed by Cundall in 1986 for use in rock mechanics problems. It has since been extended for use in soil mechanics. Static, seepage, and dynamic analyses can be performed using FLAC. F L A C uses the full dynamic equations of motion in order to model both static and dynamic systems. The equations of motion are used to compute new velocities and displacements from stresses and forces. Strains are computed from displacments and new stresses are derived from strains through a constitutive or stress-strain relation. This process, as shown in Fig. 5.1, performed in one timestep in each element, is repeated using the new stresses and forces. A static problem is solved by damping the equations of motion. A damping force is applied to each node which is proportional to the unbalanced force. As the system reaches equilibrium, the unbalanced force decreases. New velocities and displacements New stresses or forces Stress/Strain Relation (Constitutive Equation) Figure 5.1 FLAC Calculation Cycle (FLAC 3.3 manual, 1995) 65 Chapter 5. Verification of Total Stress Procedure F L A C is more often used to perform static while either S H A K E or F L U S H is commonly used in dynamic analyses. S H A K E is a one-dimensional dynamic analysis program in which strain dependent damping and shear moduli are represented by equivalent constant values which are obtained by iteration. FLUSH is the two-to-three dimensional version of S H A K E . F L A C is able to handle the nonlinearity of the shear modulus better than S H A K E but the computation time in F L A C is significantly slower than that in S H A K E . In the proposed procedure, a constant equivalent shear modulus ratio, G / G m a x , is used to account for the nonlinearity, for simplicity and to decrease computation time. This constant shear modulus ratio, G / G m a x , can be used in the F L A C model to produce the similar results to S H A K E . The fraction of G m a x in S H A K E varies with depth in order to obtain moduli compatible with effective shear strain at all points of the system. In contrast, the fraction of G m a x used in F L A C is assumed constant with depth. The type of damping used in F L A C is different from that used in S H A K E . There are essentially two types of damping available in dynamic F L A C : local damping and Rayleigh damping. Unlike the damping in S H A K E , Rayleigh damping is frequency dependent. The critical damping ratio, ^ , at any angular frequency, m„ of the system can be found from the mass and stiffness proportional damping constants (FLAC 3.3 manual, 1995). From equations relating ^ with Mi, it can be shown that mass-proportional damping is dominant at lower angular frequency ranges and stiffness-proportional damping is dominant at higher angular frequency ranges. The damping ratio is a minimum at a frequency where mass and stiffness damping each supply half of the total damping force. In F L A C , the user defines this minimum or critical damping ratio £ m i n , and frequency (o m i n. The damping ratio is approximately constant over a 3:1 frequency range centering about the comin, as depicted in Fig. 5.2. Frequency independent damping, as used in S H A K E , can be approximated by choosing co m i n to be close to the center of the dominant frequency range of the model. 66 Chapter 5. Verification of Total Stress Procedure u5> 3 2 1 10 15 CO - | 1 1 r 20 25 30 — — — — — Mass Proportional Only - Stiffness Proportional Only total Figure 5. 2 Variation of Normalized Critical Damping Ratio w i th Frequency (FLAC 3.3 manual, 1995) Local damping was originally used to equilibrate static simulations in F L A C . It is attractive to dynamic analysis because unlike Rayleigh damping it is frequency independent. Local damping operates by subtracting and adding mass at a gridpoint at velocity extremes. This is performed at every cycle of calculation. The proportion of energy removed can be related to a fraction of critical damping. The local damping coefficient used in F L A C is equal to the product of pi and critical damping ratio. Local damping, mass-, or stiffness-proportional Rayleigh damping, or a combination of the latter two types of damping can be used n F L A C . Past detailed analyses to study the effect of damping in F L A C , included in the Appendix, have shown that the use of combined mass- and stiffness- proportional Rayleigh damping gives more accurate results. Mass- or stiffness-proportional Rayleigh damping only can give reasonable results if the critical damping ratio is doubled and c o m i r i is centered about the dominant frequency of the system. Local damping gives 67 Chapter 5. Verification of Total Stress Procedure acceptable results in terms of displacement although the resulting acceleration time histories are not comparable to S H A K E . 5.2.2 Case Analyzed A 100-ft (30.48-m) soil column was analyzed using both S H A K E and F L A C . A 30-sublayer or zone column was used in both S H A K E and F L A C . The column is homogeneous and elastic and its stiffness properties are stress dependent. The shear modulus and damping attenuation relationships corresponding to sand and gravel, which are available in S H A K E , were used, as shown in Fig. 5.3. The value of the maximum shear modulus number, k 2 m a x, was 49.1. The maximum damping ratio used in S H A K E was 26 based on the equations developed by Hardin and Drnevich (1972). The water table was assumed to be at ground surface. The soil density was assumed to be 0.129 kips/ft3 or 2066 kg/m 3 . The Caltechb earthquake record scaled to 0.32 g was used in the analysis. The vertical sides in the F L A C column were "attached" to simulate a continuous column, which is assumed in S H A K E . The base of the column is fixed in F L A C but is an elastic half-space in S H A K E . a) Shear Moduli Attenuation Curves b) Damping Attenuation Curves — — — — (3) ROCK and (6) TILL (SHAKE Manual, Dmax=4.6) (7) SAND (D/Dmax=1-G/Gmax, Byrne (1990)) Figure 5. 3 Shear Modul i and Damping Attenuat ion Curves 68 Chapter 5. Verification of Total Stress Procedure 5.2.3 Results The resulting fraction of k 2 max computed in the S H A K E analysis ranged from 0.211 to 0.358 at the top of the column. With the exception of the top of the column, the average G / G m a x ratio is 0.215. The average critical damping ranged from 0.133 to 0.177, with an average critical damping ratio of 0.175, not including the top of the column. F L A C tends to underpredict the response for corresponding input parameters. This is because F L A C damping is a minimum at 4 Hz and all other frequencies will have higher damping. Comparable results are obtained when a fraction of k 2 max of 0.215 and critical damping ratio of about 10 % about a central frequency of 4 Hz are used in the F L A C model. This central frequency corresponds to the dominant frequency of the Caltechb earthquake record. The maximum shear stress ratio, cyclic stress ratio, and acceleration time history calculated for the top of the column using S H A K E and F L A C are compared and shown in Fig. 5.4. Although it has been shown that F L A C results are comparable to those computed by S H A K E , it is recommended that S H A K E analyses be performed to verify the results using F L A C . Analyses have shown that under a similar magnitude earthquake, acceptable results can be computed in F L A C for the fraction of k 2 m a x of 0.215 and critical damping ratio of 10 %. However, it should be noted that S H A K E uses cyclic or pre-liquefaction parameters for the entire column throughout the duration of the input motion. In reality, the response may be significantly different since liquefaction in layers or zones may be triggered at different times. Because F L A C uses a time domain approach, it gives the user the ability to assign post-liquefaction parameters in zones when liquefaction is predicted to occur in a particular zone. The present comparison of the L F A C model with S H A K E aimed to verify the triggering procedure. To validate the post-triggering response of the model, the results can be compared with the empirical equations developed by Bartlett and Youd. 69 Chapter 5. Verification of Total Stress Procedure Comparison Between F L A C and S H A K E Acceleration Histories at Top of Column S o 2 CD -2 H - --4 S H A K E FLAC i1 (.1 • • 10 Cyclic Stress Ratio Profile Cyclic Stress Ratio 0 0.5 1 10 -4 20 H -30 15 I 20 Time (s) 25 Shear Stress Profile Shear Stress (kPa) 0 20 40 60 80 100 10 20 30 30 35 40 Shear Modulus Profile Shear M odulus (kPa) 0 20000 40000 _ ^ | i I 10 20 30 Figure 5. 4 Comparison of Dynamic Response Between FLAC and SHAKE 5.3 Post-Triggering Verification Using Bartlett and Youd Method The earthquake induced deformations predicted by the proposed procedure are compared to those predicted by Bartlett and Youd's (1995) relationship for ground displacements. The equation requires a number of topographical and geotechnical parameters and is based on case histories from a number of earthquakes. 70 Chapter 5. Verification of Total Stress Procedure 5.3.1 Case Analyzed A 15-zone, 15-m high column was analyzed in F L A C . A smaller column than that used in the trigger check was used to reduce computation time. Liquefaction triggering was compared to S H A K E for this column as well. The geometry and soil parameters used in the analysis are within the limits for which the Bartlett and Youd's equation are valid. The water table was assumed to be 1 m below ground surface and 14 m of the column was allowed to liquefy. An S P T (N^o value of 10 was assumed for the entire column. A G / G m a x and damping ratios of 0.22 and 0.08, respectively, were used based on comparisons with S H A K E analyses. The Caltechb earthquake record, scaled to a peak ground acceleration value of 0.32 g, was used. Three values of residual shear strength ratio varying from 0.1 to 0.3 were used in order to assess the effect of the strength on the results. Similarly, the residual or limiting strain was varied, ranging from 0.025 to 0.2. The slope of the column was varied from 0.0° to 5.0°. The corresponding input parameters for Bartlett and Youd's model are cumulative thickness of liquefied layer, T 1 5 , of 14 m and average fines percentage of 0 %. Rather than using peak ground acceleration, the input earthquake parameters for the empirical equation are magnitude and epicentral distance. Using the attenuation relationship developed by Idriss (1991) for an earthquake magnitude of 6.6 and peak ground acceleration of 0.32 g, the computed epicentral distance is about 17.5 km. This is approximately equal to that computed based on attenuation relationships developed by Joyner and Boore (1988), Campbell (1990), and Geomatrix (1992). Like the F L A C analysis, the average mean grain size in the liquefied layers and the slope were varied. 5.3.2 Results The results were plotted in terms of displacement versus slope for residual strength ratios of 0.1, 0.2, and 0.3. The effect of varying limiting shear strains as well as a comparison to the findings of Bartlett and Youd are also plotted and shown in Fig. 5.5. 71 Chapter 5. Verification of Total Stress Procedure For all three residual strength values, the results are comparable to those computed by Bartlett and Youd's relationship. In general, the displacements decrease with increasing residual strength values. For a given value of residual strain, the displacements decrease by approximately 1/3 as the residual strength increases by 0.1. Similarly, for a given residual strength ratio, the displacements increase with increasing residual strain. This is reasonable because increasing the residual strains have the effect of decreasing the post-liquefaction shear modulus or stiffness. For a residual strength value of 0.1 and slope of 2°, the displacements increase by a factor of 3, from over 1 m to over 3 m, as the residual strains increase from 2.5 % to 20 %. Similarly, the displacements increase from about 0.7 m to over 2 m and from approximately 0.5 m to 1.5 m for residual strength values of 0.2 and 0.3 respectively. The results computed for a residual strength value of 2.5% compare most favourable with those predicted by Bartlett and Youd for all residual strength ratios. This is in agreement with the observations of Byrne (1996) and Dobry (1998). In addition, based on Fig. 5.5, the residual strength ratios appear to be associated with average grain size. For example, a residual strength ratio of 0.1 corresponds to a mean grain size of 0.07 mm which infers a fine to silty sand. This is in agreement with the findings of Koester (1992) who found that sands containing fines have lower residual strengths than clean sands. Similarly, sands of mean grain size of about 0.25 mm and 0.5 mm correspond to residual strengths of 0.2 and 0.3 respectively. The model results are comparable to those predicted by Bartlett and Youd for this specific set of input parameters and earthquake history. To ensure that the model yields reliable results under different circumstances, the sensitivity of the model was tested with respect to G / G m a x ratio, critical damping ratio, water level, and most importantly different earthquake records and magnitude. 72 Chapter 5. Verification of Total Stress Procedure b. su/s'vo = 0.2 + Bart let t -Youd, D50 = 0.07mm £ Bart let t -Youd, D50 = 0.5mm Ylim = QJ2-+ -im = 0£5i~ y l i m = 0.025 T 2 3 Slope (degrees) T c. su/s'vo = 0.3 + Bart let t -Youd, D50 = 0.07mm # Bart let t -Youd, D50 = 0.5mm Ylim = 0.2 -¥ l im-= 0.1. T T 0 1 2 3 4 5 Figure 5. 5 Comparison of Displacements Predicted by Proposed Model and Bartlett-Youd 73 Chapter 5. Verification of Total Stress Procedure 5.3.2.1 Effect of G/Gmax ratio The effect of fraction of maximum shear modulus on model prediction was analyzed by varying the G / G m a x and maintaining all other parameters constant. It was assumed that the residual strength ratio is 0.1, limiting shear strain is 0.2, and the ground slope is 0.5°. The water table was maintained at 1 m below ground surface and the (N^o value for the entire column is 10. The G / G m a x was varied from 0.2 to 0.4 to 0.8. Lateral displacements of 1.92 m, 2.01 m, and 2.15 m were predicted for each of the above G / G m a x values. Based on these results, the model prediction is not significantly sensitive to the G / G m a x value, for the Caltechb earthquake record. The effect of increasing G / G m a x is to decrease the time to triggering of liquefaction in specific zones. Hence, the predicted displacements increase slightly. 5.3.2.2 Effect of Critical Damping Ratio In the previous analyses, a critical damping ratio of 0.08 centered about 4 Hz was used throughout the analyses. Despite the fact that an average critical damping ratio of 0.15 was determined from S H A K E analyses, a value of 0.08 was used because it approximates an average of the pre-liquefaction and post-liquefaction viscous damping. As discussed in Chapter 4, once liquefaction is triggered, a portion of damping is accounted for by the stress-strain loops created by the dynamic motion in which a soil element loads and unloads under different stiffnesses. The most desireable scenario is one in which the viscous damping can be reduced once liquefaction is triggered. Unfortunately, this option was not available in F L A C 3.3 since the damping cannot be varied by element or zone. The sensitivity of the model to the critical damping ratio was tested by varying the critical damping values from 0.08 to 0.12 to 0.16 under the same conditions as specified above. The lateral displacements were computed to be 0.92 m, 0.85 m, and 0.63 m for damping ratios of 0.08, 0.12, and 0.16 respectively. Based on these results, the model is moderately sensitive to the chosen damping values. For damping ratios of 0.12 and 0.16, subsequent analyses on a column under steeper slope of 2° showed a difference In results varying from 2.3 m to 1.46 m. This is likely because as the damping increases, the variation in damping with frequency also 74 Chapter 5. Verification of Total Stress Procedure increases. As a result, damping at frequencies other than the central frequency are likely to be considerably larger. Based on the present study, a critical damping ratio of 0.08 is recommended at this time. This value may vary for different earthquake records. 5.3.2.3 Effect of water table The water table was varied to assess the effect of the level of the water table on the model prediction. As the water table was lowered, the estimated displacements decreased. This is consistent with the fact that as the water table was lowered, the number of zones which liquefied was also reduced. 5.3.2.4 Effect of Earthquake Time Histories and Magnitude Including the Caltechb record, the column was analyzed using a total of 6 different earthquake time histories of different magnitudes and peak ground accelerations, PGA. Most of the records were scaled to 0.15 g and 0.27 g and the response of the column to the different records were assessed. The epicentral distances were also computed accordingly using Idriss' method for use in the Bartlett and Youd approach. The main characteristics of the earthquake histories are summarized in Table 5.1. Table 5.1 Summary of Earthquake Histories used in Bartlett-Youd Verification Earthquake Record Mag PGA (g) Vmax (cm/s) Duration (s) Epicentral Distance (km) actual Idriss (est) Caltechb (modified) 6.6 0.32 19.6 40 17.8 Wildlife (shortened) 6.6 0.15 7.2 40 31 36.4 San Fernando (a) 6.6 0.27 30.5 54.4 32.1 21.2 San Fernando (b) 6.6 0.12 17.2 45.7 42.4 43.5 Loma Prieta - Gilroygc 7.1 0.35 28.9 40 28.7 22.8 Loma Prieta - Corralitos 7.1 0.48 47.5 40 6.9 17.2 Identical soil parameters were used for all analyses to assess the effects of using different earthquake base input records. The water table was assumed to be 1 m below ground surface, (1^)60 value of 10 was used, and the critical damping ratio was 0.08. The only difference between 75 Chapter 5. Verification of Total Stress Procedure analyses is the central frequency for Rayleigh Damping. The central frequency was set at the dominant frequency of the earthquake record, which in turn was determined from fast fourier transforms. The residual strength ratio was kept constant at 0.2, while the residual strain and slope were varied to generate plots similar to Fig. 5.5. Predictions computed from Bartlett and Youd's relationship were also plotted for comparison. It should be noted that the comparison may not be a fair one because the pre-liquefaction parameters may not be appropriate. Both the Caltechb and Gilroy records were scaled down to 0.15 g to compare the responses with those of the Wildlife site and San Fernando record (b). In all cases where the earthquake histories are scaled, the distance to fault rupture, using Idriss' relationship, was recalculated for use in Bartlett and Youd's equation. As shown in Fig. 5.6, the response predicted using the Gilroy record is comparable to that predicted by Bartlett and Youd. Although not shown, the response using the Caltechb record is similar. In contrast, the model overpredicts the response with respect to Bartlett and Youd for the Wildlife and San Fernando records. In general, liquefaction was predicted to have occurred at about 5 to 6 seconds. For the Wildlife record in particular, as shown in Fig. 5.7 and 5.8, the peak ground acceleration, PGA, occurred at a later time. As a result, large pulses are applied to the system after liquefaction has been predicted to have been triggered. In addition, the acceleration amplitudes following the occurrence of the P G A are comparable to the PGA. This is not the case for the Caltechb or Gilroy records. For these records, the PGA is reached at a relatively early point in time and, in general, the subsequent acceleration amplitudes are significantly lower than the PGA. The effect of relatively high pulses after liquefaction has been predicted to occur may have contributed to the overprediction of the response with respect to the Bartlett and Youd approach for the Wildlife and San Fernando earthquake records. However, the ability of the model to respond to high acceleration amplitudes after the P G A is the reason for which this approach is prefered to other approaches in which the entire earthquake history is not considered in the calculation of the post-liquefaction deformations. The Gilroy acceleration record was scaled down to 0.27 g and the response was compared to those predicted using the Caltechb and San Fernando (a) records. As shown in Fig. 5.9, the predicted response to the Caltechb and Gilroy records is comparable to Bartlett and 76 Chapter 5. Verification of Total Stress Procedure Youd. In contrast, like the previous example, the San Fernando record tends to overpredict the response. A close examination of Table 5.1 indicates that there is a significant difference between the actual epicentral distance and that computed by Idriss' attenuation relationship. If the epicentral distance computed from Idriss' relationship is used in Bartlett-Youd's equation, then the model response and that predicted by Bartlett and Youd are comparable. This was also observed in the proposed model's response to the Loma Prieta Corralitos record. 77 Chapter 5. Verification of Total Stress Procedure a. Wildlife Record, pga = 0.15g 1 2 3 Slope (degrees) b. San Fernando Record, pga = 0.12 g + Bartlett-Youd, D50 = 0.07mm # Bartlett-Youd, D50 = 0.5mm Ylim = 0.2 2 3 Slope (degrees) c. Gilroy Record, scaled to pga = 0.15 g + Bartlett-Youd, D50 = 0.07mm # Bartlett-Youd, D50 = 0.5mm Ylim = 0.2 Ylim = 0.1 + -"" Ylim = 0.05 " Tii'm = U.U25" 2 3 Slope (degrees) Figure 5. 6 Comparison of Response to Different Earthquake Records Scaled to PGA 0.15g 78 Chapter 5. Verification of Total Stress Procedure a) Wildlife Earthquake, P G A = 0.15 g 20 Time (s) 0.5 0.4 —\ 0.3 0.2 —\ 3 0.1 o 2 0.0 o -0.1 —\ -0.2 — -0.3 — -0.4 — -0.5 b) Cal techb Earthquake, P G A = 0.32 g 10 20 Time (s) 30 40 Figure 5. 7 Acceleration Time Histories of Wildl i fe and Caltechb Earthquakes 79 Chapter 5. Verification of Total Stress Procedure -0.4 —\ 0 10 20 30 40 Time (s) d) San Fernando Earthquake (a), P G A = 0.27 g 10 20 Time (s) 30 ~l 40 Figure 5. 8 Acceleration Time Histories of Gilroy and San Fernando (a) Earthquakes 80 Chapter 5. Verification of Total Stress Procedure a. Caltechb Record, pga = 0.32g - | - Bar t le t t -Youd, D50 = 0 .07mm A Bar t le t t -Youd, DSO = 0 .5mm ,im = 0 , 2 YMrn.: J l im = - 0 . 1 -' .+ Y i i m = 0 . 0 2 5 2 3 Slope (degrees) b. San Fernando Record, pga = 0.27g - | - Bar l l e t l -Youd , D50 = 0 .07mm 4k Sar l l e t t -Youd , D50 = 0.5mm Ylim = 0 . 2 . • ^ T t a — 0 T 4 — . . . s z S ^ m ^ & ^ r ^ -Y , i m = 0 . 0 2 5 2 3 Slope (degrees) c. Gilroy Record, pga = 0.27g - | - Bar t le t t -Youd, D50 = 0 .07mm A Bar t le t t -Youd, D50 = 0 .5mm Ylim = 0 . 2 Ylim = 0.1_ Ylim = 0 . 0 5 Ylim = 0 . 0 2 5 T 2 3 Slope (degrees) Figure 5. 9 Comparison of Response to Different Earthquake Records Scaled to 0.27g 81 Chapter 5. Verification of Total Stress Procedure 5.4 Summary The proposed F L A C model has been verified, with respect the pre-liquefaction and liquefaction triggering procedure, by comparing the results with those obtained by the computer program S H A K E . The model can predict a similar response to S H A K E , although it is recommended that S H A K E analyses be performed and compared to that predicted from the F L A C model to ensure that liquefaction is triggered at the right time if at all. The model was also compared to the predictions using Bartlett and Youd's equations in order to validate the post-liquefaction deformation procedure. Reasonable comparisons of the computed displacements were found. Sensitivity analyses were also performed and it was found that the model is most sensitive to earthquake input record. Comparable results were observed when the actual epicentral distance for a particular earthquake history is the same as that computed from Idriss' relationship and input into the Bartlett and Youd equation. In additon, the analyses indicate that the model is relatively insensitive to the fraction of maximum shear modulus for the cases studied. It is suspected that these parameters could become very significant if the traditional factor of safety against liquefaction is close to 1.0. The fact that the earthquake record plays a significant part in the prediction of the earthquake induced displacements indicate that other simple methods may be deficient under certain circumstances. The proposed total stress approach is compared to the more commonly used methods by application to dynamic analyses of Coquitlam Dam in Chapter 7. Available information on Coquitlam Dam is first summarized in Chapter 6. 82 Chapter 6. Past Analyses of Coquitlam Dam CHAPTER 6 PAST ANALYSES OF COQUITLAM DAM 6.1 Introduction Coquitlam Dam is one of 61 dams at 43 sites in British Columbia maintained by B.C. Hydro. As shown in Fig. 6.1, it is located in the Port Moody Conservation Reserve approximately 20 km east of Vancouver. The dam retains Coquitlam Lake which has a maximum operating volume capacity of 230 x 10 6 m 3 . Water in the reservoir is directed via a tunnel to Buntzen Lake, for power generation. In addition, it is a water supply source for the Greater Vancouver Regional District. The municipalities of Coquitlam and Port Coquitlam are approximately 8 to 15 km immediately downstream of the dam. Coquitlam Dam is a hydraulic fill embankment dam. The existing dam is 30 m high, 290 m long, and 200 m wide at the base. The dam has an upstream composite slope of 1V:2H with berms and downstream composite slopes ranging from 1V:2.5H, 1V:3H, to 1V:5V, and an idealized section is shown in Fig. 6.2. The dam consists of upstream and downstream rockfill berms, sand and gravel shells, and a sandy silt to silty sand core with interlayered sand. Subsurface explorations to 30 m depth below the base of the dam revealed that it is underlain by glacio-fluvial deposits comprising of layers of very stiff silt, dense sand and gravel, and dense sand underlain by bedrock. A bedrock ridge outcrops downstream of the left abutment. The dam was built between 1909 and 1913 by the Vancouver Power Company. It was constructed by building two rockfill toes first and then hydraulically placing sand and gravel against the rockfill. The sandy silt core was sluiced into the centre portion of the dam through two sets of double flumes. Since construction, the dam has been rehabilitated twice. The cross-section shown in Fig. 6.2. is the present one. 83 Chapter 6. Past Analyses of Coquitlam Dam Chapter 6. Past Analyses of Coquitlam Dam The dam was rehabilited as a result of two seismic stability studies. The first study was implemented in 1979 when the Water Management Study recommended that the stability of Coquitlam Dam be reviewed under the then current design earthquake. The results of the 1979 study indicated that the dam must be rehabilited. Subsequent studies revealed that the dam required further reinforcement. This led to the implementation of the 1984 seismic stability study and rehabilitation. In addition to the dynamic analyses, the studies included field and laboratory testing programs. 1 HYDRAULIC FILL CORE 6 RAISED CREST (1955) 2 HYDRAULICALLY PLACED SHELLS 7 SILT FOUNDATION 3 ROCKFILL TOES 3 SAND & GRAVEL STRATUM 4 ROCKFILL & GRANULAR REINFORCEMENT 9 BEDROCK 5 ADDITIONAL REINFORCEMENT (1934-1905) Figure 6. 2 Typical Cross Section of Coquitlam Dam (not to scale) 6.2 Field and Laboratory Testing 6.2.1 Field Exploration Field exploration programs were carried out at Coquitlam Dam in 1979 and 1984 as part of the seismic stability studies. In 1979, the program was carried out to obtain undisturbed samples of the dam material, perform in-situ tests, and install monitoring devices. The program 85 Chapter 6. Past Analyses of Coquitlam Dam consisted of mud rotary drilling, dynamic cone penetration testing, and pressuremeter testing. Two field programs were carried out in 1984 in order to assess the potential slides in the spillway channel and further refine the rehabilitation program. Mud rotary drilling and electric cone penetration tests were performed. 6.2.1.1 Mud Rotary Drilling A total of 14 mud rotary holes were drilled through the dam and/or foundation in order to obtain soil stratigraphy, samples for laboratory testing, and permeability tests. More than 15 instruments were installed in the process. Standard penetration tests were performed in the upstream and downstream shell material. A summary of the drill holes is shown in Table 6.1. Drilling in the hydraulic-filled portion of the dam was said to be easy and rapid indicating that the core material is soft. The core material was found to comprise of very soft sandy silt to silty sand with occasional fine to medium sand seams. Full recovery of undisturbed samples was difficult with an average recovery of only 69 % for 50 thin-walled Shelby tubes. Four holes were drilled through the upstream shell. The upstream shell consists of relatively dense sands and gravels which progressively become finer and less dense toward the hydraulic fill core of the dam. The S P T blow counts corrected to a standard overburden pressure of 1 ton/ft2, decrease toward the core of the dam from an average value of 40 to 29 to finally 13. If a safety hammer is assumed for energy corrections then the average ( N ^ o values would reduce to 36 to 26 to 11 for drill holes DH84-8, DH84-9, and DH84-10 respectively. Using Skempton's relationship (1986), which correlates relative density with blow count, the corresponding relative density of the upstream shell would vary from about 75 to 60 to about 45 to 50 %. 86 Chapter 6. Past Analyses of Coquitlam Dam Table 6.1 Summary of Mud Rotary Drill Holes DRILL HOLE GROUND ELEV. (m) ZONE FIELD TESTING or SAMPLING LABORATORY TESTING DH79-1A 160.10 core Shelby sample index tests DH79-1B 160.10 core Shelby sample resonant column, cyclic triaxial DH 79-2 158.48 core Shelby sample resonant column, cyclic triaxial DH 79-3 160.11 core, Shelby sample resonant column, cyclic triaxial foundation DH84-1 160.11 left abutment DH 84-2 160.06 left abutment DH 84-3 159.89 core Shelby sample, S P T C U triaxial DH 84-4 155.81 core Shelby sample C U triaxial, cyclic triaxial DH 84-5 155.83 upstream S P T shell DH 84-6 133.90 downstream S P T index tests foundation DH 84-7 151.35 downstream S P T index tests shell DH 84-8 154.46 upstream S P T shell DH 84-9 154.44 upstream S P T shell DH 84-10 154.49 upstream S P T shell Only one hole was drilled through the downstream shell. Drilling was difficult and had to be terminated above the foundation level when circulation of drilling mud could not be maintained. This suggests that the gradation of the downstream shell is coarser than that in the upstream shell, thus has higher permeability. The average S P T (N^o value of the downstream shell, discounting extremely high values, is about 36. Based on Skempton's relationship, this corresponds to a relative density of 75 %. Similarly, only one hole was drilled through the downstream foundation. Drilling was also difficult due to the dense nature of the material and continuous mud loss, particularly at the higher elevations. Uncorrected standard penetration values were generally greater than 100. 87 Chapter 6. Past Analyses of Coquitlam Dam 6.2.1.2 In-Situ Testing 6.2.1.2.1 Dynamic Cone Penetration Testing In 1979, six dynamic cone penetration tests were conducted mainly in the core material but included portions of the upstream shell. The tests, carried out by Becker Drilling Ltd., consists of driving a "Becker sleeved cone" into the ground using a hydraulically operated 63.6-lb. trip hammer dropped from a height of 762 mm (30 in.). The "Becker sleeved cone" has a maximum diameter of 50.8 mm and uses a 200-mm long, 50.8 mm diameter sleeve mounted on 38 mm diameter rods. In all cases, the cone was driven to the base of the dam but could not penetrate the stiff silt foundation. The penetration resistances increase with depth, varying from about 2 to 10 blows per foot with an average of approximately 4 blows per foot in the core material. In 1979, it was interpreted to correspond to a C R R value of only 0.05. It should be noted that the resistances were not corrected for depth. The difference in penetration resistance between the core and part of the upstream shell material which is closest to the core is not distinct, with an average value of just less than 10. It is evident that the penetration resistance in the downstream shell material is greater. 6.2.1.2.2 Pressuremeter Testing In 1979, 11 self-boring pressuremeter tests through drillhole DH 79-1B were performed in the core material. The purpose of the testing program was to determine the stiffness characteristics of the hydraulic core material under slow cyclic loading and to examine the dilation characteristics of the core material. Eight of the tests were simple expansion tests and the remaining 3 tests were slow cyclic tests. It was found that the core is locally heterogeneous with zones of compact and loose soils. The behaviour of the soil varied from slightly dilative to slightly contractive. No increase in pore pressure was observed albeit the tests were performed above the phreatic surface. Based on the results in 1979 and on studies by Vaid et al. (1981), the C R R was interpreted to be about 0.11. 88 Chapter 6. Past Analyses of Coquitlam Dam A number of difficulties were encountered during the field testing program. It was difficult to maintain water in the casing which indicates the core material has a high permeability. This is unexpected of silty material. In addition, towards the end of the testing program, the hole collapsed when the pressuremeter was pulled out. 6.2.1.2.3 Electric Cone Penetration Testing Four electric cone penetration tests were performed in the hydraulic fill core material in 1984. The pore water pressures were measured behind the cone tip in holes C P T 84-1 and C P T 84-3. In contrast, the pore pressures were measured on the cone face for holes C P T 84-2 and C P T 84-4. The typical cone penetration resistance profile is shown in Fig. 6.3. Interpretation of the cone shows that the core consists of two types of materials. Material A has a cone tip resistance of about 1 MPa, friction ratio of over 4 % and does not generate significant pore pressure. This is interpreted as silty clay to silty sand. Material B has cone tip resistance of about 2.5 MPa, friction ratio averaging 2 %, and generates pore pressures between 300 to 500 kPa, corresponding to a silty sand type material. Both materials are interlayered with sand to silty sand and are considered susceptible to liquefaction. Based on the empirical charts developed by Stark and Olson (1995), the C R R is about 0.10 for both of these materials. The pore pressure profile indicates a heterogeneous material which exhibits dilatant to contractant behaviour. 89 Chapter 6. Past Analyses of Coquitlam Dam 90 Chapter 6. Past Analyses of Coquitlam Dam 6.2.2 Laboratory Testing A summary of laboratory tests performed on Coquitlam Dam material is shown in Table 6.2. Table 6. 2 Summary of Laboratory Tests TEST MATERIAL YEAR NO. TESTS Index Tests Core 1979, 1984 38 Downstream Shell 1984 7 Downstream Fndn 1979, 1984 9 Resonant Column Core 1980 6 C U Triaxial on undisturbed samples Core 1984 3 C U Triaxial on remoulded samples Core 1984 4 C U Triaxial on liquefied undisturbed samples Core 1984 4 C U Triaxial on liquefied remoulded samples Core 1984 2 Cyclic Triaxial on undisturbed samples Core 1979, 1984 14 6.2.2.1 Index Tests Index tests were performed on hydraulic fill material, downstream shell and foundation material. They were carried out by all the laboratories that tested Coquitlam Dam samples and consisted of moisture content, grain size, Atterberg limits, wet and dry densities, void ratio, and specific gravity. In general, the recovered hydraulic fill core material consists of soft, saturated, uniform, non-plastic, gray silt with some fine to medium sand lenses and trace organics. The material is easily disturbed. The moisture content is about 30 % and the dry and wet densities are about 1500 kg/m 3 and 2000 kg/m 3 respectively. The plasticity index averages about 5 % and the void ratio is 0.95. Tests on the downstream shells showed that it consists of well graded sand and gravel with some to trace silt. The moisture content is about 14 % and the dry and wet densities are 1800 and 2070 kg/m 3 respectively. 91 Chapter 6. Past Analyses of Coquitlam Dam 6.2.2.2 Resonant Column Tests In 1980, 6 resonant column tests were performed on the core material in order to determine the variation of shear modulus and damping with shear strains. The average value of the shear modulus at low strain, G m a x , is about 90000 kPa. The corresponding dimensionless shear modulus number, K 2 m a x . has been calculated to be about 27. 6.2.2.3 Consolidated Undrained Triaxial Test Three and four consolidated undrained (CU) triaxial tests were perfomed on undisturbed and reconstituted core samples respectively. The tests indicate that the silty core exhibits dilative behaviour under compression for both undisturbed and reconstituted samples. However, the reconstituted samples assumes a relatively more contractive behaviour than the undisturbed samples which can be categorized as a limited liquefaction type of behaviour. The strength at the phase transformation state as a fraction of the confining stress is conservatively estimated to be 0.3. This strength is developed at strains of about 1 % to 1.5 % . For undisturbed samples, the material dilates rapidly upon further straining and increases in strength. In contrast, the reconstituted samples dilate slightly upon further straining with a corresponding small increase in strength. Typical test results are shown in Fig. 6.4. As mentioned previously, the recovery of the Shelby tube samples was relatively poor at about 69 %. There is a possibility that the core material sampled correspond to the stiffer layers in the soil. Therefore, it would be reasonable and conservative to assume that the core material behaves more like the weaker reconstituted sample than the undisturbed sample. 92 Chapter 6. Past Analyses of Coquitlam Dam Figure 6 . 4 CU Triaxial Test Results of Core Material 93 / Chapter 6. Past Analyses of Coquitlam Dam 6.2.2.4 Cyclic Triaxial Tests Fourteen cyclic triaxial tests were performed in 1979 and 1984 to determine the cyclic resistance of the hydraulic fill core material. The tests indicate that the core material will liquefy in 15 cycles when the cyclic stress ratio is about 0.15 under cyclic triaxial conditions. Under simple shear conditions, the corresponding cyclic stress ratio is 0.10. It is important to note that liquefaction was considered to have triggered when a double amplitude shear strain of 5 % rather than 100 % pore pressure ratio has been attained. A plot of the cyclic stress ratio versus number of cycles to liquefaction is shown in Fig. 6.5. I 2 3 4 6 6 7 89I0 20 30 40 5060708090KX) 200 300 NUMBER OF CYCLES, N Figure 6. 5 Cyclic Resistance of Core Material of Coquitlam Dam 6.2.2.5 Post Cyclic Monotonic Tests In 1984, monotonic triaxial tests were performed in compression on samples after cyclic loading to determine the strength properties of the core material after it had been liquefied. In general, the material exhibited limited liquefaction behaviour with the remoulded samples exhibiting less dilative behaviour than the undisturbed samples. The tests showed that the samples gained strength to the value before liquefaction. This corresponds to a residual strength ratio of 0.3. However, the strains required to reach the strength was 25 % rather than only 1 %. It 94 Chapter 6. Past Analyses of Coquitlam Dam is important to note that in most cases, the samples were cyclically loaded until double amplitude strains of 5 % have been attained rather than when 100 % pore pressure ratio has been developed. It has been shown that the post-liquefaction strength and stiffness decrease with increasing pore pressure ratio. In addition, as summarized in Chapter 2, recent studies have shown that the liquefied strength of material is strongest in compression loading. In the field, soil elements are subjected to simple shear, and extension loading as well, which have weaker undrained strengths. The typical results are shown in Fig. 6.6. 6.3 Seismic Stability Analyses 6.3.1 1979 Analyses Preliminary seismic analyses were carried out on Coquitlam Dam in 1979/1980 in which the triggering of liquefaction only was assessed. Liquefaction potential was evaluated, in terms of factor of safety, by comparing the cyclic resistance ratio (CRR) with the cyclic stress ratio (CSR) induced by the design earthquake. The assessment was performed using a simplifed procedure as well as a comprehensive one. In the simplified procedure, the liquefaction resistance was obtained empirically from blow count from Becker Cone Penetration testing, and was equal to 0.05. The cyclic stress ratio was obtained empirically from the equation (Eq. 4.2) developed by Seed and Idriss. The peak ground acceleration (PGA) was assumed to be 0.12 g, corresponding to the Maximum Credible Earthquake (MCE) evaluated in 1979/1980. The results of the analysis indicate that the upstream shell and upper portion of the core will liquefy. The computed cyclic stress ratio contours on the 1979 Coquitlam Dam cross section is shown in Fig. 6.7. 95 Chapter 6. Past Analyses of Coquitlam Dam 8 12 16 A X I A L STRAIN % CE iii 8 12 16 A X I A L STRAIN % 200 U l cn z> cc Q. 2 o X -400 -600 900 600 TEST ru \^ -TEST TL TESTS , rL3. TL4,1 UP a TL 2R 300 600 900 1200 1500 1800 P ( k P a ) Figure 6. 6 CU Triaxial Test Results of Liquefied Core Material 96 Chapter 6. Past Analyses of Coquitlam Dam In the more comprehensive procedure, the liquefaction resistance was obtained from laboratory testing. A s discussed earlier, the C R R was determined to be 0.10. The cyclic stress ratio was determined from a two-dimensional analysis of the deepest cross-section of the dam. Static stresses were determined using a finite element program called ISBILD, in which the construction sequence of the dam was simulated. The dynamic stresses were computed using the program QUAD4 which consideres both variable stiffness and damping. The San Fernando record scaled to a peak acceleration of 0.12 g was used as input base motion. The results indicate that the upper upstream shell and core material would liquefy under the M C E . The cyclic stress ratio profile is shown in Fig. 6.8. Because liquefaction may cause slumping in the upper part of the dam, rockfill berms were added to the slopes. The rehabilitation was carried out in 1980. SIMPLIFIED SEED - HAND CALCULATIONS Figure 6. 7 1979 Computed Cyclic Stress Ratio Contours Based on Hand Calculations 97 Chapter 6. Past Analyses of Coquitlam Dam QUAP 4 Note: Liquefaction predicted to have triggered in shaded zones Figure 6. 8 1979 Computed Cyclic Stress Ratio Contours Based on Dynamic Analyses 6.3.2 1984 Analyses Between 1980 and 1984, further studies indicated that the dam needed to be further reinforced. In result, in 1984, the liquefaction potential of proposed remediated dam sections were evaluated under the operating base earthquake, O B E , which has a P G A of 0.20 g and the M C E which has a P G A of 0.35 g. In addition, the flow slide and deformation potential considering the loss in stiffness when the core liquefies were assessed. The simplified approach was used in the triggering analysis. The cyclic resistance of the core material and upstream shell material were obtained using data from cyclic triaxial tests and in-situ testing. The C R R of the upstream shell was interpreted to be 0.285, based on S P T blow count information corrected from overburden stress. As shown in Fig. 6.9, the results indicated that under the M C E , all of the core material and a portion of the upstream shell will liquefy. 98 Chapter 6. Past Analyses of Coquitlam Dam Uote: Liquefaction predicted in shaded zone Figure 6. 9 1984 Computed Cyclic Stress Contours Based on Hand Calculations A limit equilibrium analysis was conducted to assess the flow slide potential of the dam. A number of failure surfaces were generated and a factor of safety against flow slide was computed. Friction angles of 2° and 20° for the liquefied core and upstream shell respectively were assumed. These strengths were considered to be equivalent to using full effective strength angle of 35° and pore pressure ratios, ru, of 95% in the core and 50 % in the shell. The results indicated that a flow slide above the upstream rockfill berm may occur under the M C E . However, a flow slide was not predicted to occur under the O B E . Earthquake deformations due to inertia force and to loss in stiffness when the core liquefies were analyzed using Newmark analysis and finite element post-liquefaction stress analysis respectively. The yield acceleration under the O B E , determined from limit equilibrium analysis, was computed to be 0.14g and 0.17g for the upstream and downstream slopes respectively. The displacements due to inertia force is predicted to about 7 cm in the upstream direction and 1 cm in the downstream direction. For the M C E , the critical acceleration is zero and therefore a flow slide is predicted. Under the proposed remediated design, it was considered that the top of the dam would unlikely drop below El. 155.5 due to the geometry of the dam. 99 Chapter 6. Past Analyses of Coquitlam Dam Finite element post earthquake static analyses were performed using the computer program S O I L S T R E S S . The post-cyclic stiffness parameters, obtained from laboratory tests, were reduced to 1/25 of the pre-cyclic parameters. The results indicate that under the O B E , the dam crest may settle vertically about 0.6 m and the upper portion of the downstream slope may move outward 0.3 m. Under the M C E , only an additional crest settlement of 0.3m was predicted. The dam deformation pattern is shown in Fig. 6.10. Figure 6.10 SOILSTRESS Deformation Analysis Based on Loss in Stiffness In addition, the deformation due to dissipation of excess pore water pressures was predicted. Based on a bulk modulus of compressiblity of 4 x 1 0 s m 2 /kN and an estimated average excess pore pressure generated during liquefaction of 300 kPa, the volumetric strains were computed to be about 1.2 %. Assuming that volumetric strains are due to vertical strains only, the crest settlement was conservatively estimated to be 0.4 m. The total estimated displacements of the reinforced dam was computed to be the sum of the displacements due to inertia force, loss in stiffness, and dissipation of pore pressure. Under the O B E , the total vertical displacement was calculated to be 1 m. A vertical displacement of 6.5 m, corresponding to a loss in freeboard, was estimated under the M C E . Because the normal operating reservoir level is less than the level of the rockfill berm, this amount of deformation was considered acceptable. In 1986, the dam was rehabilitated to the section shown in Figure 6.9. 100 Chapter 6. Past Analyses of Coquitlam Dam 6.4 Other Studies Dam breach analyses have been performed for Coquitlam Dam. The studies showed that should the unlikely failure of the dam occur under any circumstance, very little warning time is available to the downstream residents and the number of people potentially at risk is estimated to be about 32,000. 6.5 Summary Coquitlam Dam has been rehabilited twice since construction as a result of seismic stability analyses in 1979 and 1984. The 1979 seismic analyses comprised of a liquefaction triggering assessment only. The core and upper portion of the upstream shell were predicted to liquefy under the 1979 M C E . Rockfill berms were added to the slopes due to the possibility of slumping of the upper dam. Subsequent analyses in 1984 indicated that additional rehabilitation was necessary. Liquefaction triggering, flow slide, and deformation analyses were performed on the proposed remediated section. Under the M C E , which had P G A of 0.35 g, the core and a portion of the upstream shell were predicted to liquefy, resulting in the possibility of a flow slide. However, it was considered unlikely that overtopping would occur under the normal operating reservoir level if the dam were remediated to the proposed geometry. Consequently, the dam was remediated to the proposed geometry. Since 1984, the state-of-practice in earthquake induced deformation analysis has progressed. Although failure by overtopping under the normal operating reservoir level was not predicted in 1984, the most current analysis procedures were applied to ensure that the seismic performance of the dam is still acceptable. 101 Chapter 7. Total Stress Analyses of Coquitlam Dam CHAPTER 7 TOTAL STRESS DYNAMIC ANALYSES OF COQUITLAM DAM 7.1 Introduction Since the last (1984) seismic analysis of Coquitlam Dam, the state-of-practice procedure has changed due to a better understanding of materials under cyclic loading conditions. Although the earthquake-induced deformations that were estimated previously accounted for the effects of inertia forces and loss in stiffness, better procedures are now available. As discussed in Chapter 3, the previously calculated deformation due to the inertia forces is only one-dimensional and does not account for the deformation due to two-dimensional deformation of the structure. For a dam with high consequences of failure, the seismic stability should be re-evaluated using the most current procedures. Three different total stress deformation procedures were performed. The proposed method is a total stress dynamic analysis. The other two methods include the calculations of the deformations due to loss in stiffness only and those due both to loss in stiffness and earthquake velocity pulse, similar to Jitno's method. The results obtained by these three methods were compared. Static and seepage analyses, in which construction sequence and reservoir filling were simulated, were also performed to obtain the pre-earthquake initial stresses. All static and dynamic analyses of the dam were performed using F L A C version 3.3. 7.2 Pre-Earthquake Analysis Prior to performing the dynamic analyses on the dam, the pre-earthquake initial stresses are calculated. This involved simulating construction of the dam and performing an analysis due to the loading of the reservoir water. The pore water pressures in the dam and hence the effective stresses can be calculated. 102 Chapter 7. Total Stress Analyses of Coquitlam Dam 7.2.1 Construction Analysis A two-dimensional plane strain static analysis was performed on the idealized cross-section of the Coquitlam Dam using the computer code F L A C . The computer program allows multiple stage loading in which the construction sequence can be simulated by adding the soil materials in layers. Equilibrium is obtained by "stepping" through the calculation sequence prior to the addition of a new layer. Once equilibrium has been established, the gridpoint displacements of the top layer are set to zero, since displacements are incurred during the stepping process. The cross section of the grid used for the analysis is shown in Fig. 7.1. This cross-section is the same as that used in previous seismic analyses and is considered to represent the deepest cross section of the dam. The grid is composed of about 3610 zones. The base of the cross-section is fixed in the horizontal and vertical directions and the upstream and downstream end boundaries are fixed in the horizontal direction only. The dam was constructed in 10 layers each of over 2 m thick. Although it is likely that pore pressures were developed during construction due to the hydraulic fill placement method, the analysis was performed using drained stiffness parameters. It is judged that by now the pore pressures would have had dissipated. Because the stresses of interest are the "current" stresses and not those representing the stresses immediately after construction, this procedure is considered acceptable. 7.2.1.1 Model Parameters The Mohr Coulomb model, which is a built-in feature of the F L A C computer program, was used in the analysis. This model was modified so that the shear and bulk modulus parameters are stress dependent. The main input parameters are the shear modulus number, kg, bulk modulus number, kb, shear modulus exponent, n, bulk modulus exponent, m, friction angle, <(>, dilation angle, u, and cohesion, c. The parameters used to represent the dam materials were based on results of laboratory testing, in-situ testing, and published values available from literature. As summarized in Chapter 6, the dam and its foundation are made up of 6 different materials. The dam is comprised of rockfill toes, rockfill berms, upstream and downstream sand 103 Chapter 7. Total Stress Analyses of Coquitlam Dam and gravel shells, and a hydraulic fill sandy silt to silty sand core. The dam is underlain by a relatively thin stiff silt layer which overlies a dense sand layer. In these analyses, 8 sets of parameters were assigned to represent the dam materials. The upstream sand and gravel shells were further subdivided into 3 different zones based on the in-situ test results. The relative density of the upstream shell material decreases as it approaches the hydraulic fill core. The different material zones are shown in Fig. 7.2. The soil parameters for the hydraulic fill core were determined based on laboratory test results. The bulk modulus value was computed directly from laboratory test results. Numerical single element tests were performed in which the shear modulus, friction angle, and dilation angle were varied to produce a best fit with the laboratory results. The parameters obtained were compared with the recommended values proposed by Stark et al. (1994) to assess whether or not they are reasonable. The upstream and downstream shell parameters were based on empirical correlation with relative density values developed by Byrne et al. (1987) for sands. Using the relative density which was estimated from in-situ test results, parameters are selected for the hyperbolic model as suggested by Byrne et al. (1987). The parameters were converted to correspond to the Mohr-Coulomb model by simulating a single element test using both the Mohr-Coulomb and the hyperbolic models. Best-fit Mohr-Coulomb parameters were determined by trial and error. A plot comparing the Mohr-Coulomb response with the hyperbolic response is shown in Fig. 7.3. Similar plots were generated for each of the sand and gravel and dense sand foundation material properties. Properties for the stiff silt were determined in a similar manner from a data base of silty material developed by Stark et al. (1994). The upstream and downstream rockfill parameters were obtained from values recommended by Saboya and Byrne (1991) based on laboratory testing and back-analysis of rockfill dams. The Mohr-Coulomb parameters used for static analyses are shown in Table 7.1. 104 Chapter 7. Total Stress Analyses of Coquitlam Dam Figure 7.1 FLAC Grid Used For the Analysis of Coquitlam Dam (a) Dense Sand Foundation (b) Stiff Silt Foundation (c) Upstream and Downstream Rockfill (d) Upstream Sand and Gravel Shell (a) (e) Upstream Sand and Gravel Shell (b) (f) Upstream Sand and Gravel Shell (c) (g) Downtream Sand and Gravel Shell (h) Hydraulic Fill Core Figure 7. 2 Material Zones Used For the Analysis of Coquitlam Dam 105 Chapter 7. Total Stress Analyses of Coquitlam Dam Table 7.1 Summary of Dam Material Parameters used in Static Analysis Dam Material kg n kb m (j) de u de Pdry g g kg/m3 Dense Sand Foundation 425 0.5 650 0.25 60 1800 Stiff Silt Foundation 80 1.0 115 1.0 45 1750 Upstream and Downstream Rockfill 125 0.25 250 0.35 48 2000 Upstream Sand and Gravel Shell (a) 350 0.5 450 0.25 47 1800 (b) 200 0.5 400 0.25 43 1800 (c) 125 0.5 300 0.25 38 1800 Downstream Sand and Gravel Shell 350 0.5 450 0.25 47 1800 Hydraulic Fill Sandy Silt Core 33 0.8 45 0.8 37 8 1500 400 300 GO CO £ 200 co 100 Mohr-Coulomb Hyperbolic ! , ! , ! , ! o . o 2 0.04 0.06 0.08 Strain Figure 7. 3 Mohr-Coulomb Parameter Fit to Hyperbolic Parameter 0.1 7.2.1.2 Results The results of the construction were plotted in terms of displacement vectors and vertical or principal stress contours. The displacement vectors, shown in Fig. 7.4, are oriented almost exclusively vertically downward with some horizontal components along the upstream and downstream boundaries of the dam. The largest magnitude vectors appear to be in the middle of the dam. A plot of the vertical displacement profile along the centerline of the dam, depicted in Fig. 7.5., confirms that the largest vertical displacements occur in the middle of the dam. These 106 Chapter 7. Total Stress Analyses of Coquitlam Dam observations are in agreement with both case histories and expected for a dam with a relatively stiff foundation. The stiffness of dam material decreases with increase in dam height but is countered by the decrease in overburden stress with increase in dam height so that the maximum deformation occurs at about midheight of the dam. The vertical stress contours shown in Fig. 7.6 show low stresses in the core and slightly higher stresses in the stiffer shells, indicative of arching. This is not unreasonable because the difference in stiffness between the shell and core material is significant. Although the 1979 static analysis indicated that there are no areas of arching it should be noted that the relative stiffness between the core and upstream and downstream shells assumed in that analysis was comparable. However, subsequent laboratory and in-situ test results indicate that the core material is very soft. 7.2.2 End of Reservoir Filling Analysis The end of reservoir filling state can be modeled by performing either a coupled or an uncoupled stress-flow analysis. A coupled stress-flow analysis is performed by solving for stresses, displacements, and pore pressures together. An uncoupled analysis is performed by solving the flow separately from the stress to obtain steady state pore pressures or seepage forces. Having computed the pore pressures, the resulting stresses and deformations can be computed from two methods as follows: Method A 1. Specify total unit weights - saturated where appropriate. 2. Specify boundary water pressures. 3. Specify internal pore water pressures. Method B 1. Specify submerged unit weights below the phreatic surface. 2. Specify seepage forces below the phreatic surface. Method A was selected for Coquitlam Dam. 107 Chapter 7. Total Stress Analyses of Coquitlam Dam 150 cr o CO > CD 100-1 J I I !_ _1 I , I L _ I I I I I l_ 0 50 100 150 200 250 Figure 7. 4 End of Construction Displacement Vectors 300 _i i u 350 150-o ca > JS L U 100 0 0.4 J I . I L 0 50 100 150 200 250 300 350 Figure 7. 5 Vertical Displacements Due to Construction Along Centreline of Dam in metres 03 > CD 1 50— 100-I I I I , I I I J , I 1 1 L_ 50 100 150 200 250 300 350 Figure 7. 6 E n d of Construct ion Vertical Stress Contours in kPa 108 Chapter 7. Total Stress Analyses of Coquitlam Dam Both stress and seepage analyses were performed using F L A C . The boundary water pressure under full reservoir elevation of 155.0 m was applied on the upstream boundary. The pore pressures and saturation were fixed for the upstream boundary and the pore pressures only were fixed on the top and downstream boundary of the dam. The steady state pore pressures are computed by performing a seepage analysis in which the mechanical computation option in F L A C is turned off and the flow computation option is turned on. The analysis was carried out until flow equilibrium is reached in which the inflow is approximately equal to the outflow. Having reached flow equilibrium, the flow computation option in F L A C is turned off and the mechanical computation option is turned on. When the flow computation option is turned off, F L A C only subtracts porewater pressures from the total stresses to obtain effective stresses and does not consider porewater pressures as loads. To overcome this problem, the pore pressures were added to the total stresses for at least one time step. The boundary water pressures were applied incrementally to avoid "shocking" the system. Computation was performed until force equilibrium of the system was reached. 7.2.2.1 Flow Model Parameters The parameters required in seepage analysis in F L A C are the permeability and void ratio or porosity. The seepage parameters for the dam materials were based on laboratory testing and published typical values (Craig, 1992). The vertical permeability for the hydraulic fill core was determined from the consolidation portion of triaxial testing. The horizontal permeability was based on Hazen's empirical formula (Craig, 1992) relating permeability with effective grain size as follows: k = lO'2-Df0 Eq. 7. 1 where k = permeability in m/s D-io = effective grain size in mm. The permeability for the upstream and downstream sand and gravel were also determined using Hazen's formula. The permeability parameters for the remaining materials were 109 Chapter 7. Total Stress Analyses of Coquitlam Dam based on typical values. Table 7.2 summarizes the seepage properties for Coquitlam Dam materials. Table 7. 2 Summary of Dam Material Properties for Seepage Analysis Dam Material k porosity (m/s) Dense Sand Foundation S e 3 0.3 Stiff Silt Foundation 5 e 6 0.4 Upstream and Downstream Rockfill 5 e 2 0.43 Upstream Sand and Gravel 10"4 0.35 Downstream Sand and Gravel 10"3 0.3 Hydraulic Fill Core k x - 5e"5 0.39 k v - 1.5e'6 Since the permeability constant, K, required in F L A C is the proportionality constant in Darcy's Law and is expressed in terms of pressure rather than head, the above permeability values were divided by the unit weight of water for input into F L A C . 7.2.2.2 Results The results of the seepage analysis were plotted in terms of flow vectors, pore pressure and total head contours and are shown in Figs. 7.7, 7.8, and 7.9 respectively. As shown in Fig. 7.7, the largest flow vectors are in the downstream rockfill toe. Although the inflow and outflow are equal, the output flow is concentrated over a smaller area than the input flow, hence the outflow velocity vectors are greater. Both the pore pressure and total head contours appear to reasonable. As expected, the pore pressures increase with depth below the phreatic surface. The total head contours are similar to those obtained from SEEPA/V analyses performed in 1995. The results of the end of reservoir filling stress analysis, which is intended to account for the boundary loads, are shown in Fig. 7.10 and 7.11 in terms of vertical effective stress contours and displacements. The vertical effective stresses appear to be reasonable. The displacement vectors, shown in Fig. 7.10, indicate that the crest moves upwards and downstream and are reasonable. However, the magnitude of the displacements are likely too large. 110 Chapter 7. Total Stress Analyses of Coquitlam Dam Chapter 7. Total Stress Analyses of Coquitlam Dam This is due to the assumption that the water level rises from the bottom of the grid to the reservoir water level. A more accurate assumption would be to model the water level to rise from Elevation 150 m to the reservoir water level. Also, the reload moduli should be used in order to be compativle with the stress path of the type of soil element. However, the stresses remain more or less the same if the unloading moduli are used. 200-> 150 100-> 4 1 | 4111, i T 5 0 100 150 200 250 300 Figure 7.10 End of Reservoir Filling Displacement Vectors 350 150-ro > UJ 100-50 100 1000 150 200 250 300 350 Figure 7.11 End of Reservoir Filling Vertical Stress Contours in kPa 112 Chapter 7. Total Stress Analyses of Coquitlam Dam 7.3 Total Stress Dynamic Analysis 7.3.1 Introduction Three total stress deformation analyses were carried out, as follows: 1. Gravity Only 2. Gravity Plus Velocity Pulse 3. Proposed Total Stress Dynamic Analysis (Self Triggering and Gravity Plus Base Acceleration) The first two methods are considered to be traditional procedures in which liquefaction triggering, flow slide, and deformation analyses are performed separately. The proposed method considers liquefaction triggering and the liquefaction induced displacements in the same process and will be presented in the next section. The traditional methods were performed to compare with the proposed method. The first deformation analysis procedure (Gravity Only) is similar to Lee's Modified Modulus approach in that the deformations are caused by gravity loading only as a result of the loss of strength and stiffness in the liquefied soil. The liquefaction effect, in which the liquefied soil behaves like a liquid, is handled by adjusting the stress distribution so that the horizontal stresses in the liquefied zones are equal to the vertical stresses and the shear stresses are equal to zero. In addition, the liquefied zones are assigned reduced strength and stiffness parameters. This stress state and reduced stiffness result in an imbalance in forces within the structure and cause deformation to occur until equilibrium is achieved by stress redistribution. Although the inertia forces due to seismic loading are not considered, analyses using the dynamic option in F L A C allow consideration of inertia forces induced by the initial acceleration of the soil elements under the out of balance forces. The second deformation analysis (Gravity Plus Velocity Pulse) is comparable to the extended Newmark approach proposed by Jitno and Byrne (1995). The deformations are caused by gravity loading and a velocity pulse equal to the maximum ground velocity pulse. Only one 113 Chapter 7. Total Stress Analyses of Coquitlam Dam velocity pulse is considered because the deformations prior to triggering are negligible. The earthquake pulses which occur after liquefaction are accounted for using the residual or limiting strain values which have been correlated to factor of safety against liquefaction. The post-liquefaction parameters suggested by Byrne (1991) was shown in Table 4.1. In the above method the entire dam was considered to have a velocity at the time of liquefaction equal to the maximum ground velocity at the time of liquefaction. The liquefied zones were then assigned their reduced stiffness and strength parameters and the initial stresses were redistributed. All liquefied elements are assumed to liquefy at the same time. 7.3.2 Liquefaction Triggering Analysis The liquefaction potential was assessed by comparing the cyclic resistance ratio with the developed cyclic stress ratio. The cyclic resistance ratio, C R R , for the various dam materials were evaluated from laboratory test results and various in-situ test results such as cone penetration tip resistance and standard penetration blow count. The C R R was corrected for overburden stress using the chart shown in Fig. 2.4. The corrected cyclic resistance ratios assigned for the dam materials are summarized in Table 7.3. Table 7. 3 Corrected Cyclic Resistance Ratio for Coquit lam Dam Materials Dam Material CRRcORR Source Dense Sand Foundation 0.5 S P T blow count Stiff Silt Foundation 0.5 S P T blow count Upstream and Downstream Rockfill will not liquefy (free draining) Upstream Sand and Gravel Shell (a) 0.4 (avg) (b) 0.3 (avg) (c) 0.15 (avg) Downstream Sand and Gravel Shell will not liquefy (not saturated) Hydraulic Fill Sandy Silt Core 0.1 C P T tip, cyclic triaxial tests The cyclic stress ratios induced by the earthquake base motion were computed using a one-dimensional analysis using the computer code S H A K E . Because the Coquitlam Dam slopes are relatively flat and wide, boundary effects are judged to be minimal. In addition, as mentioned 114 Chapter 7. Total Stress Analyses of Coquitlam Dam in the preceding chapters, it has been shown that no significant difference in results were observed between one- and two-dimensional analyses. Several soil columns representing the section between the upstream to the downstream sand and gravel shells were analyzed in S H A K E . Columns representing the rockfill zones of the dam were not analyzed because the rockfill berms are not considered to liquefy. The soil parameters used in the analysis are the unit weights, maximum shear modulus number and the maximum damping ratio in percent. The maximum damping ratios for both the sand and gravel shell and core material were evaluated from relationships developed by Hardin and Drnevich (1972). The maximum shear modulus numbers for the sand and gravel shell and hydraulic fill core were determined from correlations with S P T blow count and resonant column tests respectively. The shear modulus number, (k2)m a x, is determined based on the S P T blow count value through the following relationship, developed by Tokimatsu and Seed (1987): ( * 2 L x = 2 0 - ( A r , ) J Eq.7.2 The S H A K E input parameters are summarized in Table 7.4 Table 7. 4 SHAKE input parameters Dam Material y s a t (kips/ft3)1 (k2)max D m a x (%) Upstream Sand and Gravel Shell (a) 0.129 (2.067) 65 26.2 (b) 0.129 (2.067) 55 26.2 (c) 0.129 (2.067) 45 26.2 Downstream Sand and Gravel Shell 0.129 (2.067) 65 31.2 Hydraulic Fill Core 0.1217 (1.95) 27 20 1. Note: Values in parentheses in metric units of kg/m x 1000. The base input motion used in the S H A K E analysis corresponds to the Caltechb earthquake record scaled to 0.32 g. As shown in Fig.7.12, for an earthquake with pga of 0.32 g, almost the entire hydraulic fill core and a portion of the upstream sand and gravel shell will liquefy. The factor of safety against liquefaction triggering in the core material is relatively low ranging from 0.5 to almost 1.0. 115 Chapter 7. Total Stress Analyses of Coquitlam Dam 7.3.3 Flow Slide Analysis The rise in porewater pressure associated with liquefaction causes a loss in strength and stiffness in the liquefiable zones. Although the dam slopes are stable under static conditions, the development of residual strength of the weakened liquefied core and portion of the upstream shell material may lead to the occurrence of a flow slide. The flow slide potential of Coquitlam Dam was assessed by using a limit equilibrium analysis computer program X S T A B L version 5.1. As presented in the preceding section, nearly the entire dam core and a portion of the upstream sand and gravel shell will liquefy under the Caltechb earthquake scaled to a P G A of 0.32g. Undrained strength ratios of 0.1 and 0.2 were assigned to the liquefied core and upstream sand and gravel shells respectively. Although laboratory test results indicate that the residual strength ratio of the core material is about 0.3, a value of 0.1 was assigned to the core material because the laboratory samples were tested in compression, which tends to give a higher residual strength. In addition, the tested samples were not liquefied to 100% pore pressure ratio, which also results in higher strengths than cases where sample are liquefied to 100 % pore pressure. Back-calculated values from cases histories of flow failure of silty sand to sandy silt materials have been shown to range from 0.1 to 0.2 (Ishihara, 1993). The strengths assigned to the liquefied zones are considered to be conservative. The results, depicted in Fig. 7.13, show that a flow slide is not likely to occur. The factor of safety in the upstream slope is 1.49, and 1.87 in the downstream slope. The difference in result from the 1984 analyses, which predicted that a flow slide in the upstream direction is likely to occur is due to the strength values used. In 1984, friction angles of 2° and 20° were assumed for the liquefied core and upstream sand and gravel shell. To produce the same factors of safety in X S T A B L , these friction angles correspond to undrained residual strengths of 5.7° and 11.3°, respectively. 116 Chapter 7. Total Stress Analyses of Coquitlam Dam 7.3.4 Deformation Analysis 7.3.4.1 Model Parameters Although a flow slide is not predicted to occur, the deformations associated with liquefaction may still be significant and were assessed. Deformation analyses performed presently should consider the effects of liquefaction as well as the inertia forces induced by the base motion. Two commonly used deformation analyses methods were used. The post-liquefaction stress strain curve for both analyses was assumed to be bilinear. Similar to the limit equilibrium analysis, the residual strength ratio was assumed to be 0.1 and 0.2 for the liquefied core and the upstream sand and gravel materials, respectively. The residual strain values were based on Byrne's recommended values as shown in Table 4.1. The limiting strain values were assumed to be 15% for the upstream shell material and 25% for the core material. A maximum ground velocity of 19.6 cm/s was used for the second deformation procedure. This velocity corresponds to the Caltechb earthquake record scaled to 0.32 g. 117 Chapter 7. Total Stress Analyses of Coquitlam Dam Note: Factor of safety decreases with lighter fill tone I i i 1 1 1 1 ' 1 ' 1 ' , 100 120 140 160 180 200 220 240 260 280 Figure 7.12 Liquefaction Triggering of Coquitlam Dam Post-Liquefaction Limit Equilibrium 10 most critical sur faces, Minimum Bishop F O S = 1.867 170 r 145 h 120 I 1 1 1 1 1 1 ' 1 ' 1 ' 1 ' 1 ' 1 8p 105 130 155 180 205 230 255 280 Post-Liquefaction Limit Equilibrium 10 most critical sur faces, Minimum Bishop FOS = 1.486 1 7 0 L 1 4 5 h 1 2 0 I ' 1 ' 1 ' ' ' 1 ' 1 ' 1 ' 1 1 1 8 0 1 0 5 1 3 0 1 5 5 1 8 0 2 0 5 2 3 0 2 5 5 2 8 0 Figure 7.13 Limit Equilibrium Flow Slide Analyses 118 Chapter 7. Total Stress Analyses of Coquitlam Dam The analyses were performed using the dynamic mode in F L A C in order to account for the inertia forces due to the motion of the deforming mass. Rayleigh damping of about 10% about a central frequency of 4 Hz was used. The frequency of 4 Hz corresponds approximately to the predominant undamped frequency, based on analysis of the acceleration history of the crest of the dam, 7.3.4.2 Results The velocity pulse was applied in both the upstream and downstream horizontal directions as well as in a downward vertical direction to assess the effect of direction of velocity pulse on the response of the structure. The results for both approaches are compared in Table 7.5. Table 7. 5 Results of Liquefaction Induced Displacment of Coquitlam Dam Crest Displacement (m) PROCEDURE Horizontal (+ve d/s) Vertical (+ve up) Gravity Only -0.0836 -1.09 Gravity Plus Horizontal Velocity Pulse (d/s) -0.0103 -1.102 Gravity Plus Horizontal Velocity Pulse (u/s) -0.0341 -1.1075 Gravity Plus Vertical Velocity Pulse (down) -0.0348 -1.075 As shown in the above table, the displacements due to inertial loads represented by single velocity pulses result in an increase in crest deformation by about 0.1 m larger than those due the loss in stiffness only. In fact, for the case where a velocity pulse is applied in a vertical downward direction, the predicted displacements are actually less than those computed based on gravity loading or loss in stiffness only. However, the difference is less than 0.015 m can be considered more or less negligible. The computed vertical crest displacement due exclusively to a loss in stiffness of 1.09 m is comparable to that found in 1984 of about 0.9 m under the M C E . The deformation vector pattern, shown in Fig. 7.14, is also similar. It is interesting to note that the patterns in the 119 Chapter 7. Total Stress Analyses of Coquitlam Dam upstream and downstream shells are similar in shape to the failure surfaces predicted by X S T A B L , shown earlier in Fig. 7.10. 0 50 100 150 200 250 300 Figure 7.14 Typical Post-Liquefaction Displacement Vector Pattern In contrast, the displacements due to inertia force cannot be compared to those found in 1984. In 1984, the displacements due to inertia force were computed for a one-dimensional potential failure surface as opposed to a two-dimensional system as performed in the present analyses. In addition, the peak horizontal velocities for both the O B E and M C E of 0.36 m/s and 0.40 m/s respectively, estimated in 1984, are greater than that assumed in the present analyses which is only 0.19 m/s. Despite this fact, the computed displacements are slightly greater than those predicted in 1984 under the O B E of about 0.06 m. Displacements of over 6.5 m, corresponding to a flow slide, were predicted under the M C E in 1984. Based on these results, a velocity pulse of 0.19 m/s does not appear to affect the predicted displacements significantly. A plot comparing the horizontal velocity history of the crest of the dam under the cases of gravity loading only and gravity loading plus horizontal downstream velocity pulse shows that horizontal velocities calculated using the first method (gravity only) are slightly greater than that calculated using the second method (Fig. 7.15). If a velocity pulse of about 0.3 m/s were applied to the dam, it is likely that a greater response will be observed. However, because the Caltechb input motion corresponding to P G A of 0.32 g was chosen, a 120 Chapter 7. Total Stress Analyses of Coquitlam Dam velocity pulse corresponding to the peak horizontal velocity should be used for comparative purposes. 0.40 - i 0.20 o CD > 0.00 H -0.20 L o s s in St i f fness P lus Inertia Force (Hor izonta l Veloc i ty Pu lse) — — L o s s in St i f fness Only 1 I 1 I 1 I 1 I 1 I 1 I 115000 120000 125000 130000 135000 140000 145000 N o . S t e p s Figure 7.15 Comparison of Velocity History of Crest of Dam 7.4 Proposed Total Stress Dynamic Analysis Procedure 7.4.1 Model Parameters The same material parameters used in the above triggering analysis were used in the pre-liquefaction of the dam. The shear modulus values were based on S P T blow count values where available. A Fish subroutine function, written in the F L A C data file, can be used to compute the shear modulus values based on S P T blow count input. The shear modulus values assigned for the rockfill material was slightly stiffer than those computed for the dense portion of the sand and gravel shells. Comparison with the results of the S H A K E analyses showed that the average fraction of maximum shear modulus, G / G m a x , should be about 0.3 for all zones in the dam. Post-liquefaction parameters for all material zones in the dam were assigned in the input data file. All materials were assigned a residual strength ratio of 0.3 with the exception of the 121 Chapter 7. Total Stress Analyses of Coquitlam Dam softer portion of the upstream sand and gravel shell and the hydraulic fill core material. Residual strength ratios of 0.1 and 0.2 were assumed for the hydraulic fill and soft portion of the upstream shell material respectively. A liquefaction strain value of 2.5 %, which was in agreement with the findings of the verification analyses, was used for all dam materials. A summary of the input parameters can be found in Table 7.6. Table 7 . 6 Summary of Parameters Used in Proposed Model Dam Material ( N i ) 6 0 Sf/a'vo Dense Sand Foundation 120 0.3 Stiff Silt Foundation 100 0.3 Upstream and Downstream Rockfill 70 0.3 Upstream Sand and Gravel Shell (a) 35 0.3 (b) 25 0.3 (c) 12 0.2 Downstream Sand and Gravel Shell 35 0.3 Hydraulic Fill Sandy Silt Core 27 0.1 Similar to the deformation analyses, the present analyses were performed using the dynamic mode option in F L A C . Combined stiffness- and mass-proportional Rayleigh damping assuming a minimum critical damping ratio of 4 % about a central frequency of about 4 Hz was used. The central frequency corresponds to the predominant frequency of the input base motion. A minimum damping ratio of only 4 % was assumed because in the event that the core liquefies, additional damping due to the cyclic or 'ratchet' effect of the post-liquefaction model will occur. Because the core is subjected to low or negligible static shear stresses, cyclic loading may cause large cyclic shear stress-strain loops. As a result, damping will take place due to these loops. As presented in Chapter 4, this cyclic 'ratchet' effect is caused by the different loading and unloading moduli. Free field boundaries were applied at the upstream and downstream boundaries in order to minimize wave reflections and achieve free field conditions. The Caltechb earthquake history was input as a velocity history at the base of the geometry because, as described in the Appendix, F L A C version 3.3 does not integrate acceleration histories accurately. 122 Chapter 7. Total Stress Analyses of Coquitlam Dam 7.4.2 Results The analysis predicts that the entire core and a large portion of the upstream sand and gravel shell will liquefy. A plot showing the zones predicted to liquefy by this analysis indicates that the results are somewhat comparable to those found from the S H A K E analyses (Fig. 7.12). However, a larger portion of the upstream shell is predicted to liquefy in the proposed procedure than that predicted by S H A K E . This discrepancy in results may be due to either two-dimensional effects or to the damping values used. It is likely that the difference is due to the latter because the fraction of critical damping used in the S H A K E analyses were approximately 17 %, which is considerably greater than the 4 % used in these analyses. In result, the system is significantly underdamped in the pre-liquefaction phase of the response. In this case, the the error is on the conservative side. However, it should be noted that 4 % damping is applied to frequencies of about 4 Hz. For frequencies other that 4 Hz, the damping ratio is greater than 4 %. 0 50 100 150 200 250 300 Figure 7.16 Zones Predicted to Liquefy by Proposed Method A plot of the shear stress versus shear strain curve in the core material, although not shown here, indicated that liquefaction occurred at approximately 6 s. This corresponds approximately to the time when the computed maximum shear stress has already occurred in the S H A K E analyses. In S H A K E , the maximum shear stress in the core was estimated to occur at about 5.7 s to 6.1 s. 123 Chapter 7. Total Stress Analyses of Coquitlam Dam Using the proposed procedure, the crest is predicted to deform 0.6 m and 1.1 m in the horizontally upstream and vertically downward directions, respectively. The displacement in the vertical direction is comparable to those predicted to occur due to the loss in stiffness only as well as those due to loss in stiffness and inertia force. In the proposed procedure, however, the horizontal displacement is significantly greater than those predicted by the previous methods. This is likely due to the larger upstream sand and gravel zone predicted to liquefy in this method than in the other method. In fact, as shown in Fig. 7.17, the largest displacements occur in the upstream sand and gravel berm. In this area, the predicted maximum displacements are about 1.51 m in the upstream and upwards direction. Although the displacements are significant, they are considered unlikely to cause overtopping because in fact the highest elevation of the berm and the core never drop below the reservoir water level. In an additional run, the critical damping ratio was increased to 8 % in order to assess the effect of the damping on the response. Although not shown, the upstream liquefied zone is only very slightly reduced. The horizontal and vertical displacements of the crest are 0.35 m in the upstream direction and 0.65 m downwards respectively. These values are nearly half of the displacements predicted when using a critical damping ratio of 4 %. Another different run was carried out to study the effect of time of liquefaction triggering. The entire liquefiable upstream sand and gravel shell and hydraulic fill core material were predicted to liquefy at a dynamic time of 6 s. The predicted displacements of the crest are 0.083 m upstream and 0.935 m downwards in the horizontal and vertical directions. The vertical displacement increases by nearly 0.3 m from that predicted when liquefaction in different zones are triggered at different times. In contrast, the horizontal displacement is approximately 0.3 m less than in the previous run. A plot of the displacement vectors, shown in Fig. 7.18, indicates that the displacement pattern is similar to those obtained in the gravity only and gravity plus velocity pulse deformations analyses. This is expected because liquefaction in liquefiable zones in those analyses are assumed to all occur at the same time. 124 Chapter 7. Total Stress Analyses of Coquitlam Dam m 1 0 ( H 0 5 0 1 0 0 1 5 0 2 0 0 2 5 0 3 0 0 Figure 7.17 Post-Liquefaction Displacement Pattern Predicted by Proposed Procedure Figure 7.18 Post-Liquefaction Displacement Pattern Assuming Liquefiable Zones Triggered at the Same Time 125 Chapter 7. Total Stress Analyses of Coquitlam Dam The cyclic shear stress-strain plot of a zone in the core material is shown in Fig. 7.15. The curves indicate that large loops are traced during post liquefaction response of the core material. As a result, damping is accounted for directly in model for this material. This indicates that the use of a critical damping ratio of only 4 % is reasonable, although the preliquefaction triggering may be higher. However, the error would be on the conservative side and would be acceptable. Assuming a direct relationship between damping ratio and displacement for a critical damping ratio of 4 %, horizontal and vertical displacements of about 0.16 m and 1.8 m are estimated if all liquefiable materials are assumed to liquefy at the same time. By comparing these results with those of the previous run of 0.58 m and 1.07 m in the horizontal and vertical directions, whereby liquefaction in different zones is triggered at different times, it appears that the time of liquefaction triggering is very important. 100 — | CO -100 1 Figure 7.19 Cyclic Shear-Strain Plot Predicted in Dam Core Material 126 Chapter 7. Total Stress Analyses of Coquitlam Dam 7.5 Deformation Due to Consolidation The deformation caused by the dissipation of excess porewater pressures will add to the total displacement due to liquefaction. A semi-empirical approach was used to evaluate these deformations induced by consolidation. The volumetric deformations can be estimated from a relationship developed by Ishihara and Yoshimine (1992). As summarized in Chapter 2, materials liquefied to 100 % pore pressure ratio behave like a liquid and cannot resist shear stresses. Upon drainage, the material consolidates under the applied stresses causing volumetric strains. This results in surface settlement. As described previously, the in-plane (Gx), out-of-plane (az) horizontal stresses, and the vertical stresses (ay) are set equal to the pore water pressure to simulate soil liquefaction. The bulk and shear moduli in the liquefied zones are adjusted such that volumetric strains only will occur due to the change in effective stresses. The magnitude of the strains was estimated from Ishihara and Yoshimine's chart. The change in stresses is equal to the effective stresses under gravity loading. This method was verified in FLAC; using both a single element and a column geometry, prior to its application on the dam. According to Ishihara and Yoshimine's relationship, a volumetric strain of about 3 % is predicted. It should be noted, however, that laboratory tests have shown that the volumetric strains associated with silty sands are generally about 1.2 % to 1.5 %, somewhat less than that obtained from the chart for clean sands. In result, the computed deformations are likely to be very conservative. For a volumetric strain of 3 % in the liquefied elements due to consolidation, the crest setllement was computed to be approximately 0.6 m. The displacement pattern due to volumetric deformation only is shown in Fig. 7.16. The resulting total crest movement including that due to inertia forces is about 1.6 m down in the vertical direction and 0.7 m in the upstream horizontal direction. A plot of the magnified deformed grid for comparison with the original grid is shown in Fig. 7.17. 127 Chapter 7. Total Stress Analyses of Coquitlam Dam 128 Chapter 7. Total Stress Analyses of Coquitlam Dam 7.6 Summary The proposed Total Stress Dynamic Analysis procedure was applied to Coquitlam Dam. The proposed method captures both liquefaction triggering and post liquefaction displacement prediction in a relatively simple procedure. This method allows liquefaction of elements to be triggered at different times, which can signifantly affect the predicted deformation for a two-dimensional geometry with different material properties. The results were discussed and compared with those obtained using traditional triggering methods and deformation analyses considering loss in stiffness only and loss in stiffness plus velocity pulse. The results of the proposed procedure were compared to those using S H A K E and predicted a larger liquefied zone in the upstream shell material. The reason for this is likely the small value of critical damping ratio. A low value of critical damping ratio is required because in the event that the core liquefies, damping is accounted for in the stress-strain loops. This value has a significant effect on the predicted deformations. For the Coquitlam Dam applied under a Caltechb earthquake motion, crest displacement results estimated using the proposed procedure are comparable to those using the more traditional deformation methods although the pattern of displacement is different. Significant deformation is predicted in the upstream sand and gravel shell berm, which can be attributed to the larger zone of liquefaction in the upstream shell material. However, overall, the proposed procedure compares favourably with the traditional methods and most importantly illustrates the importance of trigger time in the liquefaction-induced deformation analyses. The results of all of the present analyses, in terms of maximum displacement and deformation pattern, are comparable to those obtained in 1984 under the O B E . This is reasonable considering that the peak ground earthquake acceleration used in these analyses is similar to the 1984 O B E . A larger magnitude design earthquake may cause significantly larger displacements in the crest, however, due to the width and stiffness of the rockfill berms, overtopping of the dam may not occur. This should be verified if the design earthquake is significantly larger than that used in the present analyses and if more up to date laboratory or field 129 Chapter 7. Total Stress Analyses of Coquitlam Dam testing indicates that the portion of the upstream shell closest to the rockfill berms is less dense than assumed in these analyses. 130 Chapter 8. Conclusions and Recommendations CHAPTER 8 CONCLUSIONS AND RECOMMENDATIONS Earthquake induced liquefaction of loose saturated soils can cause severe damage to earth structures. One well-known example is the failure of the Lower San Fernando Dam in which the driving stresses exceed the residual strengths of the liquefied soil. However, the development of significant deformations is not restricted to situations in which the driving stresses exceed the residual strengths. Laboratory tests and field examples indicate that liquefied soils can still transmit some finite amount of cyclic shear stresses. As a result, inertia forces are generated even in the post liquefaction phase. This can cause significant displacements particularly if there is an initial static shear stress acting on the structure. A total stress procedure was presented in which the deformations due to the reduced stiffness and strength of liquefied soil as well as those due to the earthquake induced inertia forces are computed. Unlike other total stress procedures, such as Jitno and Byrne's in which inertia forces after liquefaction are accounted for indirectly in the residual strains based on a factor of safety, the present method includes the inertia forces in the analysis. In addition, liquefaction is triggered in each element at different times. Triggering of liquefaction is predicted by weighting the cyclic shear stresses induced by a prescribed motion and converting each half cycle into an equivalent number of cycles. Liquefaction is triggered when a specified number of cycles has been exceeded. Upon liquefaction, post liquefaction properties are assigned to the liquefied element and the stresses are reset. A FISH subroutine which performs the procedure has been written for use in the computer code F L A C . The ability of the procedure to predict liquefaction triggering was verified against the computer code S H A K E . The results indicate that the method yields results comparable to S H A K E . However, the results are sensitive to the model parameters chosen for the analysis. For 131 Chapter 8. Conclusions and Recommendations the present it is recommended that S H A K E analyses be performed in order to compare the results with the proposed procedure. The predicted liquefaction induced displacements were also compared to those predicted by Bartlett and Youd's empirical equation for sloping ground. The results indicate that the displacements are comparable. Sensitivity analyses on the input parameters of the model show that the pre-liquefaction parameters do not significantly affect the results. This may be because the high acceleration amplitudes in the earthquake record used in these analyses only occur early on in the time history. As a result, the predicted post liquefaction response would not be significantly affected by the time of liquefaction triggering. However, if relatively high acceleration amplitudes are maintained throughout the duration of the earthquake, then triggering time may become important. The procedure was applied to seismic analyses of Coquitlam Dam. Coquitlam Dam is a 30 m high hydraulic fill dam with high consequences of failure, the dam has been classified as an extrerne hazard. The seismic stability of the dam has not been investigated since 1984. Since then the state of practice in seismic assessment of earth structures has progressed. The results of the proposed procedure were compared to displacements predicted using variations of the more common methods such as the modified modulus method and Jitno and Byrne's extended Newmark method. The results are comparable also illustrate the importance of taking into account that liquefaction of soil elements could occur at different times. Based on the above results obtained using the proposed method, the following recommendations are made for future research: 1. Perform a sufficient number of analyses comparing the results of liquefaction triggering with those of S H A K E in order to develop a data base of recommended fraction of shear modulus, G / G m a x and damping values accounting for earthquake base motion. 2. Perform additional one-dimensional analyses using different earthquake base imput motions to better understand the effect of different earthquakes on the proposed procedure. 132 Chapter 8. Conclusions and Recommendations 3. Compare the results of the liquefaction triggering with the results using a two-dimensional method such as FLUSH 4. Verify the procedure with a number of actual case histories. 5. Fine-tune the liquefaction triggering procedure so that the non-linear stress strain behaviour of soil is directly incorporated by using different shear moduli compatible with strain level at a particular time during dynamic motion. As a result, damping is directly incorporated into the procedure and will not adversely affect the post-liquefaction triggering phase of the procedure. 6. Develop a non-linear shear stress-strain relationship for use in F L A C , which is not possible in S H A K E or F L U S H , and verify with case histories and compare with other procedures. 133 Bibliography BIBLIOGRAPHY 1. Bartlett, S.F. and Youd, T.L. , 1992. "Empirical Prediction of Lateral Spread Displacement," Proc , 4th US-Japan Workshop on Earthquake Resistant Design, Lifeline Facilities and Countermeasures for Soil Liquefaction, Honolulu, Hawaii, N C E E R , State University of New York, Buffalo. 2. B C Hydro, 1980. "Coquitlam Dam and Reservoir - Memorandum on Seismic Stability," Report No. H1171, March 1980. 3. B C Hydro, 1984. "Coquitlam Dam Rehabilitation - Memorandum on Design," Report No. H1757, September 1984. 4. Byrne, P.M., 1990. "A Cyclic Shear-Volume Induced Coupling and Pore-Pressure Model for Sand," Soil Mechanics Series No. 145, Dept. of Civil Engineering, University of British Columbia. 5. Byrne, P.M., 1991. "A Model for Predicting Liquefaction Induced Displacements," Proceedings, 2nd Int. Conf. on Recent Advances in Geotechnical Earthquake Engineering and Soil Dynamics, St. Louis, Missouri, Vol. 2, pp. 1027-1035. 6. Byrne, P.M., 1996. "Modelling Liquefaction," F L A C Seminar, Vancouver, B C , 7. Byrne, P.M., Cheung, H., Yan, L., 1987. "Soil Parameters for Deformation Analysis of Sand Masses," Canadian Geotechnical Journal, 24(3), 366-376. 8. Byrne, P.M. and Mclntyre, J.D., 1994. "Effective Stress Liquefaction Analysis at the Wildlife Site," 3rd International Conference on Recent Advances in Geotechnical Earthquake Engineering and Soil Dynamics, St. Louis, Missouri, Vol. 2, pp. 1864-1898. 9. Campbell, K.W., 1990. "Empirical Prediction of Near-Source Soil and Soft-Rock Ground Motion for the Diablo Canyon Power Plant Site, San Luis Obispo Country, California," Report Prepared for Lawrence Livermore National Laboratory, Dames and Moore, Evergreen, Colorado. 10. Castro, G. , 1969. "Liquefaction of Sands," Ph.D. Thesis, Harvard Soil Mechanics Series, No. 81, Harvard University, Cambridge, Mass. 11. Chern, J . C , 1985. "Undrained Response of Saturated Sands with Emphasis on Liquefaction and Cyclic Mobility," Ph.D. Thesis, The University of British Columbia, Vancouver, B C . 12. Craig, R.F., 1992. Soil Mechanics, 5th Edition, Chapman and Hall, London, 427 pp. 134 Bibliography 13. De Alba, P., Seed, H.B., and Chan, C.K., 1976. "Sand Liquefaction in Large-Scale Simple Shear Tests," J . Geotech. Engrg., A S C E , 102(9), pp. 909-927. 14. Duncan, J.M., 1992. "State-of-the-Art: Static Stability and Deformation Analysis," Proc. Specialty Conference on Stability and Performance of Slopes and Embankments II, A S C E , New York, Vol. 1, 222-266. 15. Erten, D. and Mahler, M.H., 1995. "Cyclic Undrained Behaviour of Silty Sand," Soil Dynamics and Earthquake Engineering, 14, pp. 115-123. 16. Finn, W.D., 1996. "Evaluation of Liquefaction Potential for Different Earthquake Magnitudes and Site Conditions," Proc. Symposium on Recent Developments in Seismic Liquefaction Assessment, Vancouver, B C . 17. Geomatrix Consultants, 1992. "Seismic Ground Motion Study for West Sand Francisco Bay Bridge," Draft Report to Caltrans, Division of Structures, Sacramento, California. 18. Haile, J.P., Brouwer, K.J., Manning, I.B., Byrne, P.M., 1996. "Design and Seismic Displacement Analyses for the Kensington Tailings Dam, Alaska," ICOLD, Chile, 1996. 19. Hamada, M., Towhata, I., Yasuda, S., and Isoyama, R., 1987. "Study of Permanent Ground Displacements Induced by Seismic Liquefaction," Computers and Geotechnics, 4(4), pp. 197-220. 20. Holzer, T.L., Bennet, M.J., and Youd, T.L., 1989. "Lateral Spreading Field Experiments by the US Geological Survey," 2nd US-Japan Workshop on Liquefaction, Large Ground Deformation and Their Effects on Lifelines," Edited by T.D. O'Rourke and M. Hamada, Technical Report NCEER-89-0032, State University of New York at Buffalo, NY, pp. 82-101. 21. Ishihara, K., 1982. "Evaluation of Soil Properties for Use in Earthquake Response Analysis," International Symposium on Numerical Models in Geomechanics, Zurich, Switzerland, 13-17 September, 1982 22. Ishihara, K., 1993. "Liquefaction and Flow Failure During Earthquakes," Geotechnique, 43(3), 351-415. 23. Ishihara, K. and Kosecki, J . , 1989. "Cyclic Shear Strength of Fines-Containing Sands," Earthquake and Geotechnical Engrg., Japanese Society, SM & F E , pp. 101-106. 24. Ishihara, K., Tatsuoka, F, and Yasuda, S., 1975. "Undrained Deforamtion and Liquefaction of Sand Under Cyclic Stresses," Soils and Foundations, 15(1), pp. 29-44. 25. Ishihara, K., Yasuda, S., and Yoshida, Y., 1989. "Liquefaction-Induced Flow Failure of Embankment and Residual Strength of Silty Sands," Soils and Foundations, 26(3), pp. 69-80. 135 Bibliography 26. Ishihara, K. and Yoshimine, M., 1992. "Evaluation of Settlements in Sand Deposits folloiwng Liquefaction During Earthquakes," Soils and Foundations, 32(1), pp. 173-188. 27. Jitno, H., 1995. "Liquefaction Induced Deformations of Earth Structures," Ph.D. Thesis, Department of Civil Engineering, University of British Columbia, Vancouver, B C . 28. Joyner, W.B. and Boore, D.M., 1988. "Measurement, Characterization, and Prediction of Strong Ground Motion," Earthquake Engineering and Soil Dynamics II - Recent Advances in Ground Motion Evaluation, Geotechnical Special Publication 20, A S C E , pp. 43-102. 29. Koester, J.P. , 1994. "The Influence of Fines Type and Content on Cyclic Strength," Ground Failures Under Seismic Conditions, Geotechnical Special Report No. 44, A S C E , pp. 17-33. 30. Kuerbis, R.H., 1989. "Effect of Gradation on Fines Content on the Undrained Response of Sands," M.A.Sc. Thesis, The University of British Columbia, Vancouver, B C . 31. Kuerbis, R.H., Negussey, D., and Vaid, Y.P. , 1988. "Effect of Gradation and Fines Content on the Undrained Response of Sand," A S C E Conference of Hydraulic Fill Structures, Geotechnical Special Publication 21, pp. 330-345. 32. Lee, K.L., 1974. "Seismic Permanent Deformations in Earth Dams," Report No. U C L A - E N G -7497, School of Engineering and Applied Science, University of California at Los Angeles. 33. Lee, K.L. and Albaisa, A., 1974. "Earthquake Induced Settlements in Saturated Sands," J . Soil Mech. and Fndn. Engrg. A S C E , 100(GT4), pp. ? 34. Lee, K.L. and Seed, H.B., 1970. "Undrained Strength of Anisotropically Consolidated Sand," Journal of Soil Mechanics and Foundations Division, A S C E , Vol. 96, No. SM2, pp. 411-428. 35. Lysmer, J . , Udaka, T., Seed, H.B., and Hwang, R., 1975. "LUSH: A Computer Program for Complex Response of Analysis of Soil-Structure Systems," Report E E R C 74-4, Earthquake Engineering Research Center, University of California, Berkeley. 36. Newmark, N., 1965. "Effects of Earthquake on Dams and Embankments," Geotechnique, 15(2), pp. 139-160. 37. Pillai, V .S. and Byrne, P.M., 1994. "Effect of Overburden Pressure on Liquefaction Resistance of Sand," Canadian Geotechnical Journal, 31, pp. 53-60. 38. Pillai, V .S. and Stewart, R.A., 1994. "Evaluation of Liquefaction Potential of Foundation Soils at Duncan Dam," Canadian Geotechnical Journal, 31, pp. 951-966. 39. Prakash, S. and Sandoval, J.A., 1992. "Liquefaction of Low Plasticity Silts," Soil Dynamics and Earthquake Engineering, 11, pp. 373-379. 136 Bibliography 40. Saboya, F. and Byrne, P.M., 1991. "Rockfill Parameters for Analysis of Concrete Face Rockfill Dams," Department of Civil Engineering, University of British Columbia, Vancouver, B C . 41. Salgado, F.M. and Pillai, V.S. , 1993. "Seismic Stability and Deformation Analyses of Duncan Dam," Proc. 46th Annual Canadian Geotechnical Conference, Saskatoon, pp. 259-270. 42. Sasaki, Y., Towhata, I., Tokida, K., Yamada, K., Matsuomo, H, Tamari, Y., and Saya, S. 1992. "Mechanism of Permanent Displacemnet of Ground Caused by Seismic Liquefaction," Soils and Foundations, 32(3), pp. 79-96. 43. Schnabel, P.B., Lysmer, J . , Seed, H.B., 1972. "SHAKE, A Computer Program for Earthquake Response Analysis of Horizontally Layered Sites," Earthquake Engineering Research Center, Report No. E E R C 72-12. 44. Seed, H.B., 1979. "Soil Liquefaction and Cyclic Mobility Evaluation for Level Ground during Earthquakes," J . Geotech. Engrg., A S C E , 105(2), pp. 201-255. 45. Seed, R.B. and Harder, L.F., 1990. "SPT-Based Analysis of Cyclic Pore Pressure Generation and Undrained Residual Strength," Proc. H. Bolton Seed Memorial Symposium, Vancouver, BC , 2, 351-377. 46. Seed, H.B. and Idriss, I.M., 1967. "Analysis of Soil Liquefaction: Niigata Earthquake," J . Soil Mech. Found. Div., A S C E , 93(SM3), pp. 83-108. 47. Seed, H.B. and Idriss, I.M., 1971. "Simplified Procedures for Evaluating Soil Liquefaction Potential," J . Soil Mech. Found. Div., A S C E , 97(SM9), pp. 1249-1273. 48. Seed, H.B. and Lee, K.L., 1966. "Liquefaction of Saturated Sands during Cyclic Loading," J . of the Soil Mech. and Found. Div., A S C E , 92(SM6), pp. 105-134. 49. Seed, H.B., Lee, K.L., Idriss, I.M., and Makdisi, F.I., 1975. "The Slides in the San Fernando Dams During the Earthquake of February 9, 1971," J . Geotech. Engrg., A S C E , 101(GT7), pp. 651-688. 50. Seed, H.B., Tokimatsu, K., Harder, L.F., and Cjung, R., 1985. "Influence of S P T Procedures in Soil Liquefaction Resistance Evaluations," J . Geotech. Engrg., A S C E , 111(12), pp. 51. Singh, S., 1994. "Liquefaction Characteristics of Silts," Ground Failures Under Seismic Conditions, Geotechnical Special Report No. 44, A S C E , pp. 105-115. 52. Stark, T.D., Ebeling, R.M., and Vettel, J . J . , 1994. "Hyperbolic Stress-Strain Parameters for Silts," J . Geotech. Engrg., A S C E , 120(2), pp. 420-441. 137 Bibliography 53. Stark, T.D. and Mesh, G. , 1992. "Undrained Shear Strength of Sands for Stability Analysis," J . Geotech. Engrg., A S C E , 118(11), pp. 1727-1747. 54. Stark, T.D. and Olson, S.M., 1995. "Liquefaction Resistance Using C P T and Field Case Histories," J . Geotech. Engrg., 121(12), A S C E , 856-1216. 55. Taboada-Urtuzuastegui, V.M. and Dobry, R., 1998. "Centrifuge Modeling of Earthquake-Induced Lateral Spreading in Sand," J . Geotech. Engrg., A S C E , 124(12), pp. 1195-1206. 56. Tatsuoka, F., Sasaki, T., and Yamada, S., 1984. "Settlement in Saturated Sand Induced by Cyclic Undrained Simple Shear," Proc , 8th World Conference on Earthquake Engineering, San Fransisco, Vol. 3, pp. 255-262. 57. Thevanayagam, S., Ravishankar, K., Mohan, S., 1994. "Steady-State Strength, Relative Density, and Fines Content Relationships for Sands," Transportation Research Record 1547, 61-67. 58. Thomas, J . , 1992. "Static, Cyclic and Post Liquefaction Undrained Behaviour of Fraser River Sand," M.A.Sc. Thesis, Department of Civil Engineering, University of British Columbia, Vancouver, B C . 59. Tokimatsu, K. and Seed, H.B., 1987. "Evaluation of settlements in Sand due to Earthquake Shaking," J . Geotech. Engrg, A S C E , 113(8), pp. 861-878. 60. Towhata, I. and Ishihara, K., 1985. "Undrained Strength of Sand Undergoing Cyclic Rotation of Principal Stress Axes," Soils and Foundations, Japanese Soc. S M & F E , 25(2), pp. 135-147. 61. Vaid, Y.P. , Byrne, P.M., and Hughes, J .M. , 1981. "Dilation Angle and Liquefaction Potential," J . Geotech. Engrg., A S C E , 107(GT7), pp. 1003-1008. 62. Vaid, Y .P. and Chern, J . C , 1983. "Mechanism of Deformation During Cyclic Undrained Loading of Saturated Sands," International Journal of Soil dynamics and Earthquake Engineering, 2(3), pp. 171-177. 63. Vaid, Y.P. and Chem, J . C , 1985. "Cyclic and Monotonic Undrained Response of Saturated Sands," Advances in the art of testing soils under cyclic conditions, A S C E Convention, Detroit, Michigan. 64. Vaid, Y.P. and Finn, W.D.L., 1979. "Static Shear and Liquefaction Potential," J . Geotech. Engrg., A S C E , 105(10). 65. Vaid, Y.P. and Thomas, K., 1995. "Liquefaction and Postliquefaction Behaviour of Sand," J . Geotech. Engrg., A S C E , 121(2), pp. 103-173. 138 Bibliography 66. Vaid, Y.P. , . "Liquefaction of Silty Soils," Ground Failures Under Seismic Conditions, Geotechnical Special Report No. 44, A S C E , pp. 1-15. 67. Vaid, Y.P. , Chung, E.K.F., and Kuerbis, R., 1989. "Preshearing and Undrained Response of Sand," Soils and Foundations, Tokyo, Japan, 29(4), pp. 49-61. 68. Vaid, Y.P. , Sivathayalan, S., Uthayakumar, M., Eliadorani, A., Konrad, J .M. , Hofman, B.A., and Robertson, P.K., 1998. "Static and Cyclic Undrained Behaviour of Soils at the C A N L E X Sites," Draft Version of C A N L E X report. 69. Yasuda, Y.P. , Masada, T., Yoshida, N., Nagase, H., Kiku, H., Itafuji, S., Mine, K., and Sato, K., 1994. "Torsional Shear and Triaxial Compression Tests on Deformation Characters of Sands Before and After Liquefaction," Proc. 5th US-Japan Worshop on Liquefaction Large Ground Deformation, N C E E R , State University of New York at Buffalo. 70. Yasuda, N., Matsuomo, N., 1993. "Dynamic Deformation characteristices of Sands and Rockfill Materials," Canadian Geotechnical Journal, 30, pp. 747-757 71. Yasuda, Y.P. , Nagase, H., Kiku, H., and Uchida, Y., 1991. "A Simplified Procedure for the Analysis of the Permanent Ground Displacement," Proc. 3rd Japan-US Workshop on Earthquake Resistant Design of Lifeline Facilities and Countermeasures for Soil Liquefaction, edited by T.D. O'Rourke and M.Hamada, Report No. NCEER-91-0001, N C E E R , State University of New York at Buffalo, pp. 225-236. 72. Yasuda, Y.P. , Nagase, H., Kiku, H., and Uchida, Y., 1992. "A Mechanism and a Simplified Procedure for the Analysis of Permanent Ground Displacement due to liquefaction," Soils and Foundations, 32(1), pp. 149-160. 139 Appendix I APPENDIX I 140 Appendix I TO: Dr. Peter Byrne DATE: 13 November, 1997 FROM: Charissa Dharmasetia SUBJECT: Damping and Integration in Dynamic F L A C As requested a package was prepared summarizing the results of a dynamic analysis on a column using the programs F L A C and S H A K E in order to assess the effect of damping in F L A C on the system. A total of 20 simple cases in F L A C were examined using two different input base accelerations, four types of damping for each input motion, and three values of central frequencies for cases where Rayleigh damping was used. Two cases were analyzed in S H A K E using the same two input base accelerations as in F L A C . Analyses in F L A C and S H A K E were carried out to at least four seconds after the end of dynamic motion. In general, it was found that F L A C and S H A K E results compare favorably, with the same damping factor, when Rayleigh damping is used. In addition, the resulting base displacement histories obtained from F L A C were examined and compared to base displacements manually calculated by numerically double integrating the acceleration histories. The formulations developed by Wilson and Clough (1962) were used to integrate the acceleration histories. A total of eight cases in F L A C , varying material properties and damping were performed using two different input base acceleration. The results show that F L A C does not calculate the displacements correctly and that the base displacement histories, which should only be dependent on input accelerations, vary with damping factors and material properties. DAMPING CASES ANALYZED A 100-ft (or 30.48-m) simple column was analyzed in both S H A K E and F L A C . The column is homogeneous and elastic and its properties are stress independent and have no modulus or damping attenuation. The damping factor used was 0.1 and the shear modulus input was 1840 141 Appendix I kips/ft2 in S H A K E corresponding to 88100 kPa in F L A C . The water table was assumed to be at ground surface. The soil density was assumed to be 0.1185 kips/ft3 in S H A K E and 1900 kg/m 3 in F L A C . The column was analyzed under two different input motions: a simple harmonic input with 0.2g amplitude and a frequency of 4 Hz lasting 10 seconds, and the Caltechb earthquake record reduced to last 37 seconds. The analysis was carried out to a time of at least 4 seconds after the end of earthquake motion in order to observe the effects of damping. In addition, in F L A C , the column was analyzed using the four different types of damping available: 1. Rayleigh Damping (both mass and stiffness proportional); • Central Frequency, fm i n, = dominant input frequency; • fm i n < dominant input frequency; • fmin > dominant input frequency; 2. Rayleigh Damping (mass proportional damping only); • Central Frequency, fmin, = dominant input frequency; • fm i n < dominant input frequency; • fm i n > dominant input frequency; 3. Rayleigh Damping (stiffness proportional damping only); • Central Frequency, f m i n , = dominant input frequency; • fm i n < dominant input frequency; • fm i n > dominant input frequency; 4. Local Damping. Rayleigh Damping was originally used in structural analysis to damp the natural oscillation modes of the system. As a result, it is frequency dependent. Its equations are expressed in matrix form 142 Appendix I with components proportional to mass and stiffness. Bathe and Wilson developed an expression whereby the critical damping ratio, X, at any angular frequency, CQJ, of the system can be found from the mass and stiffness proportional damping constants. From the equation, it can be shown that mass-proportional damping is dominant at lower angular frequency ranges and stiffness-proportional damping is dominant at higher angular frequency ranges. The critical damping ratio reaches a minimum at a frequency where mass and stiffness damping each supplies half of the total damping force. In F L A C , the user defines this minimum critical damping ratio, \ m m , and frequency, fmin. The damping ratio is constant over approximately a 3:1 frequency range about the minimum frequency, fm i n. As a result, frequency independent damping can be approximated by choosing 2~ to lie in the centre of the dominant range of frequencies present in the numerical model. This 'minimum frequency', fm i n is termed central frequency in F L A C . In F L A C , combined mass- and stiffness-proportional damping, mass-proportional damping only, and stiffness-proportional damping only can be specified. Theoretically, if combined mass and stiffness proportional damping is used, F L A C model response to frequencies outside the 3:1 frequency range about fm i n, should be over damped. If mass-proportional damping only is used, model response to frequencies, f, less than f~ should be 'overdamped' and response to f greater than f should be 'underdamped'. Conversely, if stiffness-proportional damping only is used, model response to f less than fm i n should be 'underdamped' and response to f greater than f m i n should be 'overdamped'. Local damping was originally used to equilibrate static simulations in F L A C . It is reported to be attractive to dynamic analysis because it is frequency independent and as a result, unlike the Rayleigh damping option, frequency does not have to be specified. Local damping operates by subtracting and adding mass at a gridpoint at velocity extremes, thus is performed twice per oscillation cycle. The proportion of energy removed can be related to fraction of critical damping. The local damping coefficient used in F L A C is equal to the product of pi and critical damping ratio. 143 Appendix I According to the manual, local damping appears to give good results for a simple case but is untested for more complicated situations. It was tested here, on a supposedly simple case. Another one of the more attractive features of local damping is that it uses larger time steps than Rayleigh damping. Unfortunately, it was recommended that in order to be able to compare these results to those with Rayleigh damping, the same time step as that used in Rayleigh damping should be used. The maximum shear stress profile and the acceleration history of the top of the column were used to compare the results of S H A K E and F L A C . The acceleration history of the base of the column was also plotted in order to check whether the input motion is the same as that in S H A K E . RESULTS Simple Harmonic Wave Rayleigh Damping (both mass and stiffness-proportional) Central Frequency, fmin = 4 Hz The acceleration histories of the top of the column, analyzed in F L A C and S H A K E , compare favorably. In both S H A K E and F L A C , the initial response is greater than the 'steady state' acceleration amplitude. However, these maximum initial amplitudes are greater in S H A K E (2.62 m/s) than in F L A C (2.4 m/s). The 'steady state' amplitudes are identical. After shaking, the residual oscillations in both F L A C and S H A K E match although the amplitudes in S H A K E are slightly greater than in F L A C . The acceleration at the top of the column comes to rest at the same time in S H A K E as in F L A C . The shear stress profiles determined from F L A C and from S H A K E also match up nearly perfectly. This is shown on the attached plots. Central Frequency, fm i n = 1.5 Hz 144 Appendix I The initial acceleration history response is significantly greater in the S H A K E model than in the F L A C model. Similarly, the 'steady state' and the post earthquake oscillation responses are also greater in magnitude in the S H A K E model than in F L A C . This indicates that the F L A C system is more damped than the S H A K E system. This response corresponds to what is expected from theory. However, uncharacteristic of overdamped systems, the F L A C model comes to rest at a later time than the S H A K E model. The F L A C shear stress profile is similar in shape to the S H A K E profile although seems shifted down. The magnitude of the F L A C profile is slightly less than that of the S H A K E profile. Central Frequency, fmin = 8 Hz The magnitudes of the initial as well as the 'steady state' acceleration responses of the top of the column are greater in the S H A K E model than in the F L A C model. Similarly, the response after base motion has ended is greater in magnitude in S H A K E than in F L A C . This response corresponds to theory. Unlike the previous case, however, the residual oscillations come to a rest earlier in F L A C than in S H A K E . The shear stress profiles predicted by S H A K E and F L A C models are similar in shape. However, the magnitude of the profile in F L A C is greater than in S H A K E . Once again, this shows that the F L A C system is more damped than the S H A K E system. Rayleigh Damping (mass-proportional damping only) Central Frequency, fmin = 4 Hz The steady state acceleration history responses of the top of the column from F L A C and S H A K E are similar. However, unlike the both mass- and stiffness-proportional damping case, the initial peak in acceleration in F L A C (3 m/s2) is greater than in S H A K E (2.62 m/s). In addition, after shaking, the residual oscillations do not match in phase or magnitude. The amplitudes of the oscillations in F L A C are greater than those in S H A K E and come to rest at a later time. This 145 Appendix I indicates that the F L A C system is less damped than the S H A K E system. This is consistent with theory because for Xm\n = 0.1, fm i n = 4 Hz, mass-proportional damping only supplies half of the specified damping at 4 Hz (in other words, damping ratio is actually only about 0.05). The shear stress profile in F L A C and S H A K E match in shape but not in amplitude. The amplitude in F L A C is greater than in S H A K E . If the critical damping ratio operating at 4 Hz is doubled to 0.2, the acceleration history response in F L A C matches to that in S H A K E . Oscillations after the end of dynamic motion come to a rest at the same time. This corresponds to theory since the actual damping ratio supplied by mass proportional damping operating at 4 Hz is 0.1. The shear stress profile in F L A C and S H A K E match in both shape and amplitude. Central Frequency, fmin = 1. 5 Hz The acceleration history response in F L A C is more transient that in previous runs. A 'steady state' response takes longer to achieve than previously. The initial response is significantly greater in amplitude in the F L A C model than in the S H A K E model. The amplitude gradually decreases to a steady state amplitude slightly greater than in S H A K E . After dynamic motions has stopped, residual accelerations continue to oscillate but decrease at a slower rate than the S H A K E response. Ten seconds after the end of dynamic motions, the accelerations in F L A C still have not come to a rest. This behaviour shows that the F L A C model is less damped than the S H A K E model. The response is even less damped than the previous run. These observations are consistent with theory, since even less damping is applied at the system frequency of 4 Hz than at 1.5 Hz. The shear stress profile response in the F L A C model is similar in shape but greater in magnitude than the S H A K E response. 146 Appendix I Central Frequency, fmin = 8 Hz The acceleration history response of the F L A C model compares very well with that obtained from S H A K E . This response is also consistent with theory since for Nin = 8 Hz, damping for system frequencies less than 8 Hz is greater than those for 8 Hz. It appears that the critical damping ratio of the system at 4 Hz is about 0.1 as opposed to 0.05 at frequency 8 Hz. The maximum shear stress profile response of F L A C also compares favorably with that of the S H A K E model. Rayleigh Damping (stiffness-proportional damping only) Central frequency, f,, = 4 Hz The acceleration history response of the top of the column in F L A C is slightly greater than the response in the S H A K E model. The F L A C model also comes to a rest at a later time after dynamic motion than the S H A K E model. This indicates that the F L A C model is less damped than the S H A K E model and is consistent with theory. According to theory, a critical damping ratio of 0.05 is applied at frequency of 4 Hz because stiffness damping contributes to only half of the specified critical damping ratio at the central frequency. The magnitude of the maximum shear stress profile is also greater in the F L A C model than in the S H A K E model. If the fraction of critical damping operating at the center of frequency 4 Hz is set to 0.2, then the stiffness-proportional contribution is 0.1 at 4 Hz. Theoretically, the F L A C response should be comparable to the S H A K E response. This response is somewhat observed, however, in general, the magnitude of the top acceleration history response is slightly less in F L A C than in S H A K E . However, the F L A C model comes to a rest at a later time than the S H A K E model. This may be because not enough damping is applied at the other system frequencies. 147 Appendix I The shear stress profile in the FLAC.model is comparable to that in S H A K E model. However, the response is F L A C is slightly less damped than that in the S H A K E model at depths greater than 15 m. Another significant drawback to using the stiffness-proportional damping only option in Rayleigh damping, besides the poor results is that the time step used is significantly smaller than that used in both mass- and stiffness-proportional damping or mass-proportional damping only. The time step used in the stiffness-proportional damping only option is on the order of 10:5 rather than 10-4 used in the previous options. As a result, the time it takes to run analyses using this option is significantly longer than in other options. Central frequency, fmin= 1. 5 Hz The initial, steady state, and post earthquake magnitude of the acceleration response of the top of the column is slightly greater in S H A K E than in the F L A C model. In general, the F L A C model is more damped than the S H A K E model. However, the F L A C model comes to a rest at a later time than the S H A K E model. This is likely because not enough damping is applied at the other system frequencies. This response is reasonable because the damping ratio at a frequency of 1.5 Hz should be greater than about 0.05. Because there is a linear relationship between damping ratio and frequency in stiffness-proportional damping, the damping at 3 Hz should be 0.1. Therefore, at a frequency of 4 Hz, the damping is greater than 0.1. The magnitude of the shear stress profile for depths less than about 16 m is greater in S H A K E than in F L A C . However, for greater depths, the reverse is true. Central frequency, fmin = 8 Hz 148 Appendix I The acceleration history response in the F L A C model is somewhat erratic but is in general greater in magnitude than the S H A K E response. The post-dynamic response is also greater in magnitude in F L A C than in S H A K E . In F L A C , the model does not come to a rest even 10 s after the end of dynamic motion. This indicates that the F L A C model is less damped and is consistent with theory. The shear stress profiles of the F L A C and S H A K E models are similar in shape. However, the magnitude in the F L A C model is greater than in the S H A K E model. Local Damping The damping coefficient was set equal to 0.314, which is equal to the product of the fraction of critical damping and pi. The time step was also set to be equal to that used in Rayleigh damping of 8.274e-4 in order to be able to compare the response. The top acceleration history response in F L A C does not compare in amplitude or shape with the response in S H A K E . The steady harmonic shape of constant amplitude is somewhat observed in F L A C , however, there are periodic waves of amplitudes of about 7.0 m/s2. In addition, the amplitude of the response drops suddenly in comparison to S H A K E but continues to oscillate at a very small amplitude even 10 s after dynamic motion has stopped. Similarly, the shear stress profiles in F L A C and S H A K E do not match in shape. However, the magnitudes of the shear stresses in F L A C are somewhat comparable those in S H A K E . This is shown on the attached plots. Caltechb earthquake motion Prior to analyzing the column in F L A C under different types of damping, the column was analyzed with the Caltechb input motion undamped in order to determine the frequency components that contained the most energy. These frequency components are required to specify the damping in Rayleigh damping to ensure damping is centered about the most dominant frequencies. The 149 Appendix I significant frequencies depend on the natural as well as on the forcing frequencies. In order to determine the significant frequencies, a fast fourier transform was performed on the acceleration history of the top of the undamped column. As shown on the attached spectral density plots, the dominant frequencies of the column are about 2 Hz followed by about 5 Hz followed by 9 Hz. The centre of the most significant frequencies is about 4 Hz, or 3.72 Hz to be more precise. This central frequency of about 4 Hz is coincidental with the dominant frequency of the Caltechb acceleration record. Incidentally, plots of the variation of critical damping ratio with frequency at central frequencies of 3.72 Hz, 1 Hz, and 9 Hz were also plotted on the Fourier spectrum plots in order to show the damping ratios used at the significant frequencies. A fast fourier transform was also performed on the Caltechb record in order to check for the frequency content to ensure that it will not affect the numerical accuracy of wave transmission. According to Kuhiemeyer and Lysmer, the following condition must be met for accurate representation of wave transmission through a model: where Al is the spatial element size, X is the wavelength associated with the highest frequency component that contains appreciable energy, f is highest frequency component in hertz, C is the speed of propagation giving the lowest natural mode Cp (p-wave) or C , (s-wave)). A/ 10 C 150 Appendix I where G is the shear modulus, K is the bulk modulus, and p is the mass density. For G = 88100 kPa, K = 234933 kPa, p = 1900 kg/m 3 C s = 215 m/s, C p = 431 m/s, therefore, C = C s . From the spectral density plot of the Caltechb acceleration record, the highest appreciable frequency component is less than 10 Hz. The associated wavelength is 22 m. From the above equation, Al = 2.2 m. All of the grid elements are 1.01 m. Therefore, the Caltechb acceleration input does not have to be filtered for high frequencies to ensure numerical accuracy for wave transmission. Rayleigh Damping (Mass- and stiffness-proportional damping) Central Frequency, fmin, = 3.72 Hz The acceleration histories of the top of the column of F L A C and S H A K E are comparable. However, for dynamic time less than about 10 s, when the acceleration amplitude is the most significant, the response in S H A K E is slightly greater in amplitude than the response in F L A C . This response is consistent with theory because the significant frequencies during this time is likely about 2.0 Hz, in which case, the F L A C response would be more damped. However, after the time of 10 s, the responses of F L A C and S H A K E are almost identical, including the time at which the top acceleration comes to a rest. The shear stress profile of F L A C matches with that of S H A K E in magnitude and shape down to a depth of 15 m. Beyond a depth of 15 m, the shear stress magnitude of F L A C is slightly greater than that of S H A K E . It is reasonable that the response in F L A C and S H A K E are not identical because the maximum shear stresses with depth may not necessarily occur under the same mode. However, it is expected that the F L A C response be less than the S H A K E response in magnitude because at all frequencies other than the central frequency, the damping ratio is greater than 0.1. 151 Appendix I Central Frequency, fmirt = 1 Hz The acceleration history response at the top of the column in F L A C is lesser in magnitude than in S H A K E at all times. The F L A C response also comes to a rest earlier than S H A K E after the end of dynamic motion. This response is consistent with theory because the F L A C model is more damped than the S H A K E model. Similarly, the shear stress profile is lesser in magnitude in the F L A C model than in the S H A K E model. However, at a depth greater than about 20 m, the F L A C response becomes greater than the S H A K E response. Central Frequency, fmin = .9 Hz For dynamic time less than 15 s, the acceleration response in F L A C is less in magnitude than in S H A K E . The response is in fact more damped than in the previous case. However, after a time of 15 s, the response in F L A C is comparable to the response in S H A K E . This response is reasonable because for the system frequency of 2 Hz, the critical damping ratio is significantly greaterthanOA. However, for frequencies greater than about 5 Hz, the critical damping ratio is about 0.1. The shear stress profiles of the F L A C and S H A K E models are comparable. For depths less than about 17 m, the S H A K E profile is slightly greater than the F L A C profile. However, below 17 m, the F L A C profile is greater than the S H A K E profile. Rayleigh Damping (Mass-proportional damping only) Central Frequency, fmin = 3.72 Hz It is difficult to compare the acceleration history response of the top of the column in F L A C with that in S H A K E because they are out of phase. In general, however, the F L A C response is slightly greater in amplitude than the S H A K E response. This indicates that the F L A C response is less damped than the S H A K E response. This response agrees with theory. 152 Appendix I The shear stress profile in F L A C is also greater in magnitude than that in the S H A K E model. In addition, the shape of the F L A C profile is slightly different than that of the S H A K E model. If the critical damping ratio is doubled to 0.2, the acceleration response of the F L A C model becomes more comparable to the S H A K E model. For dynamic times less than 10 s, there are only occasional times where the F L A C response is less than the S H A K E response. For times greater than 10 s, the F L A C and S H A K E responses are almost identical. This response is reasonable. Similarly, shear stress profile of F L A C is comparable to that of S H A K E . However, the F L A C profile is slightly greater in magnitude than that of S H A K E . Central Frequency, fmin = 1 Hz At all times, the acceleration response in F L A C is significantly greater in amplitude than in S H A K E . At a time of 40 s, 3 s after the end of dynamic motion, the F L A C model has not yet come to a rest and continues to oscillate at a high amplitude relative to the S H A K E model. This response is consistent with theory because at all frequencies greater than 1 Hz, the damping ratio is less than 0.05. The magnitude of the shear stress profile in F L A C is significantly greater than in S H A K E . This supports the theory that the F L A C model is appreciably less damped than the S H A K E model. Central Frequency, fmin = 9 Hz For dynamic times less than 15 s, the acceleration response in the F L A C model is slightly less in amplitude than the response in the S H A K E model. However, for times greater than 15 s, the response in F L A C and S H A K E are nearly identical. This response is reasonable because for frequencies less than 9 Hz, the damping ratio is greater than 0.05. For frequencies less about 5 153 Appendix I Hz, the damping ratio is less than 0.1. Because the most significant frequencies are about 5 Hz or less, the F L A C model is more damped than the S H A K E model. The shear stress profile of the F L A C model is comparable to that of the S H A K E model. However, the magnitude of the F L A C profile is slightly greater than that of the S H A K E model. Although the most significant frequencies are more damped in F L A C than in S H A K E , it is possible that the frequencies causing the greatest shear stresses are less damped than in S H A K E . Rayleigh Damping (stiffness-proportional damping only) Central Frequency, fmin = 3.72 Hz The acceleration response in the F L A C model is significantly greater in amplitude than in the S H A K E model. For times less than 10 s, the F L A C response is in phase with the S H A K E response. However, for times greater than 10 s, the F L A C response becomes out of phase. The magnitude of the F L A C response is only slightly greater than the S H A K E response after a time of 20 s. At time 40 s, when the S H A K E model has already come to a rest, the F L A C model continues to oscillate at a very small amplitude. This response is reasonable because for frequencies less than about 8 Hz (~2*fmin), the damping ratio is less than 0.1. The shear stress profile in F L A C is significantly greater than the profile in S H A K E . When the critical damping ratio of 0.1 operating at central frequency of 3.72 Hz is doubled to 0.2, the acceleration history response in the F L A C model is only slightly greater in amplitude than in the S H A K E model. For times greater than 20 s, the F L A C response is more comparable in magnitude. Although the critical damping ratio is doubled, the F L A C response is still less damped than the S H A K E response, this is probably because at the most significant frequency of the system of about 2 Hz, F L A C damping is less than 0.1. 154 Appendix I The shear stress profile in the F L A C model is greater in magnitude than in the S H A K E model although to a lesser extent than in the case where ^ m i n = 0.1. Central Frequency, fmin = 1 Hz With the exception of the response between times 10 and 15 s, the acceleration response in F L A C is less than the response in S H A K E . Between the times of 10 and 15 s, the F L A C response is comparable to the S H A K E response. This response is reasonable because for system frequencies of greater than 2 Hz, the damping ration is equal to or greater than 0.1. Since most of the significant frequencies are greater than 2 Hz, F L A C response theoretically should be equally damped or only slightly more damped than the S H A K E response. The shear stress profiles in the F L A C and S H A K E models are comparable. For depths less tan 15 m, the S H A K E response is slightly greater than the F L A C response. For depths greater than 15 m, the reverse is true. Central Frequency, fmin = 9 Hz The acceleration history response in F L A C is significantly greater than the S H A K E response at all times. The F L A C response gradually decreases in magnitude with time, since the amplitude of the forcing motion is relatively small at times greater than 10 s. At a time of 40 s, the F L A C model has not yet come to a rest. This response is consistent with theory since at all frequencies less than 9 Hz, the critical damping ratio is less than or equal to only 0.05. Similarly, the shear stress profile in the F L A C model is significantly less than in the S H A K E model. Local Damping Similar to the harmonic case, the dynamic time step was set to match the step used in the Rayleigh analysis. In this case, the time step used was 1.058e-4. 155 Appendix I For a damping factor of 0.314, the top acceleration history of the F L A C model is greater in magnitude than the response in S H A K E . In fact, the F L A C accelerations show no signs of coming to rest towards the end of the history. The shear stress profile obtained in F L A C is greater in magnitude than that obtained from S H A K E . In addition, as opposed to the S H A K E results, the F L A C profile is jagged not smooth. BASE DISPLACEMENT ANALYSIS Recently, concerns have been expressed over whether F L A C calculated base displacements correctly in dynamic F L A C . Previously, it had been found that the calculated base displacement histories in F L A C varied with material properties. In order to address these concerns, the following cases were analyzed: 1. Simple Harmonic Base Motion - Acceleration Input: • Rayleigh Damping (mass and stiffness proportional), f m j n = 4 Hz; • Rayleigh Damping (mass-proportional damping only), fm i n, = 4 Hz; • Rayleigh Damping (mass-proportional damping only), fm i n = 8 Hz; • G = 22,000 kPa, B = 60,000 kPa. 2. Caltechb Earthquake Record - Acceleration Input: • Rayleigh Damping (mass and stiffness proportional), fm i n = 3.72 Hz; • Rayleigh Damping (mass and stiffness proportional), fm i n = 1 Hz; • Rayleigh Damping (mass-proportional damping only), fm i n = 3.72 Hz; • G = 22,000 kPa, B = 60,000 kPa. 156 Appendix I In all cases, the base displacement and velocity histories were plotted and compared to the displacements and velocities integrated from the input acceleration records. Both acceleration records were integrated using formulations developed by Wilson and Clough(l 962). All cases using the simple harmonic base motion were plotted on one plot and all cases using the Caltechb Earthquake Record were plotted on another. INTEGRATION Wilson and Clough presented a step-by-step integration procedure to calculate the response of lumped-mass systems to arbitrary dynamic loads. This approach is an alternative to the mode-superposition method of solving the differential equations representing the equilibrium of the lumped-mass system. The step-by-step method involves the direct numerical integration or the equilibrium equations. It can be assumed that the acceleration varies linearly (or higher order) within a time increment. This method was applied to both the simple harmonic acceleration record and the Caltechb acceleration record. It was assumed that the acceleration varies linearly with time. The integration assuming linear variation was compared to the results of the direct integration of the simple harmonic acceleration record and was found to match. The following equations were used: {*}, ={a} + y - { J t } , {x}t={b} + ~ { x } , o where At2 157 Appendix I RESULTS Simple Harmonic Motion Acceleration Input As shown on the attached plot, the base displacement histories computed in FLAC for any of the above cases do not match the history integrated from the acceleration input motion. Although the amplitude of the displacements match, none of the histories in FLAC match each other. There seems to be a tendency for all histories to move towards one direction to different degrees varying with material properties and damping. This tendency is carried out even after earthquake motion has ceased with the exception of cases where mass-proportional only damping is used. For these cases, the relative displacements with time remain constant. This indicates that F L A C is not integrating properly. Because displacements are integrated directly from velocities, it is likely that F L A C is not integrating the velocities properly. A plot comparing the base velocity histories computed for the various cases examined is also attached. Contrary to previous beliefs, the plot shows that for all cases, the F L A C results match the velocity history integrated from the acceleration input motion. Strangely, it appears that F L A C does calculate the base velocities properly. The shear strain history response in F L A C for the case where Rayleigh Damping, with central frequency of 4 Hz, at a depth of 55 ft or 16.8 m was compared to that in S H A K E . As shown in the attached plot, the response in F L A C and S H A K E match nearly perfectly. A plot of the maximum shear strain profile is also attached. It shows that the shear strain profile also match perfectly. These results indicate that F L A C computes the relative displacements correctly which would explain the reason why in the damping analyses the results in F L A C compare well with those in S H A K E under certain conditions. This is reasonable because the stresses are computed from strain rates which in turn are computed from velocities. The resulting displacements, however, are wrong. 158 Appendix I Velocity Input All of the above four cases were re-run using the first integral of the acceleration history, the velocity history, as the input motion. The base displacement and base acceleration histories of all cases were plotted and compared with the histories determined by integration. The shear strain history at a depth of 55 ft as well as the maximum shear strain profile response in F L A C were plotted and compared to the results in S H A K E . As shown on the attached plot, the base displacement histories of all cases match perfectly with the displacement input history integrated from velocity. Similarly, the base acceleration history also matches perfectly. Like the case where the acceleration base input was used, the shear strain history at depth 55 ft and the maximum shear strain profile in F L A C matches nearly perfectly with that in S H A K E . The acceleration histories of the base as well as the top of the column in the F L A C model using combined mass and stiffness Rayleigh damping, fm i n = 4 Hz, were plotted with the corresponding histories in S H A K E . These plots show that the responses, like the case where an acceleration input was used, are comparable to the response in S H A K E . The maximum shear stress profile in F L A C is also comparable to the profile in S H A K E . Caltechb Input Motion Acceleration Input Similar to the results using simple harmonic wave input, the attached plot shows that the base displacements are not computed correctly in F L A C . In these runs, however, the differences in displacement histories between the various cases are more subtle. The reason is likely because 159 Appendix I the differences are small compared to the total magnitude of displacements. In general, F L A C tends to overpredict the displacements, in this case it is approximately by 10 percent. The plot of the comparison of the velocity histories with the integrated velocity history shows that the velocity is computed correctly in FLAC. The shear strain history response at a depth of 55 ft in F L A C also compares favorably with the response in S H A K E . The maximum shear strain profiles in FLAC and S H A K E also compare at a depth above 15 m. Below 15 m, the F L A C response is slightly greater in magnitude than the S H A K E response. Velocity Input If the velocity history, integrated from the Caltechb acceleration history, is used as the base input motion, the base displacement history used in F L A C matches perfectly with the displacements obtained from direct numerical integration of the velocity history. Similarly, the base acceleration histories of all the cases match with the original Caltechb acceleration history. The shear strain history maximum shear strain profile in F L A C and S H A K E also match comparably. The response is the same as when an acceleration input is used. The acceleration histories of the base and top of the column in the F L A C model using a velocity base motion are compared with the S H A K E responses. The responses are the same as when an acceleration base input is used. Comments Although F L A C appears to calculate the velocities and not the displacements correctly, R is likely the reverse is true because it is unreasonable that the displacements are improperly calculated when they are only simply integrated from velocities only. This was examined by integrating the velocity history output by F L A C by using the Wilson and Clough method to obtain displacement history and comparing with the displacement history directly output by FLAC. It was found that the integrated displacements match perfectly with F L A C . Thus F L A C does integrate displacements 160 Appendix I correctly. The problem is likely because a negligible error is introduced in the integration of the velocities from the acceleration histories which is magnified in the displacement histories when they are integrated from the velocity history. Although F L A C does not integrate the velocities correctly when an acceleration base input motion is used, the errors are negligible. As a result, the stresses and strains are not appreciably affected. The errors become apparent only when displacements are required. CONCLUSIONS Maximum shear stress profile and the surface or top acceleration histories were used to compare the dynamic results of a 100-ft homogeneous and elastic column modeled in F L A C and S H A K E in order to assess the effects of the different types of F L A C damping on the system. The types of damping compared in the analyses are: Rayleigh damping with combined mass and stiffness contribution, Rayleigh damping with mass-proportional damping only, Rayleigh damping with stiffness-proportional damping only, and local damping. For cases where Rayleigh damping was used, the central frequency was varied about the significant frequencies of the system. All cases were analyzed by applying a simple harmonic base motion and the Caltechb earthquake motion. The results indicate the following: • For identical critical damping ratios in F L A C and S H A K E , the F L A C response is comparable to the S H A K E response provided that Rayleigh damping with combined mass- and stiffness- proportional damping contribution, centered about the significant frequencies of the system, is used. The significant frequencies depends on the natural frequencies of the system as well as the input motion and can be determined by performing a fast fourier transform on the undamped history response of interest. 161 Appendix I • If the central frequency is specified outside the 3:1 frequency range of the significant frequencies of the system, the response is 'overdamped'. • Comparable results in F L A C using Rayleigh damping with mass-proportional damping only or stiffness-proportional damping only can be achieved by varying the central frequency as well as the damping ratios operating at the central frequency, provided that the system has only one significant frequency, which is unlikely. • In Rayleigh damping, mass-proportional damping is dominant over the lower frequency ranges, whereas stiffness-proportional damping is dominant over the higher frequency ranges. • Although local damping has the advantage of having a larger time-step than Rayleigh Damping and is frequency independent, it has been found to be completely useless. Base displacement histories were also compared in order to assess whether they are computed correctly. The displacements were compared for cases where damping and material properties were varied. All cases were analyzed by applying a simple harmonic acceleration and velocity history inputs as well as the Caltechb acceleration and velocity history dynamic inputs. The results indicate the following: • When an acceleration history input is applied, negligible errors are introduced in the integration of the history to obtain the velocity history. The errors, however, are magnified when the velocities are further integrated to obtain displacements. Displacements are, thus, overestimated. • Because the errors appear to be negligible, the stresses and strains computed are not significantly affected since they are computed from velocities. 162 Appendix I • It is recommended that acceleration histories are integrated and dynamic input is applied as a velocity history. If you have any question or comments, please do not hesitate to call me at 822-5974. Charissa Dharmasetia 163 Appendix I 7 (ui) gidaa c E 3 O o O o II ase nping CQ <5 II O 1 1 >» **-tor t3 O co c: tl X "§ c O) . c ion 1 Q . E <5 . c Q CO -C L— <u w> g> d) W u u C < •2 ts HI • g AK e » •Q 6 . 0 I .5 :a CO Q . O e 0 T J • § -c C 0 CO 0 O «> < •8 - 1 0 CQ u. E «*-o 0 c 6 c o CO s CO a. E o o (ZS/Ul) U0!)BJ3|30DV o o S o " a.?£ ,2 l« • co 1— 11 k. 0 o r; « c 11 C 8 §• 1 t I •S S ° 8-f o H UJ 1 £ ^ 'S -o" < S g x to e-s "° I i c ? o <0 c3 S o i 0 CO "S S i CO w a . E o o 3 ¥ 3 S o 5 (zs/ui) uCMiEJeiaoov 164 Appendix I o in (ui) M»d»Q c E 3 O o -» - o O ll a g> £ s-« c CQ | _ O w >•! i 0 St S I i X ? Q _ • c C Q. ,Oi ° •§ f 8 » s < S g y ^ § 1 t l •a -8 J ™ » 1 O a • 3 ? 3 u . 6 C o O •= c Ul o * JO s <5 s o. E o o (ZS/Ul) U0!JBJ0|903V O o 5 a. s .2 t O "5 c i y II U . C o </> X c o o O o o < UJ < X tn T J C CO O < _ l u . o Q. .cn •E £ v> Ct to ~_ o -9 O § ° ? 3 c o in _ 5 CO VI CL E o O o 3 (ZS/ui) uo|ieiS|S03V 165 Appendix I LO (ui) mdoci o ro E o o -6 to ° 3 -6 >> ¥ ~-6 *.« s I i x ? s ~ S & C CL E O -8 is £ » S, © o o < LU < X co T3 o l ! .5 £ o S o c o ro < a. E o o o oo E • o O • 2 lu (ZS/Ul) UO||EJ3|30QV o o 5 a . s >a o g-s i - s "> ° « J : » § -fc_ — o in o o If t l E -s. to qj » L ? < S i Ul I f 1^ Ui •o 1% c a O - 2 ro co a . E o O o c o in (ZS/Ul) U011BJ3|333V 466 Appendix I 168 Appendix I LO (ui) uidao. E 3 O ll " *> Q. m s £ m 1 ° - 2 " ° .8 ? >•€ § ° " S = -S -a X c £ C I § 2 a t <u » •§ O - tn O £ .o aj < § E LU a < o> 3. T . S E yZ CL W £ • Q •o 3 §. O a t ? < 3s _ l o 5 U . E * o .tJ O c t5 #- i1 iS g gjUl » J o co « 5 Q. u -E o o (ZS/Ul) UO|lBJ3|000V c E o a: 0 S H I c i t "° u S ° c » c § 1 | | C ® c £ to s g> • s i < i ; ui < x T J C CO i l l O 8 9" DJ 5. •S E 9- 5 E Q •H - c i cj CO Co ^ I 5 o O o 3 (ZS/Ul) U0|IBJ3|300V 1691 Appendix I in c E 3 O o o 0) CL in a x: CO o -~ in (ui) uidsa I fi O to O - -* - ° 'I O n .C © ? ! rag*-: CD | o ° 8 S >•!§ ° u- S SL § a c S- o •2 * S, ra 0) ! 6 o ™- » O E <*> Ml j i * I «• < O l 5> I .c S CO | E c ? ° ra a •§> O l | U. E «•- o £ O ? u -s m « O SC ! • M S Ul a *o 1 1 3 E C o o (zs/ui) uo|;eje|aootf c E N o "> " 5 ; o " s °> E ° » S >> S c S: s 0 u c ' * S i 1 -g •§ •2 t 8 ra S g 11 ^ o **» 1 3 « c 8 < s s UJ S I X > J ? £ ^ ! ? « 3 <° -g, E £ s i ra -CL E o o o 3 (zs/ui) uoiiejaiaoov 170 Appendix I in (UJ) uidaa c i «! O ll O II 0 » 55 CQ a £ M- 2 E 0 g 8 s | *« 1 c w C I I .2 -8 " •12. 8 * | < Is UJ I § ^ Cn X s a: co | . e O j a g O IL E i o = 5 C L i S O * « » 5 eg CO Q . E o o (zs/ui) uo!iej3|aoov c E CN, 3 0 O i , ° i s I— E O •= O £-1 S. o 5 3 I. « c I X * * « S. cu 0 a £ ~ ? 3= co -s 2. 1 °> g « g. < I a UJ § E < .c Ct co a s I '5co o I 5 | 5 u. I i o §1 CO w o. E o o CQ 5 (ZS/U|) U0!)eJ3|333V 1 7 1 -172 173 Appendix I in i--o _ in (ui) mdaa o CM c E 3 o o -« - a O ll s t §. CU m ui •jr 5 5 ° 8 a. g* 5 S Sit X c o 9- n I f 9. E 6— m CC 111 < X CO TJ C ^ £ CO o s O S 0-211 U. E •*- o ' O c O C m o * « « a E o o (ZS/UI) U0!ieJ9|30DV c E o <-> 5 o 11 Dl CL S O 9->< 1 s O O -Q X -8 a 1.? CD i= 5. o » 6 o c - a < 9 ° x CO TJ c CB " I S ^ a • "-6 3 "S 9 K ^ tt '5. S s i CO Ul CL E o O '3 (zs/ui) UO|lBia|0D0V 174 Appendix I LO E 3 O O o 0) o Q-CO U) 0) i _ 55 .c CO E 3 o a -^ _ o O II CD g> £ & Co p m | 0 g i >>:? -o j ; . a l l •2 •§ " S s S 0) CD E — TO (D W Q O U < <^ 8 <t < «£ 1 S 4 CO | . u -•a •§ 1 O 3 o < 3 3 J ou. LL E H— O o = C Lu o * tn s co * a E o o (ZS/Ul) UOI1BJ9I300V o " 5 o 11 Q. S 6 5 II « S | x ff-CD ~ c t «) 3 CD « c — ( D C S 5 Q < I f a) Ct < ^ E i CO CL . E y ' s » C 3^ CO LU O Scj 5 33 LU o 8 in 5 a. E o o o 3 (ZS/Ui) uoueJ3|a33v 1 7 5 ' Appendix I 10 c E 3 o o o CD cu c CO o _ co o O -^_ o O ll CD ? •2 I I c 1 ^ O -8 ° ca 0> OJ o o < Ul < x CO T J C CO O < w> CO I f 2 cc .g\£ I* o 3 e c o c o •£5 CO CL E o o (ZS/ui) u o j i e j a i a o o v O CL S o o cn 6 1 R! cn c c X « | 5 ® o CD 11 WJ to • i t I* 3 <J B 3 E o c o at CO C L E o o (zs /ui ) u o g i e i a i c o o v 176' Appendix! E 3 O o CD n o CO 10 CO CD CO CO CD CO o co (ui) utctea c E 3 o £ »- ° II O II c S s l £ B- -Cu ;:. • -CQ m o _ 2 n ° S S 0 . c S | i .2 f 8 * - CD CO m S 0) j j CD ~ o S S " E ^ 0) D) Zi S < fjj CO |g 1 ! s IS o | Lu E i O 8 # » -3 O < ui i £ CO « 3 Q. C E o o (zs /u i ) u o i i e i a i s s o v E o i O ^ -a .sM= I- I o O s i c5 ~ CD 5 c P CD p O » ID O c" » < I s ui < X co •o c CO O < _l Ii. o a .E ft S- a E o E g 6 i O CO CO O) < o o (ZS/Ul) U O J I B J 3 I 3 3 3 V 178?" Appendix I co (ui) u-dsa tt 5. a E O 0 - « •I- ° ^ O n e CD at Cu to - 2 § t l § FBI 1 ? « _ <D W C ct. o a s 2 a » • I f g « | O - rtj O § Q III £ f < o> cc 35 II •a •§ o C L. g « »s o 3 8 <. •§ " - J o > IL E o o f S C u, E o o (zs /ui ) uogie ja t soov I 3 ° s " t r u e ° tt ° ° - -S g O g - ' 9 H E =J *5 .-9 O fa x | | C U V ) O (f n CO -S ^-!« s tt - 3 -S 0) 4= Q . " w I < I 2 • •• 9, s> tt cc f § E £ TO * E < X CO TJ C CO O < E 3 o 5 O § s i at ^ co *- X • O 100 CL <t E o o (zs /ui ) u o i i e j o i a o o v 17m Appendix I to r-(iu) gjdaa c E 3 o fi ° * « o ° £ S 5- -CQ TO O »- 2 » ° 8 S o * S s i i X c § C a . S ° « 8 • - c o , « ? ! 8 s i ^ s? CO ? ! s O l s U. E S O c ~ c u s o S U J £ S o Q - u. E o o CM >-( j S / U l ) U 0 | I B J 3 | 3 0 0 V C E O o .c; 0 » F CL -S £ O g- ~-o o c * ^ 5 | § CO S g> ai 8 I a> § I O " ' a O C" » < I S ui § t < s.? x t | CO B. S •D § ° o . 2 < a E i 0 g | ^ -S> <2 * • 1 J o CZ V) % E 5 o o (zs/iu) uoitejaiaoav 480 ! Append/*/ LO (ui) wdaa c E o 1 « 0 ; O II -p a> g-E CQ •» II «•— • -. O <" * > . -S o o S i o -8 S. 0) a> o o < < X CA u c ca O < _ l LU o S a fc to «> a H o tn | ^ co « S Q. 5 E o o (zs/ui) uooejaiaosv 2 5 ^ o s._-CO 5 9 X ^ e < L U < X co CO O I E .9 1= l l < 3 Ct LU g j LU a 5 ui ™ a: u CO (o E * o o ( Z S / U l ) U 0 ! ) B J 3 | 3 0 0 V 181 Appendix I E O o o 0) o 0. CO CO CD w o _ LO o co (ui) M l d a a E - £ CD CO ca o V> x c o co < UJ < C ca O < s a JJQ o .91 tn * I f co o CL p E o o 0 0 9 T T (zs /u i ) uojiejaiaoov c E - * o 1 o ^ °> Z. 0 « w cn £ Q. S •£ H | cs O ^ t i c s co X o I S 8 6 * S g k_ -a> U 1 0 tn O c" S < 9 € Ul 1 t < s.g> to 9.1 •° 12 s - §> CD o q, Si S t o 8 .a-§ ! « " CO CO a. E o o 6 5 < o 6 (ZS/Ul) U0! l8J3|3O3V 163? Appendix I Appendix I c o o >. E s a. E " 55 CD Ir TT CD < I H- <U O E - s <n S a> 5 cn § 0) !E E c m <D O E » s s-l a> < cn - i ra u-E o O (iu) luduiaaeidsja TD O s § oi O O Xi c c o to o o to 2 § o o> 'c to o 3 f£ >» O CO « X c a g a. E 55 8 o c CO o E o CO "D <1) 2 cu Ol o 2 o c < .9 o o o > Itl T3 CU C CC CO V) & o < to >• x: >. 'o "o o o CD > > O CD < cn U-CO c m CD 1 CD C o U) CO CL E o O I I (S/Ul) /UpoiOA O CL >. i _ O CO X c CO f y ^ < m X "> CO To o f <J to —I i Li-ra CL E o o (%) "IBJ'S 185 Appendix I -186 Appendix I . 187
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Total stress dynamic analysis of Coquitlam Dam
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Total stress dynamic analysis of Coquitlam Dam Dharmasetia, Charissa W. 2000
pdf
Page Metadata
Item Metadata
Title | Total stress dynamic analysis of Coquitlam Dam |
Creator |
Dharmasetia, Charissa W. |
Date Issued | 2000 |
Description | Earthquake-induced liquefaction of loose saturated sands can cause large deformations to occur resulting in flow slides. The near catastrophic failure of the Lower San Fernando Dam due to the 1971 earthquake is one of the most well known examples of a flow slide. However, the concern is not limited to flow slide situations in which the driving stresses are greater than the residual strength after liquefaction. Significant deformations can occur for a non-flow slide condition. These are caused by a reduction in stiffness of liquefied materials as well as from inertia forces caused by the earthquake motion. Limit equilibrium type analyses have been used reliably in the prediction of the occurrence of a flow slide. However, difficulty has been encountered in the prediction of earthquake-induced deformation in the event that a flow slide is not predicted to occur. The problems have arisen from the modeling of the behaviour of liquefied soils. Soil exhibits stiff behaviour under cyclic loading and the strains associated up to the point of triggering of liquefaction are small. Upon liquefaction triggering, soil behaves like a liquid and initially strains under very small shear stresses. However, upon further deformation, the soil dilates and regains stiffness and strength. Most deformation analysis procedures overly simplify the effects of earthquake inertia forces and do not adequately model the stress-strain behaviour of liquefied soil. The proposed total stress procedure attempts to take account of both of the above effects. The procedure is separated into three main phases of cyclic-induced liquefaction behaviour of sands: the pretriggering, triggering, and post-triggering response of soils. The triggering of liquefaction in each element of a soil structure is predicted by weighting the cyclic shear stresses induced by a prescribed base motion. Upon triggering of liquefaction, liquefied stress-strain parameters are assigned to zones predicted to liquefy as they occur. The pre-triggering and triggering phases of the procedure were verified using SHAKE analyses. Similarly, the post-triggering phase of the procedure was compared with results obtained from Bartlett and Youd's empirical equation. In both cases, reasonable agreement was found. The method was finally applied to the Coquitlam Dam and the results compared with two of the more commonly used deformation analyses such as variations of the Modified Modulus method and Jitno and Byrne's extended Newmark method. The predicted results from all three methods are in reasonable agreement. |
Extent | 11497679 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-07-10 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0063512 |
URI | http://hdl.handle.net/2429/10592 |
Degree |
Master of Applied Science - MASc |
Program |
Civil Engineering |
Affiliation |
Applied Science, Faculty of Civil Engineering, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2000-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_2000-0382.pdf [ 10.97MB ]
- Metadata
- JSON: 831-1.0063512.json
- JSON-LD: 831-1.0063512-ld.json
- RDF/XML (Pretty): 831-1.0063512-rdf.xml
- RDF/JSON: 831-1.0063512-rdf.json
- Turtle: 831-1.0063512-turtle.txt
- N-Triples: 831-1.0063512-rdf-ntriples.txt
- Original Record: 831-1.0063512-source.json
- Full Text
- 831-1.0063512-fulltext.txt
- Citation
- 831-1.0063512.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0063512/manifest