A TWO-DIMENSIONAL NON-LINEAR STATIC AND DYNAMIC RESPONSE ANALYSIS OF SOIL STRUCTURES By RAJARATNAM SIDDHARTHAN B.Sc. The University of Sri-Lanka, Peradeniya Campus, 1977 M.A.Sc. The University of British Columbia, 1981 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY IN THE DEPARTMENT OF CIVIL ENGINEERING We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA MARCH 1984 (c) Rajaratnam Siddharthan, 1984 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 Cl y\/^ <C^v-^) \ r^Lfi-fi The University of British Columbia 1956 Main Mall Vancouver, Canada V6T 1Y3 Date ii ABSTRACT A method of analysis in two-dimensions to predict static and dynamic response of soil structures, including soil-structure interaction has been presented herein. The static and dynamic analyses can be performed in either effective or total stress mode or a combination of both modes. Non-linear stress-strain behaviour of soil has been modelled by using an incrementally elastic approach in which tangent shear modulus and tangent bulk modulus were taken as the two "elastic" parameters. The material response in shear was assumed to be hyperbolic coupled with Masing behaviour during unloading and reloading. Responses to changes in mean normal stress was assumed to be non-linear, elastic and stress dependent. Slip or contact elements have been incorporated in the analysis to represent the interface characteristics between soil and structural elements. The properties of the slip elements were assumed to be elastic, perfectly plastic, with failure at the Interface given by the Mohr-Coulomb failure criterion. In the static analysis proposed here, gravity may be switched on at once for the completed soil structure or the construction sequence can be modelled by layer analysis. The stress-strain conditions determined by the static analysis give the in-situ stress condition before the dynamic analysis. In the dynamic effective stress analysis, the residual porewater pressures are calculated using a modification of the model proposed by Martin, et al. (1975). The parameters, Gmax and _ax are modified for the effects of residual porewater pressure. The dynamic iii response study includes the prediction of post earthquake deformations. The predictive capability of the new method of analysis has been verified by comparing the recorded porewater pressure and accelerations of two centrifuged models subjected to simulated earthquakes, to those computed by the new method. This method has also been used to compute response of an offshore drilling island supporting a tanker mounted drilling rig. Results suggest that the common practice of neglecting soil-structure interaction may not be appropriate for islands which support heavy tanker type of structures. At present one-dimensional methods are used for computing the response of these islands. Comparative studies are also reported to asses the validity of this procedure. iv TABLE OF CONTENTS PAGE ABSTRACT ii TABLE OF CONTENTS v LIST OF FIGURES viiLIST OF TABLES xiNOMENCLATURE xv ACKNOWLEDGEMENTS ix CHAPTER 1 - INTRODUCTION 1 1.1 - SCOPE OF THIS THESIS 4 1.2 - ORGANIZATION OF THE THESIS 6 CHAPTER 2 - CRITICAL REVIEW OF SEED, ET AL. METHOD FOR COMPUTING DYNAMIC DEFORMATIONS 7 CHAPTER 3 - GENERAL GUIDELINES FOR DYNAMIC ANALYSES 14 3.1 - ONE-DIMENSIONAL RESPONSE ANALYSIS BY FINN, ET AL. (1977) 16 3.1.1 - Shear Stress-Strain Relationship 17 3.1.2 - Porewater Pressure Model 21 3.1.3 - Modification of Properties for Residual Porewater Pressure 22 3.1.4 - Influence of Strain Hardening 25 3.1.5 - Dissipation of Porewater Pressure 26 3.2 - LABORATORY VERIFICATION OF EFFECTIVE STRESS RESPONSE ANALYSIS 27 3.3 - FIELD VERIFICATION OF EFFECTIVE STRESS REPONSE ANALYSIS 30 V PAGE 3.4 - POREWATER PRESSURE MODEL IN PRACTICE 33 3.5 - DISCUSSION 36 CHAPTER 4 - TWO-DIMENSIONAL STATIC ANALYSIS OF SOIL STRUCTURES 38 4.1 - STRESS-STRAIN RELATIONSHIP 34.1.1 - Reasons for Selecting Gt and Bt 39 4.1.2 - Hyperbolic Shear Stress-Strain Relationship 40 4.1.2.1 - Estimation of Gmav 44.1.2.2 - Estimation of xmax 3 4.1.2.3 - Influence of Over-Consolidation 45 4.1.2.4 - Effects of Unloading 44.1.3 - Tangent Bulk Modulus Bt 7 4.2 - PHYSICAL MODELLING 48 4.3 - SIMULATION OF CONSTRUCTION SEQUENCE 49 4.3.1 - Method of Analysis 50 4.3.2 - Incremental Porewater Pressure 50 4.3.3 - Computation of Incremental Stresses and Strains 52 4.4 - SHEAR-VOLUME COUPLING 55 4.5 - INTERFACE REPRESENTATION 63 4.5.1 - Method of Analysis of Slip Element 68 4.6 - SELECTION OF SOIL PARAMETERS 70 4.6.1 - Obtaining Shear Stress-Strain Parameters 71 4.6.2 - Obtaining Bulk Modulus Parameters 72 vi PAGE CHAPTER 5 - TWO-DIMENSIONAL DYNAMIC ANALYSIS 73 5.1 - FORMULATION OF THE PROBLEM 75.1.1 - Incremental Equilibrium Equations 74 5.1.2 - Correction Forces 76 5.2 - DYNAMIC STRESS-STRAIN RELATIONSHIP 78 5.2.1 - Volume Change Behaviour 80 5.2.2 - Dynamic Shear Stress-Strain Behaviour 82 5.2.2.1 - Skeleton Curve for Dynamic Loading 82 5.2.2.2 - Unloading and Reloading 83 5.2.3 - Modelling the Effects of Residual Porewater Pressure 85 5.2.3.1 - The Behaviour of Samples with Ts 86 5.2.3.2 - Residual Porewater Pressure Generation Model 95 5.3 - DAMPING MATRIX [Cj 9 5.4 - BOUNDARY CONDITIONS 102 5.5 - ANALYSIS OF SOIL STRUCTURE SYSTEMS 104 5.5.1 - Slip Element in Dynamic Analysis 104 5.6 - SOLUTION SCHEME 105 5.7 - COMPUTATION OF POST EARTHQUAKE DEFORMATION 106 CHAPTER 6 - VERIFICATION OF THE METHOD OF ANALYSIS 108 6.1 - CAMBRIDGE CENTRIFUGE TESTS 106.2 - COMPARATIVE STUDY Ill 6.2.1 - Results of Test 1 115 6.2.2 - Results of Test 2 123 vii PAGE 6.3 - APPLICABILITY OF THE METHOD OF ANALYSIS 132 CHAPTER 7 - APPLICATION OF THE METHOD OF ANALYSIS: TANKER ISLAND RESPONSE 137 7.1 - INTRODUCTION 137.2 - ANALYSIS OF A TYPICAL TANKER ISLAND 138 7.2.1 - Results of Tanker Island Problem 142 7.3 - SOME PRACTICAL CONCLUSIONS 151 CHAPTER 8 - SUMMARY AND CONCLUSIONS 153 8.1 - SUMMARY 158.2 - CONCLUSIONS 4 8.3 - SUGGESTIONS FOR FURTHER STUDY 156 REFERENCES ' 157 APPENDIX I 16APPENDIX II 178 APPENDIX III 185 viii LIST OF FIGURES PAGE FIGURE 1.1 - Caisson-Retained Island 3 (De Jong and Bruce, 1978) FIGURE 2.1 - Conversion of Shear Strain Potential to Equivalent Shear Forces 10 FIGURE 3.1(a) - Initial Loading Curve 9 FIGURE 3.1(b) - Masing Stress-Strain Curves for Unloading and Reloading 19 FIGURE 3.2 - Hysteretic Characteristics 20 FIGURE 3.3(a) - Initial Effective Stress Condition in a Simple Shear Apparatus 23 FIGURE 3.3(b) - Intermediate Effective Stress Condition in a Simple Shear Apparatus 23 FIGURE 3.4 - Relationship Between Volumetric Strains and Porewater Pressures in Constant Strain Cyclic Simple Shear Tests, D_ = 45% 29 FIGURE 3.5 - Cyclic Stress Ratio Vs. Number of Cycle for Initial Liquefaction for Various Sands 31 FIGURE 3.6 - Predicted and Measured Porewater Pressures in Constant Stress Cyclic Simple Shear Tests, Dr = 45% 32 FIGURE 3.7 - Measured and Computed Ground Accelerations 34 FIGURE 3.8 - Measured and Computed Porewater Pressures at a Depth of 6m 34 FIGURE 4.1 - Element Stresses 4ix LIST OF FIGURES PAGE FIGURE 4.2 - Mohr Circle Diagram 44 FIGURE 4.3 - Mohr Envelope 6 FIGURE 4.4 - Effects of Loading and Unloading 46 FIGURE 4.5 - Sequential Procedure in Dam Construction 51 FIGURE 4.6 - Schematic Diagram Showing the Incremental Analytical Procedure 53 FIGURE 4.7(a) - Typical Plots for vs for Dense and Loose Sands 56 FIGURE 4.7(b) - Typical Plots for y vs for Dense and .Loose Sands 56 FIGURE 4.8 - Stress-Strain Behaviour of Ottawa Sand in Drained Simple Shear (After Vaid et al. 1981).... 58 FIGURE 4.9 - Typical Drained Triaxial Test Results on (a),(b) Dense Sacramento River Sand a) Principal Stress Ratio vs Axial Strain 59 b) Volumetric Strain vs Axial Strain 59 (After Lee, 1965) FIGURE 4.10 - Variation of Dilation Angle with Mean Normal Stress for Various Sands (After Robertson, 1982) 60 FIGURE 4.11 - Idealised Relationship Between ^ and 62 FIGURE 4.12 - Definition of Slip Element 65 FIGURE 4.13 - Plot of Typical Shear/Normal Stress vs Shear/ Normal Displacement 65 X LIST OF FIGURES PAGE FIGURE 5.1 - Dynamic Shear Stress-Strain Relationship: Skeleton Curve 84 FIGURE 5.2 - Dynamic Shear Stress-Strain Relationship: Unloading and Reloading 84 FIGURE 5.3 - Permanent and Cyclic Strains or Permanent and Cyclic Porewater Pressures 88 FIGURE 5.4 - Cyclic and Residual Behaviour of Pore Pressure and Axial Strain for ICT and ACT Tests 88 FIGURE 5.5 - q vs p' Plot for ICT and ACT Tests 91 FIGURE 5.6 - Normalized Curves for Pore Pressure Build-up Equation 94 FIGURE 5.7 - Stregth Curves for (a) Ekofisk Sands with D„ = 85%; (b) Sand with D_ = 77% 94 (After Rahman, et al. 1977) FIGURE 6.1 - Submerged Island Showing Transducer Positions ... 109 FIGURE 6.2(a) - Input Base Motion in Test 1 112 FIGURE 6.2(b) - Input Base Motion in Test 2 11FIGURE 6.3 - Liquefaction Resistance Curves of Medium Dense Leighton-Buzzard Sand 114 FIGURE 6.4(a) - Recorded Acceleration of ACC1244 in Test 1 116 (With and Without Slip Elements) FIGURE 6.4(b) - Computed Acceleration of ACC1244 in Test 1 116 (With and Without Slip Elements) FIGURE 6.5(a) - Recorded Acceleration of ACC1225 in Test 1 117 xi LIST OF FIGURES PAGE FIGURE 6.5(b) - Computed Acceleration of ACC1225 in Test 1 117 (With and Without Slip Elements) FIGURE 6.6(a) - Recorded acceleration of ACC734 in Test 1 118 FIGURE 6.6(b) - Computed Acceleration of ACC734 in Test 1 118 (With and Without Slip Elements) FIGURE 6.7(a) - Recorded and Computed Porewater Pressure of PPT2330 in Test 1 120 FIGURE 6.7(b) - Recorded and Computed Porewater Pressure of PPT68 in Test 1 120 FIGURE 6.7(c) - Recorded and Computed Porewater Pressure of PPT2338 in Test 1 121 FIGURE 6.7(d) - recorded and Computed Porewater Pressure of PPT2342 in Test 1 121 FIGURE 6.8(a) - Recorded Acceleration of ACC1244 in Test 2 125 FIGURE 6.8(b) - Computed Acceleration of ACC1244 in Test 2 125 (With Slip Elements) FIGURE 6.9(a) - Recorded Acceleration of ACC1225 in Test 2 126 FIGURE 6.9(b) - Computed Acceleration of ACC1225 in Test 2 126 (With Slip Elements) FIGURE 6.10(a)- Recorded Acceleration of ACC734 in Test 2 127 FIGURE 6.10(b)- Computed Acceleration of ACC734 in Test 2 127 (With Slip Elements) FIGURE 6.11(a)- Recorded and Computed Porewater Pressure of PPT2330 in Test 2 130 xii LIST OF FIGURES PAGE FIGURE 6.11(b)- Recorded and Computed Porewater Pressure of PPT68 in Test 2 130 FIGURE 6.11(c)- Recorded and Computed Porewater Pressure of PPT2338 in Test 2 131 FIGURE 6.11(d)- Recorded and Computed Porewater Pressure of PPT2342 in Test 2 131 FIGURE 6.12(a)- Effective Stress Paths in Test 1 135 FIGURE 6.12(b)- Effective Stress Paths in Test 2 136 FIGURE 7.1 - Schematic of Tanker Island 139 FIGURE 7.2 - Distribution of Maximum Dynamic Shear Ratio in Sand 143 FIGURE 7.3 - Distribution of Maximum Dynamic Shear Stress Ratio in Sand 143 FIGURE 7.4 - Distribution of Residual Porewater Pressure Ratio in Sand 145 FIGURE 7.5 - Distribution of Maximum Dynamic Shear Strain in Sand 145 FIGURE 7.6 - Distribution of Maximum Dynamic Displacement .... 147 FIGURE 7.7 - Post Earthquake X and Y Displacements 148 FIGURE 7.8 - Acceleration Response Spectrum for the Motion at Berm Surface 150 FIGURE 7.9 - Distribution of Residual Porewater Pressure Ratio in Sand 150 FIGURE Al.l - Iso-Parametric Element 167 xiii LIST OF FIGURES PAGE FIGURE A2.1 - Slip Element 179 xiv LIST OF TABLES PAGE TABLE 4.1 Variation of Exponent r\, with Plastic Index PI 42 TABLE 6.1 Soil Properties 113 TABLE 6.2 Recorded and Computed Maximum Accelerations 119 TABLE 6.3 Recorded and Computed Maximum Residual Porewater Pressures 124 TABLE 6.4 Recorded and Computed Maximum Accelerations 129 TABLE 6.5. Recorded and Computed Maximum Residual Porewater Pressures 133 TABLE 7.1 Static Soil Properties 138 TABLE 7.2 Dynamic Soil Properties 140 TABLE 7.3 ^s^avo and Kr Values 141 NOMENCLATURE xv A = Variables defined in slip element stiffness matrix formulation in Appendix II. a = Constant used to compute damping matrix, a^ - a^ = Constants defined in text. B = Variable defined in slip element stiffness matrix formulation in Appendix II. B = Strain-displacement matrix. Bt = Tangent bulk modulus. b = Constant used to compute damping matrix, b^ - b^ = Constants defined in text. [C] = Damping matrix. C_ = Interpolation matrix in slip element stiffness matrix formulation. C^ - C^ = Volume change constants. Cg = Cohesion is slip element. C' = Effective Cohesion. I) = Elasticity matrix. Dr = Relative density of soil. d = Depth of centre of gravity of an element below the top surface. {E} = Matrix defined in Appendix III. Et = Young's modulus. E_ = One-dimensional rebound modulus, e = Void ratio. NOMENCLATURE Maximum void ratio. Minimum void ratio. Matrix defined in Appendix III. Column vector of global inertia forces at time t Column vector of global spring forces at time t. Column vector of element spring forces at time t Shear strength of a slip element. Shear force per unit area in a slip element. Shear force per unit area in a slip element. Tangent shear modulus for loading and unloading. Tangent shear modulus. Maximum shear modulus. Hardening constants. Unit vector. Matrix defined in Appendix I. Jacobian matrix. Non-dimensional shear modulus constant. Bulk modulus constant. Anistropic consolidation ratio. Normal joint stiffness/unit area. Experimental constant used in porewater pressure model. Shear joint stiffness unit area. Slip element stiffness matrix in element axis. Global tangent stiffness matrix. xvii NOMENCLATURE [Kxy] = Slip element stiffness matrix in xy axis. KQ = Co-efficient of lateral earth pressure at rest. (K2)max = Shear modulus constant used to calculate G_ax. K = Constant which depends on the type of material. * [K ] = Porewater pressure matrix. K = Shear modulus constant for clay. K = Permiability of soil, z J L = Length of a slip element. H = Distance of a point of slip element in the element direction. [M] = Mass matrix. m = Experimental constant used in porewater pressure model. N = Scale factor of the centrifugal acceleration. [N] = Interpolation matrix. N^ = Number of cycles to cause initial liquefaction. N5Q = Number of cycles to cause a porewater pressure ratio of 50%. n = Experimental constant used in porewater pressure model. n = Bulk modulus exponent. OCR = Over consolidation ratio. {P} = Global column vector of applied loads. Pfl = Atomospheric pressure. {P } = Correction force vector. 1 corr-" PQ = Transformation matrix in slip element, p' = Average effective principle stress. xviii NOMENCLATURE 0^,0^ = Constant matrices defined in the text. q = Principle stress difference. R^ ,R_2 = Constant matrices defined in the text. T_ = Transformation matrix. t = time. U = Residual porewater pressure. U = Maximum residual porewater pressure, max r {U} = Column vector of incremental porewater pressures in elements. {TJ} = Displacement vector. u^ = Tangential displacement of node i. {uQ} = Element porewater pressure vector. V = Volume of an element. v^ = Normal displacement of node i. WgX = External work done. Wjjy = Internal work done. wn = Normal displacement in a slip element. w = Shear dispalcement in a slip element. s {X},{X},{X} = Column vectors whose components are relative displacement, velocity and acceleration of nodes. X^ = Base acceleration, x = Horizontal distance, z = Vertical distance. ix X ACKNOWLEDGEMENTS I would like to extend my sincere gratitude and appreciation to my research supervisor, Professor W.D. Liam Finn, for his guidance and keen interest and for making many valuable suggestions to improve the presentation of this thesis. I am greatly indebted to Professors, P.M. Byrne, V.P. Vaid, R.G. Campanella, and A.T. Bui for reading the thesis and for making constructive suggestions. Thanks are also due to Professor D.L. Anderson, for being always available to discuss any problem and for his useful and stimulating discussions. Special thanks are also due to Mr. F.H. Lee, for use of his test data and to Mr. Scott Steedman and Mr. Richard Dean for converting the test data to a form compatible with the computer at the University of British Columbia. The financial support of the Natural Sciences and Engineering Research Council of Canada, Department of Energy Mines and Resources, Government of Canada, Earth Technology Corporation, Long Beach California and Exxon Production Research Company, Houstan, Texas, is gratefully acknowledged. Finally, I wish to thank my wife Visha, for her patience, sacrifices and for the excellent job in typing this thesis. This thesis would not have seen this day without her whole-hearted support. I also remain indebted to my parents for their good wishes and continued moral support. It was due to their blessings and encouragement that I have completed this work. 1 CHAPTER 1 INTRODUCTION In engineering practice, it is generally agreed that the performance of soil structures subjected to seismic loading should be evaluated in terms of deformations rather than in terms of factors of safety. The allowable displacements may vary from a few inches to many feet depending on the functional aspects of the structures considered. Since the middle sixties many analytical methods for assessing earthquake induced deformations in soil structures have been proposed. These methods may be classified into two broad categories: one-dimensional methods and two-dimensional methods. The one-dimensional methods assume that deformations occur in parallel planes and that the material properties are either constant or vary normal to the planes only. The methods proposed by Newmark (1965) and Goodman, et al. (1966), assume that failure develops along well-defined failure planes and compute displacements of rigid blocks of soils. More recently the method proposed by Iai and Finn (1982) for long slopes accounts for the flexibility of the soil deposit. Often soil deposits cannot be characterized adequately as one-dimensional deposits, and the variability of properties in two or even three-dimensions must be considered. For example, in the response analyses of zoned dams, two or three-dimensional analyses are essential. When the third dimension is very much larger than the other two-dimensions and the properties do not vary significantly in this direction, a 2 two-dimensional response analysis is usually adequate. The geometry of the modes of deformations may also dictate two or three-dimensional analysis. For example in the analysis of embedded structures rocking may be an important deformation mode in addition to translation and therefore, at least a two-dimensional analysis is necessary. There are a number of two-dimensional methods available to compute seismic deformations. Some of these methods are based on elastic-plastic soil behaviour (Finn, et al. 1973; Mroz, et al. 1979; Prevost, 1979). These methods are complicated to use and have had very limited validation. The method proposed by Seed, et al. (1973), to compute seismic deformations of earth dams has found wide application in practice. This is a semi-analytical method in which the results of stress analysis and data from cyclic triaxial tests are used to estimate potential displacements in dams. Non-linearity of the soil is taken into account using an iterative elastic approach to achieve soil properties compatible with the computed strains. In recent years new types of structures have emerged, for which it is important to determine deformations under earthquake loading. Examples are the man-made sand islands which support drilling platforms for gas and oil explorations in the Beaufort Sea. These islands carry drilling equipment on either sand-filled caissons or tankers (Fig. 1.1). The deformations of these structures and islands during earthquakes are an important design consideration. A general satisfactory method for computing deformations of these structures is not available at present. Indeed the state-of-the-art for analysing deformations was recently assessed in a report on earthquake engineering research by the Fig. 1.1. Caisson-Retained Island (De Jong and Bruce, 1978). 4 National Research Council of the United States (USNRC, 1982; Finn, 1983) in the following terms: "Many problems in soil mechanics, such as safety studies of earth dams, require that the possible permanent deformations that would be produced by earthquake shaking of prescribed intensity and duration be evaluated. Where failure develops along well - defined failure planes, relatively simple elastoplastic models may suffice to calculate displacements. However, if the permanent deformations are distributed throughout the soil, the problem is much more complex, and practical, reliable methods of analysis are not available. Future progress will depend on development of suitable plasticity models for soil undergoing repetitive loading. This is currently an important area of research". For realistic predictions of stresses and displacements in soil structures the stress-strain behaviour of soils should be modelled as closely as possible. This is of course a difficult task since the stress-strain behaviour of soils is extremely complex. Using simple stressv strain relationships, a two-dimensional method for computing transient and permanent deformations in soil structures is presented in this thesis. 1.1 SCOPE OF THIS THESIS This thesis presents a method for two-dimensional static and seismic response analysis of soil structures. The earthquake loading occurs after static equilibrium or steady state conditions have been established. Soil properties such as strength and stiffness which control the response to earthquake loading depend on the effective stresses in the soil structure. Therefore, it is important to evaluate in-situ static stresses. These stresses are determined by a static analysis which takes into account non-linear stress-dependent 5 response of soil to load. Because soil behaviour depends on the loading path, the construction sequence of the soil structure is carefully modelled. A number of satisfactory methods of static analyses are already available (Kulhawy, et al. 1969; Duncan, et al. 1978). Nevertheless an independent method is presented here which uses a consistent set of material parameters in both static and dynamic analyses. This results In a much more efficient, cost effective solution to the problem of dynamic response analysis. The method for dynamic analysis takes into account the non linear hysteretic stress-strain behaviour of soils. The analysis may be carried out in either an effective stress or total stress mode, using an appropriate stress-strain relation. For effective stress response analysis, residual porewater pressures must be known. Therefore, a porewater pressure generation model has been developed for predicting seismically induced porewater pressures. The porewater pressure model is a generalization of the one-dimensional model of Martin, et al. (1975). The predictive capability of the new method for dynamic analysis has been verified by comparing the recorded porewater pressures and accelerations of a centrifuged model subjected to simulated earthquakes to those computed by the new method. This method has also been used to compute response of an offshore drilling island supporting a tanker-mounted drilling rig. The porewater pressures, stresses, accelerations and displacements in the island have been determined. At present one-dimensional methods are used for computing the response of these islands. Comparative studies were conducted to assess the validity of this procedure. Development of a non-linear method of analysis is very 6 difficult. A number of approximations have to be made to achieve a practical useful program. These approximations have been examined at length in this thesis and suggestions have been made for future research. 1.2 ORGANIZATION OF THESIS A critical review of the method proposed by Seed, et al. (1973) for computing seismic deformations is presented in Chapter 2. The formulations, basic assumptions and limitations of the model developed for the analysis of static response are presented in Chapter 3. The complete treatment of soil-structure interaction has also been included. The proposed two-dimensional dynamic response analysis is an extension of the one-dimensional response analysis of Finn, et al. (1977). Therefore, a detailed description of their model, its application in practice and in the laboratory is given in Chapter 4. The proposed method for dynamic response analysis is presented in Chapter 5. Details of the verification of the method is presented in Chapter 6. The method is used to compute response of a typical drilling tanker island subjected to seismic loading. The results of the analysis including implications for engineering design are discussed in Chapter 7. A brief summary, suggestions for future work and conclusions are given in Chapter 8. 7 CHAPTER 2 CRITICAL REVIEW OF SEED, ET AL. METHOD FOR COMPUTING DYNAMIC DEFORMATIONS The state-of-the art report on analysis of permanent deformations in earth structures (USNRC, 1982), suggests that the use of simple elastoplastic models for computing deformations may be adequate if deformations develop along well defined slip planes. When a well defined slip surface does not occur and deformations are distributed through out the soil structure, an analysis at least in two-dimensions is necessary. The most widely used two-dimensional method is the one that was proposed by Seed, et al. (1973, 1979). Detailed description and limitations of their method are presented below. 2.1 SEED, ET AL. METHOD (1973, 1979) The basic steps in the Seed, et al. method can be summarized as follows: a) Determine pre-earthquake or steady state condition that exists in the soil structure before the earthquake. b) Determine the design time history of base acceleration for the site where the earth structure is situated. c) Compute the time history of dynamic shear stresses throughout the soil structure using a two-dimensional dynamic response 8 analysis. Appropriate dynamic stress-strain relationship and damping should be used. d) Apply these stresses to undisturbed samples of soil consolidated to the initial static stresses in the soil structure to determine the strains and residual porewater pressures. e) Based on the porewater pressure data, determine the minimum factor of safety against total failure by limiting equilibrium methods after reducing the strength of elements which have developed significant seismic porewater pressures. f) If the soil structure is found to be safe against total failure, assess the overall deformation of the soil structure from the strains induced by the combined effects of static and dynamic loads as determined from the laboratory test data. Seed, et al. (1973) proposed an equivalent linear elastic method to model the dynamic non-linear, hysteretic behaviour of soils. The fundamental assumption in this type of approach is that the dynamic response of a non-linear hysteretic material may be approximated satisfactorily by a damped, elastic model if the properties of that model are chosen appropriately. The appropriate properties are obtained by an iterative process. In the dynamic finite element analysis, the stress-strain properties of the soil are defined in each finite element by the Poisson's ratio, v, and shear strain dependent shear moduli and equivalent viscous damping ratios. An average or effective shear strain (usually assumed to be 65% of the maximum shear strain) is computed in each finite element and shear moduli and damping ratios are selected compatible with these average 9 strains and are used in the next iteration. The procedure is repeated until no significant changes in moduli or damping ratios are necessary. The response determined during the last iteration is considered to be a reasonable approximation of the non-linear response. Since the final analysis with strain compatible soil properties is elastic the permanent deformation in the soil structures cannot be computed by this method. Deformations are estimated from the static and dynamic stresses with aid of strain data from cyclic triaxial tests. It is assumed that, when the static and dynamic stresses in a given finite element are simulated as closely as possible on a sample in a cyclic triaxial test, the resulting axial strain is the strain potential of the finite element in the dam. In practice this strain is converted to a shear strain potential by multiplying by a factor (1+v). The strain potential is the strain that develops in an unconstrained soil element under the specified loading. Since the finite elements in the soil structure are interconnected, the strains obtained by the above procedure are not the strains that will develop in the soil structure but are an indication of its potential for straining under the given seismic excitation. Serff, et al. (1976), have proposed a procedure for converting the strain potentials to a set of compatible deformations. The shear stress corresponding to the shear strain potential in a finite element is determined from the stress-strain curve (Fig. 2.1a). The shear stresses are converted to shear force and applied to the nodes of the finite element (Fig. 2.1b). The deformations under these nodal forces are then determined by a static analysis and are assumed to be seismically induced permanent deformations. This technique of computing compatible Fig. 2.1. Conversion of Shear Strain Potential to Equivalent Shear Forces. 11 deformation is sometimes referred to as a strain harmonizing technique. One of the serious limitations of any iterative elastic method, used to model non-linear behaviour, is that the solutions given by the method may not be unique (Desai, et al. 1977). This is because the solutions obtained in the last iteration may depend on the assumed soil properties of the first iteration. Equivalent linear method of analysis may overestimate the seismic response of non-linear hysteretic materials (Finn, et al. 1978a). The overestimation in linear methods occurs because of resonance. Resonance occurs when the fundamental period of the input motion corresponds to the fundamental period of the deposit as defined by the final set of compatible properties in the iterative equivalent linear method of analysis. Since the analysis is carried out with the final set of constant stiffnesses for the entire duration of the input motion, there is time for resonant response to build up. The stiffness properties in non-linear materials change constantly for every time step. When resonant response is a function primarily of the method of analysis, it is called pseudo-resonance. The Seed method is a total stress method and it does not take into account the effects of increasing porewater pressure on soil stiffness. Since response of soils is controlled by effective stresses, the validity of a total stress response analysis is questionable. Finn, et al. (1978) compared responses predicted by total stress and effective stress methods. They used two one-dimensional computer programs, SHAKE (Schnabel, et al. 1972) and DESRA1 (Lee, et al. 1975) for this purpose. 12 SHAKE is a total stress program which models soil as a damped equivalent linear elastic material and DESRA1 is an effective stress program which models soil as a non-linear hysteretic material. Finn, et al, (1978) concluded that the total stress analysis tends to overestimate the dynamic response when porewater pressures exceeded about 30% of the effective overburden pressure. The major difficulty in practice with the equivalent linear method is that a direct computation of permanent deformations is not possible. The concept of strain potentials has to be used to estimate permanent deformations (Serff, et al. 1976). There are two inconsisten cies in this procedure. First, the computed strains in the final iteration obtained with strain compatible soil properties are ignored, as not being correct, whereas the computed stresses are assumed to be representative stresses in the ground. Knowing that stresses and strains have a one to one relationship for a given loading, the arbitrary decision to ignore the computed strains is somewhat inconsistent. Second, the strains from the last iteration are ignored, whereas the strains computed in intermediate iterations were used to obtain compatible moduli and damping ratios. This type of inconsistent assumptions make the final estimated deformations somewhat arbitrary. When porewater pressures are allowed to dissipate in samples subjected to cyclic undrained tests, deformations occur. This plastic deformation is not accounted for in the approach proposed by Seed. Noting the limitations with the Seed method, which has been in use for about 10 years, the state-of-the-art report on analysis of permanent deformations in earth structures recommends this topic should be the subject of active research for the next ten years (USNRC, 1982). In 13 this thesis, an attempt which will allow direct computation of transient and permanent deformations of soil structures in a consistent manner, is presented. Procedures have been developed to model non-linear hysteretic behaviour of soil, taking into account the effect of porewater pressure generation on soil properties. 14 CHAPTER 3 GENERAL GUIDELINES FOR DYNAMIC ANALYSES The performance of an earth dam under earthquake induced ground motion Is an important concern in seismically active areas. The ground accelerations induced by the earthquake can cause large inertia forces throughout the dam. These forces which reverse in direction many times during an earthquake induce alternating stresses and strains in the dam. If these strains and associated displacements are large enough, large slumping and slope instability may result, leading to over topping and eventual failure of the dam. Newmark (1965), in his pioneering work on effects of earthquakes on dams, recommended that the performance of a soil structure during seismic loading should be assessed in terms of displacement and not in terms of the "factor of safety against failure" along an assumed failure surface. The allowable or satisfactory displacements depend mainly on the functional role of the soil deposit or the structure founded on the deposit. For example, for critical structures such as nuclear reactors and gravity platforms, the allowable displacement may be only a few inches; however for earth dams many feet may be acceptable. The main object of this thesis is to present a two-dimensional dynamic response analysis of soil structures, taking into account all important factors that influence the behaviour of soil deposits. The analysis predicts displacements, stresses, strains and acceleration fields etc, during and after the earthquake. 15 The determination of stresses, strains and displacements induced in a dam by an earthquake is a complex analytical problem and a number of simplifying assumptions must be made. Foremost is the assumption that the performance of an earth dam which is essentially a three-dimensional structure, can be interpreted from the performance of transverse cross-sections subjected to the same seismic loading. The plane strain condition is assumed to prevail in the proposed two-dimensional analysis. The main reasons for this assumption are the high cost and high computer storage requirements needed for a three-dimensional analysis. In the dynamic response analysis of continuous systems such as earth dams, non-uniform mass and stiffness distributions are present. The finite element approach, which can model the variation in stiffness and mass with extreme ease, has been adopted. In the dynamic response analysis of saturated soils, a decision must be made initially as to whether the analysis shall be carried out in terms of total stress or effective stress. Saturated loose cohesionless soils subjected to repetitive loading generate residual porewater pressures, and if sufficient drainage does not occur, reduction in effective stresses will result. Since deformations are controlled by effective stresses and soil properties such as moduli and strength are functions of effective stresses, an effective stress response analysis is always preferable for those type of soils. Effective stress response analysis is more difficult to perform. It requires porewater pressure generation and dissipation models and additional computations are necessary to estimate current effective stresses. Studies by Finn, et al. (1978a) on the response of level saturated sandy sites to seismic excitation showed that effective stress analyses are not generally 16 required unless the porewater pressures are likely to exceed 30% - 40% of the effective overburden pressure. In dynamic analyses, it is assumed that the loading imposed by seismic excitation is superimposed on long term equilibrium conditions. The dynamic soil properties, such as strength and stiffness, which control the response of soil structures to seismic loading, depend on initial insitu effective stress condition. Furthermore in effective stress response analysis, the computation of current effective stresses as porewater pressure develops is important. To do this, a static analysis, which uses material parameters applicable for both static and subsequent dynamic analysis has been developed. Finn, Lee and Martin (1977), presented a one-dimensional dynamic response analysis taking into account all important factors that affect the soil behaviour. The two-dimensional response analysis proposed in this thesis is an extension of their one-dimensional analysis. Therefore, a review of their method of analysis, including specifically how the non linear, hysteretic soil behaviour has been modelled is presented below. The predictive capability of their method has been verified in the laboratory and in the field. Some details on these verification procedures are also included. 3.1 ONE-DIMENSIONAL RESPONSE ANALYSIS BY FINN, ET AL. (1977) In horizontally layered deposits the assumption that the shear waves propagate vertically leads to a shear beam type of deformation pattern in the deposit. Then, only the stress-strain relationship in shear is required in the analysis. 17 The important factors that must be considered when computing the dynamic response of soils are; a) The nonlinear stress dependent stress-strain behaviour. b) The modelling of unloading-reloading. c) Contemporaneous generation and dissipation of porewater pressures. d) Hysteretic and viscous damping. e) Strain hardening. All these factors have been taken into account in the stress-The seismic loading imposes irregular loading pulses which consists of loading, unloading and reloading. The soils exhibit different behaviour in each of the above phases. The relationship between shear stress, x, and shear strain, y, for the initial loading phase under either drained or undrained loading conditions is assumed to be hyperbolic and given by, strain relations presented by Finn, et al. (1977). 3.1.1. Shear Stress-Strain Relationship G x = f(y) = max y (3.1) fl + G max IYI) T max 18 in which, G = maximum shear modulus and = shear max max strength. This initial loading or skeleton curve is shown in Fig. 3.1a. The unloading - reloading has been modelled using Masing behaviour (Masing, 1926). This implies that the equation for the unloading curve, if unloading occurs from (xr,Yr), is given by, Gmax " Vr) 1 + G v-Y (3-2) max I' 'rI 2 T max which is simply, T - Tr f(y - Yr) (3.3) The shape of the unloading - reloading curve is shown in Fig. 3.1b. Lee (1975) proposed rules for extending the Masing concept to irregular loading. He suggested that the unloading and reloading curves should follow the skeleton loading curve if the magnitude of the previous maximum shear strain is exceeded. In Fig. 3.2a, the unloading curve, beyond B, becomes the extension of the initial loading in the negative direction, i.e, BC. In the case of general loading history, further assumptions have to be made. If the current loading curve intersects the curve described by the previous loading curve, the stress strain curve should follow the previous loading curve. The above rules should apply also to unloading. Two typical examples are provided in Fig. 3.2b;(a) if loading along path BC is continued, the loading path is assumed to be 3.1(a). Initial Loading Curve. pig. 3.1(b). Masing Stress Strain Curves for Unloading and Reloading. I " (a) first unloading (b) general reloading Fig. 3.2. Hysteretic Characteristics. O 21 BCAM, (b) if unloading along path CPB is continued, the unloading path will be ABP' . Newmark and Rosenblueth (1971) have suggested a similar procedure. 3.1.2 Porewater Pressure Model Consider a sample of saturated sand under a vertical effective i stress, oy. During a drained cyclic simple shear test, a cycle of shear strain, y» causes an increment in volumetric compaction strain, Aey<j, due to grain slip. During an undrained shear test starting with the same effective stress system, the cyclic shear strain, y, causes an increase in porewater pressure, AU. It was shown by Martin, et al. (1975) that for fully saturated sands and assuming water to be incompressible, AU - Er Aevd (3.4) in which E„ = one-dimensional rebound modulus of sand at an t effective stress av« Martin, et al. (1975), also showed that under simple shear conditions the volumetric strain increment, AeV(j, is a function of the total accumulated volumetric strain, £vcp anc* tne amplitude of the shear strain cycle, y, and is given by, C3evd Ae = CAy-C e ) + —— (3.5) vd 1 2 vd Y + C, e , ' 4 vd 22 in which , C2, Cg and are volume change constants that depend on the sand type and relative density. An analytical expression for the rebound i modulus Er, at any effective stress level av is given by Martin, et al. (1975), as, - do' . E = ~— = (a*) /{m K (a ) } (3.6) r de v' 1 r v vo' ' v ' vr i in which, aVQ is the initial value of the effective stress and Kr, m and n are experimental constants for sand. The increment in porewater pressure, AU, during a given loading cycle with maximum shear strain, y, may now be computed using equations (3.4), (3.5) and (3.6). 3.1.3 Modification of Properties for Residual Porewater Pressure The residual porewater pressure, TJ reduces G„„„ and r r > max Tmax* These values should be updated as residual porewater pressure develops. Hardin, et al. (1972) assumed that Gmax is independent of stress history and suggested, * 1/2 G = K (a') ' (3.7) max m in which K is a constant, depends on soil type and relative density. The initial and current effective stress conditions in a simple shear apparatus with zero initial porewater pressure are shown in Fig. 3.3. Here it is assumed that the ratio between horizontal and vertical Fig. 3.3(a). Initial Effective Stress Condition in a Simple Shear Apparatus. m i«2K»«r- -u) 3 Fig. 3.3(b). Intermediate Effective Stress Condition in a Simple Shear Apparatus. Co 24 effective stresses is a constant KQ, where KQ is the coefficient of lateral pressure at rest. For the initial effective stress condition, ^max^ ls 8iven bv> 1/2 1 +9 V ' (G ) = K* fV>l (Ko)1/2 max 3 v ; o For the current stress conditions, 1/2 * 1+2K (G ) = K (if*) r;D-u)i/2 <3-9> max v 3 ; ; n on dividing the equation (3.9) by (3.8) one obtains, (Gmax) - U . n _ r vo il/2 .) L a' J (G ) o-' J (3.10) max vo o Therefore, knowing (^max^o' a'vo anc* ^> tne maximum shear modulus at the current effective stress condition can be calculated using equation 3.10. The shear strength (Tmax^o ^or tne Initial effective stress condition is given by (Finn, et al. 1977). f ^ 1 + K 9 o i 1 - K 1/2 (W = {(—-—°) sin2 0 - (^-^)2} 0' = C a' (3.11) o 1K 2 J K 2 ' ' vo vo 25 * in which 0' is the angle of internal friction and C is a constant which depends on soil properties. For the current stress condition, (W) = C*(o-' - U) (3.12) n vo Dividing equation (3.12) by (3.11) one obtains, (t ) max * ^ • (3.13) (T ) a* max vo o ^Tmax^o' "vo and U are known, the maximum shear strength at current effective stress condition can be calculated from equation (3.13). 3.1.4 Influence of Strain Hardening During seismic loading of dry sand or saturated sand under drained conditions, the sand structure hardens due to grain slip. Finn, et al. (1977), used following equations to modify Gmax and T • max evd (G ) = (G ) {1 + — ~~ } (3.14a) max max 1 H, + H„E ,J nn n 1 2 vd and (3.14b) 26 in which, (Gmax)nn and (T_ax)nn are the modified maximum shear modulus and shear strength in the n th cycle and , » H3 and are hardening constants. The stress-strain behaviour for one-dimensional analysis is now completely defined by equations, 3.1, 3.2, 3.4, 3.5, 3.6, 3.10, 3.13 and 3.14. In laboratory cyclic simple shear tests most of the volume changes in dry sands and the increases in porewater pressure in undrained saturated sands occur during the unloading portion of the load cycle. Therefore, Finn, et al. (1977), used modifications to the stress-strain curve to take account of strain hardening and porewater pressure only during the unloading phases of the loading. 3.1.5 Dissipation of Porewater Pressure If the saturated sand deposit can drain during shaking there will be simultaneous generation and dissipation of porewater pressure. The rate of increase of porewater pressure will be less than that in completely undrained sand. The amount of drainage depends on the permea bility and compressibility of the sand, drainage path and duration of shaking. The distribution of porewater pressure at time t is given by, _. _ E j_ <_z_ _J-\ E 8£vd ot 5z W 5zJ ot 'w (3.15) 27 in which, kz is the permeability and yw tne unit weight of water. The term containing ey(j represents the internal generation of porewater pressure (Finn, et al. 1977). The stress-strain relationship outlined above can be very easily extended to non-uniform loading using an incrementally elastic analysis in time domain. The dynamic response can be computed for each time step by numerically solving equation (3.15) and the equation of motion, as explained by Finn, et al. (1977). 3.2 LABORATORY VERIFICATION OF EFFECTIVE STRESS RESPONSE ANALYSIS The basic assumptions made in the formulation of the stress-strain relationship presented above can be broadly categorized into two groups: Those made in the formulation of porewater pressure model and those made in modelling loading, unloading and reloading. The fundamental assumption that was made in the formulation of the porewater pressure model, was that the porewater pressures in an un-drained test can be obtained from volumetric strains measured in a drained test on a similar sample with same history of shear strain loading. This means that there should be a unique relationship between volumetric strains in drained tests and porewater pressures in undrained tests for a given sand at corresponding strain histories. Finn (1981) reported results of an extensive laboratory program to investigate this basic assumption. Volumetric strains were measured in drained Ottawa sand samples at relative densities Dr = 45% and 60% when subjected to constant strain cycles in a simple shear apparatus. Porewater pressures were also measured in undrained cyclic tests at the 28 same relative densities and initial effective confining pressures. Volumetric strains EV(J are shown plotted against porewater pressure ratios U/o^. in Fig. 3.4 for D„ = 45%. Each point on the curve represents a set of values of ey(j and U/avo for a given number of cycles with equal cyclic strain amplitudes. It can be noticed that there is a slight difference in the applied shear strain amplitudes. But these small deviations are not important. The data indicate a unique relationship between volumetric strain and porewater pressure ratios. The slope of this curve normalized with respect to confining pressure will give the rebound modulus E_. Martin, et al. (1975) suggested that E_ can be evaluated from the unloading curve in an oedometer, under static conditions. But Finn (1981) showed that the rebound modulus measured in the oedometer is higher than the modulus computed from the slope of the curve shown in Fig. 3.4. He used the E„ values computed from the slope of the curve shown in Fig. 3.4 to verify the porewater pressure model. Finn (1981) maintained that the the strain hardening effect (equation. 3.14) should not be included in the stress-strain relationship when predicting the behaviour of sands under undrained conditions. This is because net volumetric strains do not occur during undrained conditions. If drainage is allowed to occur, then the effects of strain hardening should be included. The porewater pressure model coupled with the stress-strain relationship can be employed to predict liquefaction strength curves. The strength curve plots of the cyclic shear stress ratio t/ovo versus the number of cycles to cause initial liquefaction, for normally consolidated (0CR=1) and over consolidated sands, obtained analytically o b D •> M IO 0) V-1.00 -•> o» c >nfin 0.75 -u o |c 050 -0 ^_ 3 to • 025 -CL • fc. O 0. 0 L T T T T T Sand type : Ottawa sand (C-109) a' s 200 kN/m2 . Relative density * 45% vo Legend Shear strain amplitudes Drained Undrained 0.056 % 0.100 % 0.210 % 0300 % o 0056% • 0.100% * 0.200% A 0314% ) 02 0.4 0.6 08 1.0 12 1.4 Volumetric strain in percent, € y(j % Fig. 3.4. Relationship Between Volumetric Strains and Porewater Pressures in Constant Strain Cyclic Simple Shear Tests, D_ = 45%. 1.6 18 ho 30 and experimentally, are shown in Fig. 3.5. The experimental curve was obtained from undrained constant volume cyclic simple shear tests. The comparison between the computed and measured liquefaction strength curves is very good. This means that the assumptions made in the formulation of the non-linear hysteretic effective stress-strain relation ship are valid. But, liquefaction resistance curve prediction is an extreme case. Finn (1981) used this effective stress model to predict porewater pressure development during undrained tests when subjected to constant cyclic shear stress in a simple shear apparatus. This constant cyclic shear stress loading results in an irregular strain history as the porewater pressures develop. Further the model parameters (i = 1,4) used were obtained from constant cyclic shear strain tests. Typical results obtained from two undrained tests are shown in Fig. 3.6. The agreement between measured and computed porewater pressures is remarkably good, indicating that all the assumptions made in the porewater pressure model and stress-strain relationship are reasonable. 3.3 FIELD VERIFICATION OF EFFECTIVE STRESS RESPONSE ANALYSIS A unique opportunity to investigate the capability of the one-dimensional effective stress response analysis was provided recently when data became available on the dynamic response of an artificial island in Tokyo Bay to the Mid-Chiba earthquake of 1980. Owi Island No.1 is an artificial island located on the west side of Tokyo Bay. It was constructed with materials dredged from the nearby sea. A test site at the south end of the island is instrumented to record porewater pressures and ground accelerations during earthquakes. 0.40 o - > b 0.30 o o w S 0.20 w •> O *•= 0.10 o Sand type ; Ottawa sand (C-109) orv'0 «200 kN/m2 , Relative density » 45-47% Analytical curve O,«,A,A Experimental data 3 6 10 30 60 100 Number of cycles for initial liquefaction, NL Fig. 3.5. Cyclic Stress Ratio vs. Number of Cycles for Initial Liquefaction for Various Sands. Sand type : Ottawa sand (C-109) o-v'0 « 200 kN/m2 , Relative density = 45% T/cr' =0.074 T/crV*a065 vo / / o Experimental curve Analytical curve 6 10 30 Number of cycles , N 60 100 200 Fig. 3.6. Predicted and Measured Porewater Pressures in Constant Stress Cyclic Simple Shear Tests, Dr = 45%. 33 Porewater pressures are recorded by piezometers installed at depths of 6m and 14m. The Mid-Chiba earthquake, with magnitude M = 6.1, shook the Tokyo Bay area on September 25, 1980. Finn, et al. (1982) computed the response of Owi Island No.1 to the Mid-Chiba earthquake using a one-dimensional effective stress response analysis. The first 10 sees. of the recorded ground accelerations are shown in Fig. 3.7(a). During the first 4 sees, very low accelerations occurred. Significant accelerations developed between 4 and 6 sees., and thereafter only low level excitation was recorded. The ground accelerations computed by Finn, et al. (1982) are shown in Fig. 3.7(b). Except for some minor differences between 8 -10 sees, range the computed recording was very similar to the recorded motions. The porewater pressures recorded at the 6m depth on Owi Island No. 1 are shown in Fig. 3.8(a). The recorded porewater pressure has two components: transient and residual. The transient porewater pressures are instantaneous response of porewater to changes in total applied stresses and residual porewater pressures occur due to plastic volume changes. The one-dimensional response analysis used by Finn, et al. (1982) computes the residual porewater pressure component and is shown in Fig. 3.8(b). Comparison between recorded and computed porewater pressures is very good. 3.4. POREWATER PRESSURE MODEL IN PRACTICE To apply the porewater pressure model in dynamic effective analyses, 7 constants must be known; four C. (i = 1,4) constants to CJ a.: s (a) o.o i.o i i— 4.0 1.0 TIME 0.0 10.0 (b) i i— 4.0 (.0 TIME 0.0 —I 10.0 Fig. 3.7. Measured (a) and Computed (b) Ground Accelerations (Acc. in ft/sec2. Time in Sees). 8-;HWv(M*w r-TIME '•• o.o 10.0 Fig. 3.8. Measured (a) and Computed (b) Porewater Pressure at a Depth of 6m. (Porewater Pressures in lb/ft2, Time in Sees). LO 35 compute incremental volumetric strain and 3 constants Kr, m and n to represent rebound characteristics. Cyclic simple shear apparatus has to be used to obtain these constants. A number of laboratories still do not have simple shear apparatus to do these tests. Over the years a procedure has evolved from a number of practical experiences by which the direct measurement of these constants can be avoided (Finn, et al. 1982). This is done by modifying the model parameters such that it will match with the experimental liquefaction strength curve and give the right rate of porewater pressure generation. The liquefaction strength curve and the rate of porewater pressure development can be experimentally obtained by doing cyclic triaxial tests or cyclic simple shear tests on field samples. A study of a number of trial analyses to predict the undrained behaviour of samples in simple shear has revealed the following: a) The shape of the liquefaction resistance curve is sensitive to the constants C^, especially C^. b) The variation of Kr shifts the liquefaction potential up or down without changing the shape appreciably. In reality the shape of the liquefaction resistance curve for a number of sands is similar and the values of C^ given in the literature give the shape of typical liquefaction potential curves. In practice a trial and error procedure is adopted to get values for the model constants. The procedure is outlined below: 36 1) Performing a number of analyses by varying K„, and select the value for K_ such that the computed liquefaction resistance curve matches the experimental liquefaction resistance curve. 2) For this selected K_ value, calculate the development of porewater pressure with number of cycles and compare with the laboratory porewater pressure curve. 3) If these porewater pressure curves are not similar, alter and repeat the analysis. It should be noted that is the only parameter that is used in the calculation of incremental volumetric strain in the first cycle. Therefore, estimates of C^, can be interpreted from the residual porewater pressure recorded in the first cycle. This type of trial and error procedure can be employed to obtain relevant model constants such that the corresponding porewater pressure development and liquefaction resistance curves are sufficiently close to the ones observed in the laboratory. 3.5 DISCUSSION In the response analysis of horizontally layered deposits subjected to horizontal accelerations, a shear beam type of deformation pattern is assumed in the deposit. Therefore, only the stress-strain relationship in shear is required. The tangent shear modulus is used as the elastic parameter in the incrementally elastic response analysis proposed by Finn, et al. (1977). To extend their model to two-dimensions, two elastic parameters are required. A detailed 37 description of the extension of the one-dimensional stress-strain relationship to two-dimensions is discussed in Chapter 4. It has been observed in the laboratory that the presence of static shear stress affects the porewater pressure response of samples subjected to cyclic loading (Finn, et al. 1978; Vaid, et al. 1979). The porewater pressure model of Finn, et al. (1977) is strictly applicable to one-dimensional deposits, where static shear stress is zero. Therefore, in extending their model to two-dimensions, the influence of static shear must be accounted for. 38 CHAPTER 4 TWO DIMENSIONAL STATIC ANALYSIS OF SOIL STRUCTURES A static response analysis to compute in-situ effective stresses is necessary because the dynamic soil properties, such as strength and stiffness, depend on in-situ effective streses. A number of satisfactory incremental elastic methods, which model the construction sequence of the soil structures, are already available (Ozawa, et al, 1973; Duncan, et al, 1978, Byrne, et al, 1982). The static analysis presented in this theis is based mainly on the methods proposed by these authors. The method proposed in this thesis uses a consistent set of material parameters in both the static and dynamic analyses; procedures also have been incorporated to apply correction forces during the application of the load increments. 4.1 STRESS-STRAIN RELATIONSHIP A number of stress-strain relations have been proposed in the computation of in-situ static stresses in soil deposits. They can be divided broadly into linear, bilinear, elasto-plastic, visco-plastic and non-linear models. Some of these models are very complex and even for simple monotonic types of loading are expensive to use in computational schemes. 39 Two elastic constants are required for any two-dimensional isotropic, elastic or incrementally elastic analysis. Tangent shear modulus Gt and tangent bulk modulus Bt were selected as the elastic constants. In selecting these elastic constants, considerations were also given to the fact that the stress-strain formulation proposed here has to be extended to model dynamic loading conditions. It will be shown in Chapter 5, that the selection of these parameters greatly reduces the amount of computation time in the dynamic analysis. The stress-strain model proposed here, like almost all other static stress-strain models for soils, can model only saturated soils under fully drained or undrained conditions, and dry soils. The parameters selected to model the stress-strain behaviour should be based on test results which represent as closely as possible the loading conditions that exist in the field. For example, in the analysis of long term stability of earth dams, one should chose parameters from drained test results. A description of the stress-strain model, and the selection of relevant parameters for the model are discussed in detail, in this chapter. 4.1.1 Reasons for Selecting Gt and B_t In general, strain in an isotropic, homogeneous, linear elastic medium can be divided into two components: volumetric strain and deviatoric strain. The volumetric strain is related to mean normal stress through the bulk modulus. The deviatoric strain is related to deviatoric stress through the shear modulus. These two independent material moduli can be evaluated independently by applying uniform changes in AO corresponding stresses. Therefore, by selecting tangent bulk modulus and tangent shear modulus as two independent elastic constants, better controls on stresses and strains can be imposed. 4.1.2 Hyperbolic Shear Stress-Strain Relationship A number of researchers have used a hyperbolic stress-strain relationship to predict the behaviour of a soil deposit (Konder, et al. 1963). The hyperbolic relationship between i and y is given in terms of Gmax and xmax as> G y max ' T = Gmax Ml max (l + W m) in which, T,y = are the shear stress and shear strain G_ax = tangent shear modulus as y+0 Tmax = ultlmate shear strength 4.1.2.1 Estimation of Gmax Experimental data have shown that for sands and silts under drained conditions, Gmax " fK» e» 0CR) <4-2) in which, 41 o"m = current effective mean normal stress e = void ratio OCR = over consolidation ratio here OCR is defined as: _ Maximum past ma.jor principal stress Current major principal stress The following non-dimensional equation is widely used for Gmax' Gmax * KG Pa t^j"2 (4-3) in which, KQ = a non-dimensional constant for a given soil. Pfl = atmospheric pressure. The value of KQ depends mainly on void ratio or relative density of the soil, grain contact characteristics such as angularity and roughness of the soil particles etc., and also on previous loading history. An equation similar to (4.3) has been proposed by Hardin, et al. (1972) and Seed, et al. (1970) for the computation of GMO„ for sandy soils for dynamic analyses. The equation given by Hardin, et al, (1972) includes the effect of previous stress history. They proposed, G = K_ P fm/P ) (OCR) (4.4) max G a ^ a; v ' 42 in which the exponent n depends on the plastic index of the soil (Hardin, et al. 1972). Values of n are given in Table 4.1. Table 4.1 Variation of Exponent, r\ with Plastic Index, PI PI% 0 0 20 0.18 40 0.30 60 0.41 80 0.48 >100 0.5 For normally consolidated non-plastic soils under drained condition typical values for KQ varies between 200 and 800 (Byrne, 1979). For clayey soils under undrained conditions, GMAX can be related to the undrained strength SU through an equation, GMAX=KSU (4.5) where K is a constant for a given clay. 43 4.1.2.2 Estimation of Tmax For soils under drained conditions, tmax> is the maximum i shear stress that can be applied by keeping axial stresses ox, and i av at their respective values, where the x-axis is taken horizontal and the y-axis vertical for convenience. In principle, this is similar to the estimation of maximum deviatoric stress (o"j)mQv in the analysis presented by Kulhawy, et al. (1969). They estimate (a(j)max assuming that the minor principle stress remains constant. Let us consider a case where the initial vertical, horizontal i i and shear stresses are av, av, and -r respectively (Fig. y x xy 4.1). Fig. 4.2, shows the corresponding Mohr circle diagram and the Mohr envelope. The Mohr envelope is defined by c' and 0' . The points L and M represent the initial stress state. The application of shear stress will increase the size of the Mohr circle and the largest Mohr circle is the one that is tangent to the Mohr envelope. i OA = OP radius of the largest circle keeping ay and i ax constant and it is given by, i i OP = r c fgx + qy)i sin0' (4.6) "•tan0' 2 ' From triangle ABO, AB - Sax plane = maximum shear stresses on the horizontal Fig. 4.2. Mohr Circle Diagram . 45 = /(OA2 - OB2) and so, f-*}2sin20' - H^}2] 1/2 (4.7) x + max This equation reduces to the equation presented by Hardin and Drnevich (1972) if av and a are replaced by KQa ou vo and a, vo* The estimation of T max for a soil element under undrained conditions can be made based on standard field tests, laboratory tests or may be based on estimation of insitu effective stress conditions. 4.1.2.3 Influence of Over-Consolidation Compaction is generally used to obtain a certain density in dam construction. So some parts of a dam are over consolidated due to the compaction pressure. The effect of over consolidation on G__„ is Illcl A. already shown in the equation (4.4). Over consolidation has an influence on the value of xmax also. A typical Mohr envelope for plastic soils would look like that in Fig. 4.3. Different c1 and 0' values should be used depending on whether the soil is in NC state or OC state. 4.1.2.4 Effects of Unloading In geotechnical problems which involve excavation or reduction 46 co CO Q> CO u cd u CO O-C Region N«C Region CO CO 0) tH *J CO M Rt CO Normal Stress Fig. 4.3. Mohr Envelope. Slope G2-(Loadin Slope G4 (Unloading^ Gx> G2 Shear Strain Fig. 4*4 . Effects of Loading and Unloading, 47 in applied load, some soil elements will experience unloading. The modulus, G^, during unloading from a strain, y> as shown in Fig. 4.4, is higher than the modulus, G2, corresponding to loading from the same strain level if the mean normal stress remains a constant. Unloading and reloading can be modelled using the procedures presented in Chapter 3. However, if the strain ranges of interest are small, the difference between G2 and is not large. Under these circumstances changes in modulus need not be modelled. 4.1.3 Tangent Bulk Modulus Bt It is assumed that the tangent bulk modulus (Bt) is elastic for any soil under drained conditions and it is a function of the current t mean normal stress am only. Duncan, et al. (1978), suggested that, B = K, P A}n (4*8) t b a LP ' where, = Bulk Modulus constant n* = exponent Typical values of bulk modulus constant vary between 300-1000 and the exponent varies between 0.3 and 0.6. depends mainly on relative density of the soil, soil type and previous loading history. In elastic (or incrementally elastic) analysis in isotropic and homogeneous materials, a change in shear stress with constant mean normal stress will not result in any volume change. But, soils under 48 constant mean-normal stress exhibit volume change when subjected to shearing stress. In the analysis presented herein, allowance has been made to include the volumetric strains that occur due to shear stresses. A detailed description of how this is done is given in Section 4.4. It should be noted that the bulk modulus constant defined here includes only the effect of mean normal stress. Therefore, for the estimation of K^, isotropic consolidation tests performed in triaxial test equipment must be used. Duncan, et al. (1978) described a procedure for determining from conventional triaxial test data in which the mean-normal stress is not held constant. Values of determined in this manner must be considered approximate. 4.2 PHYSICAL MODELLING The domain of interest is assumed to be an assembly of a finite number of elements, connected at the nodal points. The formulation of the finite element equations including the effect of porewater pressure (Christian, et al, 1970) is presented in Appendix I. The equations are, {P} = [KT] {A} + [K*] {U} (4.9) where, {p} = global column vector of incremental applied loads [KFC] = global tangent stiffness matrix {A} = global column vector of incremental displacements [K*] = global stiffness matrix defined in the Appendix I 49 {u} = column vector of incremental porewater pressures in the elements. The tangent element stiffness matrix [kt] depends on two factors; the tangent moduli and the shape functions adopted in the finite element formulation. The shape functions give the variation of displacements within an element in terms of the nodal displacements. The simplest shape function, which assumes a linear variation of displacements, gives constant strain within a triangular element. But experience has shown that the results obtained from such elements do not predict stresses and strains accurately. Therefore, quadrilateral elements which have a linear strain variation within an element are used. For soil structures such as dams, layered deposits etc., elements of arbitrary quadrilateral shape are very appropriate because they are fairly simple and can be used to model the geometry of these soil structures quite accurately. The element stiffness matrix, [kt] for an isoparametric quadrilateral element is given in Appendix 1. 4.3. SIMULATION OF CONSTRUCTION SEQUENCE Dams are constructed sequentially. Since the behaviour of dam materials are non-linear and stress path dependent, a realistic computation of stresses and strains requires that the construction sequence be modelled. An analysis based on single stage construction or gravity switch on, will give final stresses and strains different from those calculated by following the construction sequence. 50 4.3.1 Method of Analysis Fig. 4.5 shows a schematic representation of the sequential procedure involved in dam construction. There may be pre-existing elements on which subsequent layers will be placed. The construction lifts are compacted until the required density is obtained. This type of layer by layer construction procedure is carried out until the required dimensions of the dam are obtained. In modelling the construction sequence, the incremental stresses, strains and deformations are computed for every new layer added. This is done by solving equation (4.9) for the incremental loads caused by placing a fresh layer. The final stresses, strains and displacements of the dam are simply the algebraic sum of all the incremental values computed for all the layers. 4.3.2 Incremental Porewater Pressure Static analysis can be carried out in a total or effective stress mode or a combination of both. In the combination mode, some elements may be in an effective stress mode and some may be in a total stress mode. When the effective stress principle is used in the finite element formulation, the porewater pressure term is introduced as shown in the right hand side of the equation (4.9). For the elements for which the total stress mode is assumed to be applicable, the element stiffness matrix [kt] is based on a total stress-strain relationship. Furthermore, the porewater pressure, uQ, in these elements should be set to zero. However, if the effective stress 52 mode is selected for some elements, the element stiffness for these elements should be based on an effective stress-strain relationship and a value for uQ is required. It should be remembered that the reason for the static analysis is to estimate the in-situ static condition which is a "long term" condition. In view of this, the global porewater pressure vector {Ti] used in the equation (4.9) should correspond to the long term value. Estimates of uQ for the elements can be made using a number of methods such as hydraulic model tests, electrical analogy etc. Measured porewater pressures in the field also may be used. The matrix {IT} , formulated using element porewater pressures uQ can now be used in equation (4.9) to compute stresses and displacements. 4.3.3 Computation of Incremental Stresses and Strains Strains given by the finite element analysis are a measure of changes in shape of the elements from some reference state. It is assumed that the condition of newly placed elements after they have settled under their own weight is the reference state (Ozawa, et al. 1973). The total strains are obtained by adding incremental strains caused by the construction layers about this reference state. An incremental elastic analysis can be carried out in a number of ways (Desai and Abel, 1972). The approach adopted here is shown schematically in Fig. 4.6. Estimates of the increments in stresses and strains due to a load increment are determined using moduli values corresponding to stress-strain level before the load increment was applied. New moduli corresponding to the average of the strains before Initial Estimates or Computed Values ^o 't t 1—r Layer Loading (incremental) Estimated New Moduli Stresses and Strains Calculate K G'T, B[ from W - W + M —1> Average Strains M = KI + M Moduli and Stresses V. i Initial Estimates or computed Values with modified G'T and B£ Layer Loading (incremental) Final Stresses and Strains {e} = {e} + {AE} {a} = {a0} + {Aa} Next Layer Loading Fig. 4.6. Schematic Diagram Showing the Incremental Analytical Procedure 54 the load increment is applied and the strains computed after the increment are used to compute more correct incremental stresses and strains for the same load increment. These incremental stresses and strains are added to initial stresses and strains to obtain the initial condition for the next load increment. Recall from section 4.1.2, that the relationship between shear stress and shear strain is assumed to be hyperbolic. Therefore, the shear stresses are computed using the shear strains at the end of the load increment. In doing this, as pointed out by Desai and Abel (1972), equilibrium is not necessarily satisfied. Under these circumstances equilibrium correction forces may be applied to satisfy equilibrium condition. In the method adopted here, correction forces that correspond to the changes in the shear stresses computed using the procedure outlined in Appendix I, are applied to the next load increment. Before placing a fresh layer the stress condition in previously placed elements (pre-existing) are known. Therefore, in the initial analysis for the load increment, moduli for the pre-existing elements are known. However, for the freshly placed elements, moduli must be based on estimated stresses in these elements. Ozawa, et al. (1973) suggested that the stresses can be estimated using the equations, and y ' s T = 0.5 v d sin <* lxy is o (4.10) in which d is the depth of center of gravity of the element from the top surface, <=0 is the slope of the top surface and ys is the unit weight 55 of soil. For totally submerged elements Ys should be replaced by the submerged unit weight y'« 4.4. SHEAR-VOLUME COUPLING The tangent bulk modulus Bfc defined in section 4.1.3, relates an increment in volumetric strain, Ae^, to an increment in effective i mean normal stress, Aam« But in soils volumetric strains occur also due to changes in shear stresses. This additional volumetric strain must be accounted for in any realistic modelling of soil behaviour. The characteristic drained behaviour of initially loose and dense sand samples in a simple shear apparatus is shown in Fig. 4.7. Initially for small shear strains y, both the loose and the dense samples undergo volume reduction. But later, over a considerable range of strain, they exhibit volume expansion (dilation). For both samples in the dilation range, the variation of volumetric strain ev vs Y ^s linear initially and then ev approach fixed values at very high strain levels. The region of interest in the ey vs Y plot in typical geotechnical problems would be the initial compaction region and the linear dilation region. The rate of volume change in the linear dilation region, is larger for the dense sand than for the loose sand. Hansen, (1958) suggested using dilation angle vQ to characterize the dilation rate. He defined v as, sinv = o dy~ = tanPo (4.11) 56 Fig. 4.7. (b). Typical Plots for £ Vs T for Dense and Loose Sands . 57 where 8Q is the slope angle. The negative sign is introduced since compressive volumetric strain is considered to be positive. For a given type of sand the angle vQ was found to be a function of the relative density and confining pressure. The dilation angle vQ increases with the relative density of the soil. This was clearly shown by Vaid, et al. (1981). They performed drained simple shear tests with constant vertical stress ay0 = 200 kPa, (Fig. 4.8) on Ottawa sand, (C-109), at various relative densities. The dilation angle, which is the slope of the linear dilation portion of the plot ev and y is found to increase with the relative density of the soil. The dilation is also a function of the mean normal stress level. This was shown by Lee (1965) who performed drained triaxial tests on dense Sacramento River sand samples of constant Dr = 100%. Fig. 4.9 shows a series of tests with consolidation pressures varying from 0.1 MPa to 13.7 MPa. Several important features of the test data can be noted in Fig. 4.9. Firstly, dense samples at high consolidation pressures behave like loose samples; secondly, failure in terms of maximum principal stress ratio occurs at increasing strain levels as the consolidation pressure increases; and thirdly, the dilation angle decreases and becomes negative (compaction) with increasing consolidation pressure. Fig. 4.10 shows the variation of the dilation angle vQ with mean normal stress for a number of sands which were at an initial relative density of 80 percent. It is interesting to observe that the variation of vQ lies within a narrow band (Robertson, 1982) and also the variation is linear with logarithm of mean normal stress. Based on the experimental data presented above, the following approximation for analytical purposes can be made for medium dense and 58 Fig. 4.8. Stress-Strain Behaviour of Ottawa Sand in Drained Simple Shear. (After Vaid et. al., 1981) Fig. 4.9 Typical Drained Triaxial Test Results on Dense Sacramento River Sand. (a) Principal Stress Ratio Vs Axial Strain (b) Volumetric Strain Vs Axial Strain (After Lee, 1965) ~ 20 V « TJ Ul -I o z < z o 10 ^'/(degree*) Dj.- 802 32 32 35 •37 33 35 x V Chattahoochee Sand Mol Sond Monterey Sond Glacial Sand SATAF Leigh ton Buzzard Sond Vesie and Clough 1968 OeBeer 1965 ViKet ond Mitchell 1981 Hirshfteld ond Pbuloe 1963 Boldi et ol. 1981 Cole 1967 0.1 0.5 I 5 10 50 100 .1 500 1000 MEAN NORMAL STRESS, crm Kg/cm3 Fig. 4.10. Variation of Dilation Angle, with Mean Normal Stress for Various Sands. (After Robertson, 1982). ON o 61 dense sands: the volumetric strain due to shear stresses (ev) is negligibly small (Varadarajan, et al. 1980, Byrne, et al. 1982) at low shear strain levels and above this level it varies linearly with y. This means that the plot of ey vs y can be idealised as in Fig. 4.11, where yQ is the shear strain above which the volumetric strain due to shear stress is important. It should be noted that the value of vQ should be modified for the changes in mean normal stress according to some variation such as shown in Fig. 4.10. 4.4.1 Analytical Formulation There are a number of ways of modelling shear-volume coupling. One is to modify the elasticity matrix I) (Appendix I) such that Aex and Aey depend also on shear stress increment (Verruijt, 1977). But this type of approach will give rise to an unsymmetrical stiffness matrix, which unduly complicates the computations. A simpler way is to keep the I) matrix as it is and to incorporate the volume change the same way as the temperature variations are analysed in structural mechanics (Zienkiewicz, et al. 1967; Byrne, 1979, Byrne, et al, 1982). This is accomplished in the following manner. a) The incremental shear strains in all elements are computed for the increment in load, neglecting shear-volume coupling. b) vQ can be estimated for the new mean normal stress, using Fig. 4.10, and then Aey is computed using equation 4.11. c) The volumetric strain Aev then is split into Aex and ^ft. Here it is assumed that, 63 Aed = Aev = 0*5 Aev* Let us define a strain x y v T vector Ae such that its components are the estimated "0 dilational strains. Then AeT is given by; JAe^, ~0 Aey, 0}. d) The nodal forces corresponding to this strain vector Ae^ can be computed as, /// I*?. A lo dv (4.12) V (see Appendix I) Now these forces can be added to the applied incremental load in a) and new strains and stresses can be computed. For computing incremental stresses, the following equation should be used, Aa = D_(Ae - A_eQ) (4.13) where, A£ = strain vector computed for the modified applied load. e) Now steps b) ->• d) can be carried out until convergence occurs in stress and strain increments under the applied incremental load. 4.5 INTERFACE REPRESENTATION It may be necessary to allow relative displacement to occur at the interface between two finite elements to model slip surfaces in the 64 field. Slip elements, which are sometimes referred to as elements of zero thickness, can be used to model this relative displacement. Slip elements can be assumed to be placed along the boundaries between the two-dimensional elements representing soil and structural elements or wherever it is anticipated that relative movements or separation between elements may occur. The slip is assumed to occur only along this direction and this occurs when the shearing forces in the slip element exceeds the shear strength at the interface. Goodman, et al. (1968) have developed a two-dimensional slip element with eight degrees of freedom to represent joint and fault behaviour in rock mechanics problems. Fig. 4.12 shows a slip element with nodes I,J,K and L, in global and element axes. The forces at any point in a slip element are the shear force f„ and the normal force f„ expressed per unit area of the element. The force-displacement relationship is assumed to be: (4.14) n n n where Kg, Kn = joint stiffness per unit length in shear and normal directions respectively. wg, wn = Shear and normal displacement at the point of interest. The definition of unit joint stiffness needs clarification. Imagine a direct shear test being performed along an interface element of unit thickness. At first, when a normal force is applied, the element 65 X Shear/Normal Displacementa v^,*^ Fig. 4.13. Plot of Typical Shear/Normal Stress vs Shear/Normal Displacement. 66 shortens as the asperities in the joint deform. A typical plot of normal deformation at the joint and the force applied per unit length is shown in Fig. 4.13. For analytical purposes the relationship can be approximated to a straight line and the slope is given by 1^. Similar tests can be performed in the tangential direction and a plot between f„ and w„ can be obtained. The slope of this curve will give K„. Using the equation 4.14 and also assuming a linear variation of displacement within the joint element, a stiffness matrix K _ can be —-oil obtained in local or element co-ordinates. This stiffness matrix relates the nodal forces and the nodal displacements. The displacement vector here is simply, uT = {uIf V-p Uj, vj, uK, vK, uL, vL} It has been shown in Appendix II that Ksn is given by, 0 Ks 0 ~Ks 0 -2KS 0 2Kn 0 0 "Kn 0 -2K, 2Ks 0 "2Ks 0 "Ks 0 2Kn 0 -2*n 0 ~Kn 2Ks 0 Ks 0 sym 2Kn 0 2Ks Kn 0 2Kn n (4.15) 67 where, L is the length of the element. To get the stiffness matrix in global co-ordinates a simple transformation is used, where T= transformation matrix containing direction cosines and is given by; (4.16b) in which, cosec sin<* ., .r N PQ = r o o-i (4.16c) '--sinoc cos« J o o and a = angle of inclination of the slip element with the horizontal. Even though this type of formulation does not include rotation of the element directly, this is taken into account since all 8 d.o.f. have been considered. The displacement field variation assumed within the slip element is consistent with the displacement field in an isoparametric quadrilateral finite element. Furthermore, both elements have the same degrees of freedom. Thus establishing a global stiffness matrix can be evaluated simply treating the slip element like any other element. 68 4.5.1 Method of Analysis of Slip Elements It was assumed that the tangential stress-displacement relationship (fg vs wg) in the slip element is elastic-perfectly plastic and the plastic region where slip occurs is defined by a simple Mohr-Coulomb type of yield criterion. For incrementally elastic analysis, the values of K„ and K_ should be kept constant until yield occurs. After the yield, K„ is set to a very small value. This small value can be viewed as the residual shear stiffness. The value of K^ is kept at its original value. But, if the normal force on the element is negative, meaning that the separation of the joint occurs, then the values Kg and IC^ should be set to a small value. To investigate whether yielding is possible or not, the shear and normal stress in the slip element should be determined. Since a linear variation of displacement is assumed within an element, the shear and normal stresses vary from point to point within the element. The average incremental values of shear stress Afg and normal stress Afn for a load increment can be estimated from equation (4.14) as, Afg = Ks (Awg) and Afn = 1^ (Awn) (4.17) where Aw. and Aw are incremental average shear and normal 5 XI displacement in the element respectively. Now expressions for Aw and Awn are, 69 Aws = <Autop>ave " (Aubottom>ave = (AuK+AuL)/2 - (AuI+AuJ)/2 (4.18) and Awn = (AvtQp)ave - (Avbottom)ave = (AvK+AvL)/2 - (AVl+Avj)/2 (4.19) From equations 4.17 to 4.19, Af„ and Af_ can now be written as, Afs = Ks [(AuK+AuL)/2 - (AuI+AuJ)/2] and Af~n = 1^ [(AvK+AvL)/2 - (AvI+AvJ)/2] (4.20) Total shear and normal stresses fg and f can be obtained by adding up the incremental stresses for all the load steps. Mohr-Coulomb failure criterion gives the shear strength fmax in the element at any time as, fmax = cs + fn tan< (4.21) where c and 0 are the cohesion and friction angle required to define s s the failure criterion. If fmax is greater than the absolute value of fg then the slip element nodes could separate and this is modelled as mentioned above, by reducing the Kg to a small residual value. The separation of a slip element is indicated by negative fn. It should be noted that all calculations for the computation of fg and f are performed in the local axes. Since the displacements from finite element analysis are given with respect to global axes, they 70 must be transformed to get displacements in the local axes by using the inverse of the transformation matrix T. 4.5.2 Factors that Influence Joint Parameters In the analytical model presented above, three distinct joint parameters were introduced. (1) K„, the unit stiffness along the element. (2) KR, the unit stiffness across the element. (3) shear strength, fmax defined by cg and 0g. These parameters model the behaviour of slip elements adequately. The values of K„, K_ and fmav will depend on (1) the S II 111 cL 2v surface roughness of the adjacent elements (2) shape and characteristics of the asperities, and (3) contact area ratio between the joint walls. Details on how these parameters can be obtained in the laboratory are given by Goodman, et al. (1968). 4.6 SELECTION OF SOIL PARAMETERS The process of obtaining representative values for soil properties is probably one of the difficult tasks in stress analysis. It should be emphasized that individual estimation of soil properties is not important. But the stress-strain variation given by the selected soil parameters should give the best fit to the observed laboratory behaviour of soil samples in the stress (or strain) range of interest. 71 4.6.1 Obtaining Shear Stress - Strain Parameters After deciding what drainage condition is likely to occur in the field, a laboratory test can be performed simulating the drainage conditions. For example to obtain parameters for an effective stress-strain relationship, a series of drained simple shear tests can be performed to obtain plots % vs y for various constant mean normal stress levels. Simple shear tests are ideal since the mean normal stress during the test remains reasonably constant. A simple, trial and error method can be employed to obtain values for Gmov and that fit these curves in the stress or strain ranges of interest. Then knowing the stress levels corresponding to a test, the best estimates for the parameters KQ, C' and 0' can be obtained easily. The effective stresses in simple shear at the beginning of the drained loading conditions can be assumed to be o"vo and KQ t t o"vo. Then the mean normal stress am is given by (1+2KD) aVQ/3 and this is in general assumed to remain constant during the test. But if conventional triaxial tests are performed on the samples, then Ojjj varies as the axial load varies, and therefore, it is not easy to obtain these parameters. If triaxial test data only are available then again the above procedure can be carried out by considering the shear i stress and the shear strain on the failure plane. However, since o"m increases during shearing the shear stress-strain curve obtained by this i procedure may be interpreted for a constant average am- The average i i o"m can be assumed to be the mean of om, corresponding to the beginning and the end of the test. 72 4.6.2 Obtaining Bulk Modulus Parameters The tangent bulk modulus Bt was assumed to be (equation 4.8) given by, da' a' n B^ " ~A = K, P {—) (4.22) t de b a KV ' v ' vm a Integrating this equation one gets, fK p (4.23) 1-n vm b a ' taking logarithms both sides, (l-n)log(o-'m) = log {KbP^_n(l-n)} + log (effl) (4.24) i.e., log(evm) = (l-n)log(a'm) - log {Kbp(1_nh-n)} The slope of the plot log (evm) vs log (am) obtained from drained isotropic triaxial test results will give (1-n), and from this n can be determined. Using the value of n, and the intercept on log (e.yjjj) axis, can be calculated. Now knowing the bulk modulus parameters, and n, the tangent bulk modulus can be computed at any given mean normal stress. 73 CHAPTER 5 TWO-DIMENSIONAL DYNAMIC ANALYSIS 5.1 FORMULATION OF THE PROBLEM The general dynamic equilibrium equations for a linear system at any time are given by (Clough and Penzien; 1975), [M] {X} + [C] {X} + [K] {X} = {P} (5.1) in which {x}, {£}> {x} = column vectors whose components X^, X^, and X^ give the relative acceleration, velocity and displacement with respect to the base motion respectively of a node, [M] = mass matrix [c] = damping matrix [K] = stiffness matrix {p} = inertia force vector, which is defined as, -[MJ{l}xb where, {i} is a vector with all components unity and X^ is the base acceleration. In two-dimensional problems the base acceleration may have two components: horizontal and vertical. If the it'1 equation in (5.1) is written for the horizontal direction then X^ is the base acceleration in the horizontal direction, and if the equation is for the 74 vertical direction then X^ is the base acceleration in the vertical direction. 5.1.1 Incremental Equilibrium Equations for Non-Linear Systems In the analysis of non-linear systems, the material properties change with time. An incrementally elastic approach has been adopted to model the non-linear material behaviour. Incremental dynamic equilibrium equations for any time interval, At, can be obtained by replacing the variables in equation (5.1) by incremental values. This leads to, [M] {AX} + [C] {AX} + [K] {AX} = {AP} (5.2) in which, [M], [C] and [K] are the mass, damping and stiffness matrices relevant to the time interval for which the above equations are written. It is always assumed that the mass matrix is constant. The mass matrix can be obtained by two ways: lumped mass method and consistent mass method. In the lumped mass method, the mass matrix is obtained by lumping the mass of a finite element equally at its nodes. The consistent mass matrix, is obtained using the same interpolation functions used in the finite element formulation. The lumped mass matrix is very simple to compute and it has only diagonal terms, whereas the consistent mass matrix is somewhat harder to compute and has off-diagonal terms. Even though the consistent mass method is more accurate, the presence of off-diagonal terms greatly increases the computation time required to solve the dynamic 75 equilibrium equations. For the accuracy level required in typical geotechnical problems, the lumped mass method is considered appropriate. In general, the damping matrix [c] and the stiffness matrix [K] which are introduced in the incremental equilibrium equation (5.2) depend on the distribution of velocity and displacement in the structure. Appropriate values for the time interval between any time t and t+At can be determined only by an iteration procedure, because the velocity and displacement at the end of a time increment depend on the initial stiffness and damping values. This type of iterative solution scheme for every time step is very expensive. Therefore, in practice tangent damping and tangent stiffness matrices which correspond to time t are used with appropriate corrections to the results. It will be explained later how correction forces can be introduced into the solution scheme. The stress-strain law used to determine the tangent stiffness matrix at time t, [KFC] is described in Section 5.2. The [c] matrix will be assumed to be a constant throughout the dynamic analysis and the procedure to evaluate this is presented in Section 5.3. With, [K] = [Kt]t and (5.3) [C] = [C] (5.4) the dynamic incremental equilibrium equations can be re-written as, 76 [M] {AX} + [C] {AX} + [Kt]t {AX} = {AP} (5.5) When computing response for a random loading history, equations (5.5) have to be solved for every time step. During the procedure [Kt]t and {AP} are updated. The step by step integration procedure proposed by Newmark (1959) or Wilson's 0-method (Wilson, et al. 1973) have been adopted to integrate the equations. These procedures are described in Appendix III. 5.1.2. Correction Forces The numerical integration procedure is based on three significant assumptions: (1) the variation of acceleration within a time step is assumed to vary in some known fashion e.g. linear or constant (2) the damping and stiffness properties remain constant during a time step and (3) the response at time t + At can be evaluated from the known response at time t. In general neither of these assumptions is entirely correct, even though the errors may be small when the time step is short. If errors accumulate from step to step gross errors and even solution instabilities may occur. These problems can be avoided by imposing the condition of global equilibrium at each step of the analysis. The global equilibrium equations at time t in terms of all force components are, {flit + {fDlt + {fslt = {Plt (5.6) in which {^i}t» I^D^t' and representing inertia, damping, and mass system at any time t and {p}t t. Since the mass matrix and constant during dynamic analysis, by, {fg}t are the column vectors spring forces acting on the discrete is the inertia force vector at time the damping matrix were assumed to be {fj}t and {felt are simply given {flit = [M] (X}t and (5.7) {fD}t = [C] {X}t (5.8) The spring forces {^slt» can ^e computed by expressing element stresses in terms of nodal forces, applied at the nodes of the elements. The nodal forces that correspond to the dynamic stresses in an element {f^elem* ^rom Appendix I are, in which _B is the matrix that relates strain to nodal displacements. In this manner nodal forces for all the finite elements can be computed and 78 the vector sum of all these nodal forces will give the global vector {fS}f If the solutions obtained at time t are accurate then the right and left hand sides of the equation (5.6) will be identical. But, in general it will not be so. The correction force vector {pcorr} is given by, (Pcorr} = lp}t " {fllt " {fD}t " {fsk <5-10) From equations (5.6) to (5.9), equation (5.10) can be re-written as, {Pcorr> - ^t " M Wt - [C] {*} - J /// BC Oj dv all elements (5.11) To impose equilibrium, the correction force vector {PCOrr^ can ^e added to the incremental equilibrium equations formulated at time t. This is accomplished by adding {pcorr} to the right hand side of the equation (5.5). 5.2 DYNAMIC STRESS-STRAIN RELATIONSHIP In the proposed incrementally elastic analysis in the time domain, the tangent shear modulus Gt and tangent bulk modulus Bt were selected as the two "elastic parameters". Some reasons for selecting Gt and Bt for static analysis were presented in Chapter 4. There is a further very important reason for adopting these parameters in the dynamic 79 analysis. In dynamic analysis, the moduli have to be changed for every time step. This means that the element stiffness matrix has to be re-formulated each time. This time consuming procedure can be simplified somewhat if Gt and Bt are used as the "elastic" constants. The elasticity matrix D_ (Appendix I) under plane strain conditions, is given by, D = Bt + 3 Gt Bt " 3 Gt 0 = B. 1 0 1 0 Bt - 3 Gt Bt + 3 Gt 0 0 + G. 0 G 4. 3 -2 3 0 -2 3 4. 3 0 (5.12) (5.13) Bt%+ Gt*2 (5.14) where and Q2 are two constant matrices. From the formulation of stiffness matrix presented in Appendix I, the element stiffness matrix is, k=J//B D B dv V (5.15) Now substituting for D from (5.14) the equation (5.15) can be rewritten as, 80 k. = B J/J ^ £ B_ dv + G /// B1 L B dv (5.16) V V i.e. kfc - Bt Rx + Gt R2 (5>17) where R != /// I* £L B dv V x (5.18) (5.19) The constant matricies R^ and R^ have to be computed only once. The current lct matrix may now be obtained by multiplying the constant matrices R, and R9 by the current values of Bt and Gt. 5.2.1 Volume Change Behaviour With regards to volume change behaviour during dynamic loading, the soil deposits can be divided into two basic groups. The first group includes soils which can undergo volume changes under the load increments induced by the base excitation and the second group includes soils which cannot. Saturated gravels and dry deposits belong to the first group of soils. Recall equation (4.8) which relates B,. to a', 81 nun 3 J (5.20) This equation may be used to compute Bt. This means that Bt has to be modified for every time step. However, it is known that the changes in mean normal stresses in the soil elements, due to seismic excitation is small and furthermore the volume change behaviour does not influence the response of soil structures significantly. Therefore, for simplicity, Bt may be kept constant throughout the dynamic analysis. An average Bt for elements can be evaluated based on in-situ mean normal stresses using equation (5.20). This is because the load pulses during seismic loading induce stresses such that their mean values are initial in-situ i stresses. It should be noted here that even if the changes in am are considered, for typical values of n the changes in Bt will be small. Laboratory results have revealed that, as plastic volume changes occur during cyclic loading, the soil samples harden leading to higher bulk modulus. This increase in bulk modulus due to strain hardening effect can be modelled the same way as the increase in maximum shear modulus was modelled by Finn, et al. (1977). This was accomplished by introducing hardening constants (equation 3.14). Loose saturated, sandy soils and saturated clays belong to the second group of soils. In saturated soils volume change can occur only by porewater drainage. Within the short duration of typical seismic excitation, the occurrence of appreciable amount of drainage is questionable in soils of low permeability. In view of this, the dynamic analysis proposed here assumes that no drainage occurs during the dynamic 82 loading. In saturated gravels appreciable residual porewater pressure does not develop because of its high permeability. For the second type of soils, to simulate the condition of no volume change, the bulk modulus is set to a very high value during dynamic loading conditions. 5.2.2 Dynamic Shear Stress-Strain Behaviour In the formulation of a complete dynamic effective stress-strain relationship the following basic aspects should be considered; 1) soil behaviour under initial loading, unloading and reloading phases. 2) residual porewater pressure generation and its effects. 5.2.2.1 Skeleton Curve for Dynamic Loading Under dynamic loading conditions the relationship between dynamic (or cyclic) shear stress, %c> and dynamic shear strain, yc» Is assumed to be hyperbolic. The hyperbolic relationship (equation 4.1) is defined by Gmax and -cmax. Seed et. al, (1970) proposed that maximum shear modulus, ^max ^or sandy soils is a function of effective mean normal stress only, and given by, (5.21) 83 in which, (K2)max is a constant which depends on the relative density for a given soil. Seed et. al, suggested that for sands (K2)max varies between 20 and 100. Hardin et. al, (1972) suggested that the ultimate shear strength Tmax can ^e calculated using insitu effective stresses and static shear strength parameters such as c' and 0' (equation 4.7). They pointed out that for the level of dynamic strain (yc<l%) induced by seismic loading, the hyperbolic curve in terms of static TMOV is in ci x satisfactory for dynamic loading. Unlike dynamic analyses in one-dimension, an initial static shear stress, TG, is present in the analyses in two-dimensions. The presence of is causes the available shear strength to be different, depending on the direction of shearing. A typical relationship between %c and Yc Is shown in Fig. 5.1. The available shear strength in the direction of initial static shear stress is (T .„ - T„) and in the mcix s opposite direction it is (xmax + xg)• In dynamic analysis, the current practice is to neglect the influence of xg on available shear strength and to assume that the %c vs yc relationship is symmetrical about both xQ and YC axes. 5.2.2.2. Unloading and Reloading All the basic assumptions used to model unloading and reloading phases in Chapter 3 are also adopted for dynamic analyses. However, the equations have been modified to reflect the effect of static shear on available shear strength. For example the equation for the unloading curve KL from K in Fig. 5.2 is given by, 84 Fig. 5.1. Dynamic Shear Stress-Strain Relationship: Skeleton Curve. Fig. 5.2. Dynamic Shear Stress-Strain Relationship: Unloading and Reloading. 85 u tK; i + G max Y-YK|/ maxl (5.22a) where tmaxi = Tmax + Ts» and '''K and ^YL are tne dynamic shear stress and strain corresponding to reversal point K. The equation for the reloading curve LM, is given by, (T " TL} = 1 + G Gmax (? - V max IY (5.22b) max2 where TMAX2 = TMAX - -c s , and Tl, yL are the dynamic shear stress and strain corresponding to reversal point L. 5.2.3. Modelling the Effects of Residual Porewater Pressure One of the important factors in seismic response studies of soil deposits comprised of saturated cohesionless materials is the influence of residual porewater pressure generation. The presence of residual porewater pressure reduces the resistance to deformation. The presence of initial static shear stress, T , affects the s generation of residual porewater pressure in loose saturated sands significantly (Finn, et al. 1978b, Vaid, et al. 1979). The porewater pressure model developed by Martin, et al. (1975), does not account for xg. The porewater pressure model adopted here is an extension of the model of Martin, et al. (1975) to account for the effects of xg. 86 5.2.3.1 The Behaviour of Samples with Ts_ The original procedure that accounts for the presence of x s was presented by Seed and Lee (1969). They hypothesized that the behaviour of an element in the field with an initial static stress ratio, x/a' = ocr and subjected to a shear stress history on its failure plane is identical to the behaviour of a representative laboratory sample consolidated such that it has the same initial stress ratio <*r on its failure plane and is subjected to the same shear stress history. In the field, for typical soil structures, the failure plane for earthquake excitation can be assumed to be the horizontal plane (Seed, et al. 1973). Then the field stress ratio is simply «r = Ts/o*y0, in which Oy0 is the initial vertical effective stress. If simple shear equipment is used to test representative samples then the initial vertical and shear load can be applied such that the proper ratio «r = t /o' is present on the horizontal plane. However if triaxial s y *J apparatus is used, the axial a[c and radial pressures 0*3 should be such that the initial effective stress ratio on a plane (45+0'/2) to the horizontal (failure plane) has the same ratio <rr. It can be deduced from the above hypothesis of Seed and Lee (1969), that the behaviour of horizontally layered deposits can be interpreted from isotropically consolidated triaxial (ICT) tests. Anisotropically consolidated triaxial (ACT) tests have to be used if Comprehensive laboratory results are available on the behaviour of cyclic triaxial samples under various consolidation conditions. On the contrary, only very limited simple shear test data are available. This is 87 especially true when xg * 0. Therefore, in this section only triaxial test data are used to explain the difference in behaviour of samples, with and without Tg on the failure plane. A number of researchers have studied the behaviour of soil samples subjected to ICT tests and ACT tests (Seed, et al, 1969; Finn, et al. 1978b; and Selig, et al. 1981). There are a number of basic differences between the behaviour of samples subjected to ICT tests and ACT tests. In typical cyclic triaxial tests on saturated cohesionless soils the porewater pressure and axial strains develop with increasing number of deviatoric load cycles. The porewater pressure (Ut) recorded at any time is the sum of residual porewater (or permanent) pressure (U) and the cyclic (or transient) porewater pressure (Uc) • The cyclic porewater pressure is an instantaneous increment of porewater pressure which is a function of current changes in the mean normal and shear stresses and the residual value is the porewater pressure due to plastic deformation. The residual porewater pressure may be recorded when the applied cyclic load is zero. The total porewater pressure is then given by (Fig. 5.3), Ut = U + Uc (5.23a) Similarly, the total axial strain at any time also can be separated into elastic and residual or plastic components, 88 8 (SI NUMBER OF CYCLES Fig. 5.3. Permanent and Cyclic Strains or Permanent and Cyclic Pore Pressure. COMPRESSION 0 50 K» ISO 200 2D I i i i i NO. OF CYCLES . N Fig. 5.4. Cyclic and Residual Behaviour of Pore Pressure and Axial Strain for ICT and ACT Tests (after Selig, et. al, 1981). 89 et = ep+ ec (5.23b) Results of an ICT test and an ACT test are presented in Fig. 5.4 to illustrate the differences in behaviour between these two tests (Selig, et al. 1981). Two samples of Ossterschelde sand with initial porosity of about 41.5% were consolidated isotropically (sample a) and anisotropically (sample b). The cyclic deviatoric stress for sample (a) was 0.3 kg/cm2 and for sample (b) was 0.45 kg/cm2. The axial strain and porewater pressure response of these two tests are shown over 200 cycles of stress. The axial strain record shown in Fig.5.4 for the ICT test is almost symmetrical about the X axis and the average value of et a small value. The peak values of axial strain occurs when the applied deviatoric load is maximum. This means that the peak values are the cyclic components of the axial strain. The residual or permanent component of the axial strain is small. The cyclic component remained relatively small until the residual porewater pressure developed to about 60% of the consolidation pressure and thereafter increased rapidly. At the end of the test, when the applied load was zero, the axial strain was found to be small. However, in the ACT test, the mean axial strain increased with number of cycles of loading and at the end of the test the axial strain was quite large. The residual or permanent strain in this test increased with number of cycles while the cyclic strain remained relatively small. In the ICT test, the porewater pressure increased more rapidly than for the ACT test (Fig. 5.4b). After about 150 cycles of loading, the porewater pressure ratio (U/a^) for the ICT test rapidly approached 100%, while the ratio for ACT test approached a limiting value. 90 The cyclic porewater pressures vary depending on the magnitude of deviatoric stress. The deviatoric load plays two roles. Firstly, the increase in deviatoric load results in an increase in porewater pressure. Secondly the deviatoric load causes shear stress on the failure plane, which may cause dilation of the sample. This will lead to a drop in porewater pressure. Therefore, the sum of these two effects will give the net cyclic component of the pressure. It is accepted that the behaviour of saturated sands are not significantly affected by the cyclic porewater pressure, Uc. Only residual porewater pressure affects the behaviour of saturated sands significantly. It is easy to understand the difference in behaviour between ICT tests and ACT tests with regards to the limit of porewater pressure generation if both are compared in effective stress space, such as a q, p' plot, where q is the principle stress difference given by (o^ - o"3) and p' is given by (a| + a-j)/2. Points L and M in Fig. 5.5 represent the initial stress state in the (q,p*) plot for the ICT and ACT tests. As the residual porewater pressure increases, the effective stress state at any time, when no deviatoric load is present in the triaxial sample, is given by, o{ = afc - U (5.24a) 0-3 = a3c - U (5.24bthen the corresponding p' and q are given by, a! + ol al + al , 1 3 lc 3c P 2 2 " (5.25a) Fig. 5.5. q vs p Plot for ICT and ACT Tests. 92 and q = a{c - o$c (5.25b) This means that if the effective stress paths are plotted for both tests when the applied deviatoric load is zero they will be parallel to the p' axis and will move towards the failure line. If a number of load cycles of sufficiently high cyclic stress ratio are applied, the effective stress path may move very close to the failure envelope. If further load increments are applied, the sample will behave such that stress states above OA and below OB cannot occur. This means that maximum residual porewater pressure values recorded will be such that the effective stress path will be on the failure line. Therefore, the maximum residual porewater pressure for any ACT test is b2 and for any ICT test it is &2 (Fig. 5.5). Using simple geometric principles it can be shown that (Chern, 1981; Chang, 1982), al U = b = "J2" {l + K - (K - l)/sin0*} , . max 2 2 1 cc ' (5.26) where Kc is the anisotropic consolidation ratio defined as Kc = a'lc/03c* When Kc = 1 (i.e. ICT test) Umax is simply given by cr^. It should be noted that the ACT test reported here does not experience reversal of shear stresses at any time. If the total stress path is below the positive q axis then reversal in shear stress occurs. It is easy to show that reversal will occur if, 93 °-dy > °3c <Kc ~ 1) (5.27) Finn and Byrne (1976) have shown that stress reversal is required if a condition of initial liquefaction is to be achieved in anisotropic tests. However, whether there is reversal or not, the anisotropically consolidated samples strain progressively during cyclic loading in contrast to isotropically consolidated samples. The rate of porewater pressure development in a given sample subjected to cyclic loading in a triaxial apparatus is mainly governed by relative density, consolidation pressures, and applied cyclic stress ratio, °'dy/'2a,3c* compare the porewater pressure development characteristics of an ACT test and an ICT test, the tests have to be performed under similar conditions. Finn, et al, (1978b) presented a plot of the porewater pressure ratio Um/cr3c vs N/N50 (Fig.5.6) in which Um is the maximum porewater pressure recorded at any time during a cycle, N is the number load cycles and N^Q is the number of cycles to cause the porewater pressure ratio of 50%. The figure clearly shows that the rate of porewater pressure development is different in ACT and ICT tests. The interpretation of these test results should be done with care, since N^Q is also a function of KC. In presenting laboratory results on dynamic properties of saturated cohesionless soil samples it is customary to provide liquefaction potential curves. These curves are a plot of cyclic stress ratio °'dy/2o'l3c versus number of cycles to liquefaction, N^. As explained above, in some ACT tests it is not possible to reach liquefaction defined as U/o^c = i. Therefore, the definition of uAu tubas ufao—mkm~ Number of Cycles to Liquefaction, NL Fig. 5.7. Strength Curves for (a) Ekofisk Sand with Dr = 85%; (b) Sand with Dr = 77% (after Rahman et. al, 1977). 95 has to be changed to cycles to cause U = Umax or cycles to cause some specified strain (Seed, et al. 1969). Often in practice double amplitude axial strain (peak to peak) of 5% is used. It can be seen from Fig. 5.4 that for the ICT test, a single amplitude strain of 2 1/2% is equivalent to 5% double amplitude. This is because the axial strain record is fairly symmetrical about the X axis. However, this is not true in ACT tests. Therefore, in presenting liquefaction potential curves for ACT tests it may be necessary to define as number of cycles for some maximum single amplitude strain or double amplitude or the condition where U = Umax« Typical plots of liquefaction potential curves for two dense samples with and without xg are given in Fig. 5.7 (Rahman, et al. 1977). Two observations can be made from these figures. Firstly, the liquefaction potential curves are very similar in shape. Secondly, the curves for ACT tests are above the curves for ICT test indicating higher resistance when tg is present. However, the second conclusion may not be true for loose samples (Vaid, et al. 1979). In presenting the differences between samples with and without T„, only triaxial results were used. Vaid, et al. (1979) carried out a s number of simple shear tests with and without % . The conclusions drawn by their investigation and results reported by Seed, (Seed, 1983) are also very similar to those presented above. 5.2.3.2 Residual Porewater Pressure Generation Model The porewater pressure model proposed by Martin, et al. (1975) was modified to include the effects of tg. With regards to generation of 96 porewater pressure, it has been documented in the preceeding section that the presence of Tg, has three basic effects. They are: the position of the liquefaction potential curve is altered, there is a limit to the amount of residual porewater pressure, and its rate of generation is different when x * 0. The attempts made to account for these effects s in the porewater pressure generation model are presented in this section. Recall the equation (3.5) which relates the incremental volumetric strain to cyclic shear strain amplitude, Ae , = C. (y - C0 e .) + . r— (5.28) vd 1 ' 2 vd Y + C/ £ J ' 4 vd In one-dimensional response analysis only shear strain components are present. But in an analysis in two-dimensions, vertical and horizontal strains also occur. However, if it is assumed that only the cyclic component of shear strain yXy contributes to Aevd, then the equation (5.28) can still be used. This assumption is quite reasonable since the major strain that occurs during seismic excitation in typical soil structures which have adequate static factor of safety is shear strain (Serff, et al. 1976), and furthermore, only cyclic components of shear strain are responsible for grain slip. Based on a number of shaking table tests, Pyke, et al. (1975), showed that all three components of acceleration contribute to volumetric strain. Analytical studies by Seed, et al. (1975) have also revealed that under multi-directional shaking porewater pressures build up faster than under unidirectional stress conditions, and that the shear stress ratio Td^avo to cause liquefaction under multi-directional shaking 97 conditions is about 10% less than that required under unidirectional shaking conditions. They suggested reducing the stress ratio Td/'o"vo obtained from simple shear test results by about 10% to account for this effect. A two-dimensional dynamic analysis can account only for the horizontal and the vertical components of acceleration which act on the plane of the soil structure. The effects of the third acceleration component can be accounted for by modifying either the incremental volumetric strain AeV£j or the dynamic shear stress ratio V0vo-The porewater pressure model of Martin, et al. (1975) was modified such that the liquefaction resistance curve and the generation rate of the porewater pressure matches the behaviour observed in the laboratory sample with a given xs/ay0' Here the liquefaction resistance curve is obtained by defining as number of cycles required to reach the condition of residual porewater pressure U = U ^, where TJ „ r r max' max is given by equation 5.26. Recall the equation for E , (1-m) E = ,, rn-m r mK (a ) r vo (5.29) The term Kr in equation (5.29) and constants C± through may be adjusted to model the laboratory behaviour. A proposed trial and error procedure to accomplish this has been outlined in Chapter 3. It is worthwhile noting that the liquefaction potential curves generated for various Kr values are also very similar in shape to the curves obtained 98 for various T /tr^ ratios in the laboratory. This indicates that a reasonable matching of these curves is possible. In the computation of E"r using equation (5.29), ayQ and 0^ are substituted for ayo and respectively. Here Oy0 and Oy are initial vertical effective stress and current vertical effective stress respectively. o"yD is known from the static analysis performed before the dynamic analysis. o"y is obtained by computing current effective stresses. Recall equation (4.9), {P} = [Kj {A} + [K*] {U} (5.30) This equation can be used to obtain the response of the deposit to increases in residual porewater pressures by setting {p} = {o}. The porewater pressure matrix {u} is formed from the residual porewater pressures in the elements, as in the case of static analysis. The incremental displacements, stresses and strains given by adopting this procedure is the response of the deposit to softening .of the elements. Furthermore, these incremental strains can also be viewed as permanent components of the strains. The current effective stress system can now be used to modify Gmax and xmax values, and the dynamic analysis can be continued with a hyperbolic stress-strain relationship compatible with the current effective stress system. There is a limit to the amount of residual porewater pressure that can be achieved in a triaxial apparatus. This is easy to estimate, because the ultimate stress state of a sample has to be on the Mohr 99 envelope and in a triaxial apparatus the effective stress path followed by the sample is predictable. However, in two-dimensional analysis the effective stress path followed by an element cannot be predicted before hand, and therefore the maximum residual porewater pressure that can be developed cannot be known before the dynamic loading. The procedure used to calculate effective stresses can be used to impose limits on the amount of residual porewater pressure. This can be accomplished by allowing porewater pressure generation to occur during dynamic loading only until the Mohr circle drawn for the current effective stress system touches the Mohr envelope. 5.3 DAMPING MATRIX [c] Two fundamentally different damping phenomena are associated with soils, namely material damping and radiation damping. The material damping can be viewed as a measure of energy dissipation when waves travel through soils. The loss of energy due to waves travelling away from the region of interest is known as radiation damping. The material damping can be divided broadly into two groups: viscous and hysteretic damping. The energy dissipation in viscous damping depends on the velocity of motion or strain rate and it is frequency dependent. Whereas hysteretic damping involves frictional loss of energy that is largely independent of frequency but depends on the magnitude of displacement or strain. Laboratory tests on soil samples have shown that most of the energy dissipation in soils occurs through internal friction which is hysteretic. When modelling the non-linear behaviour of soil by an incrementally elastic approach, the effect of hysteretic damping has 100 been included already in the analysis. Viscous damping is still needed to take into account the effect of flow of water inside the soil structure. Furthermore, a small amount of viscous damping is necessary to control pseudo high frequency response introduced by the numerical integration procedure. The damping due to viscous effects can be accounted for through the use of Rayleigh damping. The damping matrix [c] is given by a linear combination of [M] and [K] giving, [C] = a [M] + b [K] (5.31) in which a and b are constants. The stiffness matrix [K] varies with time during the dynamic analysis and therefore [c] would have to be computed for every time step. But knowing that the amount of viscous damping is very small compared to hysteretic damping, the time consuming procedure of computing [c] at all time steps may be unnecessary. Therefore, [c] is assumed to be a constant and can be evaluated using the tangent stiffness matrix [KT] at time t=0. Then, [C] = a [M] + b [KT]T=0 (5.32) This will give a damping ratio X for the nth mode as, bo) n 2(JO 2 n (5.33) 101 where the is the nth mode frequency. Equation (5.33) shows that the mass-proportional component of the damping is inversely proportional to the frequency while the stiffness proportional component is directly proportional to the frequency. Lee (1975) proposed that only stiffness-proportional damping should be used and suggested using, a = 0 and b = 0.005 (5.34) In substituting these values for a and b in (5.33) one gets, Xn = 0.0025u>n (5.35) It should be remembered that in typical soil structures only lower modes of vibration govern the response, and therefore, it is unnecessary to include higher mode components. The equation (5.35) implies that the damping ratio increases linearly with the frequency. Therefore, the response due to high frequency components of the input motion will be damped out significantly. This is one of the advantages of using stiffness proportional damping ratios for soils. From the equation (5.35) the damping ratio at the fundamental frequency is given by, \L = 0.0025Wl (5.36) Typical periods of vibration of soil structures may vary between 0.5 to 1.5 sec. This means that the typical damping ratio for the fundamental mode at the start of the dynamic loading varies between 1% - 3%. 102 The average stiffness of a soil structure during dynamic loading is less than the stiffness at the start of the dynamic loading. Therefore, when using the damping matrix [c] computed using [K]t=0, the effective damping ratio may be higher than the range shown above. 5.4 BOUNDARY CONDITIONS Appropriate boundary conditions in terms of forces or displacements have to be specified at all boundaries. In the dynamic analyses involving earthquake excitation, two types of bottom boundary conditions are often used: 1. A fixed bottom boundary located at the top of a rigid layer. In general, base rock or a stiff soil layer can be assumed to be rigid. 2. A bottom boundary located at the top of a soil layer or soft rock with constant elastic properties. This type of boundary is generally known as transmitting boundary. For both the above boundary conditions the earthquake motion is specified at the bottom boundary. If the second type of bottom boundary conditions is used, the domain of interest need not be extended down to a rigid layer. This procedure reduces the number of degrees of freedom leading to a reduction in computational costs. In the method of analysis presented here only the first type of boundary condition is considered. Three types of lateral boundary conditions are commonly used in two-dimensional dynamic problems, involving finite element procedures: 103 1) boundaries are located sufficiently far away from a structure so that wave reflection does not occur during analysis or is minimized. On these boundaries, forces, displacements or a combination of forces and displacements can be specified. 2) viscous boundaries are used which attempt to absorb the radiating waves by a series of dashpots and springs with constant or variable properties (Lysmer and Kuhlemeyer; 1969) 3) consistent boundaries can be provided close to the foundation of structures. These boundaries attempt to reproduce the far field response in a way consistent with the finite element formulation used to model the dynamic problem. This is accomplished by formulating a frequency dependent boundary stiffness matrix which can be obtained by solving the elastic wave propagation problem in a layered half-space (Lysmer and Wass; 1972). Of the three types of lateral boundaries, the analysis with the consistent boundaries is by far superior to the others with respect to accuracy. In the analysis with consistent boundaries, only a small region needs to be considered, thus reducing the number of degrees of freedom. But unfortunately, the formulation is strictly applicable only to linear (or iterative linear) problems and for solutions in the frequency domain only. The lateral boundaries in an incrementally elastic approach in the time domain therefore should be located as far away from a structure as practicable. The usual way to model the lateral boundaries is to allow the nodes on these boundaries to move only in the horizontal direction. 104 5.5 ANALYSIS OF SOIL - STRUCTURE SYSTEMS The response of a structure founded on a soil deposit is affected by the local soil conditions. The peak acceleration, frequency content and the spatial distribution of the response characteristics may all be affected. By including the structure in the finite element domain for the response analysis, the coupled seismic response of the soil and structure may be determined. The presence of the structure has two major effects on the soil deposit. It increases the effective stresses and it also provides additional inertia forces. Therefore, for soils which exhibit non-linear stress dependent behaviour, an uncoupled analysis in which soil and structural systems are uncoupled may not be applicable. Uncoupled analysis cannot predict the response of buried structures where strong soil-structure interaction occur. 5.5.1 Slip Elements in Dynamic Analysis Relative displacement which may occur between soil and structure during strong shaking can be modelled using slip elements. However, the slip element model described in Chapter 4 has to be modified for use under dynamic loading conditions. Ideally, once the slippage stops in a slip element, the top and bottom nodes of the element should have the same acceleration and velocities (Nadim and Whitman; 1982). However, when computing accelerations and velocities by numerically integrating the equation (5.5), different values will be unavoidable for the top and bottom nodes. In order to overcome this difficulty, when no slippage 105 occurs In a slip element, the velocities and accelerations of the bottom nodes were made equal to those of the top nodes. 5.6 SOLUTION SCHEME A step by step solution scheme is carried out to obtain the dynamic response in the time domain. A brief outline of the procedure is given below: 1) based on the current values of G , T and v. max' max ' t at time t the tangent shear modulus Gt is calculated for all the elements using the current stress-strain curve, for either initial loading, unloading or reloading whichever is appropriate. 2) the global stiffness matrix [Kt]t is determined. 3) knowing the base acceleration value at (t+At), new values for {x}, {x}, {x} at time (t + At) and increments, Ay and {ACT} are computed by employing any of the direct integration methods to solve the equations (5.5) as detailed in Appendix III 4) if stress or strain reversal occurs in any element during this time interval, At, the dynamic analysis is repeated for a shorter time interval. 5) using the shear strain increment, increments in volumetric strain and then porewater pressure are computed using equation (5.28) and (5.29). 106 6) using increments in residual porewater pressure U, as virtual loads static analysis is performed to determine current effective stresses. These effective stresses are then used to update and X]nax values. A computer program TARA-2 has been developed incorporating all these basic steps. 5.7 COMPUTATION OF POST EARTHQUAKE DEFORMATION It is often required to predict the displacements at various points on the soil structure at the end of an earthquake. This is referred to as dynamic residual displacements. To compute these, an earthquake record with enough tailing zeros should be used so that the free damped vibration response can be included in the analysis. The cyclic and permanent components of displacement response for saturated sands and silts, are assumed to occur under undrained conditions. There will be additional deformations in these soils when the residual porewater pressure dissipates. In undrained simple shear cyclic tests on saturated sands and silts the potential volumetric strain ey<j, which occurs due to grain slip, is reflected as residual porewater pressure. When the residual porewater pressure dissipates, volumetric strain, eV(j occur in the sample resulting in settlement. The prediction of deformations due to the dissipation of residual porewater pressure can be accomplished as follows. 107 One method is to treat this problem as a two-dimensional consolidation problem with known initial porewater pressures, and to solve for the deformation at discrete time intervals as the drainage occurs. A second method is to compute deformations using the volumetric strain is accumulated at the end of the dynamic loading, eV(j» This can be accomplished by treating this volumetric strain the same way the volumetric strain ey in shear-volume coupling was modelled in Chapter 4. In the program TARA-2 the second method is used. The final or post earthquake deformation after the residual porewater pressure has dissipated is the sum of the deformation calculated ic from the modified ey(j and the residual dynamic deformation predicted at the end of the dynamic analysis. 108 CHAPTER 6 VERIFICATION OF THE METHOD OF ANALYSIS The validation of computational techniques requires good prototype data for comparison. The common shaking table test cannot represent the range of in-situ pressures experienced by the soil elements in the field. For application to full scale design problems we require data from centrifugal models where similarity can be achieved with the self-weight of the prototype. Here in-situ stresses are directly scaled to those of the full scale event and therefore, the stress dependent stress-strain properties of the soil (especially cohesionless soils) can be matched in model and prototype. Centrifuge modelling laws are used to deduce prototype response from model response. The principles of centrifugal modelling have been discussed in detail by Schofield (1981). 6.1 CAMBRIDGE CENTRIFUGE TESTS A number of centrifuge tests on submerged islands were conducted by Lee (1983) using the Cambridge University Centrifuge. Full details of Cambridge centrifuge equipment and test procedures have been given by Schofield (1981). Fig. 6.1 shows the model island used in the tests. The island was a 90mm high with side slopes at 3:1 and a crest width of 200mm. The centrifuge acceleration used in the tests was 40g. This means that the corresponding prototype island is of height 3.6m, with side slopes 3:1 and has a crest width of 8m. The structural loading on the LVDT 82260 90mm I 720mm • I I LVDT 82273 ACC 734 ACC 1244 • PPT 2 338 OPPT 2332 Mild steel plates 3 PPT 2331 Concrete base ACC 258^ Container Base 750 mm 776 mm • Direction of positive acceleration measured on accelerometers LVDT Pore pressure transducer *• Accelerometer Fig. 6.1. Submerged Island Showing Transducer Locations. o KO 110 island was simulated by using mild steel plates of various thicknesses. The island rests on a concrete base which in turn is bolted to the centrifuge container. The base shaking was generated when the rotary arm on wheels follows a track mounted on the wall of the centrifuge chamber. The island was instrumented by 6 DJB A23 piezoelectric accelerometers, 6 Druck PDCK81 pore pressure transducers and 2 LVDTs. The location of these instruments are also shown in Fig. 6.1. The test soil was fine Leighton-Buzzard sand, mostly of size passing between B.S. sieve size No. 120 and No. 200. The relative density was 60 - 70%, with emax = 1.03 and emin = 0.63 (Lee, 1983). In centrifugal modelling, if the model pore fluid is the same as the prototype pore fluid, excess pore pressures would be able to dissipate N2 times faster in the model than in the prototype, whilst the earthquake would only occur N times faster. Here N is the scale factor of the centrifugal acceleration given as a ratio of gravity. Therefore, to model the prototype porewater drainage condition, it is necessary to use a pore fluid of viscosity N times that of water in the model test. A special silicone oil was used as the model pore fluid in Lee's tests. In the tests carried out on the island, the theoretical input wave-form was intended to be 12 sinusoidal pulses with a constant period of 0.5secs. However, the actual input motion to the island was more complicated due to resonances and mechanical linkage clearances interfering with the input motion, especially during the initiation of the base motion. The results of two centrifuge tests were made available to the writer. The average contact prototype pressure on the islands for these two tests were 15kPa (Test 1) and 31kPa (Test 2) respectively. The input Ill accelerations measured by the accelerometer mounted on top of the concrete base for Test 1 and Test 2 are shown in Fig. 6.2(a), and (b). The recorded maximum base accelerations in Test 1 and Test 2 were O.llg and 0.17g respectively. The island responded very differently in these two tests. The comparative study carried out to predict performance of both the tests are reported in this chapter. 6.2 COMPARATIVE STUDY The sand used in the model test was at an average relative density Dr = 65%. Typical properties which are consistent for medium dense sand of Dr = 65% for static and dynamic analysis are given in Table 6.1. The liquefaction resistance curve for the sand obtained by using UBC simple shear apparatus without any static bias (tg=0) is shown in Fig. 6.3. This liquefaction resistance curve matches the predicted liquefaction resistance curve when porewater pressure model constant, Kr = 0.012. As explained in Chapter 5, the porewater pressure model constants can be selected appropriately to account for the behaviour of samples with initial static stress, TG. Since test data on samples with TG, are not available, it was decided to ignore the changes in the pore water pressure generation rates, and to account only for the changes in liquefaction resistance. The liquefaction curves which correspond to non zero ts/°"v0 were obtained by scaling the laboratory liquefaction curve. The scaling was done by using available laboratory data on medium dense Ottawa sand (Vaid, et al. 1979). The changes in the liquefaction resistance curve can be specified by associating the appropriate Kr ~i 1 1 r 3.0 4.6 5.0 Time in Sees Fig. 6.2(a). Input Base Motion in Test 1. • rc l_ tu _| —i max = 0.17g —i 1 1 1 1 1 1 1 1 i r-2.0 3.0 4.0 5.0 6.0 7.0 Time in Sees 8.0 Fig. 6.2(b). Input Base Motion in Test 2. 113 Table 6.1 Soil Properties Properties Total unit weight kN/m3 Bulk modulus constant Bulk modulus exponent constant n Shear modulus parameter (K2)max Angle of internal friction Effective cohesion Coefficient KQ a, b values used to compute [C] C^-K^ Constants Rebound modulus constants m,n Static Dynamic 18.1 800 0.4 19.3 38.0 0.0 0.45 18.1 high 55.0 38.0 0.0 0.0,0.005 0.75,0.79,0.459,0.73 0.43,0.62 04\ Number of Cycles to Initial Liquefaction, Fig. 6.3. Liquefaction Resistance Curve of Medium Dense Leighton-Buzzard Sand. 115 value with each static shear stress ratio, "^/^vo* Linear interpolation may be used to get the Kr value corresponding to an" intermediate -us/avo. For each model test two analyses were performed. One with no slip elements between soil and structure and the other with slip elements. The following slip element properties were used, Kg = 6.3 x 105 kN/m2/m, % » 6.3 x 105 kN/m2/m Cg = 0.0 0g = 35° 6.2.1 Results of Test 1 The recorded acceleration time histories of the model island have very high frequency components which contain negligible energy. This type of high frequency electrical noise is unavoidable in centrifuge testing as it may originate due to ambient sources such as the electric motor driving the centrifuge, and also due to centrifuge vibrations. Dean (1983), suggested it is necessary to filter out very high frequency components from output quantities. The computed and recorded acceleration time histories shown here have been smoothed once using a three point average scheme, suggested by Dean (1983). In using this scheme, the current value at any point in time is computed as the sum of 1/4 of the value the previous point, 1/2 the value of the current point and 1/4 of the value of the next point. Fig. 6.4 to Fig. 6.6 show the smoothed recorded and computed acceleration time histories of accelerometers A1244, A1225 and A734. During the first 1.5 seconds of shaking low accelerations with very high 116 ~i 1 1 1 r aO 1.0 2.0 3.0 4.0 5.0 Tine in Sees 6.0 r "5.0 8.0 Fig. 6.4(a). Recorded Acceleration of ACC1244 in Test 1. -i 1 r Q.D 1.0 2-0 3.0 4.0 5.0 Time In Sees r T "7.0 e.c Fig. 6.4(b). Computed Acceleration of ACC1244 in Test 1. (with and without Slip Elements). 117 !"c Tine in Sees 6.0 -r r—T \ 7.0 8.0 Fig. 6.5(a). Recorded Acceleration of ACC1225 in Test 1. Acceleration of ACC1225 * Teat 1. Fig. 6.5(b). Computed (with and without SUp Elements). Fig. 6.6(a). Recorded Acceleration of ACC734 in Test 1. Pig. 6.6(b). Computed Acceleration of ACC734 in Test 1. (with and without Slip Elements). 119 frequency were recorded in all accelerometers. After that the amplitude of acceleration steadily increased as in the case of input motion, upto 5.5 seconds and then subsided. The acceleration time histories computed by TARA-2 with and without slip elements are very similar and therefore, only one of them is shown in Fig. 6.4.(b), through Fig. 6.6(b). The frequency and magnitude of the computed acceleration response are very similar to corresponding recorded response values. Table 6.2 shows the computed and recorded maximum acceleration of all three accelerometers. Table 6.2 Recorded and Computed Maximum Accelerations Maximum Acceleration % 8 Instrument Location Computed by TARA-2 Recorded Without Slip Elements With Slip Elements A1244 13.3 11.6 11.6 A1225 15.9 12.5 12.5 A734 13.9 12.7 12.7 The comparison between recorded and computed maximum acceleration values are very good. Four porewater pressure development plots obtained experimentally and computed by TARA-2 are presented in Figures 6.7 (a), (b), (c) and (d) . In this test very low porewater pressures were developed. During the low level shaking of the first second, the response 120 Recorded Computed With and Without Slip Elements I I 1 1 1 1 1 1 1 1 1 4.0 5.0 -6.0 "7.0 6.0 at Tine In Sees Fig. 6.7(a). Recorded and Computed Porewater Pressure of PPT2330 in Test 1. x x x— Recorded Computed Without Slip Elements Computed With Slip Elements i 1 r i.O 5.0 Time In Seas Fig. 6.7(b). Recorded and Computed Porewater Pressure of PPT68 in Test 1. 121 Fig. 6.7(c). Recorded and Computed Porewater Pressure of PPT2338 in Teat 1. Recorded Computed With and Without Slip Elements i e.c Fig. 6.7(d). Recorded and Computed Porewater Pressure of PPT2342 in Test 1. 122 of the island is essentially elastic and porewater pressures recorded are the instantaneous response to changes in total stress. Such porewater pressures develop from the elastic coupling of soil and water. Later during the period of more severe shaking, plastic volumetric strains occur resulting in the development of residual porewater pressures which are independent of the instantaneous states of stress. The recorded porewater pressures during this time have both residual and instantaneous components. After 6 seconds of shaking the input motion subsides over the next two seconds to zero. During this time the magnitude of the instantaneous component of porewater pressure is small. Drainage may occur during the excitation depending on the drainage characteristics of the sand. Since generation of porewater pressure after 6 seconds of excitation is very small, changes in porewater pressure can be caused only by drainage during this time. A close examination of recorded porewater pressures after 6 seconds of excitation reveals that the porewater pressures in all four locations except at the transducer P2342, which is located at the middle of the island are more or less a constant. At this location porewater pressure increases slightly indicating movement of water towards the center of the island. However, since these changes are small it is reasonable to assume that drainage in the island is negligible during the base excitation. TARA-2 computes only residual porewater pressures, so there are no fluctuations due to changes in instantaneous stress levels in the computed curves. Furthermore, since no drainage was assumed during the excitation the computed residual porewater pressures increase consistently. However, the rate of increase in porewater pressures during 123 low level excitation which occur before 1.5 seconds and after 6 seconds is relatively small. When rigid joint connection is assumed between a heavy, stiff structural element and an adjacent soil element in an analysis, the dynamic strains developed in the adjacent soil element are very small due to compatibility requirements in displacements. However, by introducing slip elements between the structure and soil, the relative movement which may occur between soil and structure during strong base excitation can be accounted for. The results from TARA-2 analyses, with and without slip elements, indicate that computed porewater pressures are different only in the transducers located just below the structures. The predicted porewater pressures just below the structure in the analysis which incorporates slip elements, are slightly higher than the analysis that assumes rigid connection between soil and structure. The comparison between recorded and computed porewater pressures is good. Only four porewater pressure time histories from the model tests are available. However, maximum residual porewater pressures, which can be interpreted as the mean recorded values after the excitation has subsided, are available for all six transducers. Table 6.3 shows the computed and recorded maximum residual porewater pressures at all the transducer locations. 6.2.2 Results of Test 2 The smoothed recorded and computed acceleration time histories of accelerometers A1244, A1225 and A734 are shown in Figs. 6.8 through 6.10. Table 6.3 Recorded and Computed Maximum Residual Porewater Pressures Residual Porewater Pressure, kPa Transducer Location Computed by TARA-2 Recorded Without Slip Elements With Slip Elements P2330 0.4 0.4 0.4 P2331 0.4 0.3 0.4 P2332 4.0 2.9 3.7 P68 3.0 1.5 4.4 P2338 2.4 0.6 1.8 P2342 4.0 4.2 4.2 125 T I 2.0 T" 3.D 4.0 5.0 Time In Sees Fig. 6.8(a). Recorded Acceleration of ACC1244 in Test 2. Fig. 6.8(b). Computed Acceleration of ACC1244 in Test 2. (with Slip Elements). + 25% 1 1 1 1 1 1 1 1 1 1 1 1 I 1 I I 13 LC 2.0 3.0 4.0 5.0 6.0 ">.0 8.0 Time In Sees Fig. 6.9(a). Recorded Acceleration of ACC1225 in Test 2. + 2 6% "1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ' a0 1.0 2.0 3.0 4.0 5.0 6.C "7.0 8.TJ Time in Sees Fig. 6.9(b). Computed Acceleration of ACC1225 in Test 2. (with Slip Elements). 127 0.0 "i i 1 1 1 r 2.0 3.0 4.0 5.0 Time in Sees i 1 1 r 6.0 "7.0 s.c Fig. 6.10(a). Recorded Acceleration of ACC734 in Test 2. Time in Sees Fig. 6.10(b). Computed Acceleration of ACC734 in Test 2 (with Slip Elements). 128 The base excitation for Test 2, shown in Fig. 6.2(b), has the following characteristics. Low level excitation occurs in the first second and then the amplitude of acceleration increases steadily to a constant maximum value in the next two seconds. This maximum amplitude is maintained for 2 seconds and then it subsides over the next 3 seconds. The acceleration records, except the record obtained in accelerometer A734 which is located on top of the structure, show the variation of acceleration amplitude being similar to that of the input motion. The acceleration histories recorded on the top of the structure dropped to very low values after 4 seconds of excitation. Unlike Test 1, the response computed with slip elements were found to be different from those computed without slip elements. The response analysis also showed that when slip elements are used significant amount of slip occurs between soil and structure. The computed accelerations without slip elements are lower than the values computed when slip elements are used in the analysis. In the computed acceleration histories shown in Figs. 6.8 to Fig. 6.10, only accelerations computed using slip elements are presented. Lines of constant accelerations have been drawn in Figs. 6.8 to 6.10 to aid the interpretation of results. The variation of computed acceleration histories are very similar to that of the input motion. The maximum recorded accelerations which were observed to be associated with high frequencies (indicated by sharp spikes) and maximum computed accelerations for all three accelerometers in the island are shown in Table 6.4. The comparison between recorded and computed accelerations are not good. However, a closer look at the other peak values of corresponding recorded and computed acceleration histories suggests satisfactory comparison. 129 The recorded acceleration at any time has two components: the acceleration component transmitted through the soil from the base and the acceleration component transmitted to soil through side walls and top of the centrifuge container due to the container itself vibrating. The computed acceleration history accounts only for the acceleration component transmitted through the soil. The accelerations transmitted through the container have high frequency components. The presence of these high frequency components may be responsible for the discrepancies between the recorded and computed acceleration histories. Table 6.4 Recorded and Computed Maximum Accelerations Maximum Accleration, % g Accelerometer No. Computed by TARA-2 Recorded Without Slip Elements With Slip Elements A1244 24.0 15.1 18.2 A1225 42.5 15.5 23.1 A734 23.9 15.8 18.2 Four recorded and computed porewater pressure development plots are presented in Figs. 6.11 (a), (b), (c) and (d). In this test unlike Test 1, very high porewater pressures were developed. During the low level of excitation of the first second very low porewater pressures were recorded. With the onset of more severe shaking, very high porewater 130 Fig. 6.11(a). Recorded and Computed Porewater Pressure of PPT2330 in Test 2. Fig. 6.11(b). Recorded and Computed Porewater Pressure of PPT68 in Test 2. Fig. 6.11(c). Recorded and Computed Porewater Pressure of PPT2338 in Test 2. Fig. 6.11(d). Recorded and Computed Porewater Pressure of PPT2342 in Test 2. 132 pressures were developed in all transducers except in the transducer P2330. A close examination of the recorded porewater pressure plots after 6 seconds of excitation reveals that, except for the porewater pressure transducer P2342, which is located at the middle of the island, significant dissipation of porewater pressure occurred. This is because in this test very high porewater pressures were developed leading to high pressure gradients. The tranducer P2342 behaved differently because it is too far from the free draining boundaries and at this location inward flow of water occurs. The comparison between computed and recorded porewater pressures are very good for the two transducers (P2330, P2342) which are located well inside the island. At these transducer locations the analysis with and without slip elements gave very similar results. The porewater pressures predicted under the structure by the analysis without slip elements are very low. But, when slip elements were provided the comparison between predicted and computed porewater pressures improved. When slip elements were used in the response analysis high shear strains developed in the elements below the structure resulting in higher porewater pressures in those elements. Table 6.5 shows the computed and recorded maximum residual porewater pressures at all the transducer locations. 6.3 APPLICABILITY OF THE METHOD OF ANALYSIS The recorded residual porewater pressures are interpreted as the steady increase in porewater pressures. Therefore, any high frequency noise from ambient sources on recorded values do not affect the comparison Table 6 .5 Recorded and Computed Maximum Residual Porewater Pressures Maximum Residual Porewater Pressure, kPa Transducer No. Computed by TARA-2 Recorded Without Slip Elements With Slip Elements P2330 1.0 1.5 1.5 P2331 0.9 1.1 1.1 P2332 10.5 12.0 12.1 P68 38.0 6.4 38.1 P2338 18.0 2.2 . 18.9 P2342 22.0 19.8 21.3 134 between recorded and computed porewater pressures. However, since the comparison in acceleration response is made point by point, any high frequency noise affects the comparison. Therefore, caution should be exercised when comparing maximum accelerations. When a sand sample is subjected to undrained loading, an abrupt change in the direction of stress path occurs, as the stress path approaches the failure line (Ishihara, et al. 1975). The points at which various stress paths change direction abruptly are assumed to lie on a straight line, called the phase transformation line. The slope of the phase transformation line is a few degrees less than the failure line. Any cyclic (or monotonic) loading beyond the phase transformation line may result in very low effective stresses in very loose sands and increased effective stresses due to dilation in medium dense or dense sands. The hyperbolic stress-strain relationship assumed in this thesis is strictly applicable only for the region of stress space below the phase transformation line. Figs. 6.12 (a) and (b) show the stress paths in a q, p' plot that were followed by four elements which correspond to the locations for which porewater pressure time histories are available. The stress paths reported here are for the analysis in which slip elements were included. The stress paths followed by the elements in Test 1 are well below the failure line, where as in Test 2, two elements are on the failure line for sometime during the dynamic loading. Under these circumstances the validity of response computed in Test 2 after elements have reached the failure line is questionable. Fig. 6.12(a). Effective Stress Paths in Test 1 Fig. 6.12(b). Effective Stress Paths in Test 2. ON 137 CHAPTER 7 APPLICATION OF THE METHOD OF ANALYSIS: TANKER ISLAND RESPONSE 7.1 INTRODUCTION Man made islands of cohesionless soils have been used extensively as drilling platforms for oil and gas exploration in the Beaufort Sea. Recently, as exploration has moved to deep waters, more complex forms of island construction procedures have been introduced. The caisson-retained island (De Jong, et al. (1978), and steel tanker island are two typical examples. These newer type of construction procedures greatly reduce the amount of fill material required and also reduce some of the hazards of wave loading on exposed island beaches. The maximum set down water depth for the caisson-retained island and the tanker islands depth is fixed, generally around 6 to 9 metres. Therefore, in the case of deep water a underwater sand berm is constructed up to the set down water depth. Most of the sand berms are constructed by dumping sand excavated by suction dredges from an offshore and/or onshore borrow pit and pumped as a slurry through a pipeline directly onto the location of the island. Once the sand berm is ready, a series of caissons or tanker is brought to the location and ballasted onto the berm, and backfilled with sand, gravel or water. The drilling is then carried out from the upper structure. Because of the nature of the island construction, the dumped sand is often loose and therefore the deformation, stability and 138 liquefaction potential of the island berm during earthquakes are of major concern. 7.2 ANALYSIS OF A TYPICAL TANKER ISLAND Fig. 7.1, shows, schematically a tanker island. This island is provided with a cover of about 2m of rock fill. Typical properties of rock fill and sand for static analysis are given in Table 7.1. TABLE 7.1. Static Soil Properties Properties Rock Fill Sand Total unit weight kN/m3 18.7 18.1 Bulk modulus constant K^ 1000 800 * Bulk modulus exponent constant n 0.40 0.40 Shear modulus parameter (K2)max 24.0 16.0 Angle of internal friction 38.0 32.0 Effective cohesion 0.0 0.0 Coefficient KQ 0.45 0.45 The tanker is assumed to weigh 200,000 metric tons when fully ballasted with plan dimensions 170m and x 60m and 21m high. In the case presented here, it is assumed that the hyperbolic dynamic stress-strain relationship has equal shear strength in both Fig. 7.1. Schematic of Tanker Island 140 directions of shearing. Typical properties used for the dynamic analysis are give in Table 7.2 TABLE 7.2. Dynamic Soil Properties Properties Rock Fill Sand 1300 V. High 0.4 70 45 0.8, 0.79, 0.459 and 0.730 Rebound modulus Constants m,n 0.43, 0.62 During shaking the rock fill is assumed to be free draining and no drainage is assumed in dumped sand. A very high value was assigned to Bt to simulate very low compressibility imparted to the saturated sand by the water in the pores which is not allowed to drain. Liquefaction resistance curves are required for different static stress ratios in the island. These are specified by associating the appropriate Kr value with each static shear stress ratio ts/o"vo. The values used in the example are presented in Table 7.3. Bulk modulus constant ft Bulk modulus parameter n Shear modulus parameter (K2) C^ •* CK Constants TABLE 7.3. Ts/avo and Kr values T /a* K s vo r 0.0 0.004 0.10 0.015 0.20 0.05 The static shear stress ratios for the example problem considered here vary between 0.0 to 0.13 and therefore, the values of Kr corresponding to ratios Tg/oyo, above 0.15 are not necessary. The input motion used for the analyses is the S00E acceleration component of the Imperial Valley Earthquake of May 18, 1940 scaled to O.lg. The input motion was applied at the bottom boundary of the island. Three dynamic analyses were performed: island alone, island plus structure with and without slip elements. The properties of slip elements were selected so that some slip could occur between the structure and the island. The slip element properties were assumed to be, Ks = 6.3 x 105 kN/m2/m, KN = 6.3 x 105, kN/m2/m Cs = 0 and 0g = 30° A complete response study of the tanker island could be carried out by discretizing the entire domain into finite elements. However, since the stiffness of tanker wall is very much higher than soil, the tanker and its contents would respond like a rigid box. In view of 142 this, the tanker and its contents were modelled as a uniform rigid box. The stiffness of structural elements were selected as 103 of the rock fill elements. 7.2.1. Results for Tanker Island Problem One of the factors which influence the development of residual porewater pressure is cyclic shear stress (or cyclic strain) . Since the generation of residual porewater pressure is possible only in the sand, the maximum dynamic shear stresses induced in the dumped sand along section T-T which runs through the centre of the island, are shown in Fig. 7.2 for all three cases. This figure indicates that higher dynamic shear stresses are developed when the tanker is in place due to the inertia forces on the tanker. The induced shear stresses in the dumped sand when slip is allowed to occur between structure and adjacent soil are slightly higher than when no slip elements were provided. When slip occurs, the magnitudes of shear stress that can be transmitted to the structure is limited, dictated by the shear strength of the slip elements. Therefore, island responses with and without slip occurring between the structure and the island may be expected to differ. Despite very high cyclic shear stresses generated in the sand when tanker is in place, the greatest porewater pressure ratios are developed in the unloaded island (Fig. 7.4). This is because the vertical over-burden pressures are very much greater when the tanker is present, so that, the cyclic shear stress ratios, Tcv/Oy0, which is the most important parameter controlling the development of porewater pressure in a given sand, are actually smaller (Fig. 7.3). It can be readily seen that Maximum Dynamic Shear Stress (kN/m2) 0 20 40 60 I 1 1 1 Fig. 7.2. Distribution of Maximum Dynamic Shear Stress in Sand. Fig. 7.3. Distribution of Maximum Dynamic Shear Stress Ratio in Sand. 144 the distribution of maximum cyclic shear stress ratios are proportional to the distribution of residual porewater pressure ratios. Fig. 7.4, shows the residual porewater ratio H/aLQ distribution for all three analyses. The residual porewater pressure ratios obtained for the unloaded island are higher than for the loaded island. The results obtained in the analysis without the tanker can be viewed as the solutions at a section where the influence of the structure is negligible. Therefore, the distribution of residual porewater pressure ratios when the structure is in place will vary from lower values at the middle to higher values as one moves away from the structure. Same conclusions were drawn by Yoshimi and Tokimatsu (1977) who studied the response of a rigid structure subjected to base excitation on a shaking table. The residual porewater pressure distribution given in the analysis with slip elements is consistently higher than the analysis without slip elements. This is because lower shear stresses are induced in the latter case. The distribution of maximum dynamic shear strains for section 1-1 is shown in Fig. 7.5. Even though the shear stresses induced in the unloaded island are smaller than the loaded island, higher shear strains have developed in the unloaded island. This is because of two factors. Firstly, the in-situ overburden stresses in the unloaded island are very much smaller and therefore the stress-strain curve for a given element is softer. Secondly, when an effective stress-strain relationship is used, any increase in residual porewater pressure will soften the stress-strain curve. The generation of higher residual porewater pressures and low over burden pressures have contributed to high shear strains in the unloaded island. Fig. 7.4. Distribution of Residual Porewater Pressure Ratio in Sand. Fig. 7.5. Distribution of Maximum Dynamic Shear Strain in Sand. 146 The maximum horizontal displacements which occur during the earthquake for the section 1-1 are shown in Fig. 7.6. Much smaller dynamic displacements are computed when the tanker is in place. When slip elements were provided, slip occurred between the soil and structure and the displacements are about twice the results obtained without slip elements. Fig. 7.7 a and b show the post earthquake deformations in the X and Y directions. The post earthquake displacement is the sum of the dynamic residual displacement and the displacement due to volumetric * strain component ey^ • Two observations can be made from the results presented in this figure. Firstly, the amount of the post earthquake deformations in the X-direction are proportional to the maximum dynamic deformations. Secondly, the X-component of the displacement is of the same order of magnitude as of the Y-component of the displacement (settlement). The main contribution to the X-component of the displacement comes from the dynamic residual displacement and for the settlement the main contribution is from the volumetric strain component ft evd* In the dynamic response of structures the maximum induced acceleration in the structure is one of the main design concerns. The maximum induced accelerations given by TARA-2 in the structure with and without slip elements are 0.15g and 0.17g. This means that, if slip is prevented, the acceleration induced may be higher by as much as 15% of the acceleration when slip is allowed. The maximum acceleration computed on top of the unloaded island is 0.15g. Fig. 7.6. Distribution of Maximum Dynamic Displacement. Post Earthquake X Disp.. (cm). Post Earthquake Y Disp., (cm) Fig. 7.7. Post Earthquake X and Y Displacements. 00 149 One of the ways of presenting dynamic response is to present it in terms of response spectra. Response spectra for displacement, velocity and acceleration are often presented for the motions at the base of the structure. These results are then used by the engineers to predict the behaviour of structures and also to compute design forces, such as base shear. Fig. 7.8 shows the response spectrum for acceleration of the motion at the berm surface for all three cases considered. The damping ratio used in the computation was 3%. Inspection of this figure suggests that for the example problem considered here, the acceleration response predicted using the response spectrum of the unloaded island will be higher for structures with a very low period. However, for the structures with a period greater than 0.5 sec, the response predictions will be similar. The predominant motion of a tanker during excitation are sliding and rocking. The relative importance of these two modes can be studied by comparing results obtained by a two-dimensional and one-dimensional response analysis. Fig. 7.9, shows the computed distribution of residual porewater pressure ratio u/oyo from three two-dimensional analyses which were reported earlier and also the distribution from a one-dimensional response analysis. The one-dimensional case considered was the island with tanker, without any slip elements. The results clearly show that the maximum U/Oy0 of some elements may be predicted as low as 30% of those predicted by a two-dimensional response analysis. This means a response analysis which neglects the rocking mode of vibration is non-conservative. However, it should be mentioned that the tanker considered in this example is very tall (21m) and, rigid, and therefore the rocking mode of vibration may have been more important than usual. 600 1 r- r —i 500 Island Alone — 400 |t "Jl II •*—. Island + Structure (Slip) -300 \f\ Island + Structure ?A (No Slip) -200 -100 0 10 2-0 3-0 4-0 5-0 Period, (sec) Fig. 7.8. Acceleration Response Spectrum for the Motion at Berm Surface. 7.9. Distribution of Residual Porewater Pressure Ratio in Sand. 151 7.3 SOME PRACTICAL CONCLUSIONS In practice a number of simplifying assumptions are made when computing the response of structures founded on soil deposits. The procedure outlined by the National Building Code of Canada has two basic steps. The first step is to compute the response of the soil deposit alone for the given excitation. The second step is to compute the response of the structure, to the base accelerations obtained in step one. In predicting the performance of the structure, the results such as porewater pressure, induced strain level etc., which are obtained from step one are also considered. This means that the code in essence suggests the soil- structure systems be uncoupled and analysed independently. Figures 7.2 to 7.8 clearly show that the response of structures computed using the procedures outlined in the National Building Code of Canada may be in error. The presence of the structure has two basic influences on the soil deposit. It increases the effective stresses and it also provides additional inertia forces. Therefore, for soils which exhibit non-linear stress dependent behaviour, the uncoupled analysis proposed by the code may not be applicable. From this typical example, three basic conclusions can be drawn. First of all it raises questions about the merit of any response analysis based on uncoupled soil-structure systems. Secondly, one-dimensional representation of the domain which neglects the rocking degrees of freedom may not be applicable to tall, heavy and rigid structures. Thirdly it demonstrates the importance of incorporating slip elements in the analysis. Because of the great weight of the caissons or tankers, and 152 their large lateral dimensions, soil-structure interaction effects will always be important. In these type of problems a coupled analysis of the island and structure is required. 153 CHAPTER 8 SUMMARY AND CONCLUSIONS 8.1 SUMMARY The main purpose of this research was to develop a two-dimensional static and seismic response analysis of soil deposits including soil-structure interaction. The new method for static and dynamic analyses can be performed in either effective or total stress modes or a combination of both. Non linear stress-strain behaviour of soil was modelled by using an incrementally elastic approach in which tangent shear modulus and tangent bulk modulus were taken as the two "elastic" parameters. The material response in shear was assumed to be hyperbolic with Masing behaviour during unloading and reloading. Response to changes in mean normal stress was assumed to be non-linear, elastic and stress dependent. When a static analysis is performed in the total stress mode, the shear strength, T„,0„, Gmav, and tangent bulk modulus, Bt, of an element are kept constant throughout the analysis. The tangent shear modulus, Gt, is modified for corresponding shear strains developed during the analysis. When effective stress mode is used, the parameters, tmax> ^max anc* Bt are computed from the effective stresses. The effect of dilation during shear on volume change is taken into account. In the static analysis proposed here, gravity may be switched on 154 at once for the completed soil structure or the construction sequence can be modelled by layer analysis. The stress-strain conditions determined by the static analysis give the in-situ stress conditions before the dynamic analysis. Slip or contact elements have been incorporated in the analysis to represent the interface characteristics between soil and structural elements. The properties of the slip element were assumed to be elastic perfectly plastic, with failure at the interface given by the Mohr Coulomb failure criterion. In the dynamic effective stress response analysis, residual porewater pressures are calculated using a modification of the model proposed by Martin et.al, (1975). The parameters, Gmax> and Tmax» are m°dified for the effects of residual porewater pressure. The dynamic response study includes the prediction of post earthquake deformations. An extensive study carried out to verify the proposed method of analysis using centrifuge test data suggests that the proposed method can be successfully used to predict seismic response of structures. Seismic response of a typical tanker island computed by this method is presented. 8.2 CONCLUSIONS The work that has been presented in this thesis leads to the following conclusions. 1. A consistent and reliable method for computing transient and 155 permanent deformations in two-dimensional soil structures is needed. 2. A two-dimensional dynamic response analysis which takes into account the non-linear hysteretic stress-dependent properties of soils, has been developed in terms of both total and effective stresses. 3. The method has been verified by comparing data from centrifuged models with predictions of the method. Comparison between predicted and measured response parameters is generally very good. 4. Allowing for slip to occur between soil and structural elements is very important. Analyses which allow for slip have consistently lead to higher displacements in the structure and higher porewater pressures in the soil deposit. 5. The method has been applied to compute seismic response of a typical tanker island. The results of this study suggests that it is important that the response of structures founded on soil deposits be analysed as a coupled soil-structure systems. 6. The validity of a one-dimensional response analysis instead of a two-dimensional analysis for tanker type of structures is questionable. The porewater pressures predicted by using a one-dimensional response analysis model may be as low as 30% of those predicted by a two-dimensional response analysis. 156 SUGGESTIONS FOR FURTHER STUDY 1. Additional comparative studies should be carried out so that greater confidence could be placed on the validity of this method. Comparative studies may be performed with data from centrifuge tests or field studies. 2. Sandy materials exhibit partial stabilization at low confining pressures due to dilation. Therefore, in the response evaluation near liquefaction, it is important that the method of analysis Include the dilatant behaviour of sands. 3. In the response evaluation of more permeable soils, drainage during the seismic loading may be significant and procedures should be developed to take this into account. 157 REFERENCES 1. Byrne, P.M., (1979), Class Notes: Numerical Methods in Soil Mechanics (CE573). University of British Columbia, Vancouver, B.C., Canada. 2. Byrne, P.M., and Eldridge, T.L., (1982), "A Three Parameter Dilatant Elastic Stress-Strain Model for Sand", International Symposium on Numerical Models in Geomechanics, Zurich, Sept., pp73-79. 3. Chang, C.S., (1982), "Residual Undrained Deformation from Cyclic Loading", Journal of the Geotechnical Engineering, ASCE, Vol. 108, GT4, April, pp637-646. 4. Chern, J.C, (1981), "Effect of Static Shear on Resistance to Liquefaction", M.A.Sc. Thesis, The University of British Columbia, Vancouver, May. 5. Christian, J.T. and Boehmer, J.W., (1970), "Plane Strain Consolidation by Finite Elements", Journal of the Soil Mechanics and Foundation Engineering Division, ASCE, Vol. 96, No. SM4, July, ppl435-1455. 6. Clough, R.W., and Penzien, J., (1975), "Dynamics of Structures", McGraw Hill Book Co., New York. 7. De Jong, J.J.A., and Bruce, J.C, (1978), "Design and Construction of a Caisson-Retained-Island Drilling Platform for the Beaufort Sea", Proceedings, 10th Annual Offshore Technology Conference, Houston, Texas, May 8-11, OTC Paper No. 3294. 8. Dean, E.T.R., (1983), "FLY-14 Program Suite: in Flight Data Handling and Analysis Manual", Report, Cambridge Unversity Engineering Department, Cambridge, U.K. 158 9. Desai, C.S., and Abel, J.F., (1972), "Introduction to the Finite Element Method", Litton Educational Publishing, Inc. 10. Desai, C.S., and Christian, J.T. (1977), "Numerical Methods in Geomechanics", McGraw Hill Book Co., New York. 11. Duncan, J.M., Byrne, P.M., Wong, K.S., and Mabry, P., (1978), "Strength, Stress-Strain and Bulk Modulus Parameters for Finite Element Analyses of Stresses and Movements in Soil Masses", Report No. WCB/GT/78-02, to the National Science Foundation, April. 12. Finn, W.D.L. and Miller, R.I.S., (1973), "Dynamic Analyses of Plane Non-Linear Earth Structures", 5th World Conference on Earthquake Engineering, Rome, Session ID, Paper No. 42, pp360-367. 13. Finn, W.D.L., and Byrne, P.M., (1976), "Estimating Settlements in Dry Sands During Earthquakes", Canadian Geotechnical Journal, Vol. 13, No. 4, pp355-363. 14. Finn, W.D.L., Lee, K.W., and Martin, G.R., (1977), "An Effective Stress Model for Liquefaction", Journal of the Geotechnical Engineering Division, ASCE, Vol. 103, No. GT6, Proc. Paper 13008, June, pp517-533. 15. Finn, W.D.L., Martin, G.R., and Lee, K.W., (1978a), "Comparison of Dynamic Analyses of Saturated Sands", Proc. ASCE Geotechnical Engineering Division, Specialty Conference on Earthquake Engineering and Soil Dynamics, Pasadena, California, June 19-21, pp472-491. 16. Finn, W.D.L., Lee, K.W., Maartman, C.H., and Lo, R., (1978b), "Cyclic Pore Pressures Under Anisotropic Conditions", Proc. ASCE Geotechnical Engineering Division, Specialty Conference on Earthquake Engineering and Soil Dynamics, Pasadena California, June 19-21, pp457-468. 159 17. Finn, W.D.L., (1981), "Liquefaction Potential Development Since 1976", Proceedings, International Conference on Recent Advances in Geotechnical Earthquake Engineering and Soil Dynamics, St. Louis, Missouri, April 26-May 2, pp655-681. 18. Finn, W.D.L., Iai, S., and Ishihara, K., (1982), "Performance of Artificial Offshore Islands under Wave and Earthquake Loading", Offshore Technology Conference, Vol. 1, OTC, Paper 4220, Houston, Texas, May, pp661-672. 19. Finn, W.D.L., (1983), "Analyses of Cumulative Deformations Under Cyclic Loading", Invited Challenge Contributed to National Science Foundation Workshop on Expected Research in Soil Engineering, Blackberg, Virginia, U.S.A., Aug. 22-24. 20. Goodman, R.E., and Seed, H.B., (1966), "Earthquake Induced Displacements in Sand Embankments", Journal of the Soil Mechanics and Foundations Division, ASCE, Vol. 92, SM2, March, ppl25-146. 21. Goodman, R.E., Taylor, R.L., and Brekke, TL.L, (1968), "A Model for the Mechanics of Jointed Rock", Journal of the Soil Mechanics and Foundations Division, ASCE, Vol. 94, SM3, May, pp637-659. 22. Hansen, B., (1958), "Line Ruptures Regarded as Narrow Rupture Zones, Basic Equation Based on Kinematic Considerations', Proceedings, Brussels Conference on Earth Pressure Problems, Brussels, Vol. 1, pp. 39-48. 23. Hardin, B.O., and Drenevich, V.P., (1972), "Shear Modulus and Damping in Soils, Design Equations and Curves", Journal of the Soil Mechanics and Foundation Division, ASCE, Vol. 98, SM7, Proc. Paper 9006, July, pp667-692. 160 24. Iai, S., and Finn. W.D.L., (1982), "DONAL-2, A Computer Program for Dynamic One Dimensional Analyses of Slope Layers with Energy Transmitting Boundary", Soil Dynamics Group, University of British Columbia, Vancouver, B.C., Feb. 25. Ishihara, K., Tatsuoka, F. and Yasuda, S., (1975), "Undrained Deformation and Liquefaction of Sand Under Cyclic Stresses", Soils and Foundations, Vol. 15, No. 1, pp29-44. 26. Konder, R.L., and Zelasko, J.S.,, (1963), "A Hyperbolic Stress-Strain Formulation of Sands", Proceedings of 2nd Pan American Conference on Soil Mechanics and Foundation Engineering, Vol. 1, Brazil, 1963, pp289. 27. Kulhawy, F.H., Duncan, J.M., and Seed, H.B., (1969), "Finite Element Analysis of Stresses and Movements in Embankments During Construction", Geotechnical Engineering Research Report No. TE-69-4, Department of Civil Engineering, Unviersity of California, Berkeley, Nov. 28. Lee, F.H., (1983), "Partial Liquefaction in Centrifuge Model Embankment in an Earthquake", M.Phil. Thesis, Engineering Department Cambridge University, Cambridge, July. 29. Lee, K.W., (1965), "Triaxial Compressive Strength of Saturated Sands Under Seismic Loading Conditions", Ph.D Thesis, University of California, Berkeley, ppl-521. 30. Lee, K.W., (1975), "Mechanical Model for the Analysis of Liquefaction of Horizontal Soil Deposits", Ph.D Thesis, Department of Civil Engineering, University of British Columbia, Vancouver, B.C., Canada, Sept. 161 31. Lee, K.W., and Finn, W.D.L., (1975), "DESRA-1: Program for Dynamic Effective Stress Response Analysis of Soil Deposits Including Liquefaction Evaluation", Soil Mechanics Series, No.36, Dept. of Civil Engineering, University of British Columbia, Vancouver, B.C., Canada. 32. Lysmer, J., and Kuhlemeyer, R.L., (1969), "Finite Dynamic Model for Infinite Media", Journal of the Engineering Mechanics, Division, ASCE, Vol. 95, EM4, Sept. pp859-877. 33. Lysmer, J., and Waas, G., (1972), "Shear Waves in Plane Infinite Structures", Journal of Engineering Mechanics Division, ASCE, Vol. 98, EMI, Feb., pp85-105. 34. Martin,. G.R., Finn, W.D.L., and Seed, H.B., (1975), "Fundementals of Liquefaction Under Cyclic Loading", Journal of the Geotechnical Engineering Division, ASCE, Vol. 101, GT5, May, pp423-438. 35. Masing, G., (1926), "Eigenspannungen and Verfestigung beim Messing", Proceedings, 2nd International Congress of Applied Mechanics, Zurich, Switzerland. 36. Mroz, Z., Norris, U.A., and Zienkieicz, O.C, (1979), "Application of an Anisotropic Hardening Model in the Analysis of Elasto-Plastic Deformations of Soils", Geotechnique, 29(1), ppl-34. 37. Nadim, F., and Whitman, R.V. (1982), "Slip Elements for Earthquake Response", Proceedings, 4th International Conference on Numerical Methods in Geomechanics, Edmonton, Alberta, Canada, pp61-67. 38. Newmark, N.M., (1959), "A Method of Computation for Structural Dynamics", Journal of the Engineering Mechanics Division, ASCE, Vol. 85, EM3, July. 162 39. Newmark, N.M., (1965), "Effects of Earthquakes on Dams and Embankments", 5th Rankine Lecture, Geotechnique 15, No.2, ppl39-160. 40. Newmark, N.M., and 'Rosenblueth, E., (1971), "Fundementals of Earthquake Engineering", Prentice-Hall Inc., Englewood, Cliff, N.J., ppl62-163. 41. Ozawa, Y., and Duncan, J.M., (1973), "ISBILD: A Computer Program for Analysis of Static Stresses and Movements in Embankments", Geotechnical Engineering Research Report No. TE-73-4, Department of Civil Engineering, University of California, Berkeley, Dec. 42. Prevost, J.H., (1978), "Mathematical Modelling of Soil Stress-Strain Strength Behaviour", Proceedings, 3rd International Conference on Numerical Methods in Geomechanics, Aachen, Germany, April 6, pp347-361. 43. Pyke, R., Seed, H.B., and Chan, C.K., (1975), "Settlement of Sands Under Multidirectional Shaking", Journal of the Geotechnical Engineering ASCE, Vol. 101, No. GT4, April, pp379-398. 44. Rahman, M.S., Seed, H.B., and Booker, J.R., (1977), "Pore Pressure Development Under Offshore Gravity Structures", Journal of the Geotechnical Engineering Division, ASCE, Vol. 103, GT12, Dec, ppl419-1436. 45. Robertson, P.K., (1982), "In-situ Testing of Soil with Emphasis on its Application to Liquefaction Assessment", Ph.D Thesis, submitted to Department of Civil Engineering, University of British Columbia, Vancouver, B.C., Canada, Dec. 46. Schnabel, P.B., Lysmer, J., and Seed, H.B., (1972), "SHAKE: A Computer Program for Earthquake Response Analysis of Horizontally 163 Layered Sites", Report No. EERC 72-12, Earthquake Engineering Reseach Center, University of California, Berkeley, Dec. 47. Schofield, A.N., (1981), "Dynamic and Earthquake Geotechnical Centrifuge Modelling", International Conference on Recent Development in Geotechnical Earthquake Engineering and Soil Dynamics, Missouri, U.S.A., April 28-May 2. 48. Seed, H.B., and Lee, K.L., (1969), "Porewater Pressures in Earth Slopes Under Cyclic Loading Conditions", Proceeding, 4th World Conference on Earthquake Engineering, Chile, Vol. Ill, pplll. 49. Seed, H.B., and Idriss, I.M., (1970), "Soil Moduli and Damping Factors for Dynamic Response Analyses", Report No. EERC 70-10, Earthquake Engineering Research Center, University of California, Berkeley, December. 50. Seed, H.B., Lee, K.L., Idriss, I.M., and Makdisi, F.I., (1973), "Analysis of the Slides in the San Fernando Dams During the Earthquake of Feb. 9, 1971", Report No. EERC 73-2, Earthquake Engineering Research Center, University of California, Berkeley, June. 51. Seed, H.B., Pyke, R. and Martin, G.R. (1975), "Analysis of the Effect of Multi-directional Shaking on the Liquefaction Characteristics of Sands", Report NO. EERC 75-41, Earthquake Engineering Research Center, University of California, Berkeley, Calif. Dec. 52. Seed, H.B., (1979), "Considerations in the Earthquake-Resistant Design of Earth and Rockfill Dams", 19th Rankine Lecture, Geotechnique 29, No.3, pp215-263. 53. Seed, H.B., (1983), "Earthquake-Resistant Design of Earth Dams", Proceedings of Seismic Design of Embankments and Caverns, pp41-64. 164 54. Selig, E.T., and Chang, S.C. (1981), "Soil Failure Modes in Undrained Cyclic Loading", Journal of the Geotechnical Engineering Division, ASCE, Vol. 107, GT5, May, pp539-551. 55. Serff, N., Seed, H.B., Makdisi, F.I., and Chang, C.Y., (1976), "Earthquake Induced Deformations of Earth Dams", Report No. EERC 76-4_, Earthquake Engineering Research Center, University of California, Berkeley, Sept. 56. United States National Research Council (1982), "Earthquake Engineering Research-1982", Report by Committee on Earthquake Engineering Research, National Acadamy Press, Washington, D.C. 57. Vaid, Y.P., and Finn, W.D.L., (1979), "Effect of Static Shear on Liquefaction Potential", Journal of the Geotechnical Engineering Division, ASCE, Vol. 105, GT10, OCT., ppl233-1246. 58. Vaid, Y.P., Byrne, P.M., and Hughes, J.M.O., (1981), "Dilation Angle and Liquefaction Potential", Journal of the Geotechnical Engineering Division, ASCE, Vol. 107, GT7, July. 59. Varadarajan, A., and Mishra, S.S., (1980), "Stress-Path Dependent Stress-Strain Volume Change Behaviour of Granular Soil", International Symposium on Soils under Cyclic and Transient Loading, Swansea, ppl09-119. 60. Verruijt, A., (1977), "Generation and Dissipation of Porewater Pressures", Finite Elements in Geomechanics, Edited by, Gudehus, G., John Wiley and Sons, pp293-319. 61. Wilson, E.L., Farhomand, I., and Bathe, K.J., (1973), "Non-Linear Dynamic Analysis of Complex Structures", International Journal of Earthquake Engineering and Structural Dynamics, Vol. 1, pp241-252. 165 62. Wedge, N.E. (1977), "Problems In Non-Linear Analysis of Movements in Soils", M.A.Sc. Thesis, Submitted to the Dept. of Civil Engineering, University of British Columbia, Vancouver, B.C., Canada. 63. Yoshimi, Y., and Tokimatsu, K. (1977), "Settlement of Buildings on Saturated Sand During Earthquakes", Soils and Foundations, Japanese Society of Soil Mechanics and Foundation Engineering, Vol. 17, No. 1, March, pp23-38. 64. Zienkiewicz, O.C., and Cheung, Y.K., (1967), Finite Element Method in Structural and Continuum Mechanics", McGraw Hill Book Company. 166 APPENDIX I FINITE ELEMENT FORMULATION In the finite element analysis the entire domain of interest is divided into a finite number of elements. The variation of displacement within a finite element is assumed to be given by, {tf} = [N] {6} where (Al.l) {tf} = displacement vector, giving x and y displacements at any point within an element, here TJ*" = {u,v} {6} = displacement vector, giving x and y displacements of the nodes, and [N] = interpolation function. The type of element used in the analysis is a 4 node isoparametric element. The term isoparametric implies common (iso-) parametric description of the unknown displacements and the geometry of the element. The same interpolation functions N^ are used to express both the displacement and the geometry of the element. Isoparametric element formulation has a number of advantages; it offers efficient integrations, and differentiations and it can handle curved and arbitrary geometrical shapes. The interpolation function [N] can be selected such that it can be expressed in natural coordinates (s,t) which is a system of co ordinates intrinsic to an element (Fig. Al.l). 168 In the isoparametric concept, the coordinates of a point in the element is given by, {X} = [N] {X±} (A1.2) If we consider a four-node quadrilateral isoparametric element, the matrix [N] is composed of, „ . d-s(l-t) _ (l+s)(l-t) 1 4 W2 4 N = (i+s)(i+t) _ q-s^ (i+t) 3 4 4 4 Now equations (Al.l) and (A1.2) can be rewritten as, ul vl u2 jU, rN, 0 N2 0 N3 0 N^ 0 -I v2 M L 0 Nx 0 N2 0 N3 0 N^ u3 (A1.3) v3 u4 v4 and, (x} = fNl 0 N2 0 N3 0 Slf 0 , y, V L 0 Nx 0 N2 0 N3 0 N^-1 x3 Xl yi x2 ?3 X4 169 It can be seen that the transformation shown in (A1.4), maps the quadrilateral element into a square as shown in Fig. Al.l. If plane strain conditions are assumed the strain vector e is simply given by _et = {ex, ey, Yxy} where ex and ey are normal strains and yxy 1S tne shear strain. These strains are given in terms of displacements as, x .e. e X _ M ox e y av ay Y ' xy 8u ay ax strain matrix e from A1.3 and A1.5, ex ax 0 ax 0 ax 0 wk ax 0 = 0 ay 0 aN2 ay 0 3N3 ay 0 ay aN1 aN1 aN2 8N2 a\ at^ Y 'xy ey ax ay ax ay ax ay ax > {e} = [B] {6} But = f(s,t) and also x,y are functions of s,t. should be computed using the following relationship, u. u„ (A1.5) 2 (A1.6) (A1.7) So any derivative aN± aN± ax aN± ay as ax as ay as (A1.8) In matrix form, the global and local derivatives can be written as, 170 5N± ax ax aN± 8N± —— as as -— -— as x r-i ax (AI.9) BN. 3x ay. 3N, LJJ oN. lAi.y; at at 8t ay ay where, [j] = Jacobian matrix The derivatives of with respect to x, y can now be obtained from A1.9, as, 8N± aN± 5± • [J]_1 3N, ay at where, [j]-1 = Inverse of Jacobian matrix J, which is simply, at as Now from A1.4, the Jacobian matrix can be written as, ax By. 1^1 t^± R -I _ as as _ i as xi • i as yi ... LJJ " . 8N u 8N (A1.12) ax ay. _i fc f\ at at ^ 5t xi A ot y± 171 The components of matrix [j] can be evaluated since oN^/os, oN^/Qt and (x^, y^) are known. So, knowing the matrix J, [j]-1 also can be computed. Say, = [hi Ti2) •21 (A1.13) 22 It should be noted that these I^j are f(s,t). Therefore, substituting this in equation (ALIO), yields, oN± oNjL By" tl11 L21 22 8N. as" aN. i at (A.14) The matrix [B] (equation A1.6) which relates the strain vector to nodal displacements, has derivatives of interpolation function with respect to x and y. Now knowing these derivatives from equation (A1.4), [B] can be rewritten as, [B] = 3x8 aN, aN •11 0S + X12 9t | o i aN I aN 2i as + 122 at | An bi I aN, aN, + i. 2i as 22 at aN. aN, + i 3N2 3N2 Xn + xi2 a~T aN. 9N„ 12 at | x2i a£ + i 22 at | 172 0 3N3 aN3 0 hi as + at aN2 3N3 3N3 as + I22 at 0 X21 as + I22 at aN2 aN3 m3 aN3 BN3 as + 112 at I21 as + I22 at in as + h2 at aN^ ai^ | [n os + Ti2 at | aN„ aN, o A2i as + x22 at aN4 mh xn as- + xi2 IT I21 as~ + *22 IT | (A1.15) The stresses and strains are connected through elasticity matrix given by, {a'} = [D] {e} (A1.16) where, [a'] = effective stress vector. For 2D plane strain conditions, it is given by, io'}*- = {0X, Oy, Txy} [D] = elasticity matrix for 2D plane strain condition, [D] -B + ^ G Sym B -f G 4 B + j G 0 0 (A1.17) 173 Using the virtual work principle, the internal work (Wj^) done by applying infinitesimal virtual nodal displacement {£} is, WIN = /// M dv (A1.18) where {e} = virtual strains due to virtual displacement {£} and {a} = total stress vector (A1.19) The total stress vector can be split into effective stress and porepressure vectors. i.e.: {a} = {a'} + {uQ} Effective stress Porepressure (A1.19) Vector Vector Here {a'= {ax, ay, Txy} and {xxj*- = {uD, uD, o} in which uQ is the porewater pressure. Now, substituting {a} from equations (Al.19) and (A1.18), one gets, WIN = /// {e}* [{o'} + {u }] dv (A1.20) V Substituting for {a'} from A1.16 this reduces to, WIN = /// M + K) DV (AI.21) v 174 using equation (A1.7), {e} can be replaced by [B] {£}, then, WIN - J/J [B]' [D] [B] {6} + JJJ {6}fc [B]t {uJ dv (A1.21) V V But, external work done by the load vector {p} riding through the virtual displacement {?>}, is simply, WEX = tpJ (A1.22) The virtual work principle gives, WIN WEX i .e. {*}*{*} = {*}'/// [B]t[D][B] dv {6} v (A1.23) + {«}'/// [B]' {uQ} dv V Noting that, dv = dxdydz and also = |j| dsdtdz where dz = thickness of the element. After substituting this in (A1.23), and dividing both sides by {?}t one gets, 175 {P} = [k]{6} + [k*] {uD} (A1.24) where, [k] = element stiffness matrix lit = [/ J [B] [D][B] lJl dsdt] (A1.25) (assuming unit thickness) and, [k ] = porewater pressure matrix 11 t = [/ / [B] |j| dsdt] (A1.26) -1-1 The integrations shown above have to be evaluated numerically. The Gauss integration technique has been employed and the number of points used are 2x2. The formulation presented here is for any linear elastic material. For incrementally elastic analysis, the displacements, stresses and moduli values should be simply replaced by incremental displacements, incremental stresses and tangent moduli. After evaluating the incremental load vector j), element tangent stiffness matrix [kt], porewater pressure matrix [k*], and also estimating the incremental porewater pressure uQ, for all the elements the global incremental load-displacement relationship can be formed. This will lead to, {P} = [Kt]{A} + [K*]{U} (A1.27) 176 in which {P}, [KT], {A}, [K ] and {u} are relevant variables in global axes. By solving this equation the displacement field {A} can be obtained, and it can be used to calculate element strains and stresses using equations (A1.6.) and (A1.16) respectively. Since the shape function gives linear strain variation within an element, the strains and therefore, stresses vary within an element. For convenience, average stress and strain of an element are computed at the centre of gravity of the element. In Chapter 4, and Chapter 5, it is required to express strains and stresses in terms of nodal forces. Recall from equation (A1.24) the nodal forces are given by, {P} = J/J !>]' [D] [B] dv {£} (A1.28) V But, strains are connected to the matrix [BJ in equation (A1.7), {e} = [B] {6} (A1.29) Therefore, from (A1.28) and (A1.29), the nodal forces can be written in terms of strains as, {P} " JJJ [B]1 [D] {e} dv (A1.30) V Now, from equation (A1.16), k'l = [D] {e} (A1.31) 177 and from equations, (A1.30), and (A1.31) the nodal forces can be written in terms of stresses as, {P} = /// [Bf {a'} dv (A1.32) V The equation (A1.30) and (A1.32) can now be used to express strains and stresses in terms of element nodal forces. 178 APPENDIX II STIFFNESS MATRIX FORMULATION FOR THE SLIP ELEMENT As outlined in Chapter 4, that the force-displacement relationship at any point within a slip element has been assumed to be given by, f K 0 w (fS} = US K ] [WS] (A2.1) i.e., f_ = ICQ W where, fg and fn = shear and normal stresses Kg, ^ = joint stiffness in shear and normal directions wg, wn = shear and normal displacements The elastic stored energy, 0E in a slip element due (Fig, (A2.1) to applied forces can be obtained by, 0E = 2 t f dZ (A2* 2) in which L is the total length of the slip element. A factor half is included because the relationship between f_ and w is assumed to be linear. From (A2.1), 0g now can be written as, Fig. A2.1. Slip Element. s = tangential direction Uj = tangential displacement of node i n = normal direction Vj = normal displacement of node i 180 0E = 2 /o W' ko W dA (A2*3) Since the variation of displacements (u, v) within an element is linear the displacement at any point which is at a distance, A, from node I on the bottom edge IJ of the element is, u U) = £ Uj + (1 - L> UI (A2-4> bottom In a similar manner the following equations can be written for, u (A), v (1) and v U), top bottom top l .e. and u U) = L °K + (1 ~ L? "L (A2.5) top v(A) = L VJ + (1 + L} VI (A2.6) bottom v (x) = L VK + (1 " L? VL (A2'7) top where, u^, v^ refer to displacement in tangential and normal direction of the nodes I, J, K and L. Shear and normal displacements at any point are, w = u (A) - u (A) (A2.8) top bottom and, 181 wn = v (A) - v U) (A2.9) top bottom Now substituting for utop> ubottom, vtop and vbottom from equations (A2.4) to (A2.7), in equations (A2.8)-and (A2.9), w = H1 -i> ~i i u - ft] (A.2.10) and, w_ = [-(1 - f) f f (1 (A2.ll) From equation (A2.1), w is, Ul v, w = [ S] = — Lw J -A 0 -B 0 B 0 A 0 [ ] 0 -A 0 -B 0 B 0 A J \ VK \ (A2.12) in which, A = 1 " L AND B - L In matrix form the equation (A2.12) is, 182 w = C 6 2x1 2x§ 8X1 (A2.13) where, = the interpolation matrix and 6^ = nodal displacement matrix Substitution of (A2.13) in (A2.3) gives, K = 9 t sT cl K CJ d* E 2 •'o ooo (A2.14) Performing the matrix multiplication, T C k C ooo -A 0 0 -A -B 0 0 -B K 0 B 0 [ 8 0 B 0 K A 0 n 0 A r-A 0 -B 0 B 0 A On -A 0 -B 0 B 0 A-l K A2 0 ABK 0 -ABK 0 -A2K 0 s s s s 0 A2K 0 ABK 0 -ABK 0 -A2K n B2K n n i ABK 0 0 -B2K 0 -ABK 0 s s s s 0 ABK 0 B2K 0 -B2K 0 -ABK n n n -ABK 0 -B2K 0 B2K 0 ABK 0 s s s B2K s 0 -ABK 0 -B2K 0 0 ABK n n n n -A2K 0 -ABK 0 ABK 0 A2K 0 s s s s 0 -A2K 0 -ABK 0 ABK 0 A2K (A2.15) To perform the integration shown in equation A2.14, one should know integrations of, 183 /L A2 dA, fL B2 dA and fL AB dA Jo J o J o These are simply, JL J o A2 dA = 3 o (1 - ft2 dA L 3 /L J 0 B2 dA = = J o (f)2 dA L 3 (A2.16) and, Jo AB dA = • JL ; o (1-ftf dA 6 Now the equation (A2.14) can be written as, 0E " 2 °T Ksn 6 (A2'17> where, sn 6 2K s 0 K s 0 -K s 0 -2K s 0 2K n 0 K n 0 -K n 0 -2K n 2K s 0 -2K s 0 -K s 0 2K n 0 -2K n 0 -K n Sym 2K s 0 2K n K s 0 2K s 0 K n 0 2K Recall that elastic stored energy 0g in the formulation of a linear elastic finite element is, 0E = j 6T K 6 (A2.18) 184 Where K is the stiffness matrix of the finite element. Now, comparing (A2.17) and (A2.18), the stiffness matrix for slip element can be deduced as K„_. —sn 185 APPENDIX III STEP BY STEP INTEGRATION For the proposed incrementally elastic dynamic analysis in the time domain the "elastic" properties have to be modified for every time step. The elastic parameters depend on the level of strain in the deposit. Therefore, the displacement field in the deposit should be evaluated at every time step. This requires that the incremental dynamic equilibrium equations (equation 5.5) have to be solved numerically for every time step. Newmark's method (1959) of step by step integration is very popular and extensively used in dynamic analyses. This method basically provides numerical solution in time domain, where the solution is advanced by one discrete step at a time. In this method, two parameters <* and p are used so that the velocity and displacement at time t+At can be expressed in terms of acceleration, velocity and displacement at time t, and of the known acceleration at time t+At. For convenience let us define that, T = t+At Then the relationship in terms of <* and 8 are, {X}T = {X}t + (1 - «) At {X}t + ocAt {X}T (A3.1) and, 186 {X}T = {X}t + At {X}t + (± - 6) (At)2 {X}t + B At2 {X}T (A3.2) Newmark (1959) proposed that « = 1/2 and B = 1/4 be used for an unconditionally stable integration procedure, which incidentally corresponds to a constant average acceleration method of integration. If = = 1/2 and 8 = 1/6 are used then this method gives a linear variation of acceleration within the time step. Re-writing the incremental equilibrium equations from Chapter 5, [M] {AX} + [C] {AX} + [Kt]t {AX} = {AP} (A3.3) Now substituting for, {AX} {AX} {AX} (A3.4) and, in equation (A3.3), one gets, [M] {^ - Xj + [C] {Xj, - Xt} + [Kt]t {^ - Xt} = {AP} (A3.5) From equations (A3.1) and (A3.2), {x}^, and {^}^ can be expressed in terms of other variables as follows, 187 [|£2+[c]p£+[*t]t] (AX} = {AP} + [M] {E}t + [C] {F}t WT= m2 [{AX} " At{x}t" (i ~ P> At2 W J (A3*6) and, WT - (X}t + (!--) At {X} (A3.7) + ^ [{AX} - At{x}t - B) At* {X}J Substituting for {x}-£ and {x},p in equation (A3.5), W2 [ {AX} - At{X}t - (j - 6) At* {X}t - BAt2 {x}J + [C] [(1—) At {X}fc + ^ [{AX} - At {X}t (A3.8) - £ - B) At2 {X}J + [Kt]t {AX} = {AP} Collecting terms, and defining following simplifying symbols, and, . Wt = f Wt +'(2p " 1} At ^t (A3'10) the equation (A3.8) can be reduced to, (A3.11) 188 Recall from Chapter 5, Section 5.1.2, correction forces should be applied to restore total equilibrium at time T. The {Pcorr} evaluated using equation (5.11) can be added to right hand side of the equation (A3.11). Then the equation (A3.11) is, ^+[c]ptr+[KT]T] w = {AP} + [M] {E}t + [C] {F}t (A3.12) where, {AP} = {AP} + {Pcorr} (A3.13) The only unknown in the above equation is {AX} and therefore, {AX} can be obtained as, {AX} = [D]"1 [{AP} + [M] {E}t + [C] {F}t] (A3.14) in which, [»]-[{S.]*WJE+["tU (A3-15> Now knowing {AX}, the unknowns {x}^ and {x}T and {x}T can be evaluated. {x}T is simply, {X}T = {AX} + {X}t (A3.16) 189 From equation (A3.2), an expression for {x}T is, WT • "p^2 [{AX} - At {X}t " (| - 8) At2 {X}t] (A3.17) Substituting for {E}t from (A3.9), the equation (A3.17) can be simplified as, WT = "pit"2 {AX} - {E}t + {X}t (A3.18) From equations (A3.1) and (A3.2) an expression for {x}^, after rearranging terms is, {X}T = {X}t + (I-) At {X}t + ccAt {X}T (A3.19) Knowing {AX} by solving equation (A3.14), the response at time T can be computed using equations, (A3.18) and (A3.19). In the numerical step by step integration, the following sequence of calculations have to be performed for every time step. 1. Initial velocity {x}t and displacements {x}t are known either from values at the end of the preceding increment or as initial conditions of the problem. Based on these values, {E}t, {F}t, and {pCOrr}» are evaluated using equations (A3.9), (A3.10) and (5.11). 2. With these values and the known non-linear properties of the soil deposit, the damping matrix [c] and [Kt]t are evaluated according to appropriate equations in Chapter 5. 190 3. The matrix [D] is then calculated using equation (A3.15). 4. Using the increment in base acceleration value at time t, it is possible to evaluate the right hand side of the equation (A3. 12). 5. The equation (A3.12) can now be solved for {AX} and then displacement, acceleration and velocity vectors can be evaluated from equations (A3.16), (A3.18) and (A3.19) respectively. When step 5 has been completed, the analysis for this time increment is finished and the entire process may be repeated for the next time step. Obviously this process can be carried out consecutively for any desired number of time increments; thus the complete response history can be computed. Two important aspects have to be considered in any numerical integration procedure. They are the accuracy and stability of the procedures. Accuracy refers to how well the numerical solution matches the exact continuous solution. Stability refers to whether extraneous solutions are introduced in such a way that they increase rather than decay, and thus come to dominate the results. Usually there is an upper limit to At that is necessary to guarantee stability, and the value of that limit depends on the type of element stiffness and mass matrix as well as on « and 8. With a linear acceleration assumption (<= = 1/2, 6 = 1/6) the analysis will give good accuracy if the shortest period of the deposit is 5 to 10 times greater than At (Clough and Penzien 1975). The linear acceleration method is only conditionally stable, and it will blow up if it is applied to soil structures with the shortest period less than about 1.8 times the integration interval. Thus the time increments must be made 191 short relative to the least period of vibration contained in the system regardless of whether the higher modes contribute significantly or not. In the analysis using finite element procedures the shortest periods of vibration in general may be several orders of magnitude less than the periods associated with the significant response. In these cases, the linear acceleration method cannot be used because of the very short time step required to avoid instability; instead, an unconditionally stable method is required which will not blow up regardless of the time step. Several unconditionally stable methods are available. The constant average acceleration method (« = 1/2 , 8 = 1/4) is one of the simplest of these methods. But this assumption has been reported not to give good results than the methods with linear acceleration assumption. The Wilson 9-method (Wilson, et al. 1973), is also an unconditionally stable method. This method is a modification of linear acceleration method and is reported to be the best of all unconditionally stable methods (Clough, et al. 1975). The Wilson 8-method is based on the assumption that the acceleration varies linearly over an extended computation interval, %, such that, x = 9 At where 0 > 1.37 (A3.20) When 0=1, this method reverts to the standard linear acceleration method. The analysis procedure is exactly the same as the procedure presented above except that in the equations the time step At has to be replaced by T and also the equations to evaluate {x}-j,, {&} ^, and {x}T have to be modified. 192 Since this is essentially a linear acceleration method, here = 1/2 and 8=1/6. By inspection, the required equations can be rewritten as follows. The equation (A3.12) can be rewritten as, f^M] + 1[C] + [Kfc]t] {AX} m {~p} + [M] [c] {p}t (A3>21) in which, {E*t = x ^t + 3 Mt (A3'22) and, {F}t = 3 {X}t + \ {X}t (A3.23) After evaluating {AX}, which is over an extended time increment the displacement, velocity and acceleration values should be computed at time t = T. The acceleration at t = T can be evaluated using, WT = [*2 {AX} + {E}J {X}t (A3.24) {x},p, from equation (A3.19) with « = 1/2 and 8 = 1/6 is, {X}T = Wt +f [Wt + {X}T] (A3.25) 193 and also the displacement {x}^, from (A3.2) is, {X}T = (X}t + At (X}t + *f - [{x}T + 2 {X}J (A3.26) f It must be remembered that stability in numerical integration does not guarantee accuracy or vice versa. The Wilson 0-method imposes artificial damping in higher modes. But knowing that the response due to higher modes of vibration contributes very little to the true response of structures, this method in a way filters out the high frequency response. Therefore, this method has been found to yield realistic results in a number of dynamic analyses.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- A two-dimensional non-linear static and dynamic response...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
A two-dimensional non-linear static and dynamic response analysis of soil structures Siddharthan, Rajaratnam 1984-06-13
pdf
Page Metadata
Item Metadata
Title | A two-dimensional non-linear static and dynamic response analysis of soil structures |
Creator |
Siddharthan, Rajaratnam |
Publisher | University of British Columbia |
Date Issued | 1984 |
Description | A method of analysis in two-dimensions to predict static and dynamic response of soil structures, including soil-structure interaction has been presented herein. The static and dynamic analyses can be performed in either effective or total stress mode or a combination of both modes. Non-linear stress-strain behaviour of soil has been modelled by using an incrementally elastic approach in which tangent shear modulus and tangent bulk modulus were taken as the two "elastic" parameters. The material response in shear was assumed to be hyperbolic coupled with Masing behaviour during unloading and reloading. Responses to changes in mean normal stress was assumed to be non-linear, elastic and stress dependent. Slip or contact elements have been incorporated in the analysis to represent the interface characteristics between soil and structural elements. The properties of the slip elements were assumed to be elastic, perfectly plastic, with failure at the Interface given by the Mohr-Coulomb failure criterion. In the static analysis proposed here, gravity may be switched on at once for the completed soil structure or the construction sequence can be modelled by layer analysis. The stress-strain conditions determined by the static analysis give the in-situ stress condition before the dynamic analysis. In the dynamic effective stress analysis, the residual porewater pressures are calculated using a modification of the model proposed by Martin, et al. (1975). The parameters, G[sub max] and[sub max] are modified for the effects of residual porewater pressure. The dynamic response study includes the prediction of post earthquake deformations. The predictive capability of the new method of analysis has been verified by comparing the recorded porewater pressure and accelerations of two centrifuged models subjected to simulated earthquakes, to those computed by the new method. This method has also been used to compute response of an offshore drilling island supporting a tanker mounted drilling rig. Results suggest that the common practice of neglecting soil-structure interaction may not be appropriate for islands which support heavy tanker type of structures. At present one-dimensional methods are used for computing the response of these islands. Comparative studies are also reported to asses the validity of this procedure. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-06-13 |
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.0062642 |
URI | http://hdl.handle.net/2429/25677 |
Degree |
Doctor of Philosophy - PhD |
Program |
Civil Engineering |
Affiliation |
Applied Science, Faculty of Civil Engineering, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-UBC_1984_A1 S56.pdf [ 7.41MB ]
- Metadata
- JSON: 831-1.0062642.json
- JSON-LD: 831-1.0062642-ld.json
- RDF/XML (Pretty): 831-1.0062642-rdf.xml
- RDF/JSON: 831-1.0062642-rdf.json
- Turtle: 831-1.0062642-turtle.txt
- N-Triples: 831-1.0062642-rdf-ntriples.txt
- Original Record: 831-1.0062642-source.json
- Full Text
- 831-1.0062642-fulltext.txt
- Citation
- 831-1.0062642.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-0062642/manifest