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 B r i t i s h 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 t h i s thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA MARCH 1984 (c) Rajaratnam Siddharthan, 1984 In presenting t h i s thesis i n p a r t i a l f u l f i l m e n t of the requirements for an advanced degree at the University of B r i t i s h Columbia, I agree that the Library s h a l l make i t f r e e l y a v a i l a b l e for reference and study. I further agree that permission f o r extensive copying of t h i s t h e s i s f o r scholarly purposes may be granted by the head of my department or by h i s or her representatives. I t i s understood that copying or p u b l i c a t i o n of t h i s t h e s i s f o r f i n a n c i a l gain s h a l l not be allowed without my written permission. Department of C l y\/^ <C^v-^) \ r^Lfi-fi The University of B r i t i s h Columbia 1956 Main Mall Vancouver, Canada V6T 1Y3 Date i i ABSTRACT A method of analysis i n two-dimensions to predict s t a t i c and dynamic response of s o i l structures, including s o i l - s t r u c t u r e i n t e r a c t i o n has been presented herein. The s t a t i c and dynamic analyses can be performed i n e i t h e r e f f e c t i v e or t o t a l stress mode or a combination of both modes. Non-linear s t r e s s - s t r a i n behaviour of s o i l has been modelled by u s i n g an incrementally e l a s t i c approach i n which tangent shear modulus and tangent bulk modulus were taken as the two " e l a s t i c " parameters. The material response i n shear was assumed to be hyperbolic coupled with Masing behaviour during unloading and reloading. Responses to changes i n mean normal stress was assumed to be non-linear, e l a s t i c and stress dependent. S l i p or contact elements have been incorporated i n the analysis to represent the i n t e r f a c e c h a r a c t e r i s t i c s between s o i l and s t r u c t u r a l elements. The properties of the s l i p elements were assumed to be e l a s t i c , p e r f e c t l y p l a s t i c , with f a i l u r e at the Interface given by the Mohr-Coulomb f a i l u r e c r i t e r i o n . In the s t a t i c analysis proposed here, gravity may be switched on at once for the completed s o i l structure or the construction sequence can be modelled by layer a n a l y s i s . The s t r e s s - s t r a i n conditions determined by the s t a t i c analysis give the i n - s i t u stress condition before the dynamic an a l y s i s . In the dynamic e f f e c t i v e stress a n a l y s i s , the r e s i d u a l porewater pressures are c a l c u l a t e d using a modification of the model proposed by M a r t i n , et a l . (1975). The parameters, G m a x and _ a x are modified for the e f f e c t s of r e s i d u a l porewater pressure. The dynamic i i i response study includes the p r e d i c t i o n of post earthquake deformations. The p r e d i c t i v e c a p a b i l i t y of the new method of analysis has been v e r i f i e d 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 d r i l l i n g i s l a n d supporting a tanker mounted d r i l l i n g r i g . Results suggest that the common pract i c e of neglecting s o i l - s t r u c t u r e i n t e r a c t i o n 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 i s l a n d s . Comparative studies are also reported to asses the v a l i d i t y of t h i s procedure. i v TABLE OF CONTENTS PAGE ABSTRACT i i TABLE OF CONTENTS i v LIST OF FIGURES v i i i LIST OF TABLES x i v NOMENCLATURE xv ACKNOWLEDGEMENTS i x 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 Stre s s - S t r a i n 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 Str a i n Hardening 25 3.1.5 - D i s s i p a t i o n 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 38 4.1.1 - Reasons for Selecting G t and B t 39 4.1.2 - Hyperbolic Shear Stre s s - S t r a i n Relationship 40 4.1.2.1 - Estimation of G m a v 40 4.1.2.2 - Estimation of x m a x 43 4.1.2.3 - Influence of Over-Consolidation 45 4.1.2.4 - E f f e c t s of Unloading 45 4.1.3 - Tangent Bulk Modulus B t 47 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 S l i p Element 68 4.6 - SELECTION OF SOIL PARAMETERS 70 4.6.1 - Obtaining Shear Stre s s - S t r a i n Parameters 71 4.6.2 - Obtaining Bulk Modulus Parameters 72 v i PAGE CHAPTER 5 - TWO-DIMENSIONAL DYNAMIC ANALYSIS 73 5.1 - FORMULATION OF THE PROBLEM 73 5.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 E f f e c t s of Residual Porewater Pressure 85 5.2.3.1 - The Behaviour of Samples with T s 86 5.2.3.2 - Residual Porewater Pressure Generation Model 95 5.3 - DAMPING MATRIX [Cj 99 5.4 - BOUNDARY CONDITIONS 102 5.5 - ANALYSIS OF SOIL STRUCTURE SYSTEMS 104 5.5.1 - S l i p Element i n 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 108 6.2 - COMPARATIVE STUDY I l l 6.2.1 - Results of Test 1 115 6.2.2 - Results of Test 2 123 v i i 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 7.2 - ANALYSIS OF A TYPICAL TANKER ISLAND 138 7.2.1 - R e s u l t s of Tanker I s l a n d P roblem 142 7.3 - SOME PRACTICAL CONCLUSIONS 151 CHAPTER 8 - SUMMARY AND CONCLUSIONS 153 8.1 - SUMMARY 153 8.2 - CONCLUSIONS 154 8.3 - SUGGESTIONS FOR FURTHER STUDY 156 REFERENCES ' 157 APPENDIX I 166 APPENDIX I I 178 APPENDIX I I I 185 v i i i LIST OF FIGURES PAGE FIGURE 1.1 - Caisson-Retained Island 3 (De Jong and Bruce, 1978) FIGURE 2.1 - Conversion of Shear St r a i n P o t e n t i a l to Equivalent Shear Forces 10 FIGURE 3.1(a) - I n i t i a l Loading Curve 19 FIGURE 3.1(b) - Masing St r e s s - S t r a i n Curves for Unloading and Reloading 19 FIGURE 3.2 - Hysteretic C h a r a c t e r i s t i c s 20 FIGURE 3.3(a) - I n i t i a l E f f e c t i v e Stress Condition i n a Simple Shear Apparatus 23 FIGURE 3.3(b) - Intermediate E f f e c t i v e Stress Condition i n a Simple Shear Apparatus 23 FIGURE 3.4 - Relationship Between Volumetric Strains and Porewater Pressures i n Constant S t r a i n C y c l i c Simple Shear Tests, D_ = 45% 29 FIGURE 3.5 - C y c l i c Stress Ratio Vs. Number of Cycle for I n i t i a l Liquefaction for Various Sands 31 FIGURE 3.6 - Predicted and Measured Porewater Pressures i n Constant Stress C y c l i c Simple Shear Tests, D r = 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 44 i x LIST OF FIGURES PAGE FIGURE 4.2 - Mohr C i r c l e Diagram 44 FIGURE 4.3 - Mohr Envelope 46 FIGURE 4.4 - E f fec t s of Loading and Unloading 46 FIGURE 4.5 - Sequentia l Procedure i n Dam Construct ion 51 FIGURE 4.6 - Schematic Diagram Showing the Incremental A n a l y t i c a l Procedure 53 FIGURE 4.7(a) - T y p i c a l P lot s for vs for Dense and Loose Sands 56 FIGURE 4.7(b) - T y p i c a l P lo t s for y vs for Dense and .Loose Sands 56 FIGURE 4.8 - S t r e s s - S t r a i n Behaviour of Ottawa Sand i n Drained Simple Shear (After Vaid et a l . 1 9 8 1 ) . . . . 58 FIGURE 4.9 - T y p i c a l Drained T r i a x i a l Test Results on (a ) , (b) Dense Sacramento River Sand a) P r i n c i p a l Stress Ratio vs A x i a l S t r a i n 59 b) Volumetric S t ra in vs A x i a l S t r a in 59 (After Lee, 1965) FIGURE 4.10 - V a r i a t i o n of D i l a t i o n Angle with Mean Normal Stress for Various Sands (After Robertson, 1982) 60 FIGURE 4.11 - Idea l i sed Rela t ionship Between ^ and 62 FIGURE 4.12 - D e f i n i t i o n of S l i p Element 65 FIGURE 4.13 - P lo t of T y p i c a l Shear/Normal Stress vs Shear/ Normal Displacement 65 X LIST OF FIGURES PAGE FIGURE 5.1 - Dynamic Shear S t re s s -S t ra in R e l a t i o n s h i p : Skeleton Curve 84 FIGURE 5.2 - Dynamic Shear S t re s s -S t ra in Re la t ionsh ip : Unloading and Reloading 84 FIGURE 5.3 - Permanent and C y c l i c Stra ins or Permanent and C y c l i c Porewater Pressures 88 FIGURE 5.4 - C y c l i c and Residual Behaviour of Pore Pressure and A x i a l S t r a i n for ICT and ACT Tests 88 FIGURE 5.5 - q vs p' P lot for ICT and ACT Tests 91 FIGURE 5.6 - Normalized Curves for Pore Pressure Bui ld-up Equation 94 FIGURE 5.7 - Stregth Curves for (a) Ekof i sk Sands with D„ = 85%; (b) Sand with D_ = 77% 94 (After Rahman, et a l . 1977) FIGURE 6.1 - Submerged Is land Showing Transducer Pos i t ions . . . 109 FIGURE 6.2(a) - Input Base Motion i n Test 1 112 FIGURE 6.2(b) - Input Base Motion i n Test 2 112 FIGURE 6.3 - L ique fac t ion Resistance Curves of Medium Dense Leighton-Buzzard Sand 114 FIGURE 6.4(a) - Recorded A c c e l e r a t i o n of ACC1244 i n Test 1 116 (With and Without S l i p Elements) FIGURE 6.4(b) - Computed A c c e l e r a t i o n of ACC1244 i n Test 1 116 (With and Without S l i p Elements) FIGURE 6.5(a) - Recorded A c c e l e r a t i o n of ACC1225 i n Test 1 117 x i LIST OF FIGURES PAGE FIGURE 6.5(b) - Computed A c c e l e r a t i o n of ACC1225 i n Test 1 117 (With and Without S l i p Elements) FIGURE 6.6(a) - Recorded a c c e l e r a t i o n of ACC734 i n Test 1 118 FIGURE 6.6(b) - Computed A c c e l e r a t i o n of ACC734 i n Test 1 118 (With and Without S l i p Elements) FIGURE 6.7(a) - Recorded and Computed Porewater Pressure of PPT2330 i n Test 1 120 FIGURE 6.7(b) - Recorded and Computed Porewater Pressure of PPT68 i n Test 1 120 FIGURE 6.7(c) - Recorded and Computed Porewater Pressure of PPT2338 i n Test 1 121 FIGURE 6.7(d) - recorded and Computed Porewater Pressure of PPT2342 i n Test 1 121 FIGURE 6.8(a) - Recorded A c c e l e r a t i o n of ACC1244 i n Test 2 125 FIGURE 6.8(b) - Computed A c c e l e r a t i o n of ACC1244 i n Test 2 125 (With S l i p Elements) FIGURE 6.9(a) - Recorded A c c e l e r a t i o n of ACC1225 i n Test 2 126 FIGURE 6.9(b) - Computed A c c e l e r a t i o n of ACC1225 i n Test 2 126 (With S l i p Elements) FIGURE 6.10(a)- Recorded A c c e l e r a t i o n of ACC734 i n Test 2 127 FIGURE 6.10(b)- Computed A c c e l e r a t i o n of ACC734 i n Test 2 127 (With S l i p Elements) FIGURE 6.11(a)- Recorded and Computed Porewater Pressure of PPT2330 i n Test 2 130 x i i LIST OF FIGURES PAGE FIGURE 6.11(b)- Recorded and Computed Porewater Pressure of PPT68 i n Test 2 130 FIGURE 6.11(c)- Recorded and Computed Porewater Pressure of PPT2338 i n Test 2 131 FIGURE 6.11(d)- Recorded and Computed Porewater Pressure of PPT2342 i n Test 2 131 FIGURE 6.12(a)- E f f e c t i v e Stress Paths i n Test 1 135 FIGURE 6.12(b)- E f f e c t i v e Stress Paths i n Test 2 136 FIGURE 7.1 - Schematic of Tanker Island 139 FIGURE 7.2 - D i s t r i b u t i o n of Maximum Dynamic Shear Ratio i n Sand 143 FIGURE 7.3 - D i s t r i b u t i o n of Maximum Dynamic Shear Stress Ratio i n Sand 143 FIGURE 7.4 - D i s t r i b u t i o n of Residual Porewater Pressure Ratio i n Sand 145 FIGURE 7.5 - D i s t r i b u t i o n of Maximum Dynamic Shear S t r a i n i n Sand 145 FIGURE 7.6 - D i s t r i b u t i o n 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 - D i s t r i b u t i o n of Residual Porewater Pressure Ratio i n Sand 150 FIGURE A l . l - Iso-Parametric Element 167 x i i i LIST OF FIGURES PAGE FIGURE A2.1 - S l i p Element 179 x i v LIST OF TABLES PAGE TABLE 4.1 V a r i a t i o n of Exponent r\, with P l a s t i c Index PI 42 TABLE 6.1 S o i l 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 S t a t i c S o i l Properties 138 TABLE 7.2 Dynamic S o i l Properties 140 TABLE 7.3 ^ s ^ a v o a n d K r V a l u e s 1 4 1 NOMENCLATURE xv A = Variables defined i n s l i p element s t i f f n e s s matrix formulation i n Appendix I I . a = Constant used to compute damping matrix, a^ - a^ = Constants defined i n text. B = Variable defined i n s l i p element s t i f f n e s s matrix formulation i n Appendix I I . B = Strain-displacement matrix. B t = Tangent bulk modulus. b = Constant used to compute damping matrix, b^ - b^ = Constants defined i n text. [C] = Damping matrix. C_ = Interpolation matrix i n s l i p element s t i f f n e s s matrix formulation. C^ - C^ = Volume change constants. C g = Cohesion i s s l i p element. C' = E f f e c t i v e Cohesion. I) = E l a s t i c i t y matrix. D r = Relative density of s o i l . d = Depth of centre of gravity of an element below the top surface. {E} = Matrix defined i n Appendix I I I . E t = Young's modulus. E_ = One-dimensional rebound modulus, e = Void r a t i o . NOMENCLATURE Maximum void r a t i o . Minimum void r a t i o . Matrix defined i n Appendix I I I . Column vector of global i n e r t i a 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 s l i p element. Shear force per unit area i n a s l i p element. Shear force per unit area i n a s l i p element. Tangent shear modulus for loading and unloading. Tangent shear modulus. Maximum shear modulus. Hardening constants. Unit vector. Matrix defined i n Appendix I. Jacobian matrix. Non-dimensional shear modulus constant. Bulk modulus constant. A n i s t r o p i c consolidation r a t i o . Normal j o i n t s t i f f n e s s / u n i t area. Experimental constant used i n porewater pressure model. Shear j o i n t s t i f f n e s s unit area. S l i p element s t i f f n e s s matrix i n element a x i s . Global tangent s t i f f n e s s matrix. x v i i NOMENCLATURE [K x y ] = S l i p element s t i f f n e s s matrix i n xy ax i s . K Q = C o - e f f i c i e n t of l a t e r a l earth pressure at r e s t . ( K 2 ) m a x = Shear modulus constant used to cal c u l a t e G_ a x. K = Constant which depends on the type of m a t e r i a l . * [K ] = Porewater pressure matrix. K = Shear modulus constant for clay. K = P e r m i a b i l i t y of s o i l , z J L = Length of a s l i p element. H = Distance of a point of s l i p element i n the element d i r e c t i o n . [M] = Mass matrix. m = Experimental constant used i n porewater pressure model. N = Scale factor of the c e n t r i f u g a l a c c e l e r a t i o n . [N] = Interpolation matrix. N^ = Number of cycles to cause i n i t i a l l i q u e f a c t i o n . N 5 Q = Number of cycles to cause a porewater pressure r a t i o of 50%. n = Experimental constant used i n porewater pressure model. n = Bulk modulus exponent. OCR = Over consolidation r a t i o . {P} = Global column vector of applied loads. Pfl = Atomospheric pressure. {P } = Correction force vector. 1 corr-" P Q = Transformation matrix i n s l i p element, p' = Average e f f e c t i v e p r i n c i p l e s t r e s s . x v i i i NOMENCLATURE 0^,0^ = Constant matrices defined i n the text. q = P r i n c i p l e stress d i f f e r e n c e . R^ ,R_2 = Constant matrices defined i n the text. T_ = Transformation matrix. t = time. U = Residual porewater pressure. U = Maximum r e s i d u a l porewater pressure, max r {U} = Column vector of incremental porewater pressures i n elements. {TJ} = Displacement vector. u^ = Tangential displacement of node i . {u Q} = 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 i n a s l i p element. w = Shear dispalcement i n a s l i p element. s {X},{X},{X} = Column vectors whose components are r e l a t i v e displacement, v e l o c i t y and a c c e l e r a t i o n of nodes. X^ = Base a c c e l e r a t i o n , x = Horizontal distance, z = V e r t i c a l distance. i x X ACKNOWLEDGEMENTS I would l i k e to extend my sincere gratitude and appreciation to my research supervisor, Professor W.D. Liam Finn, for his guidance and keen i n t e r e s t and for making many valuable suggestions to improve the presentation of t h i s t h e s i s . 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 a v a i l a b l e to discuss any problem and for h i s 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 B r i t i s h Columbia. The f i n a n c i a l 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 C a l i f o r n i a and Exxon Production Research Company, Houstan, Texas, i s g r a t e f u l l y acknowledged. F i n a l l y , I wish to thank my wife Visha, for her patience, s a c r i f i c e s and for the excellent job i n typing t h i s t h e s i s . This thesis would not have seen t h i s day without her whole-hearted support. I also remain indebted to my parents for t h e i r good wishes and continued moral support. It was due to t h e i r blessings and encouragement that I have completed t h i s work. 1 CHAPTER 1 INTRODUCTION In engineering p r a c t i c e , i t i s generally agreed that the performance of s o i l structures subjected to seismic loading should be evaluated i n terms of deformations rather than i n 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 s i x t i e s many a n a l y t i c a l methods for assessing earthquake induced deformations i n s o i l structures have been proposed. These methods may be c l a s s i f i e d into two broad categories: one-dimensional methods and two-dimensional methods. The one-dimensional methods assume that deformations occur i n p a r a l l e l 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 a l . (1966), assume that f a i l u r e develops along w e l l -defined f a i l u r e planes and compute displacements of r i g i d blocks of s o i l s . More recently the method proposed by I a i and Finn (1982) for long slopes accounts for the f l e x i b i l i t y of the s o i l deposit. Often s o i l deposits cannot be characterized adequately as one-dimensional deposits, and the v a r i a b i l i t y of properties i n two or even three-dimensions must be considered. For example, i n the response analyses of zoned dams, two or three-dimensional analyses are e s s e n t i a l . When the th i r d dimension i s very much larger than the other two-dimensions and the properties do not vary s i g n i f i c a n t l y i n t h i s d i r e c t i o n , a 2 two-dimensional response analysis i s usually adequate. The geometry of the modes of deformations may also d i c t a t e two or three-dimensional a n a l y s i s . For example i n the analysis of embedded structures rocking may be an important deformation mode i n addition to t r a n s l a t i o n and therefore, at least a two-dimensional analysis i s necessary. There are a number of two-dimensional methods a v a i l a b l e to compute seismic deformations. Some of these methods are based on e l a s t i c -p l a s t i c s o i l behaviour (Finn, et a l . 1973; Mroz, et a l . 1979; Prevost, 1979). These methods are complicated to use and have had very l i m i t e d v a l i d a t i o n . The method proposed by Seed, et a l . (1973), to compute seismic deformations of earth dams has found wide a p p l i c a t i o n i n p r a c t i c e . This i s a semi-analytical method i n which the r e s u l t s of stress analysis and data from c y c l i c t r i a x i a l tests are used to estimate p o t e n t i a l displacements i n dams. Non-linearity of the s o i l i s taken into account using an i t e r a t i v e e l a s t i c approach to achieve s o i l properties compatible with the computed s t r a i n s . In recent years new types of structures have emerged, for which i t i s important to determine deformations under earthquake loading. Examples are the man-made sand islands which support d r i l l i n g platforms for gas and o i l explorations i n the Beaufort Sea. These islands carry d r i l l i n g equipment on ei t h e r s a n d - f i l l e d caissons or tankers ( F i g . 1.1). The deformations of these structures and islands during earthquakes are an important design consideration. A general s a t i s f a c t o r y method for computing deformations of these structures i s not a v a i l a b l e at present. Indeed the state-of-the-art for analysing deformations was recently assessed i n a report on earthquake engineering research by the F i g . 1.1. Ca i s son-Reta ined I s land (De J o n g and Bruce, 1978). 4 National Research Council of the United States (USNRC, 1982; Finn, 1983) i n the following terms: "Many problems i n s o i l mechanics, such as safety studies of earth dams, require that the possible permanent deformations that would be produced by earthquake shaking of prescribed i n t e n s i t y and duration be evaluated. Where f a i l u r e develops along w e l l - d e f i n e d f a i l u r e p l a nes, r e l a t i v e l y simple e l a s t o p l a s t i c models may s u f f i c e to c a l c u l a t e displacements. However, i f the permanent deformations are d i s t r i b u t e d throughout the s o i l , the problem i s much more complex, and p r a c t i c a l , r e l i a b l e methods of analysis are not a v a i l a b l e . F uture progress w i l l depend on development of s u i t a b l e p l a s t i c i t y models for s o i l undergoing r e p e t i t i v e loading. This i s c urrently an important area of research". For r e a l i s t i c predictions of stresses and displacements i n s o i l structures the s t r e s s - s t r a i n behaviour of s o i l s should be modelled as c l o s e l y as p o s s i b l e . This i s of course a d i f f i c u l t task since the s t r e s s -s t r a i n behaviour of s o i l s i s extremely complex. Using simple stressv s t r a i n r e l a t i o n s h i p s , a two-dimensional method for computing transient and permanent deformations i n s o i l structures i s presented i n t h i s t h e s i s . 1.1 SCOPE OF THIS THESIS This thesis presents a method for two-dimensional s t a t i c and seismic response analysis of s o i l structures. The earthquake loading occurs a f t e r s t a t i c equilibrium or steady state conditions have been established. S o i l properties such as strength and s t i f f n e s s which control the response to earthquake loading depend on the e f f e c t i v e stresses i n the s o i l s t r u c t u r e . Therefore, i t i s important to evaluate i n - s i t u s t a t i c stresses. These stresses are determined by a s t a t i c analysis which takes into account non-linear stress-dependent 5 response of s o i l to load. Because s o i l behaviour depends on the loading path, the construction sequence of the s o i l structure i s c a r e f u l l y modelled. A number of s a t i s f a c t o r y methods of s t a t i c analyses are already a v a i l a b l e (Kulhawy, et a l . 1969; Duncan, et a l . 1978). Nevertheless an independent method i s presented here which uses a consistent set of material parameters i n both s t a t i c and dynamic analyses. This r e s u l t s In a much more e f f i c i e n t , cost e f f e c t i v e s o l u t i o n to the problem of dynamic response a n a l y s i s . The method for dynamic analysis takes into account the non-l i n e a r h y s t e r e t i c s t r e s s - s t r a i n behaviour of s o i l s . The analysis may be c a r r i e d out i n either an e f f e c t i v e stress or t o t a l stress mode, using an appropriate s t r e s s - s t r a i n r e l a t i o n . For e f f e c t i v e stress response a n a l y s i s , r e s i d u a l porewater pressures must be known. Therefore, a porewater pressure generation model has been developed for pre d i c t i n g s e i s m i c a l l y induced porewater pressures. The porewater pressure model i s a generalization of the one-dimensional model of Martin, et a l . (1975). The p r e d i c t i v e c a p a b i l i t y of the new method for dynamic analysis has been v e r i f i e d 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 d r i l l i n g i s l a n d supporting a tanker-mounted d r i l l i n g r i g . The porewater pressures, stresses, accelerations and displacements i n the i s l a n d have been determined. At present one-dimensional methods are used for computing the response of these i s l a n d s . Comparative studies were conducted to assess the v a l i d i t y of t h i s procedure. Development of a non-linear method of analysis i s very 6 d i f f i c u l t . A number of approximations have to be made to achieve a p r a c t i c a l useful program. These approximations have been examined at length i n this thesis and suggestions have been made for future research. 1.2 ORGANIZATION OF THESIS A c r i t i c a l review of the method proposed by Seed, et a l . (1973) for computing seismic deformations i s presented i n Chapter 2. The formulations, basic assumptions and l i m i t a t i o n s of the model developed for the analysis of s t a t i c response are presented i n Chapter 3. The complete treatment of s o i l - s t r u c t u r e i n t e r a c t i o n has also been included. The proposed two-dimensional dynamic response analysis i s an extension of the one-dimensional response analysis of Finn, et a l . (1977). Therefore, a d e t a i l e d d e s c r i p t i o n of t h e i r model, i t s a p p l i c a t i o n i n p r a c t i c e and i n the laboratory i s given i n Chapter 4. The proposed method for dynamic response analysis i s presented i n Chapter 5. D e t a i l s of the v e r i f i c a t i o n of the method i s presented i n Chapter 6. The method i s used to compute response of a t y p i c a l d r i l l i n g tanker i s l a n d subjected to seismic loading. The r e s u l t s of the analysis including implications for engineering design are discussed i n Chapter 7. A b r i e f summary, suggestions for future work and conclusions are given i n Chapter 8. 7 CHAPTER 2 CRITICAL REVIEW OF SEED, ET AL. METHOD FOR COMPUTING DYNAMIC DEFORMATIONS The s t a t e - o f - t h e a r t r e p o r t on a n a l y s i s of permanent deformations i n earth structures (USNRC, 1982), suggests that the use of simple e l a s t o p l a s t i c models for computing deformations may be adequate i f deformations develop along well defined s l i p planes. When a well defined s l i p surface does not occur and deformations are d i s t r i b u t e d through out the s o i l structure, an analysis at least i n two-dimensions i s necessary. The most widely used two-dimensional method i s the one that was proposed by Seed, et a l . (1973, 1979). Detailed d e s c r i p t i o n and l i m i t a t i o n s of t h e i r method are presented below. 2.1 SEED, ET AL. METHOD (1973, 1979) The basic steps i n the Seed, et a l . method can be summarized as follows: a) Determine pre-earthquake or steady state condition that exists i n the s o i l structure before the earthquake. b) Determine the design time h i s t o r y of base a c c e l e r a t i o n for the s i t e where the earth structure i s sit u a t e d . c) Compute the time h i s t o r y of dynamic shear stresses throughout the s o i l structure using a two-dimensional dynamic response 8 a n a l y s i s . Appropriate dynamic s t r e s s - s t r a i n r e l a t i o n s h i p and damping should be used. d) Apply these stresses to undisturbed samples of s o i l consolidated to the i n i t i a l s t a t i c stresses i n the s o i l structure to determine the s t r a i n s and re s i d u a l porewater pressures. e) Based on the porewater pressure data, determine the minimum factor of safety against t o t a l f a i l u r e by l i m i t i n g equilibrium methods a f t e r reducing the strength of elements which have developed s i g n i f i c a n t seismic porewater pressures. f) If the s o i l structure i s found to be safe against t o t a l f a i l u r e , assess the o v e r a l l deformation of the s o i l structure from the st r a i n s induced by the combined e f f e c t s of s t a t i c and dynamic loads as determined from the laboratory test data. Seed, et a l . (1973) proposed an equivalent l i n e a r e l a s t i c method to model the dynamic non-linear, h y s t e r e t i c behaviour of s o i l s . The fundamental assumption i n t h i s type of approach i s that the dynamic response of a non-linear h y s t e r e t i c material may be approximated s a t i s f a c t o r i l y by a damped, e l a s t i c model i f the properties of that model are chosen appropriately. The appropriate properties are obtained by an i t e r a t i v e process. In the dynamic f i n i t e element an a l y s i s , the s t r e s s - s t r a i n properties of the s o i l are defined i n each f i n i t e element by the Poisson's r a t i o , v, and shear s t r a i n dependent shear moduli and equivalent viscous damping r a t i o s . An average or e f f e c t i v e shear s t r a i n (usually assumed to be 65% of the maximum shear s t r a i n ) i s computed i n each f i n i t e element and shear moduli and damping r a t i o s are selected compatible with these average 9 s t r a i n s and are used i n the next i t e r a t i o n . The procedure i s repeated u n t i l no s i g n i f i c a n t changes i n moduli or damping r a t i o s are necessary. The response determined during the l a s t i t e r a t i o n i s considered to be a reasonable approximation of the non-linear response. Since the f i n a l analysis with s t r a i n compatible s o i l properties i s e l a s t i c the permanent deformation i n the s o i l structures cannot be computed by t h i s method. Deformations are estimated from the s t a t i c and dynamic stresses with aid of s t r a i n data from c y c l i c t r i a x i a l t e s t s . I t i s assumed that, when the s t a t i c and dynamic stresses i n a given f i n i t e element are simulated as c l o s e l y as possible on a sample i n a c y c l i c t r i a x i a l t e s t , the r e s u l t i n g a x i a l s t r a i n i s the s t r a i n p o t e n t i a l of the f i n i t e element i n the dam. In practice this s t r a i n i s converted to a shear s t r a i n p o t e n t i a l by multiplying by a factor (1+v). The s t r a i n p o t e n t i a l i s the s t r a i n that develops i n an unconstrained s o i l element under the s p e c i f i e d loading. Since the f i n i t e elements i n the s o i l structure are interconnected, the s t r a i n s obtained by the above procedure are not the s t r a i n s that w i l l develop i n the s o i l structure but are an i n d i c a t i o n of i t s p o t e n t i a l for s t r a i n i n g under the given seismic e x c i t a t i o n . S e r f f , et a l . (1976), have proposed a procedure f o r converting the s t r a i n potentials to a set of compatible deformations. The shear stress corresponding to the shear s t r a i n p o t e n t i a l i n a f i n i t e element i s determined from the s t r e s s - s t r a i n curve ( F i g . 2.1a). The shear stresses are converted to shear force and applied to the nodes of the f i n i t e element ( F i g . 2.1b). The deformations under these nodal forces are then determined by a s t a t i c analysis and are assumed to be s e i s m i c a l l y induced permanent deformations. This technique of computing compatible Fig. 2.1. Conversion of Shear Strain Potential to Equivalent Shear Forces. 11 deformation i s sometimes referred to as a s t r a i n harmonizing technique. One of the serious l i m i t a t i o n s of any i t e r a t i v e e l a s t i c method, used to model non-linear behaviour, i s that the solutions given by the method may not be unique (Desai, et a l . 1977). This i s because the solutions obtained i n the l a s t i t e r a t i o n may depend on the assumed s o i l properties of the f i r s t i t e r a t i o n . Equivalent l i n e a r method of analysis may overestimate the seismic response of non-linear h y s t e r e t i c materials (Finn, et a l . 1978a). The overestimation i n l i n e a r 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 f i n a l set of compatible properties i n the i t e r a t i v e equivalent l i n e a r method of a n a l y s i s . Since the analysis i s ca r r i e d out with the f i n a l set of constant s t i f f n e s s e s for the ent i r e duration of the input motion, there i s time for resonant response to bu i l d up. The s t i f f n e s s properties i n non-linear materials change constantly for every time step. When resonant response i s a function p r i m a r i l y of the method of a n a l y s i s , i t i s c a l l e d pseudo-resonance. The Seed method i s a t o t a l stress method and i t does not take into account the e f f e c t s of increasing porewater pressure on s o i l s t i f f n e s s . Since response of s o i l s i s con t r o l l e d by e f f e c t i v e stresses, the v a l i d i t y of a t o t a l stress response analysis i s questionable. Finn, et a l . (1978) compared responses predicted by t o t a l stress and e f f e c t i v e stress methods. They used two one-dimensional computer programs, SHAKE (Schnabel, et a l . 1972) and DESRA1 (Lee, et a l . 1975) for th i s purpose. 12 SHAKE i s a t o t a l stress program which models s o i l as a damped equivalent l i n e a r e l a s t i c material and DESRA1 i s an e f f e c t i v e stress program which models s o i l as a non-linear h y s t e r e t i c material. Finn, et a l , (1978) concluded that the t o t a l stress analysis tends to overestimate the dynamic response when porewater pressures exceeded about 30% of the e f f e c t i v e overburden pressure. The major d i f f i c u l t y i n p r a c t i c e with the equivalent l i n e a r method i s that a d i r e c t computation of permanent deformations i s not p o s s i b l e . The concept of s t r a i n potentials has to be used to estimate permanent deformations ( S e r f f , et a l . 1976). There are two i n c o n s i s t e n -cies i n t h i s procedure. F i r s t , the computed s t r a i n s i n the f i n a l i t e r a t i o n obtained with s t r a i n compatible s o i l properties are ignored, as not being correct, whereas the computed stresses are assumed to be representative stresses i n the ground. Knowing that stresses and s t r a i n s have a one to one r e l a t i o n s h i p for a given loading, the a r b i t r a r y decision to ignore the computed s t r a i n s i s somewhat in c o n s i s t e n t . Second, the s t r a i n s from the l a s t i t e r a t i o n are ignored, whereas the s t r a i n s computed i n intermediate i t e r a t i o n s were used to obtain compatible moduli and damping r a t i o s . This type of inconsistent assumptions make the f i n a l estimated deformations somewhat a r b i t r a r y . When porewater pressures are allowed to d i s s i p a t e i n samples subjected to c y c l i c undrained t e s t s , deformations occur. This p l a s t i c deformation i s not accounted for i n the approach proposed by Seed. Noting the l i m i t a t i o n s with the Seed method, which has been i n use for about 10 years, the state-of-the-art report on analysis of permanent deformations i n earth structures recommends t h i s topic should be the subject of active research for the next ten years (USNRC, 1982). In 13 t h i s t h e s i s , an attempt which w i l l allow d i r e c t computation of transient and permanent deformations of s o i l structures i n a consistent manner, i s presented. Procedures have been developed to model non-linear h y s t e r e t i c behaviour of s o i l , taking into account the e f f e c t of porewater pressure generation on s o i l p roperties. 14 CHAPTER 3 GENERAL GUIDELINES FOR DYNAMIC ANALYSES The performance of an earth dam under earthquake induced ground motion Is an important concern i n s e i s m i c a l l y active areas. The ground accelerations induced by the earthquake can cause large i n e r t i a forces throughout the dam. These forces which reverse i n d i r e c t i o n many times during an earthquake induce a l t e r n a t i n g stresses and s t r a i n s i n the dam. If these s t r a i n s and associated displacements are large enough, large slumping and slope i n s t a b i l i t y may r e s u l t , leading to over topping and eventual f a i l u r e of the dam. Newmark (1965), i n his pioneering work on e f f e c t s of earthquakes on dams, recommended that the performance of a s o i l structure during seismic loading should be assessed i n terms of displacement and not i n terms of the "factor of safety against f a i l u r e " along an assumed f a i l u r e surface. The allowable or s a t i s f a c t o r y displacements depend mainly on the fun c t i o n a l r o l e of the s o i l deposit or the structure founded on the deposit. For example, for c r i t i c a l 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 th i s thesis i s to present a two-dimensional dynamic response analysis of s o i l structures, taking into account a l l important factors that influence the behaviour of s o i l deposits. The analysis predicts displacements, stresses, s t r a i n s and ac c e l e r a t i o n f i e l d s etc, during and a f t e r the earthquake. 15 The determination of stresses, s t r a i n s and displacements induced i n a dam by an earthquake i s a complex a n a l y t i c a l problem and a number of s i m p l i f y i n g assumptions must be made. Foremost i s the assumption that the performance of an earth dam which i s e s s e n t i a l l y a three-dimensional structure, can be interpreted from the performance of transverse cross-sections subjected to the same seismic loading. The plane s t r a i n condition i s assumed to p r e v a i l i n the proposed two-dimensional a n a l y s i s . The main reasons for t h i s assumption are the high cost and high computer storage requirements needed for a three-dimensional a n a l y s i s . In the dynamic response analysis of continuous systems such as earth dams, non-uniform mass and s t i f f n e s s d i s t r i b u t i o n s are present. The f i n i t e element approach, which can model the v a r i a t i o n i n s t i f f n e s s and mass with extreme ease, has been adopted. In the dynamic response analysis of saturated s o i l s , a decision must be made i n i t i a l l y as to whether the analysis s h a l l be c a r r i e d out i n terms of t o t a l stress or e f f e c t i v e s t r e s s . Saturated loose cohesionless s o i l s subjected to r e p e t i t i v e loading generate r e s i d u a l porewater pressures, and i f s u f f i c i e n t drainage does not occur, reduction i n e f f e c t i v e stresses w i l l r e s u l t . Since deformations are c o n t r o l l e d by e f f e c t i v e stresses and s o i l properties such as moduli and strength are functions of e f f e c t i v e stresses, an e f f e c t i v e stress response analysis i s always preferable for those type of s o i l s . E f f e c t i v e stress response analysis i s more d i f f i c u l t to perform. It requires porewater pressure generation and d i s s i p a t i o n models and a d d i t i o n a l computations are necessary to estimate current e f f e c t i v e stresses. Studies by Finn, et a l . (1978a) on the response of l e v e l saturated sandy s i t e s to seismic e x c i t a t i o n showed that e f f e c t i v e stress analyses are not generally 16 required unless the porewater pressures are l i k e l y to exceed 30% - 40% of the e f f e c t i v e overburden pressure. In dynamic analyses , i t i s assumed that the loading imposed by seismic e x c i t a t i o n i s superimposed on long term equ i l ib r ium c o n d i t i o n s . The dynamic s o i l p r o p e r t i e s , such as strength and s t i f f n e s s , which c o n t r o l the response of s o i l s t ructures to seismic load ing , depend on i n i t i a l i n s i t u e f f e c t i v e s tress c o n d i t i o n . Furthermore i n e f f e c t i v e s tress response a n a l y s i s , the computation of current e f f e c t i v e s tresses as porewater pressure develops i s important. To do t h i s , a s t a t i c a n a l y s i s , which uses mater i a l parameters app l i cab le for both s t a t i c and subsequent dynamic ana lys i s has been developed. F i n n , Lee and Mart in (1977), presented a one-dimensional dynamic response ana lys i s taking in to account a l l important factors that a f fect the s o i l behaviour. The two-dimensional response ana lys i s proposed i n t h i s thes i s i s an extension of t h e i r one-dimensional a n a l y s i s . Therefore , a review of t h e i r method of a n a l y s i s , i n c l u d i n g s p e c i f i c a l l y how the non-l i n e a r , h y s t e r e t i c s o i l behaviour has been modelled i s presented below. The p r e d i c t i v e c a p a b i l i t y of t h e i r method has been v e r i f i e d i n the laboratory and i n the f i e l d . Some d e t a i l s on these v e r i f i c a t i o n procedures are a lso i n c l u d e d . 3.1 ONE-DIMENSIONAL RESPONSE ANALYSIS BY FINN, ET A L . (1977) In h o r i z o n t a l l y layered deposits the assumption that the shear waves propagate v e r t i c a l l y leads to a shear beam type of deformation pattern i n the depos i t . Then, only the s t r e s s - s t r a i n r e l a t i o n s h i p i n shear i s required i n the a n a l y s i s . 17 The important factors that must be considered when computing the dynamic response of s o i l s are; a) The nonlinear stress dependent s t r e s s - s t r a i n behaviour. b) The modelling of unloading-reloading. c) Contemporaneous generation and d i s s i p a t i o n of porewater pressures. d) Hysteretic and viscous damping. e) S t r a i n hardening. A l l these factors have been taken into account i n the s t r e s s -The seismic loading imposes i r r e g u l a r loading pulses which consists of loading, unloading and reloading. The s o i l s exhibit d i f f e r e n t behaviour i n each of the above phases. The r e l a t i o n s h i p between shear s t r e s s , x, and shear s t r a i n , y, f o r the i n i t i a l loading phase under ei t h e r drained or undrained loading conditions i s assumed to be hyperbolic and given by, s t r a i n r e l a t i o n s presented by Finn, et a l . (1977). 3.1.1. Shear S t r e s s - S t r a i n Relationship G x = f(y) = max y (3.1) f l + G max I Y I ) T max 18 i n w h i c h , G = maximum s h e a r m o d u l u s and = s h e a r max max s t reng th . This i n i t i a l loading or skeleton curve i s shown i n F i g . 3 .1a. The unloading - re load ing has been modelled using Masing behaviour (Masing, 1926). This impl ies that the equation for the unloading curve, i f unloading occurs from ( x r , Y r ) , i s given by, Gmax " Vr) 1 + G v-Y (3-2) max I' ' rI 2 T max which i s s imply, T - T r f ( y - Y r ) (3.3) The shape of the unloading - re loading curve i s shown i n F i g . 3 .1b. Lee (1975) proposed rules for extending the Masing concept to i r r e g u l a r l oad ing . He suggested that the unloading and re load ing curves should fo l low the skeleton loading curve i f the magnitude of the previous maximum shear s t r a i n i s exceeded. In F i g . 3.2a, the unloading curve, beyond B, becomes the extension of the i n i t i a l loading i n the negative d i r e c t i o n , i . e , BC. In the case of general loading h i s t o r y , fur ther assumptions have to be made. I f the current loading curve i n t e r s e c t s the curve described by the previous loading curve, the s tress s t r a i n curve should fo l low the previous loading curve . The above rules should apply a l so to unloading . Two t y p i c a l examples are provided i n F i g . 3 .2b;(a) i f loading along path BC i s cont inued, the loading path i s assumed to be 3.1(a). Initial Loading Curve. p i g . 3 . 1 ( b ) . Masing Stress Strain C u r v e s for Unloading and Reloading. I " (a) first unloading (b) general reloading F i g . 3.2. Hysteretic Characteristics. O 21 BCAM, (b) i f unloading along path CPB i s continued, the unloading path w i l l be ABP' . Newmark and Rosenblueth (1971) have suggested a s i m i l a r procedure. 3.1.2 Porewater Pressure Model Consider a sample of saturated sand under a v e r t i c a l e f f e c t i v e i s t r e s s , o y. During a drained c y c l i c simple shear t e s t , a cycle of shear s t r a i n , y» causes an increment i n volumetric compaction s t r a i n , A e y < j , due to grain s l i p . During an undrained shear test s t a r t i n g with the same e f f e c t i v e stress system, the c y c l i c shear s t r a i n , y, causes an increase i n porewater pressure, A U . It was shown by Martin, et a l . (1975) that for f u l l y saturated sands and assuming water to be incompressible, AU - E r A e v d (3.4) i n which E „ = one-dimensional rebound modulus of sand at an t e f f e c t i v e stress a v« Martin, et a l . (1975), also showed that under simple shear conditions the volumetric s t r a i n increment, A e V ( j , i s a function of the t o t a l accumulated volumetric s t r a i n , £ v c p a n c * t n e amplitude of the shear s t r a i n cycle, y, and i s given by, C 3 e v d Ae = CAy -C e ) + — — (3.5) vd 1 2 vd Y + C, e , ' 4 vd 22 i n which , C 2 , Cg and are volume change constants that depend on the sand type and r e l a t i v e dens i ty . An a n a l y t i c a l expression for the rebound i m o d u l u s E r , a t any e f f e c t i v e s t r e s s l e v e l a v i s g i v e n by M a r t i n , et a l . (1975), as, - do' . E = ~ — = (a*) /{m K (a ) } (3.6) r de v ' 1 r v vo ' ' v ' vr i i n which, a V Q i s the i n i t i a l value of the e f f e c t i v e s t res s and K r , m and n are experimental constants for sand. The increment i n porewater pressure, AU, during a given loading cyc le with maximum shear s t r a i n , y, may now be computed using equations (3 .4 ) , (3.5) and ( 3 . 6 ) . 3 .1 .3 M o d i f i c a t i o n of Proper t ie s for Res idual Porewater Pressure The r e s i d u a l porewater p r e s s u r e , TJ reduces G „ „ „ and r r > max Tmax* These values should be updated as r e s i d u a l porewater pressure develops. Hard in , et a l . (1972) assumed that G m a x i s independent of s tress h i s t o r y and suggested, * 1/2 G = K (a ' ) ' (3.7) max m i n which K i s a constant , depends on s o i l type and r e l a t i v e d e n s i t y . The i n i t i a l and current e f f ec t ive stress condi t ions i n a simple shear apparatus with zero i n i t i a l porewater pressure are shown i n F i g . 3 .3 . Here i t i s assumed that the r a t i o between h o r i z o n t a l and v e r t i c a l F i g . 3.3(a). Initial Effective Stress Condition in a Simple Shear Apparatus. m i«2K»«r- -u) 3 F i g . 3.3(b). Intermediate Effective Stress Condition i n a Simple Shear Apparatus. Co 24 e f f e c t i v e stresses i s a constant K Q, where K Q i s the c o e f f i c i e n t of l a t e r a l pressure at r e s t . For the i n i t i a l e f f e c t i v e stress condition, ^ max^ l s 8 i v e n b v> 1/2 1 +9 V ' (G ) = K* f V > 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 d i v i d i n g the equation (3.9) by (3.8) one obtains, ( Gmax ) - U . n _ r vo i l / 2 .) L a' J (G o-' J (3.10) max vo o T h e r e f o r e , knowing ( ^ m a x ^ o ' a ' v o a n c * ^> t n e maximum s h e a r modulus at the current e f f e c t i v e stress condition can be calculated using equation 3.10. The shear s t r e n g t h ( T m a x ^ o ^ o r t n e I n i t i a l e f f e c t i v e stress condition i s given by (Finn, et a l . 1977). f ^ 1 + K 9 o i 1 - K 1/2 ( W = {(— - —°) s i n 2 0 - ( ^ - ^ ) 2 } 0' = C a' (3.11) o 1K 2 J K 2 ' ' vo vo 25 * i n which 0' i s the angle of i n t e r n a l f r i c t i o n and C i s a constant which depends on s o i l 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 ^ T m a x ^ o ' "vo and U a r e known, t h e maximum s h e a r strength at current e f f e c t i v e stress condition can be calculated from equation (3.13). 3.1.4 Influence of S t r a i n Hardening During seismic loading of dry sand or saturated sand under drained conditions, the sand structure hardens due to grain s l i p . Finn, et a l . (1977), used f o l l o w i n g equations to modify G m a x and T • max e v d (G ) = (G ) {1 + — ~ ~ } (3.14a) max max 1 H, + H„E ,J nn n 1 2 vd and (3.14b) 26 i n w h i c h , ( G m a x ) n n a n d ( T _ a x ) n n a r e t h e m o d i f i e d maximum shear modulus and shear strength i n the n th cyc le and , » H3 and are hardening constants . The s t r e s s - s t r a i n behaviour for one-dimensional a n a l y s i s i s now completely defined by equations, 3 .1 , 3.2, 3.4, 3 .5, 3 .6, 3.10, 3.13 and 3.14. In laboratory c y c l i c simple shear tests most of the volume changes i n dry sands and the increases i n porewater pressure i n undrained saturated sands occur during the unloading por t ion of the load c y c l e . Therefore , F i n n , et a l . (1977), used modi f icat ions to the s t r e s s - s t r a i n curve to take account of s t r a i n hardening and porewater pressure only during the unloading phases of the l o a d i n g . 3.1.5 D i s s i p a t i o n of Porewater Pressure I f the saturated sand deposit can dra in during shaking there w i l l be simultaneous generation and d i s s i p a t i o n of porewater pressure . The rate of increase of porewater pressure w i l l be le s s than that i n completely undrained sand. The amount of drainage depends on the permea-b i l i t y and c o m p r e s s i b i l i t y of the sand, drainage path and durat ion of shaking. The d i s t r i b u t i o n of porewater pressure at time t i s given by, _ . _ E j _ <_z_ _J-\ E 8 £ v d ot 5z W 5 z J ot 'w (3.15) 27 i n which, k z i s the permeabi l i ty and y w t n e un i t weight of water. The term c o n t a i n i n g e y ( j r e p r e s e n t s the i n t e r n a l g e n e r a t i o n of porewater pressure ( F i n n , et a l . 1977). The s t r e s s - s t r a i n r e l a t i o n s h i p out l ined above can be very e a s i l y extended to non-uniform loading using an incremental ly e l a s t i c ana lys i s i n time domain. The dynamic response can be computed for each time step by numer ica l ly so lv ing equation (3.15) and the equation of motion, as explained by F i n n , et a l . (1977). 3.2 LABORATORY VERIFICATION OF EFFECTIVE STRESS RESPONSE ANALYSIS The bas ic assumptions made i n the formulat ion of the s t re s s -s t r a i n r e l a t i o n s h i p presented above can be broadly categorized in to two groups: Those made i n the formulat ion of porewater pressure model and those made i n model l ing load ing , unloading and re load ing . The fundamental assumption that was made i n the formulat ion of the porewater pressure model, was that the porewater pressures i n an un-drained test can be obtained from volumetr ic s t ra ins measured i n a drained tes t on a s i m i l a r sample with same h i s t o r y of shear s t r a i n l o a d i n g . This means that there should be a unique r e l a t i o n s h i p between volumetr ic s t r a i n s i n drained tests and porewater pressures i n undrained tests for a given sand at corresponding s t r a i n h i s t o r i e s . F inn (1981) reported re su l t s of an extensive laboratory program to inves t iga te th i s bas ic assumption. Volumetric s t r a in s were measured i n drained Ottawa sand samples at r e l a t i v e dens i t i e s D r = 45% and 60% when subjected to constant s t r a i n cycles i n a simple shear apparatus. Porewater pressures were a lso measured i n undrained c y c l i c tes t s at the 28 same r e l a t i v e d e n s i t i e s and i n i t i a l e f f e c t i v e confining pressures. Volumetric s t r a i n s E V ( J are shown plotted against porewater pressure r a t i o s U/o^. i n F i g . 3.4 for D„ = 45%. Each point on the curve represents a set of values of e y (j and U/a v o f o r a given number of cycles with equal c y c l i c s t r a i n amplitudes. I t can be noticed that there i s a s l i g h t difference i n the applied shear s t r a i n amplitudes. But these small deviations are not important. The data i n d i c a t e a unique r e l a t i o n s h i p between volumetric s t r a i n and porewater pressure r a t i o s . The slope of t h i s curve normalized with respect to confining pressure w i l l give the rebound modulus E_. Martin, et a l . (1975) suggested that E_ can be evaluated from the unloading curve i n an oedometer, under s t a t i c conditions. But Finn (1981) showed that the rebound modulus measured i n the oedometer i s higher than the modulus computed from the slope of the curve shown i n F i g . 3.4. He used the E„ values computed from the slope of the curve shown i n F i g . 3.4 to v e r i f y the porewater pressure model. Finn (1981) maintained that the the s t r a i n hardening e f f e c t (equation. 3.14) should not be included i n the s t r e s s - s t r a i n r e l a t i o n s h i p when pred i c t i n g the behaviour of sands under undrained conditions. This i s because net volumetric s t r a i n s do not occur during undrained conditions. If drainage i s allowed to occur, then the e f f e c t s of s t r a i n hardening should be included. The porewater pressure model coupled with the s t r e s s - s t r a i n r e l a t i o n s h i p can be employed to predict l i q u e f a c t i o n strength curves. The strength curve plots of the c y c l i c shear stress r a t i o t / o v o versus the number of cycles to cause i n i t i a l l i q u e f a c t i o n , for normally consolidated (0CR=1) and over consolidated sands, obtained a n a l y t i c a l l y o b D •> M IO 0) V-1.00 -•> o» c >nfin 0.75 -u o |c 0 5 0 -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/m 2 . Relative density * 45% vo Legend Shear strain amplitudes Drained Undrained 0.056 % 0.100 % 0.210 % 0 3 0 0 % o 0 0 5 6 % • 0.100% * 0.200% A 0 3 1 4 % ) 0 2 0.4 0.6 0 8 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 exper imenta l ly , are shown i n F i g . 3 .5 . The experimental curve was obtained from undrained constant volume c y c l i c simple shear t e s t s . The comparison between the computed and measured l i q u e f a c t i o n s trength curves i s very good. This means that the assumptions made i n the formulat ion of the non- l inear h y s t e r e t i c e f f ec t ive s t r e s s - s t r a i n r e l a t i o n -ship are v a l i d . But, l i q u e f a c t i o n res i s tance curve p r e d i c t i o n i s an extreme case. F inn (1981) used th i s e f f e c t i v e s tress model to pred ic t porewater pressure development during undrained tests when subjected to constant c y c l i c shear s tress i n a simple shear apparatus. This constant c y c l i c shear s tress loading re su l t s i n an i r r e g u l a r s t r a i n h i s t o r y as the porewater pressures develop. Further the model parameters ( i = 1,4) used were obtained from constant c y c l i c shear s t r a i n t e s t s . T y p i c a l r e s u l t s obtained from two undrained tests are shown i n F i g . 3 .6 . The agreement between measured and computed porewater pressures i s remarkably good, i n d i c a t i n g that a l l the assumptions made i n the porewater pressure model and s t r e s s - s t r a i n r e l a t i o n s h i p are reasonable. 3.3 FIELD VERIFICATION OF EFFECTIVE STRESS RESPONSE ANALYSIS A unique opportunity to inves t iga te the c a p a b i l i t y of the one-dimensional e f f e c t i v e s tress response ana lys i s was provided r e c e n t l y when data became a v a i l a b l e on the dynamic response of an a r t i f i c i a l i s l a n d i n Tokyo Bay to the Mid-Chiba earthquake of 1980. Owi Is land No.1 i s an a r t i f i c i a l i s l a n d located on the west s ide of Tokyo Bay. It was constructed with mater ia l s dredged from the nearby sea. A tes t s i t e at the south end of the i s l a n d i s instrumented to record porewater pressures and ground acce lera t ions 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/m 2 , Relative density » 4 5 - 4 7 % Analytical curve O , « , A , A Experimental data 3 6 10 30 60 100 Number of cycles for initial liquefaction, N L F i g . 3 . 5 . C y c l i c S t r e s s Rat io v s . N u m b e r of C y c l e s for I n i t i a l L i q u e f a c t i o n for V a r i o u s S a n d s . Sand type : Ottawa sand (C-109) o-v'0 « 200 kN/m 2 , Relative density = 45% T / c r ' =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, D r = 45%. 33 Porewater pressures are recorded by piezometers i n s t a l l e d at depths of 6m and 14m. The Mid-Chiba earthquake, with magnitude M = 6 . 1 , shook the Tokyo Bay area on September 25, 1980. F i n n , et a l . (1982) computed the response of Owi Is land No.1 to the Mid-Chiba earthquake using a one-dimensional e f f e c t i v e s tress response a n a l y s i s . The f i r s t 10 sees. of the recorded ground acce lera t ions are shown i n F i g . 3 .7 ( a ) . During the f i r s t 4 sees, very low acce lera t ions occurred . S i g n i f i c a n t acce le ra t ions developed between 4 and 6 sees . , and thereaf ter only low l e v e l e x c i t a t i o n was recorded. The ground acce lera t ions computed by F i n n , et a l . (1982) are shown i n F i g . 3 .7 (b ) . Except for some minor d i f ferences between 8 -10 sees, range the computed recording was very s i m i l a r to the recorded motions. The porewater pressures recorded at the 6m depth on Owi Is land No. 1 are shown i n F i g . 3 .8 (a ) . The recorded porewater pressure has two components: t rans ient and r e s i d u a l . The t rans ient porewater pressures are instantaneous response of porewater to changes i n t o t a l app l ied stresses and r e s i d u a l porewater pressures occur due to p l a s t i c volume changes. The one-dimensional response ana lys i s used by F i n n , et a l . (1982) computes the r e s i d u a l porewater pressure component and i s shown i n F i g . 3 . 8 (b ) . Comparison between recorded and computed porewater pressures i s very good. 3.4. POREWATER PRESSURE MODEL IN PRACTICE To apply the porewater pressure model i n dynamic e f f e c t i v e 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 s t r a i n and 3 constants K r, m and n to represent rebound c h a r a c t e r i s t i c s . C y c l i c simple shear apparatus has to be used to obtain these constants. A number of laboratories s t i l l do not have simple shear apparatus to do these t e s t s . Over the years a procedure has evolved from a number of p r a c t i c a l experiences by which the d i r e c t measurement of these constants can be avoided (Finn, et a l . 1982). This i s done by modifying the model parameters such that i t w i l l match with the experimental l i q u e f a c t i o n strength curve and give the r i g h t rate of porewater pressure generation. The l i q u e f a c t i o n strength curve and the rate of porewater pressure development can be experimentally obtained by doing c y c l i c t r i a x i a l tests or c y c l i c simple shear tests on f i e l d samples. A study of a number of t r i a l analyses to predict the undrained behaviour of samples i n simple shear has revealed the following: a) The shape of the l i q u e f a c t i o n resistance curve i s s e n s i t i v e to the constants C^, e s p e c i a l l y C^. b) The v a r i a t i o n of K r s h i f t s the l i q u e f a c t i o n p o t e n t i a l up or down without changing the shape appreciably. In r e a l i t y the shape of the l i q u e f a c t i o n resistance curve for a number of sands i s s i m i l a r and the values of C^ given i n the l i t e r a t u r e give the shape of t y p i c a l l i q u e f a c t i o n p o t e n t i a l curves. In pr a c t i c e a t r i a l and error procedure i s adopted to get values for the model constants. The procedure i s outlined below: 36 1) Performing a number of analyses by varying K„, and select the value for K_ such that the computed l i q u e f a c t i o n resistance curve matches the experimental l i q u e f a c t i o n resistance curve. 2) For t h i s selected K_ value, c a l c u l a t e 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 s i m i l a r , a l t e r and repeat the a n a l y s i s . It should be noted that i s the only parameter that i s used i n the c a l c u l a t i o n of incremental volumetric s t r a i n i n the f i r s t c y c l e . Therefore, estimates of C^, can be interpreted from the r e s i d u a l porewater pressure recorded i n the f i r s t c y c l e . This type of t r i a l and error procedure can be employed to obtain relevant model constants such that the corresponding porewater pressure development and l i q u e f a c t i o n resistance curves are s u f f i c i e n t l y close to the ones observed i n the laboratory. 3.5 DISCUSSION In the response analysis of h o r i z o n t a l l y layered deposits subjected to h o r i z o n t a l accelerations, a shear beam type of deformation pattern i s assumed i n the deposit. Therefore, only the s t r e s s - s t r a i n r e l a t i o n s h i p i n shear i s required. The tangent shear modulus i s used as the e l a s t i c parameter i n the incrementally e l a s t i c response analysis proposed by F i n n , et a l . (1977). To extend t h e i r model to two-dimensions, two e l a s t i c parameters are required. A d e t a i l e d 37 d e s c r i p t i o n of the extension of the one-dimensional s t r e s s - s t r a i n r e l a t i o n s h i p to two-dimensions i s discussed i n Chapter 4. It has been observed i n the laboratory that the presence of s t a t i c shear s tress a f fects the porewater pressure response of samples subjected to c y c l i c loading ( F i n n , et a l . 1978; V a i d , et a l . 1979). The porewater pressure model of F i n n , et a l . (1977) i s s t r i c t l y app l i cab le to one-dimensional depos i t s , where s t a t i c shear s tress i s ze ro . Therefore , i n extending t h e i r model to two-dimensions, the inf luence of s t a t i c shear must be accounted f o r . 38 CHAPTER 4 TWO DIMENSIONAL STATIC ANALYSIS OF SOIL STRUCTURES A s t a t i c response analysis to compute i n - s i t u e f f e c t i v e stresses i s necessary because the dynamic s o i l properties, such as strength and s t i f f n e s s , depend on i n - s i t u e f f e c t i v e streses. A number of s a t i s f a c t o r y incremental e l a s t i c methods, which model the construction sequence of the s o i l structures, are already a v a i l a b l e (Ozawa, et a l , 1973; Duncan, et a l , 1978, Byrne, et a l , 1982). The s t a t i c analysis presented i n t h i s theis i s based mainly on the methods proposed by these authors. The method proposed i n t h i s thesis uses a consistent set of material parameters i n both the s t a t i c and dynamic analyses; procedures also have been incorporated to apply c o r r e c t i o n forces during the a p p l i c a t i o n of the load increments. 4.1 STRESS-STRAIN RELATIONSHIP A number of s t r e s s - s t r a i n r e l a t i o n s have been proposed i n the computation of i n - s i t u s t a t i c stresses i n s o i l deposits. They can be divided broadly into l i n e a r , b i l i n e a r , e l a s t o - p l a s t i c , v i s c o - p l a s t i c and non-linear models. Some of these models are very complex and even for simple monotonic types of loading are expensive to use i n computational schemes. 39 Two e l a s t i c constants are required for any two-dimensional i s o t r o p i c , e l a s t i c or incrementally e l a s t i c a n a l y s i s . Tangent shear modulus G t and tangent bulk modulus B t were selected as the e l a s t i c constants. In s e l e c t i n g these e l a s t i c constants, considerations were also given to the f a c t that the s t r e s s - s t r a i n formulation proposed here has to be extended to model dynamic loading conditions. It w i l l be shown i n Chapter 5, that the s e l e c t i o n of these parameters g r e a t l y reduces the amount of computation time i n the dynamic a n a l y s i s . The s t r e s s - s t r a i n model proposed here, l i k e almost a l l other s t a t i c s t r e s s - s t r a i n models for s o i l s , can model only saturated s o i l s under f u l l y drained or undrained conditions, and dry s o i l s . The parameters selected to model the s t r e s s - s t r a i n behaviour should be based on test r e s u l t s which represent as c l o s e l y as possible the loading conditions that e x i s t i n the f i e l d . For example, i n the analysis of long term s t a b i l i t y of earth dams, one should chose parameters from drained test r e s u l t s . A d e s c r i p t i o n of the s t r e s s - s t r a i n model, and the s e l e c t i o n of relevant parameters for the model are discussed i n d e t a i l , i n t h i s chapter. 4.1.1 Reasons for Selecting G t and B_t In general, s t r a i n i n an i s o t r o p i c , homogeneous, l i n e a r e l a s t i c medium can be divided into two components: volumetric s t r a i n and d e v i a t o r i c s t r a i n . The volumetric s t r a i n i s related to mean normal stress through the bulk modulus. The d e v i a t o r i c s t r a i n i s related to d e v i a t o r i c stress through the shear modulus. These two independent material moduli can be e v a l u a t e d i n d e p e n d e n t l y by a p p l y i n g u n i f o r m changes i n AO corresponding stresses. Therefore, by sel e c t i n g tangent bulk modulus and tangent shear modulus as two independent e l a s t i c constants, better controls on stresses and stra i n s can be imposed. 4.1.2 Hyperbolic Shear S t r e s s - S t r a i n Relationship A number of researchers have used a hyperbolic s t r e s s - s t r a i n r e l a t i o n s h i p to predict the behaviour of a s o i l deposit (Konder, et a l . 1963). The hyperbolic r e l a t i o n s h i p between i and y i s given i n terms of Gmax a n d xmax a s> G y max ' T = Gmax M l max ( l + W m ) i n which, T,y = are the shear stress and shear s t r a i n G_ a x = tangent shear modulus as y+0 Tmax = u l t l m a t e shear strength 4.1.2.1 Estimation of Gmax Experimental data have shown that f o r sands and s i l t s under drained conditions, Gmax " f K » e» 0CR) <4-2) i n which, 41 o"m = current e f f e c t i v e mean normal stress e = void r a t i o OCR = over consolidation r a t i o here OCR i s defined as: _ Maximum past ma.jor p r i n c i p a l stress Current major p r i n c i p a l stress The following non-dimensional equation i s widely used f o r Gmax' Gmax * K G P a t ^ j " 2 (4-3) i n which, K Q = a non-dimensional constant for a given s o i l . Pfl = atmospheric pressure. The value of K Q depends mainly on void r a t i o or r e l a t i v e density of the s o i l , grain contact c h a r a c t e r i s t i c s such as angularity and roughness of the s o i l p a r t i c l e s etc., and also on previous loading h i s t o r y . An equation s i m i l a r to (4.3) has been proposed by Hardin, et a l . (1972) and Seed, et a l . (1970) for the computation of G M O „ for sandy s o i l s for dynamic analyses. The equation given by Hardin, et a l , (1972) includes the e f f e c t of previous stress h i s t o r y . They proposed, G = K _ P fm/P ) (OCR) (4.4) max G a ^ a ; v ' 42 i n which the exponent n depends on the p l a s t i c index of the s o i l (Hardin, et a l . 1 9 7 2 ) . Values of n are given i n Table 4 . 1 . Table 4 . 1 V a r i a t i o n of Exponent, r\ with P l a s t i c Index, PI PI% 0 0 20 0 . 1 8 40 0 . 3 0 60 0 . 4 1 80 0 . 4 8 >100 0 . 5 For normally consolidated non-plastic s o i l s under drained condition t y p i c a l values for KQ varies between 200 and 8 0 0 (Byrne, 1 9 7 9 ) . For clayey s o i l s under undrained conditions, G M A X can be related to the undrained strength S U through an equation, G M A X = K S U ( 4 . 5 ) where K i s a constant for a given c l a y . 43 4.1.2.2 Estimation of Tmax For s o i l s under drained conditions, t m a x > i s the maximum i shear stress that can be applied by keeping a x i a l stresses o x, and i a v at the i r respective values, where the x-axis i s taken h o r i z o n t a l and the y-axis v e r t i c a l for convenience. In p r i n c i p l e , this i s s i m i l a r to the e s t i m a t i o n of maximum d e v i a t o r i c s t r e s s ( o " j ) m Q v i n the a n a l y s i s presented by Kulhawy, et a l . (1969). They e s t i m a t e ( a ( j ) m a x assuming tha t the minor p r i n c i p l e s t r e s s remains constant. Let us consider a case where the i n i t i a l v e r t i c a l , h o r i z o n t a l i i and s h e a r s t r e s s e s a r e av, av, and -r r e s p e c t i v e l y ( F i g . y x xy 4.1). F i g . 4.2, shows the corresponding Mohr c i r c l e diagram and the Mohr envelope. The Mohr envelope i s defined by c' and 0 ' . The points L and M represent the i n i t i a l stress state. The a p p l i c a t i o n of shear stress w i l l increase the size of the Mohr c i r c l e and the largest Mohr c i r c l e i s the one that i s tangent to the Mohr envelope. i OA = OP radius of the largest c i r c l e keeping a y and i a x constant and i t i s given by, i i OP = r c f g x + q y ) i s i n 0 ' (4.6) "•tan0' 2 ' From t r i a n g l e ABO, A B - S a x plane = maximum shear s t r e s s e s on the h o r i z o n t a l F i g . 4.2. M o h r C i r c l e Diagram . 45 = /(OA 2 - OB 2) and so, f - * } 2 s i n 2 0 ' - H^} 2 ] 1/2 (4.7) x + max This equation reduces to the equation presented by Hardin and D r n e v i c h ( 1 9 7 2 ) i f a v and a a r e r e p l a c e d by K Q a o u vo and a, vo* The estimation of T max for a s o i l element under undrained conditions can be made based on standard f i e l d t e s t s , laboratory tests or may be based on estimation of i n s i t u e f f e c t i v e stress conditions. 4.1.2.3 Influence of Over-Consolidation Compaction i s generally used to obtain a c e r t a i n density i n dam construction. So some parts of a dam are over consolidated due to the compaction pressure. The e f f e c t of over consolidation on G__„ i s Illcl A. already shown i n the equation (4.4). Over consolidation has an influence on the value of x m a x a l s o . A t y p i c a l Mohr envelope for p l a s t i c s o i l s would look l i k e that i n F i g . 4.3. Dif f e r e n t c 1 and 0' values should be used depending on whether the s o i l i s i n NC state or OC state. 4.1.2.4 E f f e c t s 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. S lope G 2-(Loadin Slope G 4 (Unloading^ G x > G 2 Shear Strain Fig. 4*4 . Effects of Loading and Unloading, 47 i n applied load, some s o i l elements w i l l experience unloading. The modulus, G^, during unloading from a s t r a i n , y> a s shown i n F i g . 4.4, i s higher than the modulus, G 2 , corresponding to loading from the same s t r a i n l e v e l i f the mean normal stress remains a constant. Unloading and reloading can be modelled using the procedures presented i n Chapter 3. However, i f the s t r a i n ranges of i n t e r e s t are small, the differ e n c e between G 2 and i s not large. Under these circumstances changes i n modulus need not be modelled. 4.1.3 Tangent Bulk Modulus B t It i s assumed that the tangent bulk modulus (B t) i s e l a s t i c for any s o i l under drained conditions and i t i s a function of the current t mean normal stress a m only. Duncan, et a l . (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 r e l a t i v e density of the s o i l , s o i l type and previous loading h i s t o r y . In e l a s t i c (or incrementally e l a s t i c ) analysis i n i s o t r o p i c and homogeneous materials, a change i n shear stress with constant mean normal stress w i l l not r e s u l t i n any volume change. But, s o i l s under 48 constant mean-normal stress exhibit volume change when subjected to shearing s t r e s s . In the analysis presented herein, allowance has been made to include the volumetric s t r a i n s that occur due to shear stresses. A d e t a i l e d d e s c r i p t i o n of how th i s i s done i s given i n Section 4.4. It should be noted that the bulk modulus constant defined here includes only the e f f e c t of mean normal st r e s s . Therefore, for the estimation of K^, i s o t r o p i c consolidation tests performed i n t r i a x i a l t est equipment must be used. Duncan, et a l . (1978) described a procedure for determining from conventional t r i a x i a l test data i n which the mean-normal stress i s not held constant. Values of determined i n thi s manner must be considered approximate. 4.2 PHYSICAL MODELLING The domain of i n t e r e s t i s assumed to be an assembly of a f i n i t e number of elements, connected at the nodal points. The formulation of the f i n i t e element equations including the e f f e c t of porewater pressure ( C h r i s t i a n , et a l , 1970) i s presented i n Appendix I. The equations are, {P} = [ K T ] {A} + [K*] {U} (4.9) where, {p} = global column vector of incremental applied loads [ K F C ] = global tangent s t i f f n e s s matrix {A} = global column vector of incremental displacements [ K * ] = global s t i f f n e s s matrix defined i n the Appendix I 49 {u} = column vector of incremental porewater pressures i n the elements. The tangent element s t i f f n e s s matrix [ k t ] depends on two fac t o r s ; the tangent moduli and the shape functions adopted i n the f i n i t e element formulation. The shape functions give the v a r i a t i o n of displacements within an element i n terms of the nodal displacements. The simplest shape function, which assumes a l i n e a r v a r i a t i o n of displacements, gives constant s t r a i n within a tria n g u l a r element. But experience has shown that the r e s u l t s obtained from such elements do not predict stresses and s t r a i n s accurately. Therefore, q u a d r i l a t e r a l elements which have a l i n e a r s t r a i n v a r i a t i o n within an element are used. For s o i l structures such as dams, layered deposits etc., elements of a r b i t r a r y q u a d r i l a t e r a l shape are very appropriate because they are f a i r l y simple and can be used to model the geometry of these s o i l structures q u i t e a c c u r a t e l y . The element s t i f f n e s s m a t r i x , [ k t ] f o r an isoparametric q u a d r i l a t e r a l element i s given i n Appendix 1. 4 . 3 . SIMULATION OF CONSTRUCTION SEQUENCE Dams are constructed sequentially. Since the behaviour of dam m a t e r i a l s are n o n - l i n e a r and s t r e s s path dependent, a r e a l i s t i c computation of stresses and s t r a i n s requires that the construction sequence be modelled. An analysis based on single stage construction or gravity switch on, w i l l give f i n a l stresses and st r a i n s d i f f e r e n t from those calculated by following the construction sequence. 50 4.3.1 Method of Ana ly s i s F i g . 4.5 shows a schematic representat ion of the sequent ia l procedure involved i n dam c o n s t r u c t i o n . There may be p r e - e x i s t i n g elements on which subsequent layers w i l l be p laced . The cons t ruc t ion l i f t s are compacted u n t i l the required density i s obta ined . Thi s type of l ayer by layer cons t ruc t ion procedure i s c a r r i e d out u n t i l the required dimensions of the dam are obta ined. In m o d e l l i n g the c o n s t r u c t i o n sequence , the i n c r e m e n t a l s t re s se s , s t r a i n s and deformations are computed for every new layer added. This i s done by so lv ing equation (4.9) for the incremental loads caused by p lac ing a f resh l a y e r . The f i n a l s t res ses , s t r a i n s and displacements of the dam are simply the a lgebra ic sum of a l l the incremental values computed for a l l the l a y e r s . 4 .3.2 Incremental Porewater Pressure S t a t i c ana lys i s can be c a r r i e d out i n a t o t a l or e f f e c t i v e s tress mode or a combination of both . In the combination mode, some elements may be i n an e f f e c t i v e s tress mode and some may be i n a t o t a l s t ress mode. When the e f f e c t i v e s tress p r i n c i p l e i s used i n the f i n i t e element formula t ion , the porewater pressure term i s introduced as shown i n the r i g h t hand s ide of the equation (4 .9 ) . For the elements for which the t o t a l s tress mode i s assumed to be a p p l i c a b l e , the element s t i f fne s s matrix [ k t ] i s based on a t o t a l s t r e s s - s t r a i n r e l a t i o n s h i p . Furthermore, the porewater pressure , u Q , i n these elements should be set to zero . However, i f the e f f e c t i v e s tress 52 mode i s s e l e c t e d f o r some e l e m e n t s , t h e element s t i f f n e s s f o r t h e s e e l e m e n t s s h o u l d be based on an e f f e c t i v e s t r e s s - s t r a i n r e l a t i o n s h i p and a v a l u e f o r u Q i s r e q u i r e d . I t s h o u l d be remembered t h a t the r e a s o n f o r t h e s t a t i c a n a l y s i s i s t o e s t i m a t e the i n - s i t u s t a t i c c o n d i t i o n w h i c h i s a " l o n g t e r m " c o n d i t i o n . I n v i e w of t h i s , the g l o b a l p o r e w a t e r p r e s s u r e v e c t o r {Ti ] used i n the e q u a t i o n (4.9) s h o u l d c o r r e s p o n d t o t h e l o n g term v a l u e . E s t i m a t e s of u Q f o r t h e e lements can be made u s i n g a number of methods s u c h as h y d r a u l i c model t e s t s , e l e c t r i c a l a n a l o g y e t c . Measured p o r e w a t e r p r e s s u r e s i n t h e f i e l d a l s o may be u s e d . The m a t r i x {IT} , f o r m u l a t e d u s i n g element p o r e w a t e r p r e s s u r e s u Q can now be used i n e q u a t i o n (4.9) t o compute s t r e s s e s and d i s p l a c e m e n t s . 4.3.3 C o m p u t a t i o n o f I n c r e m e n t a l S t r e s s e s and S t r a i n s S t r a i n s g i v e n by the f i n i t e element a n a l y s i s a r e a measure o f changes i n shape of the e l e m e n t s from some r e f e r e n c e s t a t e . I t i s assumed t h a t the c o n d i t i o n of newly p l a c e d e lements a f t e r t h e y have s e t t l e d under t h e i r own w e i g h t i s the r e f e r e n c e s t a t e (Ozawa, et a l . 1 9 7 3 ) . The t o t a l s t r a i n s a r e o b t a i n e d by a d d i n g i n c r e m e n t a l s t r a i n s c aused by t h e c o n s t r u c t i o n l a y e r s about t h i s r e f e r e n c e s t a t e . An i n c r e m e n t a l e l a s t i c a n a l y s i s can be c a r r i e d out i n a number o f ways ( D e s a i and A b e l , 1 9 7 2 ) . The a p p r o a c h a d o p t e d h e r e i s shown s c h e m a t i c a l l y i n F i g . 4.6. E s t i m a t e s of the i n c r e m e n t s i n s t r e s s e s and s t r a i n s due t o a l o a d i n c r e m e n t a r e d e t e r m i n e d u s i n g m o d u l i v a l u e s c o r r e s p o n d i n g t o s t r e s s - s t r a i n l e v e l b e f o r e t h e l o a d i n c r e m e n t was a p p l i e d . New m o d u l i c o r r e s p o n d i n g t o the average of t h e s t r a i n s b e f o r e I n i t i a l 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 = K I + M Moduli and Stresses V. i I n i t i a l Estimates or computed Values with modified G'T and B£ Layer Loading (incremental) F i n a l 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 i s applied and the s t r a i n s computed a f t e r the increment are used to compute more correct incremental stresses and s t r a i n s for the same load increment. These incremental stresses and s t r a i n s are added to i n i t i a l stresses and s t r a i n s to obtain the i n i t i a l condition for the next load increment. Re c a l l from section 4.1.2, that the r e l a t i o n s h i p between shear stress and shear s t r a i n i s assumed to be hyperbolic. Therefore, the shear stresses are computed using the shear s t r a i n s at the end of the load increment. In doing t h i s , as pointed out by Desai and Abel (1972), equilibrium i s not n e c e s s a r i l y s a t i s f i e d . Under these circumstances equilibrium c o r r e c t i o n forces may be applied to s a t i s f y equilibrium condition. In the method adopted here, correction forces that correspond to the changes i n the shear stresses computed using the procedure outlined i n Appendix I, are applied to the next load increment. Before placing a fresh layer the stress condition i n previously placed elements (pre-existing) are known. Therefore, i n the i n i t i a l a nalysis for the load increment, moduli for the pre-existing elements are known. However, for the f r e s h l y placed elements, moduli must be based on estimated stresses i n these elements. Ozawa, et a l . (1973) suggested that the stresses can be estimated using the equations, and y ' s T = 0.5 v d s i n <* lxy is o (4.10) i n which d i s the depth of center of gravity of the element from the top surface, <=0 i s the slope of the top surface and y s i s the unit weight 55 of s o i l . For t o t a l l y submerged elements Ys should be replaced by the submerged unit weight y'« 4.4. SHEAR-VOLUME COUPLING The tangent bulk modulus Bfc defined i n section 4.1.3, re l a t e s an increment i n volumetric s t r a i n , A e ^ , to an increment i n e f f e c t i v e i mean normal s t r e s s , Aa m« But i n s o i l s volumetric s t r a i n s occur also due to changes i n shear stresses. This a d d i t i o n a l volumetric s t r a i n must be accounted for i n any r e a l i s t i c modelling of s o i l behaviour. The c h a r a c t e r i s t i c drained behaviour of i n i t i a l l y loose and dense sand samples i n a simple shear apparatus i s shown i n F i g . 4.7. I n i t i a l l y for small shear st r a i n s y, both the loose and the dense samples undergo volume reduction. But l a t e r , over a considerable range of s t r a i n , they e xhibit volume expansion ( d i l a t i o n ) . For both samples i n the d i l a t i o n range, the v a r i a t i o n of volumetric s t r a i n e v vs Y ^ s l i n e a r i n i t i a l l y and then e v approach fixed values at very high s t r a i n l e v e l s . The region of i n t e r e s t i n the e y vs Y p l o t i n t y p i c a l geotechnical problems would be the i n i t i a l compaction region and the l i n e a r d i l a t i o n region. The rate of volume change i n the l i n e a r d i l a t i o n region, i s larger for the dense sand than for the loose sand. Hansen, (1958) suggested using d i l a t i o n angle v Q to characterize the d i l a t i o n r a t e . He defined v as, sinv = o dy~ = t a n P o (4.11) 5 6 F i g . 4.7. (b). Typ ica l Plots for £ Vs T for Dense and Loose Sands . 57 where 8 Q i s the slope angle. The negative sign i s introduced since compressive volumetric s t r a i n i s considered to be p o s i t i v e . For a given type of sand the angle v Q was found to be a function of the r e l a t i v e density and confining pressure. The d i l a t i o n angle v Q increases with the r e l a t i v e density of the s o i l . This was c l e a r l y shown by Vaid, et a l . (1981). They performed drained simple shear tests with constant v e r t i c a l stress a y 0 = 200 kPa, ( F i g . 4.8) on Ottawa sand, (C-109), at various r e l a t i v e d e n s i t i e s . The d i l a t i o n angle, which i s the slope of the l i n e a r d i l a t i o n portion of the plo t e v and y i s found to increase with the r e l a t i v e density of the s o i l . The d i l a t i o n i s also a function of the mean normal stress l e v e l . This was shown by Lee (1965) who performed drained t r i a x i a l tests on dense Sacramento River sand samples of constant D r = 100%. F i g . 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 i n F i g . 4.9. F i r s t l y , dense samples at high consolidation pressures behave l i k e loose samples; secondly, f a i l u r e i n terms of maximum p r i n c i p a l stress r a t i o occurs at increasing s t r a i n l e v e l s as the consolidation pressure increases; and t h i r d l y , the d i l a t i o n angle decreases and becomes negative (compaction) with increasing consolidation pressure. F i g . 4.10 shows the v a r i a t i o n of the d i l a t i o n angle v Q with mean normal stress for a number of sands which were at an i n i t i a l r e l a t i v e density of 80 percent. It i s i n t e r e s t i n g to observe that the v a r i a t i o n of v Q l i e s within a narrow band (Robertson, 1982) and also the v a r i a t i o n i s li n e a r with logarithm of mean normal s t r e s s . Based on the experimental data presented above, the following approximation for a n a l y t i c a l 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) F i g . 4 .9 T y p i c a l Drained T r i a x i a l T e s t Results on Dense Sacramento Ri v e r Sand. (a) Pr i n c i p a l Stress Ratio V s A x i a l Strain (b) Volumetric Strain Vs A x i a l Strain ( A f t e r Lee, 1965) ~ 20 V « T J Ul -I o z < z o 10 ^'/(degree*) Dj.- 802 32 3 2 3 5 •37 3 3 3 5 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/cm3Fig. 4.10. Variation of Dilation Angle, with Mean Normal Stress for Various Sands. (After Robertson, 1982). ON o 61 dense sands: the volumetric s t r a i n due to shear stresses ( e v ) i s n e g l i g i b l y small (Varadarajan, et a l . 1980, Byrne, et a l . 1982) at low shear s t r a i n l e v e l s and above t h i s l e v e l i t varies l i n e a r l y with y. This means that the p l o t of e y vs y can be i d e a l i s e d as i n F i g . 4.11, where yQ i s the shear s t r a i n above which the volumetric s t r a i n due to shear stress i s important. It should be noted that the value of v Q should be modified for the changes i n mean normal stress according to some v a r i a t i o n such as shown i n F i g . 4.10. 4.4.1 A n a l y t i c a l Formulation There are a number of ways of modelling shear-volume coupling. One i s to modify the e l a s t i c i t y matrix I) (Appendix I) such that A e x and Ae y depend also on shear stress increment ( V e r r u i j t , 1977). But t h i s type of approach w i l l give r i s e to an unsymmetrical s t i f f n e s s matrix, which unduly complicates the computations. A simpler way i s to keep the I) matrix as i t i s and to incorporate the volume change the same way as the temperature v a r i a t i o n s are analysed i n s t r u c t u r a l mechanics (Zienkiewicz, et a l . 1967; Byrne, 1979, Byrne, et a l , 1982). This i s accomplished i n the following manner. a) The incremental shear st r a i n s i n a l l elements are computed for the increment i n load, neglecting shear-volume coupling. b) v Q can be estimated for the new mean normal s t r e s s , using F i g . 4.10, and then Ae y i s computed using equation 4.11. c) The v o l u m e t r i c s t r a i n A e v then i s s p l i t i n t o A e x and ^ f t . Here i t i s assumed that, 63 A e d = Aev = 0*5 Aev* L e t u s d e f i n e a s t r a i n x y v T v e c t o r Ae such that i t s components are the estimated "0 d i l a t i o n a l s t r a i n s . Then A e T i s g i v e n by; JAe^, ~ 0 Ae y, 0}. d) The nodal forces corresponding to this s t r a i n vector Ae^ can be computed as, / / / I*?. A lo d v (4.12) V (see Appendix I) Now these forces can be added to the applied incremental load i n a) and new s t r a i n s and stresses can be computed. For computing incremental stresses, the following equation should be used, Aa = D_(Ae - A_eQ) (4.13) where, A£ = s t r a i n vector computed for the modified applied load. e) Now steps b) ->• d) can be ca r r i e d out u n t i l convergence occurs i n stress and s t r a i n increments under the applied incremental load. 4.5 INTERFACE REPRESENTATION I t may be necessary to allow r e l a t i v e displacement to occur at the interface between two f i n i t e elements to model s l i p surfaces i n the 64 f i e l d . S l i p elements, which are sometimes referred to as elements of zero thickness, can be used to model t h i s r e l a t i v e displacement. S l i p elements can be assumed to be placed along the boundaries between the two-dimensional elements representing s o i l and s t r u c t u r a l elements or wherever i t i s anticipated that r e l a t i v e movements or separation between elements may occur. The s l i p i s assumed to occur only along t h i s d i r e c t i o n and thi s occurs when the shearing forces i n the s l i p element exceeds the shear strength at the i n t e r f a c e . Goodman, et a l . (1968) have developed a two-dimensional s l i p element with eight degrees of freedom to represent j o i n t and f a u l t behaviour i n rock mechanics problems. F i g . 4.12 shows a s l i p element with nodes I,J,K and L, i n global and element axes. The forces at any point i n a s l i p element are the shear force f„ and the normal force f„ expressed per unit area of the element. The force-displacement r e l a t i o n s h i p i s assumed to be: (4.14) n n n where K g, K n = j o i n t s t i f f n e s s per unit length i n shear and normal d i r e c t i o n s r e s p e c t i v e l y . wg, wn = Shear and normal displacement at the point of i n t e r e s t . The d e f i n i t i o n of unit j o i n t s t i f f n e s s needs c l a r i f i c a t i o n . Imagine a d i r e c t shear test being performed along an i n t e r f a c e element of unit thickness. At f i r s t , when a normal force i s 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 a s p e r i t i e s i n the j o i n t deform. A t y p i c a l plot of normal deformation at the j o i n t and the force applied per unit length i s shown i n F i g . 4.13. For a n a l y t i c a l purposes the r e l a t i o n s h i p can be approximated to a s t r a i g h t l i n e and the slope i s given by 1^. Similar tests can be performed i n the tangential d i r e c t i o n and a plot between f„ and w„ can be obtained. The slope of th i s curve w i l l give K„. Using the equation 4.14 and also assuming a l i n e a r v a r i a t i o n of displacement within the j o i n t element, a s t i f f n e s s matrix K _ can be — - o i l obtained i n l o c a l or element co-ordinates. This s t i f f n e s s matrix r e l a t e s the nodal forces and the nodal displacements. The displacement vector here i s simply, u T = {u I f V-p U j , v j , u K , v K , u L , v L} It has been shown i n Appendix II that K s n i s given by, 0 K s 0 ~ K s 0 -2KS 0 2K n 0 0 " K n 0 -2K, 2 K s 0 " 2 K s 0 " K s 0 2 K n 0 - 2 * n 0 ~ K n 2 K s 0 K s 0 sym 2 K n 0 2 K s Kn 0 2 K n n (4.15) 6 7 where, L i s the length of the element. To get the s t i f f n e s s matrix i n global co-ordinates a simple transformation i s used, where T= transformation matrix containing d i r e c t i o n cosines and i s given by; (4.16b) i n which, cosec sin<* ., . r N PQ = r o o-i (4.16c) '--sinoc cos« J o o and a = angle of i n c l i n a t i o n of the s l i p element with the h o r i z o n t a l . Even though t h i s type of formulation does not include r o t a t i o n of the element d i r e c t l y , this i s taken into account since a l l 8 d.o.f. have been considered. The displacement f i e l d v a r i a t i o n assumed within the s l i p element i s consistent with the displacement f i e l d i n an isoparametric q u a d r i l a t e r a l f i n i t e element. Furthermore, both elements have the same degrees of freedom. Thus e s t a b l i s h i n g a global s t i f f n e s s matrix can be evaluated simply t r e a t i n g the s l i p element l i k e any other element. 68 4.5.1 Method of Analysis of S l i p Elements I t was assumed tha t the t a n g e n t i a l s t r e s s - d i s p l a c e m e n t r e l a t i o n s h i p ( f g vs w g) i n the s l i p element i s e l a s t i c - p e r f e c t l y p l a s t i c and the p l a s t i c region where s l i p occurs i s defined by a simple Mohr-Coulomb type of y i e l d c r i t e r i o n . For incrementally e l a s t i c a n a l y s i s , the values of K„ and K_ should be kept constant u n t i l y i e l d occurs. A f t e r the y i e l d , K„ i s set to a very small value. This small value can be viewed as the r e s i d u a l shear s t i f f n e s s . The value of K^ i s kept at i t s o r i g i n a l value. But, i f the normal force on the element i s negative, meaning that the separation of the j o i n t occurs, then the v a l u e s K g and IC^ should be set to a small value. To investigate whether y i e l d i n g i s possible or not, the shear and normal stress i n the s l i p element should be determined. Since a l i n e a r v a r i a t i o n of displacement i s assumed within an element, the shear and normal stresses vary from point to point within the element. The average incremental values of shear stress A f g and normal stress A f n for a load increment can be estimated from equation (4.14) as, A f g = K s (Awg) and A f n = 1 ^ (Awn) (4.17) where Aw. and Aw are i n c r e m e n t a l average shear and normal 5 XI displacement i n the element r e s p e c t i v e l y . Now expressions for Aw and Awn are, 69 A w s = < A u t o p > a v e " ( A u b o t t o m > a v e = (Au K+Au L)/2 - (Au I+Au J)/2 (4.18) and Awn = ( A v t Q p ) a v e - ( A v b o t t o m ) a v e = (Av K+Av L)/2 - (A V l+Avj)/2 (4.19) From equations 4.17 to 4.19, Af„ and Af_ can now be w r i t t e n as, A f s = K s [(Au K+Au L)/2 - (Au I+Au J)/2] and Af~n = 1^ [(Av K+Av L)/2 - (Av I+Av J)/2] (4.20) Tot a l shear and normal stresses f g and f can be obtained by adding up the incremental stresses for a l l the load steps. Mohr-Coulomb f a i l u r e c r i t e r i o n g i v e s the shear s t r e n g t h f m a x i n the element at any time as, fmax = c s + f n t a n < (4.21) where c and 0 are the cohesion and f r i c t i o n angle required to define s s the f a i l u r e c r i t e r i o n . If f m a x i s greater than the absolute value of fg then the s l i p element nodes could separate and t h i s i s modelled as mentioned above, by reducing the K g to a small r e s i d u a l value. The separation of a s l i p element i s indicated by negative f n . I t should be noted that a l l c a l c u l a t i o n s for the computation of f g and f are performed i n the l o c a l axes. Since the displacements from f i n i t e element analysis are given with respect to global axes, they 70 must be transformed to get displacements i n the l o c a l axes by using the inverse of the transformation matrix T. 4.5.2 Factors that Influence J o i n t Parameters In the a n a l y t i c a l model presented above, three d i s t i n c t j o i n t parameters were introduced. (1) K„, the unit s t i f f n e s s along the element. (2) K R, the unit s t i f f n e s s across the element. (3) shear strength, f m a x defined by c g and 0 g. These parameters model the behaviour of s l i p elements adequately. The v a l u e s of K„, K_ and fmav w i l l depend on (1) the S I I 111 cL 2v surface roughness of the adjacent elements (2) shape and c h a r a c t e r i s t i c s of the a s p e r i t i e s , and (3) contact area r a t i o between the j o i n t w a l l s . D e t a i l s on how these parameters can be obtained i n the laboratory are given by Goodman, et a l . (1968). 4.6 SELECTION OF SOIL PARAMETERS The process of obtaining representative values for s o i l properties i s probably one of the d i f f i c u l t tasks i n stress a n a l y s i s . It should be emphasized that i n d i v i d u a l estimation of s o i l properties i s not important. But the s t r e s s - s t r a i n v a r i a t i o n given by the selected s o i l parameters should give the best f i t to the observed laboratory behaviour of s o i l samples i n the stress (or st r a i n ) range of i n t e r e s t . 71 4.6.1 Obtaining Shear Stress - S t r a i n Parameters Aft e r deciding what drainage condition i s l i k e l y to occur i n the f i e l d , a laboratory test can be performed simulating the drainage conditions. For example to obtain parameters for an e f f e c t i v e s t r e s s -s t r a i n r e l a t i o n s h i p , a serie s of drained simple shear tests can be performed to obtain plots % vs y for various constant mean normal stress l e v e l s . Simple shear tests are i d e a l since the mean normal stress during the test remains reasonably constant. A simple, t r i a l and error method can be employed to o b t a i n v a l u e s f o r G m o v and that f i t these curves i n the stress or s t r a i n ranges of i n t e r e s t . Then knowing the stress l e v e l s corresponding to a te s t , the best estimates f o r the parameters K Q , C' and 0 ' can be obtained e a s i l y . The e f f e c t i v e stresses i n simple shear at the beginning of the dr a i n e d l o a d i n g c o n d i t i o n s can be assumed to be o" v o and K Q t t o " v o . Then the mean n o r m a l s t r e s s a m i s g i v e n by (1+2K D) a V Q/3 and t h i s i s i n general assumed to remain constant during the t e s t . But i f conventional t r i a x i a l tests are performed on the samples, then Ojjj varies as the a x i a l load v a r i e s , and therefore, i t i s not easy to obtain these parameters. If t r i a x i a l test data only are a v a i l a b l e then again the above procedure can be c a r r i e d out by considering the shear i stress and the shear s t r a i n on the f a i l u r e plane. However, since o"m increases during shearing the shear s t r e s s - s t r a i n curve obtained by th i s i procedure may be interpreted for a constant average am- The average i i o"m can be assumed to be the mean of o m , c o r r e s p o n d i n g to the beginning and the end of the t e s t . 7 2 4.6.2 Obtaining Bulk Modulus Parameters The tangent bulk modulus B t 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 t h i s equation one gets, f K p (4.23) 1-n vm b a ' taking logarithms both sides, (l-n)log(o-' m) = log {K bP^ _ n(l-n)} + log (e f f l) (4.24) i . e . , l o g ( e v m ) = ( l - n ) l o g ( a ' m ) - log {K bp( 1 _ nh-n)} The slope of the p l o t l o g ( e v m ) vs l o g ( a m ) o b t a i n e d from drained i s o t r o p i c t r i a x i a l test r e s u l t s w i l l give (1-n), and from t h i s n can be determined. Using the value of n, and the inter c e p t on log (e.yjjj) axis, can be c a l c u l a t e d . Now knowing the bulk modulus parameters, and n, the tangent bulk modulus can be computed at any given mean normal s t r e s s . 73 CHAPTER 5 TWO-DIMENSIONAL DYNAMIC ANALYSIS 5.1 FORMULATION OF THE PROBLEM The general dynamic equilibrium equations for a l i n e a r system at any time are given by (Clough and Penzien; 1975), [M] {X} + [C] {X} + [K] {X} = {P} (5.1) i n which {x}, {£}> { x} = column v e c t o r s whose components X^, X^, and X^ g i v e the r e l a t i v e a c c e l e r a t i o n , v e l o c i t y and displacement with respect to the base motion re s p e c t i v e l y of a node, [M] = mass matrix [c] = damping matrix [K] = s t i f f n e s s matrix {p} = i n e r t i a force vector, which i s defined as, -[MJ{l}x b where, {i} i s a vector with a l l components unity and X^ i s the base a c c e l e r a t i o n . In two-dimensional problems the base a c c e l e r a t i o n may have two components: h o r i z o n t a l and v e r t i c a l . If the i t ' 1 equation i n (5.1) i s w r i t t e n f o r the h o r i z o n t a l d i r e c t i o n t h e n X^ i s the base acceleration i n the ho r i z o n t a l d i r e c t i o n , and i f the equation i s for the 74 v e r t i c a l d i r e c t i o n then X^ i s the base ac c e l e r a t i o n i n the v e r t i c a l d i r e c t i o n . 5.1.1 Incremental Eq u i l i b r i u m Equations f o r Non-Linear Systems In the analysis of non-linear systems, the material properties change with time. An incrementally e l a s t i c approach has been adopted to model the non-linear material behaviour. Incremental dynamic equilibrium equations for any time i n t e r v a l , At, can be obtained by replacing the variables i n equation (5.1) by incremental values. This leads to, [M] {AX} + [C] {AX} + [K] {AX} = {AP} (5.2) i n which, [M], [C] and [K] are the mass, damping and s t i f f n e s s matrices relevant to the time i n t e r v a l f o r which the above equations are written. It i s always assumed that the mass matrix i s 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 i s obtained by lumping the mass of a f i n i t e element equally at i t s nodes. The consistent mass matrix, i s obtained using the same i n t e r p o l a t i o n functions used i n the f i n i t e element formulation. The lumped mass matrix i s very simple to compute and i t has only diagonal terms, whereas the consistent mass matrix i s somewhat harder to compute and has off-diagonal terms. Even though the consistent mass method i s more accurate, the presence of off-diagonal terms greatly increases the computation time required to solve the dynamic 75 equilibrium equations. For the accuracy l e v e l required i n t y p i c a l geotechnical problems, the lumped mass method i s considered appropriate. In general, the damping matrix [ c ] and the s t i f f n e s s matrix [K] which are introduced i n the incremental equilibrium equation (5.2) depend on the d i s t r i b u t i o n of v e l o c i t y and displacement i n the structure. Appropriate values for the time i n t e r v a l between any time t and t+At can be determined only by an i t e r a t i o n procedure, because the v e l o c i t y and displacement at the end of a time increment depend on the i n i t i a l s t i f f n e s s and damping values. This type of i t e r a t i v e s o l u t i o n scheme for every time step i s very expensive. Therefore, i n practice tangent damping and tangent s t i f f n e s s matrices which correspond to time t are used with appropriate corrections to the r e s u l t s . It w i l l be explained l a t e r how correction forces can be introduced into the so l u t i o n scheme. The s t r e s s - s t r a i n law used to determine the tangent s t i f f n e s s matrix at time t, [K F C] i s described i n Section 5.2. The [ c ] matrix w i l l be assumed to be a constant throughout the dynamic analysis and the procedure to evaluate t h i s i s presented i n Section 5.3. With, [K] = [ K t ] t and (5.3) [C] = [C] (5.4) the dynamic incremental equilibrium equations can be re-written as, 76 [M] {AX} + [C] {AX} + [ K t ] t {AX} = {AP} (5.5) When computing response for a random loading h i s t o r y , equations (5.5) have to be solved for every time step. During the procedure [ K t ] t and {AP} are updated. The step by step i n t e g r a t i o n procedure proposed by Newmark (1959) or Wilson's 0-method (Wilson, et a l . 1973) have been adopted to integrate the equations. These procedures are described i n Appendix I I I . 5.1.2. Correction Forces The n u m e r i c a l i n t e g r a t i o n procedure i s based on three s i g n i f i c a n t assumptions: (1) the v a r i a t i o n of acc e l e r a t i o n within a time step i s assumed to vary i n some known fashion e.g. l i n e a r or constant (2) the damping and s t i f f n e s s 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 i s e n t i r e l y c o r r e c t , even though the errors may be small when the time step i s short. If errors accumulate from step to step gross errors and even s o l u t i o n i n s t a b i l i t i e s may occur. These problems can be avoided by imposing the condition of global equilibrium at each step of the a n a l y s i s . The global equilibrium equations at time t i n terms of a l l force components are, { f l i t + { f D l t + { fslt = { P l t ( 5 . 6 ) i n w h i c h {^i}t» I^D^t' a n d representing i n e r t i a , damping, and mass system at any time t and {p}t t . Since the mass matrix and constant during dynamic a n a l y s i s , by, { f g } t a r e the column v e c t o r s spring forces acting on the d i s c r e t e i s the i n e r t i a force vector at time the damping matrix were assumed to be { f j } t and {felt a r e simply given { f l i t = [M] ( X } t and ( 5 . 7 ) { f D } t = [C] {X} t ( 5 . 8 ) The spring forces {^slt» c a n ^ e computed by expressing element stresses i n terms of nodal forces, applied at the nodes of the elements. The nodal forces that correspond to the dynamic stresses i n an element {f^elem* ^ r o m Appendix I are, i n which _B i s the matrix that r e l a t e s s t r a i n to nodal displacements. In t h i s manner nodal forces for a l l the f i n i t e elements can be computed and 78 the vector sum of a l l these nodal forces w i l l give the global vector { f S } f If the solutions obtained at time t are accurate then the ri g h t and l e f t hand sides of the equation (5.6) w i l l be i d e n t i c a l . But, i n general i t w i l l not be so. The corr e c t i o n force vector { p c o r r} i s given by, ( Pcorr} = l p } t " { f l l t " { f D } t " { f s k <5-10) From equations (5.6) to (5.9), equation (5.10) can be re-written as, { Pcorr> - ^ t " M W t - [C] {*} - J /// B C O j dv a l l elements (5.11) To impose equilibrium, the corr e c t i o n force vector { P C O r r ^ c a n ^ e added to the incremental equilibrium equations formulated at time t . This i s accomplished by adding { p c o r r} to the rig h t hand side of the equation (5.5). 5.2 DYNAMIC STRESS-STRAIN RELATIONSHIP In the proposed incrementally e l a s t i c analysis i n the time domain, the tangent shear modulus G t and tangent bulk modulus B t were selected as the two " e l a s t i c parameters". Some reasons for s e l e c t i n g G t and B t for s t a t i c analysis were presented i n Chapter 4 . There i s a further very important reason for adopting these parameters i n the dynamic 79 a n a l y s i s . In dynamic an a l y s i s , the moduli have to be changed for every time step. This means that the element s t i f f n e s s matrix has to be re-formulated each time. This time consuming procedure can be s i m p l i f i e d somewhat i f G t and B t are used as the " e l a s t i c " constants. The e l a s t i c i t y matrix D_ (Appendix I) under plane s t r a i n conditions, i s given by, D = B t + 3 G t B t " 3 G t 0 = B. 1 0 1 0 B t - 3 G t B t + 3 G t 0 0 + G. 0 G 4. 3 -2 3 0 -2 3 4. 3 0 (5.12) (5.13) B t % + Gt*2 (5.14) where and Q 2 are two constant matrices. From the formulation of s t i f f n e s s matrix presented i n Appendix I, the element s t i f f n e s s matrix i s , k = J / / B D B dv V (5.15) Now s u b s t i t u t i n g for D from (5.14) the equation (5.15) can be rewritten as, 80 k. = B J / J ^ £ B_ dv + G /// B 1 L B dv (5.16) V V i . e . kfc - B t Rx + G t R2 ( 5 > 1 7 ) 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 l c t matrix may now be obtained by multi p l y i n g the constant matrices R, and R 9 by the current values of B t and G t. 5.2.1 Volume Change Behaviour With regards to volume change behaviour during dynamic loading, the s o i l deposits can be divided into two basic groups. The f i r s t group includes s o i l s which can undergo volume changes under the load increments induced by the base e x c i t a t i o n and the second group includes s o i l s which cannot. Saturated gravels and dry deposits belong to the f i r s t group of s o i l s . R e c a l l equation (4.8) which re l a t e s B,. to a', 81 nun 3 J (5.20) This equation may be used to compute B t. This means that B t has to be modified for every time step. However, i t i s known that the changes i n mean normal stresses i n the s o i l elements, due to seismic e x c i t a t i o n i s small and furthermore the volume change behaviour does not influence the response of s o i l structures s i g n i f i c a n t l y . Therefore, f o r s i m p l i c i t y , B t may be kept constant throughout the dynamic a n a l y s i s . An average B t for elements can be evaluated based on i n - s i t u mean normal stresses using equation (5.20). This i s because the load pulses during seismic loading induce stresses such that t h e i r mean values are i n i t i a l i n - s i t u i stresses. It should be noted here that even i f the changes i n a m are considered, f o r t y p i c a l values of n the changes i n B t w i l l be small. Laboratory r e s u l t s have revealed that, as p l a s t i c volume changes occur during c y c l i c loading, the s o i l samples harden leading to higher bulk modulus. This increase i n bulk modulus due to s t r a i n hardening e f f e c t can be modelled the same way as the increase i n maximum shear modulus was modelled by Finn, et a l . (1977). This was accomplished by introducing hardening constants (equation 3.14). Loose saturated, sandy s o i l s and saturated clays belong to the second group of s o i l s . In saturated s o i l s volume change can occur only by porewater drainage. Within the short duration of t y p i c a l seismic e x c i t a t i o n , the occurrence of a p p r e c i a b l e amount of drainage i s questionable i n s o i l s of low permeability. In view of t h i s , the dynamic analysis proposed here assumes that no drainage occurs during the dynamic 82 loading. In saturated gravels appreciable r e s i d u a l porewater pressure does not develop because of i t s high permeability. For the second type of s o i l s , to simulate the condition of no volume change, the bulk modulus i s set to a very high value during dynamic loading conditions. 5.2.2 Dynamic Shear S t r e s s - S t r a i n Behaviour In the formulation of a complete dynamic e f f e c t i v e s t r e s s - s t r a i n r e l a t i o n s h i p the following basic aspects should be considered; 1) s o i l behaviour under i n i t i a l loading, unloading and reloading phases. 2) r e s i d u a l porewater pressure generation and i t s e f f e c t s . 5.2.2.1 Skeleton Curve f o r Dynamic Loading Under dynamic loading conditions the r e l a t i o n s h i p between dynamic (or c y c l i c ) shear stress, % c > and dynamic shear s t r a i n , y c» I s assumed to be hyperbolic. The hyperbolic r e l a t i o n s h i p (equation 4.1) i s defined by G m a x and - c m a x . Seed et. a l , (1970) proposed that maximum shear modulus, ^max ^ o r s a n d y s o i l s i s a function of e f f e c t i v e mean normal stress only, and given by, (5.21) 83 i n which, ( K 2 ) m a x i s a constant which depends on the r e l a t i v e density for a given s o i l . Seed et. a l , suggested that for sands ( K 2 ) m a x varies between 20 and 100. Hardin et. a l , (1972) suggested that the ultimate shear strength Tmax c a n ^ e calculated using i n s i t u e f f e c t i v e stresses and s t a t i c shear strength parameters such as c' and 0' (equation 4.7). They pointed out that for the l e v e l of dynamic s t r a i n (y c<l%) induced by seismic l o a d i n g , the h y p e r b o l i c c u r v e i n terms of s t a t i c T M O V i s in ci x s a t i s f a c t o r y for dynamic loading. Unlike dynamic analyses i n one-dimension, an i n i t i a l s t a t i c shear stress, T G , i s present i n the analyses i n two-dimensions. The presence of i s causes the a v a i l a b l e shear strength to be d i f f e r e n t , depending on the d i r e c t i o n of shearing. A t y p i c a l r e l a t i o n s h i p between % c and Y c I s shown i n F i g . 5.1. The a v a i l a b l e shear strength i n the d i r e c t i o n of i n i t i a l s t a t i c shear stress i s ( T . „ - T„) and i n the mcix s opposite d i r e c t i o n i t i s ( x m a x + x g ) • In dynamic an a l y s i s , the current practice i s to neglect the influence of x g on a v a i l a b l e shear strength and to assume that the % c vs y c r e l a t i o n s h i p i s symmetrical about both xQ and Y C axes. 5.2.2.2. Unloading and Reloading A l l the basic assumptions used to model unloading and reloading phases i n Chapter 3 are also adopted for dynamic analyses. However, the equations have been modified to r e f l e c t the e f f e c t of s t a t i c shear on a v a i l a b l e shear strength. For example the equation for the unloading curve KL from K i n F i g . 5.2 i s 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 t K ; i + G max Y-Y K|/ maxl (5.22a) where t m a x i = T m a x + T s » a n d '''K a n d ^YL a r e t n e dynamic shear stress and s t r a i n corresponding to reversal point K. The equation for the reloading curve LM, i s given by, ( T " T L } = 1 + G Gmax (? - V max IY (5.22b) max2 where T M A X 2 = T M A X - -c s , and T l , y L a r e t h e d y n a m i c shear stress and s t r a i n corresponding to rever s a l point L. 5.2.3. Modelling the E f f e c t s of Residual Porewater Pressure One of the important factors i n seismic response studies of s o i l deposits comprised of saturated cohesionless materials i s the influence of re s i d u a l porewater pressure generation. The presence of r e s i d u a l porewater pressure reduces the resistance to deformation. The presence of i n i t i a l s t a t i c shear s t r e s s , T , a f f e c t s the s generation of r e s i d u a l porewater pressure i n loose saturated sands s i g n i f i c a n t l y (Finn, et a l . 1978b, Vaid, et a l . 1979). The porewater pressure model developed by Martin, et a l . (1975), does not account f o r x g. The porewater pressure model adopted here i s an extension of the model of Martin, et a l . (1975) to account for the e f f e c t s of x g . 86 5.2.3.1 The Behaviour of Samples with Ts_ The o r i g i n a l procedure that accounts for the presence of x s was presented by Seed and Lee (1969). They hypothesized that the behaviour of an element i n the f i e l d with an i n i t i a l s t a t i c stress r a t i o , x/a' = oc r and subjected to a shear stress h i s t o r y on i t s f a i l u r e plane i s i d e n t i c a l to the behaviour of a representative laboratory sample consolidated such that i t has the same i n i t i a l stress r a t i o <*r on i t s f a i l u r e plane and i s subjected to the same shear stress h i s t o r y . In the f i e l d , for t y p i c a l s o i l structures, the f a i l u r e plane f o r earthquake e x c i t a t i o n can be assumed to be the hori z o n t a l plane (Seed, et a l . 1973). Then the f i e l d s t r e s s r a t i o i s simply « r = T s / o * y 0 , i n which O y 0 i s the i n i t i a l v e r t i c a l e f f e c t i v e s t r e s s . I f simple shear equipment i s used to test representative samples then the i n i t i a l v e r t i c a l and shear load can be applied such that the proper r a t i o « r = t /o' i s present on the ho r i z o n t a l plane. However i f t r i a x i a l s y *J a p p a r a t u s i s u s e d , t h e a x i a l a[c and r a d i a l p r e s s u r e s 0*3 should be such that the i n i t i a l e f f e c t i v e stress r a t i o on a plane (45+0'/2) to the h o r i z o n t a l ( f a i l u r e plane) has the same r a t i o < r r . It can be deduced from the above hypothesis of Seed and Lee (1969), that the behaviour of h o r i z o n t a l l y layered deposits can be interpreted from i s o t r o p i c a l l y consolidated t r i a x i a l (ICT) t e s t s . A n i s o t r o p i c a l l y consolidated t r i a x i a l (ACT) tests have to be used i f Comprehensive laboratory r e s u l t s are ava i l a b l e on the behaviour of c y c l i c t r i a x i a l samples under various consolidation conditions. On the contrary, only very l i m i t e d simple shear test data are a v a i l a b l e . This i s 87 e s p e c i a l l y true when x g * 0. Therefore, i n this section only t r i a x i a l test data are used to explain the diff e r e n c e i n behaviour of samples, with and without T g on the f a i l u r e plane. A number of researchers have studied the behaviour of s o i l samples subjected to ICT tests and ACT tests (Seed, et a l , 1969; Finn, et a l . 1978b; and S e l i g , et a l . 1981). There are a number of basic differences between the behaviour of samples subjected to ICT tests and ACT t e s t s . In t y p i c a l c y c l i c t r i a x i a l tests on saturated cohesionless s o i l s the porewater pressure and a x i a l s t r a i n s develop with increasing number of d e v i a t o r i c load c y c l e s . The porewater pressure (U t) recorded at any time i s the sum of r e s i d u a l porewater (or permanent) pressure (U) and the c y c l i c (or transient) porewater pressure (U c) • The c y c l i c porewater pressure i s an instantaneous increment of porewater pressure which i s a function of current changes i n the mean normal and shear stresses and the r e s i d u a l value i s the porewater pressure due to p l a s t i c deformation. The r e s i d u a l porewater pressure may be recorded when the applied c y c l i c load i s zero. The t o t a l porewater pressure i s then given by ( F i g . 5.3), U t = U + U c (5.23a) S i m i l a r l y , the t o t a l a x i a l s t r a i n at any time also can be separated into e l a s t i c and r e s i d u a l or p l a s t i c 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 e t = e p + e c (5.23b) Results of an ICT test and an ACT test are presented i n F i g . 5.4 to i l l u s t r a t e the differences i n behaviour between these two tests ( S e l i g , et a l . 1981). Two samples of Ossterschelde sand with i n i t i a l porosity of about 41.5% were consolidated i s o t r o p i c a l l y (sample a) and a n i s o t r o p i c a l l y (sample b). The c y c l i c d e v i a t o r i c stress for sample (a) was 0.3 kg/cm2 and for sample (b) was 0.45 kg/cm 2. The a x i a l s t r a i n and porewater pressure response of these two tests are shown over 200 cycles of s t r e s s . The a x i a l s t r a i n record shown i n Fig.5.4 for the ICT test i s almost symmetrical about the X axis and the average value of e t a small value. The peak values of a x i a l s t r a i n occurs when the applied d e v i a t o r i c load i s maximum. This means that the peak values are the c y c l i c components of the a x i a l s t r a i n . The residual or permanent component of the a x i a l s t r a i n i s small. The c y c l i c component remained r e l a t i v e l y small u n t i l the r e s i d u a l porewater pressure developed to about 60% of the c onsolidation pressure and thereafter increased r a p i d l y . At the end of the test, when the applied load was zero, the a x i a l s t r a i n was found to be small. However, i n the ACT t e s t , the mean a x i a l s t r a i n increased with number of cycles of loading and at the end of the test the a x i a l s t r a i n was quite large. The r e s i d u a l or permanent s t r a i n i n t h i s test increased with number of cycles while the c y c l i c s t r a i n remained r e l a t i v e l y small. In the ICT t e s t , the porewater pressure increased more r a p i d l y than for the ACT test ( F i g . 5.4b). Aft e r about 150 cycles of loading, the porewater pressure r a t i o (U/a^) for the ICT test r a p i d l y approached 100%, while the r a t i o for ACT test approached a l i m i t i n g value. 9 0 The c y c l i c porewater pressures vary depending on the magnitude of d e v i a t o r i c s t r e s s . The d e v i a t o r i c load plays two r o l e s . F i r s t l y , the increase i n d e v i a t o r i c load r e s u l t s i n an increase i n porewater pressure. Secondly the d e v i a t o r i c load causes shear stress on the f a i l u r e plane, which may cause d i l a t i o n of the sample. This w i l l lead to a drop i n porewater pressure. Therefore, the sum of these two e f f e c t s w i l l give the net c y c l i c component of the pressure. I t i s accepted that the behaviour of saturated sands are not s i g n i f i c a n t l y affected by the c y c l i c porewater pressure, U c. Only r e s i d u a l porewater pressure a f f e c t s the behaviour of saturated sands s i g n i f i c a n t l y . It i s easy to understand the difference i n behaviour between ICT tests and ACT tests with regards to the l i m i t of porewater pressure generation i f both are compared i n e f f e c t i v e stress space, such as a q, p' pl o t , where q i s the p r i n c i p l e stress difference given by (o^ - o"3) and p' i s given by ( a | + a-j)/2. Points L and M i n F i g . 5.5 represent the i n i t i a l stress state i n the (q,p*) plot for the ICT and ACT t e s t s . As the r e s i d u a l porewater pressure increases, the e f f e c t i v e stress state at any time, when no d e v i a t o r i c load i s present i n the t r i a x i a l sample, i s given by, o{ = a f c - U (5.24a) 0-3 = a 3 c - U (5.24b) then the corresponding p' and q are given by, a! + ol al + al , 1 3 l c 3c P 2 2 " (5.25a) F i g . 5.5. q vs p Plot for ICT and A C T Tests . 92 and q = a{c - o$c (5.25b) This means that i f the e f f e c t i v e s tress paths are p l o t t e d for both tests when the appl ied d e v i a t o r i c load i s zero they w i l l be p a r a l l e l to the p' axis and w i l l move towards the f a i l u r e l i n e . If a number of load cycles of s u f f i c i e n t l y high c y c l i c s tress r a t i o are a p p l i e d , the e f f ec t ive s tress path may move very c lose to the f a i l u r e envelope. I f further load increments are a p p l i e d , the sample w i l l behave such that s tress states above OA and below OB cannot occur . This means that maximum r e s i d u a l porewater pressure values recorded w i l l be such that the e f f e c t i v e s tress path w i l l be on the f a i l u r e l i n e . Therefore , the maximum r e s i d u a l porewater pressure for any ACT test i s b 2 and for any ICT test i t i s &2 ( F i g . 5 .5 ) . Using simple geometric p r i n c i p l e s i t can be shown that (Chern, 1981; Chang, 1982), al U = b = " J 2 " {l + K - (K - l ) / s in0 * } , . max 2 2 1 c c ' (5.26) where K c i s the a n i s o t r o p i c conso l ida t ion r a t i o defined as K c = a ' l c / 0 3 c * When K c = 1 ( i . e . ICT t e s t ) U m a x i s s i m p l y g iven by c r ^ . It should be noted that the ACT test reported here does not experience r e v e r s a l of shear stresses at any time. I f the t o t a l s tress path i s below the p o s i t i v e q axis then rever sa l i n shear s t ress occurs . It i s easy to show that r ever sa l w i l l occur i f , 93 °-dy > °3c <Kc ~ 1) (5.27) Finn and Byrne (1976) have shown that stress r e v e r s a l i s required i f a condition of i n i t i a l l i q u e f a c t i o n i s to be achieved i n anisotropic t e s t s . However, whether there i s r e v e r s a l or not, the a n i s o t r o p i c a l l y consolidated samples s t r a i n progressively during c y c l i c loading i n contrast to i s o t r o p i c a l l y consolidated samples. The rate of porewater pressure development i n a given sample subjected to c y c l i c loading i n a t r i a x i a l apparatus i s mainly governed by r e l a t i v e density, consolidation pressures, and applied c y c l i c stress r a t i o , °'dy/'2a,3c* compare the porewater pressure development c h a r a c t e r i s t i c s of an ACT test and an ICT t e s t , the tests have to be performed under s i m i l a r conditions. Finn, et a l , (1978b) presented a plot of the porewater pre s s u r e r a t i o U m/cr3 c vs N /N 5 0 ( F i g . 5 . 6 ) i n which U m i s the maximum porewater pressure recorded at any time during a cycle, N i s the number load cycles and N^Q i s the number of cycles to cause the porewater pressure r a t i o of 50%. The figure c l e a r l y shows that the rate of porewater pressure development i s d i f f e r e n t i n ACT and ICT t e s t s . The i n t e r p r e t a t i o n of these test r e s u l t s should be done with care, since N^Q i s also a function of K C . In presenting laboratory r e s u l t s on dynamic properties of s a t u r a t e d c o h e s i o n l e s s s o i l samples i t i s customary to p r o v i d e l i q u e f a c t i o n p o t e n t i a l curves. These curves are a plot of c y c l i c stress r a t i o °'dy/2o'l3c versus number of cycles to l i q u e f a c t i o n , N^. As explained above, i n some ACT tests i t i s not possible to reach l i q u e f a c t i o n defined as U/o^ c = i . Therefore, the d e f i n i t i o n of uAu tubas ufao—mkm~ Number of Cycles to Liquefaction, N L 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 = U m a x o r cycles to cause some s p e c i f i e d s t r a i n (Seed, et a l . 1969). Often i n p r a c t i c e double amplitude a x i a l s t r a i n (peak to peak) of 5% i s used. I t can be seen from F i g . 5.4 that for the ICT t e s t , a single amplitude s t r a i n of 2 1 / 2 % i s equivalent to 5% double amplitude. This i s because the a x i a l s t r a i n record i s f a i r l y symmetrical about the X a x i s . However, th i s i s not true i n ACT t e s t s . Therefore, i n presenting l i q u e f a c t i o n p o t e n t i a l curves for ACT tests i t may be necessary to define as number of cycles for some maximum single amplitude s t r a i n or double amplitude or the condition where U = U m a x « T y p i c a l plots of l i q u e f a c t i o n p o t e n t i a l curves for two dense samples with and without x g are given i n F i g . 5.7 (Rahman, et a l . 1977). Two observations can be made from these f i g u r e s . F i r s t l y , the l i q u e f a c t i o n p o t e n t i a l curves are very s i m i l a r i n shape. Secondly, the curves for ACT tests are above the curves for ICT test i n d i c a t i n g higher resistance when t g i s present. However, the second conclusion may not be true for loose samples (Vaid, et a l . 1979). In presenting the differences between samples with and without T„, only t r i a x i a l r e s u l t s were used. Vaid, et a l . (1979) c a r r i e d out a s number of simple shear tests with and without % . The conclusions drawn by th e i r i n v e s t i g a t i o n and r e s u l t s reported by Seed, (Seed, 1983) are also very s i m i l a r to those presented above. 5.2.3.2 Residual Porewater Pressure Generation Model The porewater pressure model proposed by Martin, et a l . (1975) was modified to include the e f f e c t s of t g . With regards to generation of 96 porewater pressure, i t has been documented i n the preceeding section that the presence of T g , has three basic e f f e c t s . They are: the p o s i t i o n of the l i q u e f a c t i o n p o t e n t i a l curve i s a l t e r e d , there i s a l i m i t to the amount of r e s i d u a l porewater pressure, and i t s rate of generation i s d i f f e r e n t when x * 0. The attempts made to account for these e f f e c t s s i n the porewater pressure generation model are presented i n t h i s s e c t i o n . Recall the equation (3.5) which relates the incremental volumetric s t r a i n to c y c l i c shear s t r a i n amplitude, Ae , = C. ( y - C 0 e .) + . r — (5.28) vd 1 ' 2 vd Y + C / £ J ' 4 vd In one-dimensional response a n a l y s i s only shear s t r a i n components are present. But i n an analysis i n two-dimensions, v e r t i c a l and horizontal s t r a i n s also occur. However, i f i t i s assumed that only the c y c l i c component of shear s t r a i n y X y contributes to A e v d , then the equation (5.28) can s t i l l be used. This assumption i s quite reasonable since the major s t r a i n that occurs during seismic e x c i t a t i o n i n t y p i c a l s o i l structures which have adequate s t a t i c f actor of safety i s shear s t r a i n ( S e r f f , et a l . 1976), and furthermore, only c y c l i c components of shear s t r a i n are responsible for grain s l i p . Based on a number of shaking table t e s t s , Pyke, et a l . (1975), showed that a l l three components of acceleration contribute to volumetric s t r a i n . A n a l y t i c a l studies by Seed, et a l . (1975) have also revealed that under m u l t i - d i r e c t i o n a l shaking porewater pressures b u i l d up f a s t e r than under u n i d i r e c t i o n a l stress conditions, and that the shear stress r a t i o T d ^ a v o t o cause l i q u e f a c t i o n under m u l t i - d i r e c t i o n a l shaking 97 conditions i s about 10% le s s than that required under u n i d i r e c t i o n a l shaking c o n d i t i o n s . They suggested r e d u c i n g the s t r e s s r a t i o Td/'o"vo obtained from simple shear test r e s u l t s by about 10% to account for this e f f e c t . A two-dimensional dynamic analysis can account only for the hor i z o n t a l and the v e r t i c a l components of a c c e l e r a t i o n which act on the plane of the s o i l s tructure. The e f f e c t s of the t h i r d a c c e l e r a t i o n component can be accounted for by modifying either the incremental v o l u m e t r i c s t r a i n A e V £ j or the dynamic s h e a r s t r e s s r a t i o V 0 v o -The porewater pressure model of Martin, et a l . (1975) was modified such that the l i q u e f a c t i o n resistance curve and the generation rate of the porewater pressure matches the behaviour observed i n the laboratory sample with a given x s / a y 0 ' Here the l i q u e f a c t i o n resistance curve i s obtained by de f i n i n g as number of cycles required to reach the c o n d i t i o n of r e s i d u a l porewater pr e s s u r e U = U ^, where TJ „ r r max' max i s given by equation 5.26. R e c a l l the equation for E , (1-m) E = , , rn-m r mK (a ) r vo (5.29) The term K r i n equation (5.29) and constants C± through may be adjusted to model the laboratory behaviour. A proposed t r i a l and error procedure to accomplish t h i s has been outlined i n Chapter 3. It i s worthwhile noting that the l i q u e f a c t i o n p o t e n t i a l curves generated for various K r values are also very s i m i l a r i n shape to the curves obtained 98 for various T /tr^ r a t i o s i n the laboratory. This indicates that a reasonable matching of these curves i s possible. In the computation of E"r u s i n g e q u a t i o n (5.29), a y Q and 0 ^ are s u b s t i t u t e d f o r a y o and r e s p e c t i v e l y . Here O y 0 and Oy are i n i t i a l v e r t i c a l e f f e c t i v e s t r e s s and c u r r e n t v e r t i c a l e f f e c t i v e s t r e s s r e s p e c t i v e l y . o"yD i s known from the s t a t i c analysis performed before the dynamic a n a l y s i s . o"y i s obtained by computing current e f f e c t i v e stresses. R e c a l l equation (4.9), {P} = [ K j {A} + [K*] {U} (5.30) This equation can be used to obtain the response of the deposit to increases i n r e s i d u a l porewater pressures by s e t t i n g {p} = {o}. The porewater pressure matrix {u} i s formed from the r e s i d u a l porewater pressures i n the elements, as i n the case of s t a t i c a n a l y s i s . The incremental displacements, stresses and s t r a i n s given by adopting t h i s procedure i s the response of the deposit to softening .of the elements. Furthermore, these incremental st r a i n s can also be viewed as permanent components of the s t r a i n s . The current e f f e c t i v e stress system can now be used to modify G m a x and x m a x v a l u e s , and the dynamic a n a l y s i s can be continued with a hyperbolic s t r e s s - s t r a i n r e l a t i o n s h i p compatible with the current e f f e c t i v e stress system. There i s a l i m i t to the amount of r e s i d u a l porewater pressure that can be achieved i n a t r i a x i a l apparatus. This i s easy to estimate, because the ultimate stress state of a sample has to be on the Mohr 99 envelope and i n a t r i a x i a l apparatus the e f f e c t i v e stress path followed by the sample i s pr e d i c t a b l e . However, i n two-dimensional analysis the e f f e c t i v e stress path followed by an element cannot be predicted before hand, and therefore the maximum res i d u a l porewater pressure that can be developed cannot be known before the dynamic loading. The procedure used to calculate e f f e c t i v e stresses can be used to impose l i m i t s on the amount of r e s i d u a l porewater pressure. This can be accomplished by allowing porewater pressure generation to occur during dynamic loading only u n t i l the Mohr c i r c l e drawn for the current e f f e c t i v e stress system touches the Mohr envelope. 5.3 DAMPING MATRIX [ c ] Two fundamentally d i f f e r e n t damping phenomena are associated with s o i l s , namely material damping and r a d i a t i o n damping. The material damping can be viewed as a measure of energy d i s s i p a t i o n when waves t r a v e l through s o i l s . The loss of energy due to waves t r a v e l l i n g away from the region of i n t e r e s t i s known as r a d i a t i o n damping. The material damping can be divided broadly into two groups: viscous and h y s t e r e t i c damping. The energy d i s s i p a t i o n i n viscous damping depends on the v e l o c i t y of motion or s t r a i n rate and i t i s frequency dependent. Whereas h y s t e r e t i c damping involves f r i c t i o n a l loss of energy that i s l a r g e l y independent of frequency but depends on the magnitude of displacement or s t r a i n . Laboratory tests on s o i l samples have shown that most of the energy d i s s i p a t i o n i n s o i l s occurs through i n t e r n a l f r i c t i o n which i s h y s t e r e t i c . When modelling the non-linear behaviour of s o i l by an incrementally e l a s t i c approach, the e f f e c t of h y s t e r e t i c damping has 1 0 0 been included already i n the a n a l y s i s . Viscous damping i s s t i l l needed to take into account the e f f e c t of flow of water ins i d e the s o i l s t r u c t u r e . Furthermore, a small amount of viscous damping i s necessary to co n t r o l pseudo high frequency response introduced by the numerical i n t e g r a t i o n procedure. The damping due to viscous e f f e c t s can be accounted f o r through the use of Rayleigh damping. The damping matrix [c] i s given by a l i n e a r combination of [M] and [K] g i v i n g , [C] = a [M] + b [ K ] (5.31) i n which a and b are constants. The s t i f f n e s s 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 i s very small compared to hysteretic damping, the time consuming procedure of computing [c] at a l l time steps may be unnecessary. Therefore, [c] i s assumed to be a constant and can be evaluated using the tangent s t i f f n e s s matrix [K T] at time t=0. Then, [C] = a [ M ] + b [ K T ] T = 0 (5.32) This w i l l give a damping r a t i o X for the nth mode as, bo) n 2(JO 2 n (5.33) 101 where the i s the nth mode frequency. Equation ( 5 .33 ) shows that the mass-proportional component of the damping i s inversely proportional to the frequency while the s t i f f n e s s proportional component i s d i r e c t l y proportional to the frequency. Lee (1975) proposed that only s t i f f n e s s - p r o p o r t i o n a l damping should be used and suggested using, a = 0 and b = 0.005 (5 .34 ) In s u b s t i t u t i n g these values f o r a and b i n ( 5 .33 ) one gets, X n = 0.0025u> n ( 5 . 3 5 ) I t should be remembered that i n t y p i c a l s o i l structures only lower modes of v i b r a t i o n govern the response, and therefore, i t i s unnecessary to include higher mode components. The equation ( 5 . 35 ) implies that the damping r a t i o increases l i n e a r l y with the frequency. Therefore, the response due to high frequency components of the input motion w i l l be damped out s i g n i f i c a n t l y . This i s one of the advantages of using s t i f f n e s s proportional damping r a t i o s for s o i l s . From the equation (5 .35 ) the damping r a t i o at the fundamental frequency i s given by, \ L = 0 . 0 0 2 5 W l ( 5 . 36 ) Typical periods of v i b r a t i o n of s o i l structures may vary between 0 .5 to 1.5 sec. This means that the t y p i c a l damping r a t i o for the fundamental mode at the s t a r t of the dynamic loading varies between 1% - 3%. 102 The average s t i f f n e s s of a s o i l structure during dynamic loading i s less than the s t i f f n e s s at the s t a r t of the dynamic loading. Therefore, when u s i n g the damping m a t r i x [c] computed u s i n g [ K ] t = 0 , the e f f e c t i v e damping r a t i o may be higher than the range shown above. 5 . 4 BOUNDARY CONDITIONS A p p r o p r i a t e boundary c o n d i t i o n s i n terms of f o r c e s or displacements have to be s p e c i f i e d at a l l boundaries. In the dynamic analyses involving earthquake e x c i t a t i o n , two types of bottom boundary conditions are often used: 1. A f i x e d bottom boundary located at the top of a r i g i d l a y e r . In general, base rock or a s t i f f s o i l layer can be assumed to be r i g i d . 2. A bottom boundary located at the top of a s o i l layer or soft rock with constant e l a s t i c properties. This type of boundary i s generally known as transmitting boundary. For both the above boundary conditions the earthquake motion i s s p e c i f i e d at the bottom boundary. If the second type of bottom boundary conditions i s used, the domain of i n t e r e s t need not be extended down to a r i g i d l a y e r . This procedure reduces the number of degrees of freedom leading to a reduction i n computational costs. In the method of analysis presented here only the f i r s t type of boundary condition i s considered. Three types of l a t e r a l boundary conditions are commonly used i n two-dimensional dynamic problems, involving f i n i t e element procedures: 103 1) boundaries are located s u f f i c i e n t l y far away from a structure so that wave r e f l e c t i o n does not occur during analysis or i s minimized. On these boundaries, forces, displacements or a combination of forces and displacements can be s p e c i f i e d . 2) viscous boundaries are used which attempt to absorb the ra d i a t i n g waves by a series of dashpots and springs with constant or va r i a b l e properties (Lysmer and Kuhlemeyer; 1969) 3) consistent boundaries can be provided close to the foundation of structures. These boundaries attempt to reproduce the far f i e l d response i n a way consistent with the f i n i t e element formulation used to model the dynamic problem. This i s accomplished by formulating a frequency dependent boundary s t i f f n e s s matrix which can be obtained by solving the e l a s t i c wave propagation problem i n a layered half-space (Lysmer and Wass; 1972). Of the three types of l a t e r a l boundaries, the analysis with the consistent boundaries i s 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 i s s t r i c t l y applicable only to l i n e a r (or i t e r a t i v e l i n e a r ) problems and for solutions i n the frequency domain only. The l a t e r a l boundaries i n an incrementally e l a s t i c approach i n the time domain therefore should be located as far away from a structure as p r a c t i c a b l e . The usual way to model the l a t e r a l boundaries i s to allow the nodes on these boundaries to move only i n the h o r i z o n t a l d i r e c t i o n . 104 5.5 ANALYSIS OF SOIL - STRUCTURE SYSTEMS The response of a structure founded on a s o i l deposit i s a f f e c t e d by the l o c a l s o i l conditions. The peak ac c e l e r a t i o n , frequency content and the s p a t i a l d i s t r i b u t i o n of the response c h a r a c t e r i s t i c s may a l l be a f f e c t e d . By including the structure i n the f i n i t e element domain for the response an a l y s i s , the coupled seismic response of the s o i l and structure may be determined. The presence of the structure has two major e f f e c t s on the s o i l deposit. I t increases the e f f e c t i v e stresses and i t also provides a d d i t i o n a l i n e r t i a forces. Therefore, for s o i l s which exhibit non-linear stress dependent behaviour, an uncoupled analysis i n which s o i l and s t r u c t u r a l systems are uncoupled may not be a p p l i c a b l e . Uncoupled analysis cannot predict the response of buried structures where strong s o i l - s t r u c t u r e i n t e r a c t i o n occur. 5.5.1 S l i p Elements i n Dynamic Analysis Relative displacement which may occur between s o i l and structure during strong shaking can be modelled using s l i p elements. However, the s l i p element model described i n Chapter 4 has to be modified for use under dynamic loading conditions. Ideally, once the slippage stops i n a s l i p element, the top and bottom nodes of the element should have the same acce l e r a t i o n and v e l o c i t i e s (Nadim and Whitman; 1982). However, when computing accelerations and v e l o c i t i e s by numerically i n t e g r a t i n g the equation (5.5), d i f f e r e n t values w i l l be unavoidable for the top and bottom nodes. In order to overcome t h i s d i f f i c u l t y , when no slippage 105 occurs In a s l i p element, the v e l o c i t i e s and accelerations of the bottom nodes were made equal to those of the top nodes. 5.6 SOLUTION SCHEME A step by step s o l u t i o n scheme i s ca r r i e d out to obtain the dynamic response i n the time domain. A b r i e f outline of the procedure i s given below: 1) based on the c u r r e n t v a l u e s of G , T and v. max' max ' t at time t the tangent shear modulus G t i s calculated for a l l the elements using the current s t r e s s - s t r a i n curve, for eit h e r i n i t i a l l o a d i n g , u n l o a d i n g or r e l o a d i n g whichever i s appropriate. 2) the global s t i f f n e s s matrix [ K t ] t i s determined. 3) knowing the base a c c e l e r a t i o n value at (t+At), new values f o r {x}, {x}, {x} at time ( t + At) and increments, Ay and {ACT} are computed by employing any of the d i r e c t i n t e g r a t i o n methods to solve the equations (5.5) as de t a i l e d i n Appendix III 4) i f stress or s t r a i n r e v e r s a l occurs i n any element during this time i n t e r v a l , At, the dynamic analysis i s repeated f o r a shorter time i n t e r v a l . 5) using the shear s t r a i n increment, increments i n volumetric s t r a i n and then porewater pressure are computed using equation (5.28) and (5.29). 106 6) using increments i n residual porewater pressure U, as v i r t u a l loads s t a t i c analysis i s performed to determine current e f f e c t i v e stresses. These e f f e c t i v e stresses are then used to update and X ] n a x values. A computer program TARA-2 has been developed incorporating a l l these basic steps. 5.7 COMPUTATION OF POST EARTHQUAKE DEFORMATION I t i s often required to predict the displacements at various points on the s o i l structure at the end of an earthquake. This i s referred to as dynamic r e s i d u a l displacements. To compute these, an earthquake record with enough t a i l i n g zeros should be used so that the free damped v i b r a t i o n response can be included i n the a n a l y s i s . The c y c l i c and permanent components of displacement response for saturated sands and s i l t s , are assumed to occur under undrained conditions. There w i l l be a d d i t i o n a l deformations i n these s o i l s when the resi d u a l porewater pressure d i s s i p a t e s . In undrained simple shear c y c l i c tests on saturated sands and s i l t s the p o t e n t i a l volumetric s t r a i n e y < j , which occurs due to grain s l i p , i s r e f l e c t e d as r e s i d u a l porewater pressure. When the re s i d u a l porewater pressure d i s s i p a t e s , volumetric s t r a i n , e V (j occur i n the sample r e s u l t i n g i n settlement. The pr e d i c t i o n of deformations due to the d i s s i p a t i o n of r e s i d u a l porewater pressure can be accomplished as follows. 107 One method i s to treat t h i s problem as a two-dimensional consolidation problem with known i n i t i a l porewater pressures, and to solve for the deformation at di s c r e t e time i n t e r v a l s as the drainage occurs. A second method i s to compute deformations using the volumetric s t r a i n is accumulated at the end of the dynamic loading, eV(j» This can be accomplished by tr e a t i n g t h i s volumetric s t r a i n the same way the v o l u m e t r i c s t r a i n e y i n shear-volume c o u p l i n g was modelled i n Chapter 4. In the program TARA-2 the second method i s used. The f i n a l or post earthquake deformation a f t e r the r e s i d u a l porewater pressure has dissipated i s the sum of the deformation calculated ic from the m o d i f i e d e y ( j and the r e s i d u a l dynamic deformation predicted at the end of the dynamic a n a l y s i s . 108 CHAPTER 6 VERIFICATION OF THE METHOD OF ANALYSIS The v a l i d a t i o n of computational techniques requires good prototype data for comparison. The common shaking table test cannot represent the range of i n - s i t u pressures experienced by the s o i l elements i n the f i e l d . For a p p l i c a t i o n to f u l l scale design problems we require data from c e n t r i f u g a l models where s i m i l a r i t y can be achieved with the self-weight of the prototype. Here i n - s i t u stresses are d i r e c t l y scaled to those of the f u l l scale event and therefore, the stress dependent s t r e s s - s t r a i n properties of the s o i l ( e s p e c i a l l y cohesionless s o i l s ) can be matched i n model and prototype. Centrifuge modelling laws are used to deduce prototype response from model response. The p r i n c i p l e s of c e n t r i f u g a l modelling have been discussed i n d e t a i l 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. F u l l d e t a i l s of Cambridge centrifuge equipment and test procedures have been given by Schofield (1981). F i g . 6.1 shows the model i s l a n d used i n the t e s t s . The i s l a n d was a 90mm high with side slopes at 3:1 and a crest width of 200mm. The centrifuge a c c e l e r a t i o n used i n the tests was 40g. This means that the corresponding prototype i s l a n d i s of height 3.6m, with side slopes 3:1 and has a crest width of 8m. The s t r u c t u r a l loading on the LVDT 82260 90mm I 720mm • I I LVDT 82273 A C C 734 A C C 1244 • PPT 2 338 OPPT 2332 Mild s t e e l p l a t e s 3 PPT 2331 Concrete base ACC 2 5 8 ^ Container Base 750 mm 776 mm • D i r e c t i o n of p o s i t i v e a c c e l e r a t i o n measured on accelerometers LVDT Pore p r e s s u r e transducer * • Accelerometer Fig. 6.1. Submerged Island Showing Transducer Locations. o KO 110 i s l a n d was simulated by using mild s t e e l plates of various thicknesses. The i s l a n d rests on a concrete base which i n turn i s 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 i s l a n d was instrumented by 6 DJB A23 p i e z o e l e c t r i c accelerometers, 6 Druck PDCK81 pore pressure transducers and 2 LVDTs. The l o c a t i o n of these instruments are also shown i n F i g . 6.1. The test s o i l was f i n e Leighton-Buzzard sand, mostly of s i z e passing between B.S. sieve si z e No. 120 and No. 200. The r e l a t i v e density was 60 - 70%, with e m a x = 1.03 and e m i n = 0.63 (Lee, 1983). In c e n t r i f u g a l modelling, i f the model pore f l u i d i s the same as the prototype pore f l u i d , excess pore pressures would be able to d i s s i p a t e N 2 times faster i n the model than i n the prototype, whilst the earthquake would only occur N times f a s t e r . Here N i s the scale factor of the c e n t r i f u g a l a c c e l e r a t i o n given as a r a t i o of g r a v i t y . Therefore, to model the prototype porewater drainage condition, i t i s necessary to use a pore f l u i d of v i s c o s i t y N times that of water i n the model t e s t . A s p e c i a l s i l i c o n e o i l was used as the model pore f l u i d i n Lee's t e s t s . In the tests c a r r i e d out on the i s l a n d , the t h e o r e t i c a l input wave-form was intended to be 12 sinusoidal pulses with a constant period of 0.5secs. However, the actual input motion to the i s l a n d was more com p l i c a t e d due to resonances and mechanical l i n k a g e c l e a r a n c e s i n t e r f e r i n g with the input motion, e s p e c i a l l y during the i n i t i a t i o n of the base motion. The r e s u l t s of two centrifuge tests were made a v a i l a b l e to the w r i t e r . The average contact prototype pressure on the islands f or these two tests were 15kPa (Test 1) and 31kPa (Test 2) r e s p e c t i v e l y . The input I l l accelerations measured by the accelerometer mounted on top of the concrete base for Test 1 and Test 2 are shown i n F i g . 6.2(a), and ( b ) . The recorded maximum base accelerations i n Test 1 and Test 2 were O . l l g and 0.17g r e s p e c t i v e l y . The i s l a n d responded very d i f f e r e n t l y i n these two t e s t s . The comparative study c a r r i e d out to predict performance of both the tests are reported i n t h i s chapter. 6.2 COMPARATIVE STUDY The sand used i n the model test was at an average r e l a t i v e density D r = 65%. T y p i c a l properties which are consistent for medium dense sand of D r = 65% for s t a t i c and dynamic analysis are given i n Table 6.1. The l i q u e f a c t i o n resistance curve for the sand obtained by using UBC simple shear apparatus without any s t a t i c bias (t g=0) i s shown i n F i g . 6.3. This l i q u e f a c t i o n resistance curve matches the predicted l i q u e f a c t i o n resistance curve when porewater pressure model constant, K r = 0.012. As explained i n Chapter 5, the porewater pressure model constants can be selected appropriately to account for the behaviour of samples with i n i t i a l s t a t i c s t r e s s , T G . Since test data on samples with T G, are not a v a i l a b l e , i t was decided to ignore the changes i n the pore-water pressure generation rates, and to account only for the changes i n l i q u e f a c t i o n r e s i s t a n c e . The l i q u e f a c t i o n curves which correspond to non zero t s/°" v 0 were obtained by s c a l i n g the laboratory l i q u e f a c t i o n curve. The s c a l i n g was done by using a v a i l a b l e laboratory data on medium dense Ottawa sand (Vaid, et a l . 1979). The changes i n the l i q u e f a c t i o n resistance curve can be s p e c i f i e d by associating the appropriate K r ~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 S o i l Properties Properties Total unit weight kN/m3 Bulk modulus constant Bulk modulus exponent constant n Shear modulus parameter ( K 2 ) m a x Angle of i n t e r n a l f r i c t i o n E f f e c t i v e cohesion C o e f f i c i e n t K Q a, b values used to compute [C] C^-K^ Constants Rebound modulus constants m,n S t a t i c 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 v a l u e w i t h each s t a t i c shear s t r e s s r a t i o , "^/^vo* L i n e a r i n t e r p o l a t i o n may be used to get the K r value corresponding to an" intermediate -u s/a v o. For each model test two analyses were performed. One with no s l i p elements between s o i l and structure and the other with s l i p elements. The following s l i p element properties were used, K g = 6.3 x 10 5 kN/m2/m, % » 6.3 x 10 5 kN/m2/m C g = 0.0 0 g = 35° 6.2.1 Results of Test 1 The recorded a c c e l e r a t i o n time h i s t o r i e s of the model i s l a n d have very high frequency components which contain n e g l i g i b l e energy. This type of high frequency e l e c t r i c a l noise i s unavoidable i n centrifuge te s t i n g as i t may orig i n a t e due to ambient sources such as the e l e c t r i c motor d r i v i n g the centrifuge, and also due to centrifuge v i b r a t i o n s . Dean (1983), suggested i t i s necessary to f i l t e r out very high frequency components from output q u a n t i t i e s . The computed and recorded a c c e l e r a t i o n time h i s t o r i e s shown here have been smoothed once using a three point average scheme, suggested by Dean (1983). In using t h i s scheme, the current value at any point i n time i s 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. F i g . 6.4 to F i g . 6.6 show the smoothed recorded and computed acceleration time h i s t o r i e s of accelerometers A1244, A1225 and A734. During the f i r s t 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 I n 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 i n a l l accelerometers. A f t e r that the amplitude of a c c e l e r a t i o n s t e a d i l y increased as i n the case of input motion, upto 5.5 seconds and then subsided. The a c c e l e r a t i o n time h i s t o r i e s computed by TARA-2 with and without s l i p elements are very s i m i l a r and therefore, only one of them i s shown i n F i g . 6.4.(b), through F i g . 6.6(b). The frequency and magnitude of the computed ac c e l e r a t i o n response are very s i m i l a r to corresponding recorded response values. Table 6.2 shows the computed and recorded maximum acc e l e r a t i o n of a l l three accelerometers. Table 6.2 Recorded and Computed Maximum Accelerations Maximum Acceleration % 8 Instrument Location Computed by TARA-2 Recorded Without S l i p Elements With S l i p 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. F o u r p o r e w a t e r p r e s s u r e d e v e l o p m e n t p l o t s o b t a i n e d experimentally and computed by TARA-2 are presented i n Figures 6.7 (a), (b), (c) and (d) . In t h i s test very low porewater pressures were developed. During the low l e v e l shaking of the f i r s t second, the response 120 Recorded Computed With and Without S l i p 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 S l i p Elements Computed With S l i p 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 S l i p Elements i e.c Fig. 6.7(d). Recorded and Computed Porewater Pressure of PPT2342 in Test 1. 122 of the i s l a n d i s e s s e n t i a l l y e l a s t i c and porewater pressures recorded are the instantaneous response to changes i n t o t a l s t r e s s . Such porewater pressures develop from the e l a s t i c coupling of s o i l and water. Later during the period of more severe shaking, p l a s t i c volumetric s t r a i n s occur r e s u l t i n g i n the development of r e s i d u a l porewater pressures which are independent of the instantaneous states of s t r e s s . The recorded porewater p r e s s u r e s d u r i n g t h i s time have both r e s i d u a l and i n s t a n t a n e o u s components. A f t e r 6 seconds of shaking the input motion subsides over the next two seconds to zero. During t h i s time the magnitude of the instantaneous component of porewater pressure i s small. Drainage may occur during the e x c i t a t i o n depending on the drainage c h a r a c t e r i s t i c s of the sand. Since generation of porewater pressure a f t e r 6 seconds of e x c i t a t i o n i s very small, changes i n porewater pressure can be caused only by drainage during t h i s time. A close examination of recorded porewater pressures a f t e r 6 seconds of e x c i t a t i o n reveals that the porewater pressures i n a l l four locations except at the transducer P2342, which i s located at the middle of the i s l a n d are more or les s a constant. At t h i s l o c a t i o n porewater pressure increases s l i g h t l y i n d i c a t i n g movement of water towards the center of the i s l a n d . However, since these changes are small i t i s reasonable to assume that drainage i n the i s l a n d i s n e g l i g i b l e during the base e x c i t a t i o n . TARA-2 computes only r e s i d u a l porewater pressures, so there are no f l u c t u a t i o n s due to changes i n instantaneous stress l e v e l s i n the computed curves. Furthermore, since no drainage was assumed during the e x c i t a t i o n the computed r e s i d u a l porewater p r e s s u r e s i n c r e a s e co n s i s t e n t l y . However, the rate of increase i n porewater pressures during 123 low l e v e l e x c i t a t i o n which occur before 1.5 seconds and a f t e r 6 seconds i s r e l a t i v e l y small. When r i g i d j o i n t connection i s assumed between a heavy, s t i f f s t r u c t u r a l element and an adjacent s o i l element i n an a n a l y s i s , the dynamic str a i n s developed i n the adjacent s o i l element are very small due to compatibility requirements i n displacements. However, by introducing s l i p elements between the structure and s o i l , the r e l a t i v e movement which may occur between s o i l and structure during strong base e x c i t a t i o n can be accounted f o r . The r e s u l t s from TARA-2 analyses, with and without s l i p elements, ind i c a t e that computed porewater pressures are d i f f e r e n t only i n the transducers located just below the structures. The predicted porewater pressures just below the structure i n the analysis which incorporates s l i p elements, are s l i g h t l y higher than the analysis that assumes r i g i d connection between s o i l and structure. The comparison between recorded and computed porewater pressures i s good. Only four porewater pressure time h i s t o r i e s from the model tests are a v a i l a b l e . However, maximum re s i d u a l porewater pressures, which can be interpreted as the mean recorded values a f t e r the e x c i t a t i o n has subsided, are a v a i l a b l e for a l l s i x transducers. Table 6.3 shows the computed and recorded maximum re s i d u a l porewater pressures at a l l the t r a n s d u c e r l o c a t i o n s . 6.2.2 Results of Test 2 The smoothed recorded and computed a c c e l e r a t i o n time h i s t o r i e s of accelerometers A1244, A1225 and A734 are shown i n F i g s . 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 S l i p Elements With S l i p 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 i n 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 e x c i t a t i o n for Test 2, shown i n F i g . 6.2(b), has the following c h a r a c t e r i s t i c s . Low l e v e l e x c i t a t i o n occurs i n the f i r s t second and then the amplitude of acceleration increases s t e a d i l y to a constant maximum value i n the next two seconds. This maximum amplitude i s maintained for 2 seconds and then i t subsides over the next 3 seconds. The acceleration records, except the record obtained i n accelerometer A734 which i s located on top of the structure, show the v a r i a t i o n of acceleration amplitude being s i m i l a r to that of the input motion. The acce l e r a t i o n h i s t o r i e s recorded on the top of the structure dropped to very low values a f t e r 4 seconds of e x c i t a t i o n . Unlike Test 1, the response computed with s l i p elements were found to be d i f f e r e n t from those computed without s l i p elements. The response analysis also showed that when s l i p elements are used s i g n i f i c a n t amount of s l i p occurs between s o i l and st r u c t u r e . The computed accelerations without s l i p elements are lower than the values computed when s l i p elements are used i n the an a l y s i s . In the computed ac c e l e r a t i o n h i s t o r i e s shown i n F i g s . 6.8 to F i g . 6.10, only accelerations computed using s l i p elements are presented. Lines of constant accelerations have been drawn i n Fig s . 6.8 to 6.10 to aid the i n t e r p r e t a t i o n of r e s u l t s . The v a r i a t i o n of computed ac c e l e r a t i o n h i s t o r i e s are very s i m i l a r to that of the input motion. The maximum recorded accelerations which were observed to be associated with high frequencies (indicated by sharp s p i k e s ) and maximum computed a c c e l e r a t i o n s f o r a l l three accelerometers i n the i s l a n d are shown i n 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 ac c e l e r a t i o n h i s t o r i e s suggests s a t i s f a c t o r y comparison. 129 The recorded a c c e l e r a t i o n at any time has two components: the acc e l e r a t i o n component transmitted through the s o i l from the base and the acce l e r a t i o n component transmitted to s o i l through side walls and top of the centrifuge container due to the container i t s e l f v i b r a t i n g . The computed ac c e l e r a t i o n h i s t o r y accounts only for the a c c e l e r a t i o n component transmitted through the s o i l . 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 ac c e l e r a t i o n h i s t o r i e s . Table 6.4 Recorded and Computed Maximum Accelerations Maximum Ac c l e r a t i o n , % g Accelerometer No. Computed by TARA-2 Recorded Without S l i p Elements With S l i p 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 i n F i g s . 6.11 (a), (b), (c) and (d). In th i s test unlike Test 1, very high porewater pressures were developed. During the low l e v e l of e x c i t a t i o n of the f i r s t 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 i n a l l transducers except i n the transducer P2330. A close examination of the recorded porewater pressure plots a f t e r 6 seconds of e x c i t a t i o n reveals that, except for the porewater pressure transducer P2342, which i s located at the middle of the i s l a n d , s i g n i f i c a n t d i s s i p a t i o n of porewater pressure occurred. This i s because i n t h i s test very high porewater pressures were developed leading to high pressure gradients. The tranducer P2342 behaved d i f f e r e n t l y because i t i s too far from the free draining boundaries and at t h i s l o c a t i o n 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 i s l a n d . At these transducer locations the analysis with and without s l i p elements gave very s i m i l a r r e s u l t s . The porewater pressures predicted under the structure by the analysis without s l i p elements are very low. But, when s l i p elements were provided the comparison between predicted and computed porewater pressures improved. When s l i p elements were used i n the response analysis high shear s t r a i n s developed i n the elements below the structure r e s u l t i n g i n higher porewater pressures i n those elements. Table 6.5 shows the computed and recorded maximum r e s i d u a l porewater pressures at a l l the transducer l o c a t i o n s . 6.3 APPLICABILITY OF THE METHOD OF ANALYSIS The recorded r e s i d u a l porewater pressures are interpreted as the steady increase i n porewater pressures. Therefore, any high frequency noise from ambient sources on recorded values do not a f f e c t 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 S l i p Elements With S l i p 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 i n a c c e l e r a t i o n response i s made point by point, any high frequency noise a f f e c t s the comparison. Therefore, caution should be exercised when comparing maximum accelera t i o n s . When a sand sample i s subjected to undrained loading, an abrupt change i n the d i r e c t i o n of stress path occurs, as the stress path approaches the f a i l u r e l i n e (Ishihara, et a l . 1975). The points at which various stress paths change d i r e c t i o n abruptly are assumed to l i e on a s t r a i g h t l i n e , c a l l e d the phase transformation l i n e . The slope of the phase transformation l i n e i s a few degrees less than the f a i l u r e l i n e . Any c y c l i c (or monotonic) loading beyond the phase transformation l i n e may r e s u l t i n very low e f f e c t i v e stresses i n very loose sands and increased e f f e c t i v e stresses due to d i l a t i o n i n medium dense or dense sands. The hyperbolic s t r e s s - s t r a i n r e l a t i o n s h i p assumed i n t h i s thesis i s s t r i c t l y a pplicable only for the region of stress space below the phase transformation l i n e . F i g s . 6.12 (a) and (b) show the stress paths i n a q, p' plot that were followed by four elements which correspond to the locations for which porewater pressure time h i s t o r i e s are a v a i l a b l e . The stress paths reported here are for the analysis i n which s l i p elements were included. The stress paths followed by the elements i n Test 1 are well below the f a i l u r e l i n e , where as i n Test 2, two elements are on the f a i l u r e l i n e for sometime during the dynamic loading. Under these circumstances the v a l i d i t y of response computed i n Test 2 a f t e r elements have reached the f a i l u r e l i n e i s questionable. Fig. 6.12(a). Effective Stress Paths in Test 1 Fig. 6.12(b). Effective Stress Paths in Test 2. O N 137 CHAPTER 7 APPLICATION OF THE METHOD OF ANALYSIS: TANKER ISLAND RESPONSE 7.1 INTRODUCTION Man made i s l a n d s of c o h e s i o n l e s s s o i l s have been used extensively as d r i l l i n g platforms for o i l and gas exploration i n the Beaufort Sea. Recently, as exploration has moved to deep waters, more complex forms of i s l a n d construction procedures have been introduced. The caisson-retained i s l a n d (De Jong, et a l . (1978), and s t e e l tanker i s l a n d are two t y p i c a l examples. These newer type of construction procedures greatly reduce the amount of f i l l material required and also reduce some of the hazards of wave loading on exposed i s l a n d beaches. The maximum set down water depth for the caisson-retained i s l a n d and the tanker islands depth i s f i x e d , generally around 6 to 9 metres. Therefore, i n the case of deep water a underwater sand berm i s 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 p i t and pumped as a s l u r r y through a p i p e l i n e d i r e c t l y onto the l o c a t i o n of the i s l a n d . Once the sand berm i s ready, a series of caissons or tanker i s brought to the l o c a t i o n and b a l l a s t e d onto the berm, and b a c k f i l l e d with sand, gravel or water. The d r i l l i n g i s then c a r r i e d out from the upper structure. Because of the nature of the i s l a n d construction, the dumped sand i s often loose and therefore the deformation, s t a b i l i t y and 138 l i q u e f a c t i o n p o t e n t i a l of the i s l a n d berm during earthquakes are of major concern. 7.2 ANALYSIS OF A TYPICAL TANKER ISLAND F i g . 7.1, shows, schematically a tanker i s l a n d . This i s l a n d i s provided with a cover of about 2m of rock f i l l . T y p i c a l properties of rock f i l l and sand for s t a t i c analysis are given i n Table 7.1. TABLE 7.1. S t a t i c S o i l Properties Properties Rock F i l l 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 ( K 2 ) m a x 24.0 16.0 Angle of i n t e r n a l f r i c t i o n 38.0 32.0 E f f e c t i v e cohesion 0.0 0.0 C o e f f i c i e n t K Q 0.45 0.45 The tanker i s assumed to weigh 200,000 metric tons when f u l l y ballasted with plan dimensions 170m and x 60m and 21m high. In the case presented here, i t i s assumed that the hyperbolic dynamic s t r e s s - s t r a i n r e l a t i o n s h i p has equal shear strength i n both Fig. 7.1. Schematic of Tanker Island 140 d i r e c t i o n s of shearing. T y p i c a l properties used for the dynamic analysis are give i n Table 7.2 TABLE 7.2. Dynamic S o i l Properties Properties Rock F i l l 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 f i l l i s assumed to be free draining and no drainage i s assumed i n dumped sand. A very high value was assigned to B t to simulate very low compressibility imparted to the saturated sand by the water i n the pores which i s not allowed to drain. Liquefaction resistance curves are required for d i f f e r e n t s t a t i c stress r a t i o s i n the i s l a n d . These are s p e c i f i e d by associating the appropriate K r value with each s t a t i c shear stress r a t i o t s / o " v o . The values used i n the example are presented i n Table 7.3. Bulk modulus constant ft Bulk modulus parameter n Shear modulus parameter (K 2) C^ •* CK Constants TABLE 7.3. T s / a v o and K r values T /a* K s vo r 0.0 0.004 0.10 0.015 0.20 0.05 The s t a t i c shear stress r a t i o s for the example problem considered here vary between 0.0 to 0.13 and therefore, the values of K r corresponding to r a t i o s T g / o y o , above 0.15 are not necessary. The input motion used for the analyses i s the S00E a c c e l e r a t i o n component of the Imperial V a l l e y Earthquake of May 18, 1940 scaled to O.lg. The input motion was applied at the bottom boundary of the i s l a n d . Three dynamic analyses were performed: i s l a n d alone, i s l a n d plus structure with and without s l i p elements. The properties of s l i p elements were selected so that some s l i p could occur between the structure and the i s l a n d . The s l i p element properties were assumed to be, K s = 6.3 x 10 5 kN/m2/m, K N = 6.3 x 10 5, kN/m2/m C s = 0 and 0 g = 30° A complete response study of the tanker i s l a n d could be c a r r i e d out by d i s c r e t i z i n g the enti r e domain into f i n i t e elements. However, since the s t i f f n e s s of tanker wall i s very much higher than s o i l , the tanker and i t s contents would respond l i k e a r i g i d box. In view of 142 t h i s , the tanker and i t s contents were modelled as a uniform r i g i d box. The s t i f f n e s s of s t r u c t u r a l elements were selected as 10 3 of the rock f i l l elements. 7.2.1. Results f o r Tanker Island Problem One of the factors which influence the development of r e s i d u a l porewater pressure i s c y c l i c shear stress (or c y c l i c s t r a i n ) . Since the generation of r e s i d u a l porewater pressure i s possible only i n the sand, the maximum dynamic shear stresses induced i n the dumped sand along section T-T which runs through the centre of the i s l a n d , are shown i n F i g . 7.2 for a l l three cases. This fi g u r e indicates that higher dynamic shear stresses are developed when the tanker i s i n place due to the i n e r t i a forces on the tanker. The induced shear stresses i n the dumped sand when s l i p i s allowed to occur between structure and adjacent s o i l are s l i g h t l y higher than when no s l i p elements were provided. When s l i p occurs, the magnitudes of shear stress that can be transmitted to the structure i s l i m i t e d , dictated by the shear strength of the s l i p elements. Therefore, i s l a n d responses with and without s l i p occurring between the structure and the i s l a n d may be expected to d i f f e r . Despite very high c y c l i c shear stresses generated i n the sand when tanker i s i n place, the greatest porewater pressure r a t i o s are developed i n the unloaded i s l a n d ( F i g . 7.4). This i s because the v e r t i c a l over-burden pressures are very much greater when the tanker i s present, so that, the c y c l i c shear stress r a t i o s , T c v / O y 0 , which i s the most important parameter c o n t r o l l i n g the development of porewater pressure i n a given sand, are a c t u a l l y smaller ( F i g . 7.3). I t can be r e a d i l y 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 d i s t r i b u t i o n of maximum c y c l i c shear stress r a t i o s are proportional to the d i s t r i b u t i o n of r e s i d u a l porewater pressure r a t i o s . F i g . 7.4, shows the r e s i d u a l p o r e w a t e r r a t i o H/aLQ d i s t r i b u t i o n for a l l three analyses. The r e s i d u a l porewater pressure r a t i o s obtained for the unloaded i s l a n d are higher than for the loaded i s l a n d . The r e s u l t s obtained i n the analysis without the tanker can be viewed as the solutions at a section where the influence of the structure i s n e g l i g i b l e . Therefore, the d i s t r i b u t i o n of r e s i d u a l porewater pressure r a t i o s when the structure i s i n place w i l l vary from lower values at the middle to higher values as one moves away from the s t r u c t u r e . Same conclusions were drawn by Yoshimi and Tokimatsu (1977) who studied the response of a r i g i d structure subjected to base e x c i t a t i o n on a shaking ta b l e . The r e s i d u a l porewater pressure d i s t r i b u t i o n given i n the analysis with s l i p elements i s co n s i s t e n t l y higher than the analysis without s l i p elements. This i s because lower shear stresses are induced i n the l a t t e r case. The d i s t r i b u t i o n of maximum dynamic shear s t r a i n s f or section 1-1 i s shown i n F i g . 7.5. Even though the shear stresses induced i n the unloaded i s l a n d are smaller than the loaded i s l a n d , higher shear s t r a i n s have developed i n the unloaded i s l a n d . This i s because of two f a c t o r s . F i r s t l y , the i n - s i t u overburden stresses i n the unloaded i s l a n d are very much smaller and therefore the s t r e s s - s t r a i n curve for a given element i s s o f t e r . Secondly, when an e f f e c t i v e s t r e s s - s t r a i n r e l a t i o n s h i p i s used, any increase i n r e s i d u a l porewater pressure w i l l soften the s t r e s s - s t r a i n curve. The generation of higher r e s i d u a l porewater pressures and low over burden pressures have contributed to high shear s t r a i n s i n the unloaded i s l a n d . 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 hori z o n t a l displacements which occur during the earthquake for the section 1-1 are shown i n F i g . 7.6. Much smaller dynamic displacements are computed when the tanker i s i n place. When s l i p elements were provided, s l i p occurred between the s o i l and structure and the displacements are about twice the r e s u l t s obtained without s l i p elements. F i g . 7.7 a and b show the post earthquake deformations i n the X and Y d i r e c t i o n s . The post earthquake displacement i s the sum of the dynamic r e s i d u a l displacement and the displacement due to volumetric * s t r a i n component e y ^ • Two o b s e r v a t i o n s can be made from the r e s u l t s presented i n t h i s f i g u r e . F i r s t l y , the amount of the post earthquake deformations i n the X- d i r e c t i o n are proportional to the maximum dynamic deformations. Secondly, the X-component of the displacement i s of the same order of magnitude as of the Y-component of the displacement ( s e t t l e m e n t ) . The main c o n t r i b u t i o n to the X-component of the displacement comes from the dynamic r e s i d u a l displacement and for the settlement the main contribution i s from the volumetric s t r a i n component ft evd* In the dynamic response of structures the maximum induced acceleration i n the structure i s one of the main design concerns. The maximum induced accelerations given by TARA-2 i n the structure with and without s l i p elements are 0.15g and 0.17g. This means that, i f s l i p i s prevented, the acc e l e r a t i o n induced may be higher by as much as 15% of the acceleration when s l i p i s allowed. The maximum ac c e l e r a t i o n computed on top of the unloaded i s l a n d i s 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 1 4 9 One of the ways of presenting dynamic response i s to present i t i n terms of response spectra. Response spectra for displacement, v e l o c i t y and acceleration are often presented for the motions at the base of the stru c t u r e . These r e s u l t s are then used by the engineers to predict the behaviour of structures and also to compute design forces, such as base shear. F i g . 7.8 shows the response spectrum for a c c e l e r a t i o n of the motion at the berm surface for a l l three cases considered. The damping r a t i o used i n the computation was 3%. Inspection of t h i s f i g u r e suggests that for the example problem considered here, the a c c e l e r a t i o n response predicted using the response spectrum of the unloaded i s l a n d w i l l be higher for structures with a very low period. However, for the structures with a period greater than 0.5 s e c , the response predictions w i l l be s i m i l a r . The predominant motion of a tanker during e x c i t a t i o n are s l i d i n g and rocking. The r e l a t i v e importance of these two modes can be studied by comparing r e s u l t s obtained by a two-dimensional and one-dimensional response a n a l y s i s . F i g . 7.9, shows the computed d i s t r i b u t i o n of r e s i d u a l porewater pressure r a t i o u / o y o from three two-dimensional analyses which were reported e a r l i e r and also the d i s t r i b u t i o n from a one-dimensional response a n a l y s i s . The one-dimensional case considered was the i s l a n d with tanker, without any s l i p elements. The r e s u l t s c l e a r l y show that the maximum U/Oy0 of some elements may be predicted as low as 30% of those predicted by a two-dimensional response a n a l y s i s . This means a response analysis which neglects the rocking mode of v i b r a t i o n i s non-conservative. However, i t should be mentioned that the tanker considered i n t h i s example i s very t a l l (21m) and, r i g i d , and therefore the rocking mode of v i b r a t i o n may have been more important than usual. 600 1 r- r —i 500 Island Alone — 400 | t " J l 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 s i m p l i f y i n g assumptions are made when computing the response of structures founded on s o i l deposits. The procedure outlined by the National Building Code of Canada has two basic steps. The f i r s t step i s to compute the response of the s o i l deposit alone for the given e x c i t a t i o n . The second step i s to compute the response of the structure, to the base accelerations obtained i n step one. In predicting the performance of the structure, the r e s u l t s such as porewater pressure, induced s t r a i n l e v e l etc., which are obtained from step one are also considered. This means that the code i n essence suggests the s o i l - structure systems be uncoupled and analysed independently. Figures 7.2 to 7.8 c l e a r l y show that the response of structures computed using the procedures outlined i n the National Building Code of Canada may be i n e r r o r . The presence of the structure has two basic influences on the s o i l deposit. It increases the e f f e c t i v e stresses and i t also provides a d d i t i o n a l i n e r t i a forces. Therefore, for s o i l s which exhibit non-linear stress dependent behaviour, the uncoupled analysis proposed by the code may not be a p p l i c a b l e . From t h i s t y p i c a l example, three basic conclusions can be drawn. F i r s t of a l l i t r a i s e s questions about the merit of any response analysis based on uncoupled s o i l - s t r u c t u r e systems. Secondly, one-dimensional representation of the domain which neglects the rocking degrees of freedom may not be applicable to t a l l , heavy and r i g i d s t r u ctures. T h i r d l y i t demonstrates the importance of incorporating s l i p elements i n the a n a l y s i s . Because of the great weight of the caissons or tankers, and 152 t h e i r large l a t e r a l dimensions, s o i l - s t r u c t u r e i n t e r a c t i o n e f f e c t s w i l l always be important. In these type of problems a coupled analysis of the i s l a n d and structure i s required. 153 CHAPTER 8 SUMMARY AND CONCLUSIONS 8.1 SUMMARY The main purpose o f t h i s r e s e a r c h was t o d e v e l o p a two-d i m e n s i o n a l s t a t i c and s e i s m i c r e s p o n s e a n a l y s i s o f s o i l d e p o s i t s i n c l u d i n g s o i l - s t r u c t u r e i n t e r a c t i o n . The new method f o r s t a t i c and dynamic a n a l y s e s can be pe r f o r m e d i n e i t h e r e f f e c t i v e o r t o t a l s t r e s s modes or a c o m b i n a t i o n o f b o t h . Non-l i n e a r s t r e s s - s t r a i n b e h a v i o u r of s o i l was m o d e l l e d by u s i n g an i n c r e m e n t a l l y e l a s t i c a p p r o a c h i n w h i c h t a n g e n t shear modulus and t a n g e n t b u l k modulus were t a k e n as the two " e l a s t i c " p a r a m e t e r s . The m a t e r i a l r e s p o n s e i n s h e a r was assumed t o be h y p e r b o l i c w i t h M a s i n g b e h a v i o u r d u r i n g u n l o a d i n g and r e l o a d i n g . Response t o changes i n mean normal s t r e s s was assumed t o be n o n - l i n e a r , e l a s t i c and s t r e s s dependent. When a s t a t i c a n a l y s i s i s per f o r m e d i n the t o t a l s t r e s s mode, t h e s h e a r s t r e n g t h , T „ , 0 „ , G m a v , a n d t a n g e n t b u l k m o d u l u s , B t , o f an element a r e k e p t c o n s t a n t t h r o u g h o u t t h e a n a l y s i s . The tan g e n t shear modulus, G t, i s m o d i f i e d f o r c o r r e s p o n d i n g s h e a r s t r a i n s d e v e l o p e d d u r i n g the a n a l y s i s . When e f f e c t i v e s t r e s s mode i s u s e d , the p a r a m e t e r s , t m a x > ^max a n c * B t a r e c o m p u t e d f r o m t h e e f f e c t i v e s t r e s s e s . The e f f e c t o f d i l a t i o n d u r i n g s h e a r on volume change i s t a k e n i n t o a c c o u n t . I n t he s t a t i c a n a l y s i s p r o p o s e d h e r e , g r a v i t y may be s w i t c h e d on 154 at once for the completed s o i l structure or the construction sequence can be modelled by layer a n a l y s i s . The s t r e s s - s t r a i n conditions determined by the s t a t i c analysis give the i n - s i t u stress conditions before the dynamic a n a l y s i s . S l i p or contact elements have been incorporated i n the analysis to represent the i n t e r f a c e c h a r a c t e r i s t i c s between s o i l and s t r u c t u r a l elements. The properties of the s l i p element were assumed to be e l a s t i c p e r f e c t l y p l a s t i c , with f a i l u r e at the interface given by the Mohr Coulomb f a i l u r e c r i t e r i o n . In the dynamic e f f e c t i v e stress response a n a l y s i s , r e s i d u a l porewater pressures are calculated using a modification of the model proposed by M a r t i n e t . a l , (1975). The parameters, G m a x> and Tmax» a r e m°dified for the e f f e c t s of r e s i d u a l porewater pressure. The dynamic response study includes the p r e d i c t i o n of post earthquake deformations. An extensive study c a r r i e d out to v e r i f y 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 t y p i c a l tanker i s l a n d computed by t h i s method i s presented. 8.2 CONCLUSIONS The work that has been presented i n t h i s thesis leads to the following conclusions. 1. A consistent and r e l i a b l e method for computing transient and 155 permanent deformations i n two-dimensional s o i l structures i s needed. 2. A two-dimensional dynamic response analysis which takes into account the non-linear h y s t e r e t i c stress-dependent properties of s o i l s , has been developed i n terms of both t o t a l and e f f e c t i v e stresses. 3. The method has been v e r i f i e d by comparing data from centrifuged models with predictions of the method. Comparison between predicted and measured response parameters i s generally very good. 4. Allowing for s l i p to occur between s o i l and s t r u c t u r a l elements i s very important. Analyses which allow for s l i p have c o n s i s t e n t l y lead to higher displacements i n the structure and higher porewater pressures i n the s o i l deposit. 5. The method has been applied to compute seismic response of a t y p i c a l tanker i s l a n d . The re s u l t s of t h i s study suggests that i t i s important that the response of structures founded on s o i l deposits be analysed as a coupled s o i l - s t r u c t u r e systems. 6. The v a l i d i t y of a one-dimensional response analysis instead of a two-dimensional analysis f or tanker type of structures i s 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 a n a l y s i s . 156 SUGGESTIONS FOR FURTHER STUDY 1 . A d d i t i o n a l comparative studies should be c a r r i e d out so that greater confidence could be placed on the v a l i d i t y of th i s method. Comparative studies may be performed with data from centrifuge tests or f i e l d studies. 2. Sandy materials exhibit p a r t i a l s t a b i l i z a t i o n at low confining pressures due to d i l a t i o n . Therefore, i n the response evaluation near l i q u e f a c t i o n , i t i s important that the method of analysis Include the d i l a t a n t behaviour of sands. 3. In the response evaluation of more permeable s o i l s , drainage during the seismic loading may be s i g n i f i c a n t and procedures should be developed to take t h i s into account. 157 REFERENCES 1. Byrne, P.M., (1979), Class Notes: Numerical Methods i n S o i l Mechanics (CE573). University of B r i t i s h Columbia, Vancouver, B.C., Canada. 2. Byrne, P.M., and Eldridge, T.L., (1982), "A Three Parameter D i l a t a n t E l a s t i c S t r e s s - S t r a i n Model for Sand", International Symposium on Numerical Models i n Geomechanics, Zurich, Sept., pp73-79. 3. Chang, C.S., (1982), "Residual Undrained Deformation from C y c l i c Loading", Journal of the Geotechnical Engineering, ASCE, V o l . 108, GT4, A p r i l , pp637-646. 4. Chern, J . C , (1981), " E f f e c t of S t a t i c Shear on Resistance to Liq u e f a c t i o n " , M.A.Sc. Thesis, The U n i v e r s i t y of B r i t i s h Columbia, Vancouver, May. 5. C h r i s t i a n , J.T. and Boehmer, J.W., (1970), "Plane S t r a i n Consolidation by F i n i t e Elements", Journal of the S o i l Mechanics and Foundation Engineering D i v i s i o n , ASCE, V o l . 96, No. SM4, July, ppl435-1455. 6. Clough, R.W., and Penzien, J . , (1975), "Dynamics of Structures", McGraw H i l l Book Co., New York. 7. De Jong, J.J.A., and Bruce, J . C , (1978), "Design and Construction of a Caisson-Retained-Island D r i l l i n g 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: i n F l i g h t 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 F i n i t e Element Method", L i t t o n Educational Publishing, Inc. 10. Desai, C.S., and C h r i s t i a n , J.T. (1977), "Numerical Methods i n Geomechanics", McGraw H i l l Book Co., New York. 11. Duncan, J.M., Byrne, P.M., Wong, K.S., and Mabry, P., (1978), "Strength, S t r e s s - S t r a i n and Bulk Modulus Parameters f o r F i n i t e Element Analyses of Stresses and Movements i n S o i l Masses", Report No. WCB/GT/78-02, to the National Science Foundation, A p r i l . 12. Finn, W.D.L. and M i l l e r , 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 i n Dry Sands During Earthquakes", Canadian Geotechnical Journal, V o l . 13, No. 4, pp355-363. 14. Finn, W.D.L., Lee, K.W., and Martin, G.R., (1977), "An E f f e c t i v e S t r e s s Model f o r L i q u e f a c t i o n " , J o u r n a l of the G e o t e c h n i c a l Engineering D i v i s i o n , ASCE, V o l . 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 D i v i s i o n , S p ecialty Conference on Earthquake Engineering and S o i l Dynamics, Pasadena, C a l i f o r n i a , June 19-21, pp472-491. 16. Finn, W.D.L., Lee, K.W., Maartman, C.H., and Lo, R., (1978b), " C y c l i c Pore Pressures Under Anisotropic Conditions", Proc. ASCE Geotechnical Engineering D i v i s i o n , S p ecialty Conference on Earthquake Engineering and S o i l Dynamics, Pasadena C a l i f o r n i a , June 19-21, pp457-468. 159 17. Finn, W.D.L., (1981), "Liquefaction P o t e n t i a l Development Since 1976", Proceedings, International Conference on Recent Advances i n Geotechnical Earthquake Engineering and S o i l Dynamics, St. Louis, Missouri, A p r i l 26-May 2, pp655-681. 18. Finn, W.D.L., I a i , S., and Ishihara, K., (1982), "Performance of A r t i f i c i a l Offshore Islands under Wave and Earthquake Loading", Offshore Technology Conference, V o l . 1, OTC, Paper 4220, Houston, Texas, May, pp661-672. 19. Finn, W.D.L., (1983), "Analyses of Cumulative Deformations Under C y c l i c Loading", Invited Challenge Contributed to National Science Foundation Workshop on Expected Research i n S o i l Engineering, Blackberg, V i r g i n i a , U.S.A., Aug. 22-24. 20. Goodman, R.E., and Seed, H.B., (1966), "Earthquake Induced Displacements i n Sand Embankments", Journal of the S o i l Mechanics and Foundations D i v i s i o n , ASCE, V o l . 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 S o i l Mechanics and Foundations D i v i s i o n , 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, V o l . 1, pp. 39-48. 23. Hardin, B.O., and Drenevich, V.P., (1972), "Shear Modulus and Damping i n S o i l s , Design Equations and Curves", Journal of the S o i l Mechanics and Foundation D i v i s i o n , ASCE, V o l . 98, SM7, Proc. Paper 9006, July, pp667-692. 160 24. I a i , S., and Finn. W.D.L., (1982), "DONAL-2, A Computer Program for Dynamic One Dimensional Analyses of Slope Layers with Energy Transmitting Boundary", S o i l Dynamics Group, U n i v e r s i t y of B r i t i s h Columbia, Vancouver, B.C., Feb. 25. Ishihara, K., Tatsuoka, F. and Yasuda, S., (1975), "Undrained Deformation and Liquefaction of Sand Under C y c l i c Stresses", S o i l s and Foundations, Vol. 15, No. 1, pp29-44. 26. Konder, R.L., and Zelasko, J.S.,, (1963), "A Hyperbolic S t r e s s - S t r a i n Formulation of Sands", Proceedings of 2nd Pan American Conference on S o i l Mechanics and Foundation Engineering, V o l . 1, B r a z i l , 1963, pp289. 27. Kulhawy, F.H., Duncan, J.M., and Seed, H.B., (1969), " F i n i t e Element A n a l y s i s of S t r e s s e s and Movements i n Embankments During Construction", Geotechnical Engineering Research Report No. TE-69-4, Department of C i v i l Engineering, Unviersity of C a l i f o r n i a , Berkeley, Nov. 28. Lee, F.H., (1983), " P a r t i a l Liquefaction i n Centrifuge Model Embankment i n an Earthquake", M.Phil. Thesis, Engineering Department Cambridge Uni v e r s i t y , Cambridge, J u l y . 29. Lee, K.W., (1965), " T r i a x i a l Compressive Strength of Saturated Sands Under Seismic Loading Conditions", Ph.D Thesis, U n i v e r s i t y of C a l i f o r n i a , Berkeley, ppl-521. 30. Lee, K.W., (1975), "Mechanical Model for the Analysis of Li q u e f a c t i o n of Horizontal S o i l Deposits", Ph.D Thesis, Department of C i v i l Engineering, U n i v e r s i t y of B r i t i s h Columbia, Vancouver, B.C., Canada, Sept. 161 31. Lee, K.W., and Finn, W.D.L., (1975), "DESRA-1: Program for Dynamic E f f e c t i v e Stress Response Analysis of S o i l Deposits Including L i q u e f a c t i o n Evaluation", S o i l Mechanics Series, No.36, Dept. of C i v i l Engineering, University of B r i t i s h Columbia, Vancouver, B.C., Canada. 32. Lysmer, J . , and Kuhlemeyer, R.L., (1969), " F i n i t e Dynamic Model for I n f i n i t e Media", Journal of the Engineering Mechanics, D i v i s i o n , ASCE, V o l . 95, EM4, Sept. pp859-877. 33. Lysmer, J . , and Waas, G., (1972), "Shear Waves i n Plane I n f i n i t e Structures", Journal of Engineering Mechanics D i v i s i o n , ASCE, V o l . 98, EMI, Feb., pp85-105. 34. Martin,. G.R., Finn, W.D.L., and Seed, H.B., (1975), "Fundementals of Liquefaction Under C y c l i c Loading", Journal of the Geotechnical Engineering D i v i s i o n , ASCE, V o l . 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 Zi e n k i e i c z , O.C, (1979), "Application of an Anisotropic Hardening Model i n the Analysis of E l a s t o - P l a s t i c Deformations of S o i l s " , Geotechnique, 29(1), ppl-34. 37. Nadim, F., and Whitman, R.V. (1982), " S l i p Elements for Earthquake Response", Proceedings, 4th International Conference on Numerical Methods i n Geomechanics, Edmonton, Alberta, Canada, pp61-67. 38. Newmark, N.M., (1959), "A Method of Computation f o r St r u c t u r a l Dynamics", Journal of the Engineering Mechanics D i v i s i o n , ASCE, V o l . 85, EM3, Ju l y . 162 39. Newmark, N.M., (1965), " E f f e c t s 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, C l i f f , N.J., ppl62-163. 41. Ozawa, Y., and Duncan, J.M., (1973), "ISBILD: A Computer Program for A n a l y s i s of S t a t i c S t r e s s e s and Movements i n Embankments", Geotechnical Engineering Research Report No. TE-73-4, Department of C i v i l Engineering, University of C a l i f o r n i a , Berkeley, Dec. 42. Prevost, J.H., (1978), "Mathematical Modelling of S o i l S t r e s s - S t r a i n Strength Behaviour", Proceedings, 3rd International Conference on Numerical Methods i n Geomechanics, Aachen, Germany, A p r i l 6, pp347-361. 43. Pyke, R., Seed, H.B., and Chan, C.K., (1975), "Settlement of Sands Under M u l t i d i r e c t i o n a l Shaking", Journal of the Geotechnical Engineering ASCE, V o l . 101, No. GT4, A p r i l , 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 D i v i s i o n , ASCE, V o l . 103, GT12, D e c , ppl419-1436. 45. Robertson, P.K., (1982), " I n - s i t u Testing of S o i l with Emphasis on i t s A p p l i c a t i o n to Liquefaction Assessment", Ph.D Thesis, submitted to Department of C i v i l Engineering, University of B r i t i s h 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 H o r i z o n t a l l y 163 Layered S i t e s " , Report No. EERC 72-12, Earthquake Engineering Reseach Center, University of C a l i f o r n i a , Berkeley, Dec. 47. S c h o f i e l d , A.N., (1981), "Dynamic and Earthquake Geotechnical C e n t r i f u g e M o d e l l i n g " , I n t e r n a t i o n a l Conference on Recent Development i n Geotechnical Earthquake Engineering and S o i l Dynamics, Missouri, U.S.A., A p r i l 28-May 2. 48. Seed, H.B., and Lee, K.L., (1969), "Porewater Pressures i n Earth Slopes Under C y c l i c Loading Conditions", Proceeding, 4th World Conference on Earthquake Engineering, Chile, V ol. I l l , p p l l l . 49. Seed, H.B., and I d r i s s , I.M., (1970), " S o i l Moduli and Damping Factors for Dynamic Response Analyses", Report No. EERC 70-10, Earthquake Engineering Research Center, Un i v e r s i t y of C a l i f o r n i a , Berkeley, December. 50. Seed, H.B., Lee, K.L., I d r i s s , I.M., and Makdisi, F.I., (1973), "Analysis of the Slides i n the San Fernando Dams During the Earthquake of Feb. 9, 1971", Report No. EERC 73-2, Earthquake Engineering Research Center, University of C a l i f o r n i a , Berkeley, June. 51. Seed, H.B., Pyke, R. and Martin, G.R. (1975), "Analysis of the E f f e c t of M u l t i - d i r e c t i o n a l Shaking on the Liquefaction C h a r a c t e r i s t i c s of Sands", Report NO. EERC 75-41, Earthquake Engineering Research Center, University of C a l i f o r n i a , Berkeley, C a l i f . Dec. 52. Seed, H.B., (1979), "Considerations i n the Earthquake-Resistant Design of E a r t h and R o c k f i l l Dams", 19th Rankine L e c t u r e , 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. S e l i g , E.T., and Chang, S.C. (1981), " S o i l F a i l u r e Modes i n Undrained C y c l i c Loading", Journal of the Geotechnical Engineering D i v i s i o n , ASCE, Vol. 107, GT5, May, pp539-551. 55. S e r f f , 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, U n i v e r s i t y of C a l i f o r n i a , Berkeley, Sept. 56. U n i t e d S t a t e s N a t i o n a l Research C o u n c i l (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), " E f f e c t of S t a t i c Shear on Liquefaction P o t e n t i a l " , Journal of the Geotechnical Engineering D i v i s i o n , ASCE, Vol. 105, GT10, OCT., ppl233-1246. 58. Vaid, Y.P., Byrne, P.M., and Hughes, J.M.O., (1981), " D i l a t i o n Angle and Liquefaction P o t e n t i a l " , Journal of the Geotechnical Engineering D i v i s i o n , ASCE, V o l . 107, GT7, J u l y . 59. Varadarajan, A., and Mishra, S.S., (1980), "Stress-Path Dependent S t r e s s - S t r a i n Volume Change Behaviour of G r a n u l a r S o i l " , I nternational Symposium on S o i l s under C y c l i c and Transient Loading, Swansea, ppl09-119. 60. V e r r u i j t , A., (1977), "Generation and D i s s i p a t i o n of Porewater Pressures", F i n i t e Elements i n 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 S t r u c t u r a l Dynamics, V o l . 1, pp241-252. 165 62. Wedge, N.E. (1977), "Problems In Non-Linear Analysis of Movements i n S o i l s " , M.A.Sc. Thesis, Submitted to the Dept. of C i v i l Engineering, University of B r i t i s h Columbia, Vancouver, B.C., Canada. 63. Yoshimi, Y., and Tokimatsu, K. (1977), "Settlement of Buildings on Saturated Sand During Earthquakes", S o i l s and Foundations, Japanese Society of S o i l Mechanics and Foundation Engineering, V o l . 17, No. 1, March, pp23-38. 64. Zienkiewicz, O.C., and Cheung, Y.K., (1967), F i n i t e Element Method i n S t r u c t u r a l and Continuum Mechanics", McGraw H i l l Book Company. 166 APPENDIX I FINITE ELEMENT FORMULATION In the f i n i t e element analysis the entire domain of i n t e r e s t i s divided into a f i n i t e number of elements. The v a r i a t i o n of displacement within a f i n i t e element i s assumed to be given by, {tf} = [N] {6} where ( A l . 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] = i n t e r p o l a t i o n function. The type of element used i n the a n a l y s i s i s a 4 node isoparametric element. The term isoparametric implies common (i s o - ) parametric d e s c r i p t i o n of the unknown displacements and the geometry of the element. The same i n t e r p o l a t i o n functions N^ are used to express both the displacement and the geometry of the element. Isoparametric element formulation has a number of advantages; i t o f f e r s e f f i c i e n t i n t e g r a t i o n s , and d i f f e r e n t i a t i o n s and i t can handle curved and a r b i t r a r y geometrical shapes. The i n t e r p o l a t i o n function [N] can be selected such that i t can be expressed i n natural coordinates (s,t) which i s a system of co-ordinates i n t r i n s i c to an element ( F i g . A l . l ) . 168 In the isoparametric concept, the coordinates of a point i n the element i s given by, {X} = [ N ] {X±} (A1.2) If we consider a four-node q u a d r i l a t e r a l isoparametric element, the matrix [N] i s 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 ( A l . l ) and (A1.2) can be rewritten as, u l v l u 2 j U , r N , 0 N 2 0 N 3 0 N^ 0 -I v 2 M L 0 N x 0 N 2 0 N 3 0 N ^ u 3 (A1.3) v 3 u 4 v 4 and, ( x } = f N l 0 N 2 0 N 3 0 S l f 0 , y, V L 0 N x 0 N 2 0 N 3 0 N^-1 x 3 X l y i x 2 ?3 X 4 169 It can be seen that the transformation shown i n (A1.4), maps the q u a d r i l a t e r a l element into a square as shown i n F i g . A l . l . If plane s t r a i n conditions are assumed the s t r a i n vector e i s s i m p l y g i v e n by _et = { e x , e y , Y x y } where e x and e y a r e normal s t r a i n s and y x y 1 S t n e shear s t r a i n . These s t r a i n s are given i n terms of displacements as, x .e. e X _ M ox e y av ay Y ' xy 8u ay ax s t r a i n matrix e from A1.3 and A1.5, e x ax 0 ax 0 ax 0 wk ax 0 = 0 ay 0 a N 2 ay 0 3N 3 ay 0 ay aN 1 aN1 aN 2 8N 2 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 r e l a t i o n s h i p , u. u„ (A1.5) 2 (A1.6) (A1.7) So any d e r i v a t i v e aN ± aN ± ax aN ± ay as ax as ay as (A1.8) In matrix form, the global and l o c a l d e r ivatives can be written as, 170 5N± ax ax aN ± 8N ± — — as as - — - — as ax r - i ax (AI.9) BN. 3x ay. 3N, L J J o N . l A i . y ; at a t 8 t ay ay where, [ j ] = Jacobian matrix The de r i v a t i v e s of with respect to x, y can now be obtained from A1.9, as, 8N± aN± 5± • [ J ] _ 1 3N, a y a t where, [ j ] - 1 = Inverse of Jacobian matrix J , which i s simply, a t a s Now from A1.4, the Jacobian matrix can be written as, ax By. 1 ^ 1 t^± R -I _ as as _ i as x i • i as y i ... L J J " . 8N u 8N (A1.12) a x a y . _ i fc f \ a t a t ^ 5 t x i A o t 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, = [ h i Ti2) •21 (A1.13) 22 It should be noted that these I ^ j are f ( s , t ) . Therefore, s u b s t i t u t i n g this i n equation ( A L I O ) , y i e l d s , oN ± oNjL By" t l 1 1 L21 22 8N. as" aN. i at (A.14) The matrix [B] (equation A1.6) which relates the s t r a i n vector to nodal displacements, has derivatives of i n t e r p o l a t i o n 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 0 S + X12 9 t | o i aN I aN 2i a s + 122 at | A n bi I aN, aN, + i . 2i as 22 at aN. aN, + i 3N 2 3N 2 X n + xi2 a~T aN. 9N„ 12 at | x2i a£ + i 22 at | 172 0 3N 3 aN 3 0 hi as + at aN 2 3N 3 3N 3 as + I22 at 0 X 2 1 as + I22 a t a N 2 aN 3 m3 aN 3 B N 3 as + 112 at I21 as + I 2 2 at i n as + h2 a t aN^ a i ^ | [ n o s + T i 2 a t | aN„ aN, o A 2 i a s + x22 a t aN 4 mh x n a s - + x i 2 IT I 2 1 as~ + *22 IT | (A1.15) The s t r e s s e s and s t r a i n s a r e c o n n e c t e d t h r o u g h e l a s t i c i t y m a t r i x g i v e n by, {a'} = [ D ] {e} (A1.16) where, [ a ' ] = e f f e c t i v e s t r e s s v e c t o r . F o r 2D p l a n e s t r a i n c o n d i t i o n s , i t i s g i v e n by, i o ' } * - = { 0 X , O y , T x y } [ D ] = e l a s t i c i t y m a t r i x f o r 2D p l a n e s t r a i n c o n d i t i o n , [ D ] -B + ^ G Sym B - f G 4 B + j G 0 0 (A1.17) 173 Using the v i r t u a l work p r i n c i p l e , the i n t e r n a l work (Wj^) done by a p p l y i n g i n f i n i t e s i m a l v i r t u a l nodal displacement {£} i s , WIN = /// M dv (A1.18) where {e} = v i r t u a l s t r a i n s due to v i r t u a l displacement {£} and {a} = t o t a l stress vector (A1.19) The t o t a l stress vector can be s p l i t into e f f e c t i v e stress and porepressure vectors. i . e . : {a} = {a'} + {uQ} E f f e c t i v e stress Porepressure (A1.19) Vector Vector Here {a'= {a x, a y , T x y } and {xxj*- = {u D, u D , o} i n which u Q i s the porewater pressure. Now, s u b s t i t u t i n g {a} from equations (Al.19) and (A1.18), one gets, WIN = /// {e}* [{o'} + {u }] dv (A1.20) V Substituting for {a'} from A1.16 th i s reduces to, WIN = /// M + K) D V (AI.21) v 174 using equation (A1.7), {e} can be replaced by [B] {£}, then, W I N - J / J [B]' [ D ] [B] {6} + J J J {6} f c [ B ] t { u J dv (A1.21) V V But, external work done by the load vector {p} r i d i n g through the v i r t u a l displacement {?>}, i s simply, WEX = t p J (A1.22) The v i r t u a l work p r i n c i p l e gives, W I N 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. Aft e r s u b s t i t u t i n g t h i s i n (A1.23), and d i v i d i n g both sides by {?} t one gets, 175 {P} = [k]{6} + [k*] {uD} (A1.24) where, [k] = element s t i f f n e s s matrix l i t = [ / J [ B] [D][B] l Jl d s d t ] (A1.25) (assuming unit thickness) and, [k ] = porewater pressure matrix 1 1 t = [/ / [B] | j | dsdt] (A1.26) -1-1 The integrations shown above have to be evaluated numerically. The Gauss inte g r a t i o n technique has been employed and the number of points used are 2 x 2 . The formulation presented here i s for any l i n e a r e l a s t i c m a t e r i a l . For incrementally e l a s t i c a n a l y s i s , the displacements, stresses and moduli values should be simply replaced by incremental displacements, incremental stresses and tangent moduli. Af t e r evaluating the incremental load vector j), element tangent s t i f f n e s s m a t r i x [ k t ] , porewater p r e s s u r e m a t r i x [ k * ] , and a l s o estimating the incremental porewater pressure u Q , for a l l the elements the global incremental load-displacement r e l a t i o n s h i p can be formed. This w i l l lead to, {P} = [K t]{A} + [K*]{U} (A1.27) 176 i n which {P}, [ K T ] , {A}, [ K ] and {u} are relevant v a r i a b l e s i n glo b a l axes. By solving t h i s equation the displacement f i e l d {A} can be obtained, and i t can be used to c a l c u l a t e element s t r a i n s and stresses using equations (A1.6.) and (A1.16) r e s p e c t i v e l y . Since the shape function gives l i n e a r s t r a i n v a r i a t i o n within an element, the s t r a i n s and therefore, stresses vary within an element. For convenience, average stress and s t r a i n of an element are computed at the centre of gravity of the element. In Chapter 4, and Chapter 5, i t i s required to express s t r a i n s and stresses i n terms of nodal forces. R e c a l l from equation (A1.24) the nodal forces are given by, {P} = J / J !>]' [D] [B] dv {£} (A1.28) V But, str a i n s are connected to the matrix [B J i n equation (A1.7), {e} = [B] {6} (A1.29) Therefore, from (A1.28) and (A1.29), the nodal forces can be written i n terms of str a i n s as, {P} " J J J [ 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 i n terms of stresses as, {P} = /// [ B f {a'} dv (A1.32) V The equation (A1.30) and (A1.32) can now be used to express s t r a i n s and stresses i n terms of element nodal forces. 178 APPENDIX I I STIFFNESS MATRIX FORMULATION FOR THE SLIP ELEMENT As o u t l i n e d i n Chapter 4, that the f o r c e - d i s p l a c e m e n t r e l a t i o n s h i p at any point within a s l i p element has been assumed to be given by, f K 0 w (fS} = U S K ] [ W S ] (A2.1) i . e . , f_ = ICQ W where, f g and f n = shear and normal stresses Kg, ^ = j o i n t s t i f f n e s s i n shear and normal d i r e c t i o n s wg, wn = shear and normal displacements The e l a s t i c stored energy, 0 E i n a s l i p element due ( F i g , (A2.1) to applied forces can be obtained by, 0 E = 2 t f dZ (A2* 2) i n which L i s the t o t a l length of the s l i p element. A factor half i s included because the r e l a t i o n s h i p between f_ and w i s assumed to be l i n e a r . From (A2.1), 0g now can be written as, Fig. A2.1. Slip Element. s = tangential direction U j = tangential displacement o f node i n = normal direction V j = normal displacement of node i 180 0 E = 2 /o W ' k o W d A ( A 2 * 3 ) Since the v a r i a t i o n of displacements (u, v) within an element i s l i n e a r the displacement at any point which i s at a distance, A, from node I on the bottom edge I J of the element i s , u U) = £ U j + (1 - L> U I ( A 2 - 4 > bottom In a s i m i l a r manner the following equations can be written f o r , 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 V J + ( 1 + L } V I (A2.6) bottom v (x) = L VK + ( 1 " L? V L ( A 2 ' 7 ) top where, u^, v^ r e f e r to displacement i n tangential and normal d i r e c t i o n 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 w n = v (A) - v U ) (A2.9) top bottom Now s u b s t i t u t i n g f o r u t o p > u b o t t o m , v t o p a n d v b o t t o m f r o m equations (A2.4) to (A2.7), i n equations (A2.8)-and (A2.9), w = H 1 -i> ~i i u - ft] (A.2.10) and, w_ = [-(1 - f) f f (1 ( A 2 . l l ) From equation (A2.1), w i s , U l v, w = [ S ] = — L w J -A 0 -B 0 B 0 A 0 [ ] 0 -A 0 -B 0 B 0 A J \ V K \ (A2.12) i n which, A = 1 " L A N D B - L In matrix form the equation (A2.12) i s , 182 w = C 6 2x1 2x§ 8X1 ( A 2 . 1 3 ) where, = the i n t e r p o l a t i o n matrix and 6^ = nodal displacement matrix Substitution of ( A 2 . 1 3 ) i n ( A 2 . 3 ) gives, K = 9 t s T cl K CJ d* E 2 •'o o o o ( A 2 . 1 4 ) Performing the matrix m u l t i p l i c a t i o n , T C k C o o o - 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 A 2 0 A B K 0 - A B K 0 - A 2 K 0 s s s s 0 A 2 K 0 A B K 0 - A B K 0 - A 2 K n B 2 K n n i A B K 0 0 - B 2 K 0 - A B K 0 s s s s 0 A B K 0 B 2 K 0 - B 2 K 0 - A B K n n n - A B K 0 - B 2 K 0 B 2 K 0 A B K 0 s s s B 2 K s 0 - A B K 0 - B 2 K 0 0 A B K n n n n - A 2 K 0 - A B K 0 A B K 0 A 2 K 0 s s s s 0 - A 2 K 0 - A B K 0 A B K 0 A 2 K ( A 2 . 1 5 ) To perform the in t e g r a t i o n shown i n equation A 2 . 1 4 , one should know integrations of, 183 / L A 2 dA, f L B 2 dA and f L AB dA J o J o J o These are simply, J L J o A 2 dA = 3 o (1 - ft2 dA L 3 / L J 0 B 2 dA = = J o ( f )2 dA L 3 (A2.16) and, Jo AB dA = • J L ; o ( 1 - f t f dA 6 Now the equation (A2.14) can be written as, 0 E " 2 ° T K s n 6 ( A 2' 1 7> 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 R e c a l l that e l a s t i c stored energy 0g i n the formulation of a l i n e a r e l a s t i c f i n i t e element i s , 0E = j 6 T K 6 (A2.18) 184 Where K i s the s t i f f n e s s matrix of the f i n i t e element. Now, comparing (A2.17) and (A2.18), the s t i f f n e s s matrix for s l i p element can be deduced as K„_. —sn 185 APPENDIX III STEP BY STEP INTEGRATION For the proposed incrementally e l a s t i c dynamic analysis i n the time domain the " e l a s t i c " properties have to be modified for every time step. The e l a s t i c parameters depend on the l e v e l of s t r a i n i n the deposit. Therefore, the displacement f i e l d i n 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 i n t e g r a t i o n i s very popular and extensively used i n dynamic analyses. This method b a s i c a l l y provides numerical s o l u t i o n i n time domain, where the s o l u t i o n i s advanced by one d i s c r e t e step at a time. In t h i s method, two parameters <* and p are used so that the v e l o c i t y and displacement at time t+At can be expressed i n terms of a c c e l e r a t i o n , v e l o c i t y and displacement at time t, and of the known acc e l e r a t i o n at time t+At. For convenience l e t us define that, T = t+At Then the r e l a t i o n s h i p i n 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) ( A t ) 2 {X} t + B A t 2 {X} T (A3.2) Newmark (1959) proposed that « = 1/2 and B = 1/4 be used for an u n c o n d i t i o n a l l y s t a b l e i n t e g r a t i o n procedure, which i n c i d e n t a l l y corresponds to a constant average acceleration method of i n t e g r a t i o n . I f = = 1/2 and 8 = 1/6 are used then t h i s method gives a l i n e a r v a r i a t i o n of acc e l e r a t i o n within the time step. Re-writing the incremental equilibrium equations from Chapter 5, [M] {AX} + [C] {AX} + [ K t ] t {AX} = {AP} (A3.3) Now su b s t i t u t i n g f o r , {AX} {AX} {AX} (A3.4) and, i n equation (A3.3), one gets, [M] {^ - X j + [C] {Xj, - X t} + [ K t ] t {^ - X t} = {AP} (A3.5) From e q u a t i o n s (A3.1) and ( A 3 . 2 ) , {x}^, and {^}^ can be expressed i n terms of other variables as follows, 187 [ | £ 2 + [ c ] p £ + [ * t ] t ] (AX} = {AP} + [M] {E} t + [C] {F} t W T = m 2 [ { A X } " A t { x } t " (i ~ P> A t 2 W J ( A 3 * 6 ) and, W T - (X} t + (!--) At {X} (A3.7) + ^ [{AX} - At{x}t - B) At* {X}J Substituting for {x}-£ and {x},p i n equation (A3.5), W 2 [ {AX} - A t { X } t - (j - 6 ) At* {X} t - BAt 2 {x}J + [C] [ ( 1 — ) At {X}fc + ^ [{AX} - At {X} t (A3.8) - £ - B) At2 {X}J + [ K t ] t {AX} = {AP} Co l l e c t i n g terms, and defining following s i m p l i f y i n g symbols, and, . W t = f W t + ' ( 2 p " 1 } A t ^ t ( A 3 ' 1 0 ) the equation (A3.8) can be reduced to, (A3.11) 188 Recall from Chapter 5, Section 5.1.2, cor r e c t i o n forces should be applied to restore t o t a l equilibrium at time T. The { P c o r r } evaluated using equation (5.11) can be added to r i g h t hand side of the equation (A3.11). Then the equation (A3.11) i s , ^ + [ c ] p t r + [ K T ] T ] w = {AP} + [M] {E} t + [C] {F} t (A3.12) where, {AP} = {AP} + { P c o r r } (A3.13) The only unknown i n the above equation i s {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 i s simply, {X} T = {AX} + {X} t (A3.16) 189 From equation (A3.2), an expression for {x}T i s , W T • "p^2 [{AX} - At {X} t " (| - 8) A t 2 {X} t] (A3.17) S u b s t i t u t i n g f o r {E} t from (A3.9), the eq u a t i o n (A3.17) can be s i m p l i f i e d as, W T = "pit"2 {AX} - {E} t + {X} t (A3.18) From equations (A3.1) and (A3.2) an e x p r e s s i o n f o r {x}^ , a f t e r rearranging terms i s , {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 i n t e g r a t i o n , the following sequence of c a l c u l a t i o n s have to be performed for every time step. 1. I n i t i a l v e l o c i t y {x} t and displacements {x} t are known eith e r from values at the end of the preceding increment or as i n i t i a l conditions of the problem. Based on these values, { E } t , { F } t , and { p C O r r } » a r e e v a l u a t e d u s i n g equations (A3.9), (A3.10) and (5.11). 2. With these values and the known non-linear properties of the s o i l deposit, the damping matrix [c] and [ K t ] t are evaluated according to appropriate equations i n Chapter 5. 190 3. The matrix [D] i s then calculated using equation (A3.15). 4. Using the increment i n base acceleration value at time t, i t i s possible to evaluate the r i g h t hand side of the equation (A3. 12). 5. The equation (A3.12) can now be solved for {AX} and then displacement, a c c e l e r a t i o n and v e l o c i t y vectors can be evaluated from equations (A3.16), (A3.18) and (A3.19) r e s p e c t i v e l y . When step 5 has been completed, the analysis f or t h i s time increment i s f i n i s h e d and the entire process may be repeated for the next time step. Obviously t h i s process can be c a r r i e d out consecutively for any desired number of time increments; thus the complete response h i s t o r y can be computed. Two important aspects have to be considered i n any numerical i n t e g r a t i o n procedure. They are the accuracy and s t a b i l i t y of the procedures. Accuracy ref e r s to how well the numerical s o l u t i o n matches the exact continuous s o l u t i o n . S t a b i l i t y r e f e r s to whether extraneous solutions are introduced i n such a way that they increase rather than decay, and thus come to dominate the r e s u l t s . Usually there i s an upper l i m i t to At that i s necessary to guarantee s t a b i l i t y , and the value of that l i m i t depends on the type of element s t i f f n e s s and mass matrix as well as on « and 8. With a l i n e a r a c c e l e r a t i o n assumption (<= = 1/2, 6 = 1/6) the analysis w i l l give good accuracy i f the shortest period of the deposit i s 5 to 10 times greater than At (Clough and Penzien 1975). The l i n e a r a c c e l e r a t i o n method i s only c o n d i t i o n a l l y stable, and i t w i l l blow up i f i t i s applied to s o i l structures with the shortest period le s s than about 1.8 times the i n t e g r a t i o n i n t e r v a l . Thus the time increments must be made 191 short r e l a t i v e to the least period of v i b r a t i o n contained i n the system regardless of whether the higher modes contribute s i g n i f i c a n t l y or not. In the analysis using f i n i t e element procedures the shortest periods of v i b r a t i o n i n general may be several orders of magnitude less than the periods associated with the s i g n i f i c a n t response. In these cases, the l i n e a r a c c e l e r a t i o n method cannot be used because of the very short time step required to avoid i n s t a b i l i t y ; instead, an unconditionally stable method i s required which w i l l not blow up regardless of the time step. Several unconditionally stable methods are a v a i l a b l e . The constant average a c c e l e r a t i o n method (« = 1/2 , 8 = 1/4) i s one of the simplest of these methods. But t h i s assumption has been reported not to give good r e s u l t s than the methods with l i n e a r a c c e l e r a t i o n assumption. The Wilson 9-method (Wilson, et a l . 1973), i s also an unconditionally stable method. This method i s a modification of l i n e a r a c c e l e r a t i o n method and i s reported to be the best of a l l unconditionally stable methods (Clough, et a l . 1975). The Wilson 8-method i s based on the assumption that the acceleration varies l i n e a r l y over an extended computation i n t e r v a l , %, such that, x = 9 At where 0 > 1.37 (A3.20) When 0 = 1 , t h i s method reverts to the standard l i n e a r a c c e l e r a t i o n method. The analysis procedure i s exactly the same as the procedure presented above except that i n the equations the time step At has to be r e p l a c e d by T and a l s o the equations to e v a l u a t e {x}-j,, {&} ,^ and { x } T have to be modified. 192 Since t h i s i s e s s e n t i a l l y a l i n e a r a c c e l e r a t i o n 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] + [ K f c ] t ] { A X } m { ~ p } + [ M ] [ c ] { p } t ( A 3 > 2 1 ) i n which, { E*t = x ^ t + 3 M t ( A 3 ' 2 2 ) and, {F} t = 3 {X} t + \ {X} t (A3.23) After evaluating {AX}, which i s over an extended time increment the displacement, v e l o c i t y and a c c e l e r a t i o n values should be computed at time t = T. The a c c e l e r a t i o n at t = T can be evaluated using, W T = [*2 {AX} + { E } J {X} t (A3.24) {x},p, from equation (A3.19) with « = 1/2 and 8 = 1/6 i s , {X} T = W t + f [ W t + {X} T] (A3.25) 193 and also the displacement {x}^, from (A3.2) i s , {X} T = ( X } t + At ( X } t + *f - [ {x } T + 2 { X } J (A3.26) f I t must be remembered that s t a b i l i t y i n numerical i n t e g r a t i o n does not guarantee accuracy or vice versa. The Wilson 0-method imposes a r t i f i c i a l damping i n higher modes. But knowing that the response due to higher modes of v i b r a t i o n contributes very l i t t l e to the true response of structures, this method i n a way f i l t e r s out the high frequency response. Therefore, t h i s method has been found to y i e l d r e a l i s t i c r e s u l t s i n 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
pdf
Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.
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 |
AggregatedSourceRepository | 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}]}"
data-media="{[{embed.selectedMedia}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0062642/manifest