UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Mechanical model for the analysis of liquefaction of horizontal soil deposits Lee, Kwok Wing 1975

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

Download

Media
831-UBC_1975_A1 L43_5.pdf [ 8.54MB ]
Metadata
JSON: 831-1.0063017.json
JSON-LD: 831-1.0063017-ld.json
RDF/XML (Pretty): 831-1.0063017-rdf.xml
RDF/JSON: 831-1.0063017-rdf.json
Turtle: 831-1.0063017-turtle.txt
N-Triples: 831-1.0063017-rdf-ntriples.txt
Original Record: 831-1.0063017-source.json
Full Text
831-1.0063017-fulltext.txt
Citation
831-1.0063017.ris

Full Text

MECHANICAL OF  MODEL  FOR  LIQUEFACTION SOIL  OF  THE  ANALYSIS  HORIZONTAL  DEPOSITS  by KWOK WING LEE B.So. (Eng.), University of Hong Kong, I966 M.Sc. (Eng.), University of Hong Kong, I968  A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY in the Department of C i v i l Engineering  We accept t h i s thesis as conforming to the required standard  THE UNIVERSITY OF BRITISH COLUMBIA September, 1975  In presenting  t h i s thesis in p a r t i a l f u l f i l m e n t of the requirements f o r  an advanced degree at the University of B r i t i s h Columbia, I agree that the Library shall make i t freely a v a i l a b l e for reference and study. I further agree that permission for extensive  copying of t h i s thesis  for s c h o l a r l y purposes may be granted by the Head of my Department or by his 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 thesis for 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 i v i l E n g i n e e r i n g The University of B r i t i s h Columbia 2075 Wesbrook Place Vancouver, Canada V6T 1W5  Date  Sept.  1975  i  ABSTRACT  During the development of l i q u e f a c t i o n In a s o i l deposit  subjected  to v i b r a t i o n  there are two processes  which work i n opposite d i r e c t i o n s . tendency  under  cyclic  loading  The volume compaction  causes  the pore  water  pressure to r i s e , and the d i s s i p a t i o n of excess pore water pressure  (consolidation) decreases i t .  Recently,  Martin,  Finn and Seed(l975) studied the mechanics .of pore pressure generation  of  water  a s o i l sample subjected to c y c l i c  loading and a r e l a t i o n between shear s t r a i n cycles, volume compaction  and  established.  pore  water  pressure  thesis  f o r saturated granular  under c y c l i c simple shear conditions. hysteretlc  compaction,  was  A material model based on t h i s r e l a t i o n s h i p  i s developed i n t h i s  a  increment  stress-strain  pore water  The model includes  relationship,  pressure  soil  rise  volume  and d i s s i p a t i o n .  Using t h i s proposed comprehensive material model, a global mechanical  model  i s constructed  to  simulate  the  l i q u e f a c t i o n (including consolidation) behavior of a thick horizontal motion.  deposit  when  In t h i s way  the  subjected coupled  response, pore water pressure the deposit  under  numerical techniques discussed i n d e t a i l .  rise  to  horizontal  problems and  of  dynamic  consolidation of  seismic loading can be analysed. used  to  solve  base  such  The  problems are  The response of a t y p i c a l saturated  11  sand d e p o s i t under earthquake l o a d i n g  i s determined u s i n g  the proposed model and the r e s u l t s show t h a t t h e model can p r e d i c t the v a r i o u s phenomena t h a t s a t u r a t e d sand d e p o s i t s e x h i b i t d u r i n g earthquakes.  The g l o b a l model a l s o makes  c l e a r the i n f l u e n c e o f p e r m e a b i l i t y p o t e n t i a l o f the s o i l  deposit.  on  the  liquefaction  ill  TABLE  OF  CONTENTS Page  ABSTRACT  1  TABLE OF CONTENTS  i i i  LIST OF TABLES  vl  LIST OF FIGURES  vil  NOTATION  xii  ACKNOWLEDGEMENTS  rvii  CHAPTER 1. INTRODUCTION. 1.1 General C h a r a c t e r i s t i c s of S o i l Liquefaction.  1  1.2 Description of Some Liquefaction Case H i s t o r i e s .  2  1.3 Spontaneous Liquefaction, Cyclic Mobility and Liquefaction.  7  1. k Purpose and Scope of Research.  10  CHAPTER 2. REVIEW OF PAST INVESTIGATIONS ON LIQUEFACTION OF SATURATED SANDS. 2.1 Liquefaction P o t e n t i a l Deposit Based on Experience.  of a Field  12  2.2 Laboratory Investigation of Liquefaction Using Scaled Sand Models.  18  2.3 Laboratory Investigation of Liquefaction Using Small Sand Samples.  26  2. JJ- Some Methods of Determining Liquefaction P o t e n t i a l of a S o i l Deposit.  38  CHAPTER 3.  CHAPTER 4.  CHAPTER 5.  CHAPTER 6  MODELLING OF SATURATED GRANULAR SOIL UNDER CYCLIC SIMPLE SHEAR CONDITIONS. 3.1  Volume Change C h a r a c t e r i s t i c s During C y c l i c Drained Shear.  3.2  Stress-Strain Under Cyclic Conditions.  3.3  Modelling Of Saturated Granular Soil Under Cyclic Simple Shear.  3.4  Computation o f S o i l Behavior Using the Proposed S o i l Model.  Relationship Simple Shear  MECHANICAL MODEL FOR THE ANALYSIS OF LIQUEFACTION. 4.1  Introduction.  4.2  F o r m u l a t i o n o f the Problem.  4.3  S o l u t i o n Scheme.  APPLICATION OF THE MECHANICAL MODEL A DYNAMIC ANALYSIS. 5.1  A Hypothetical S o i l  5.2  Response Characteristics the S o i l D e p o s i t .  5*3  Magnitude o f Pore Pressure R i s e and i t s E f f e c t on the Dynamic Response of the Deposit.  5.4  Effect of Dissipation.  Pore  SUMMARY AND CONCLUSIONS. 6.1  Summary.  6.2  Conclusions.  Deposit. of  Pressure  V  6.3 Suggestions Studies.  For  Further 147  LIST OF REFERENCES  148  APPENDIX I  155  APPENDIX I I  162  APPENDIX I I I  169  APPENDIX IV  173  vi  LIST OF TABLES Table 5-1 5-2  5-3  gaffe  Data f o r Layer S o i l Deposit  Approximation  Frequencies and Phase the Simulated Accelerogram  to 118  Angles f o r Earthquake  S'l-Values f o r Different Layers  130 130  vii  LIST OF FIGURES Figure  Page  1- 1  Modes of Foundation F a i l u r e  5  2- 1  C r i t i c a l N-Value versus Depth  15  2-2  Zone of Liquefaction  17  2-3  P r o b a b i l i t y of Liquefaction  19  2-4  Critical  Acceleration  f o r Various 21  Sands 2-5  2-6  E f f e c t of Acceleration Level Resistance to Liquefaction  on  Effect  on  of  Resistance  Surcharge  2k  Pressure  25  to Liquefaction  2-7  Idealised F i e l d Loading. Conditions  28  2-8  Idealised T r i a x i a l Test Conditions  29  2-9  C y c l i c Shear Stress Wave Forms  33  2-10  Typical Cyclic  Simple  Shear  Test  Data  Jk  2-11  T y p i c a l C y c l i c T r i a x i a l Test Data  2-12  Resistance Simple Shear  to  Liquefaction  35  in 37  e  F  Line f o r Sand B  Method of Potential  Evaluating  Liquefaction  Compaction of Dry Coarse Sand Effect of Compaction  vertical  Stress  E f f e c t of Relative Density Shear S t r a i n on Compaction Incremental Curves  Volumetric  on  and  Strain  Comparison of Measured and Computed Volumetric S t r a i n Hyperbolic S t r e s s - S t r a i n Curve Equivalent Shear I n i t i a l Loading  Modulus  D e f i n i t i o n of Equivalent Modulus And Damping Symmetrical Loop  for  Shear Ratio  Hysteretic C h a r a c t e r i s t i c s One-Dimensional S o i l Model Force-Displacement Relationship S l i p Elements  of  Volume Compaction Caused by C y c l i c Loading 'Reduced' Strain  Incremental  Volumetric  Increments of Volumetric S t r a i n , €. Hyperbolic Stress-Strain Loops (Drained Test) Showing E f f e c t of Hardening Shear Modulus as a Function of C y c l i c Shear S t r a i n and Volumetric Strain Shear Modulus as a Function C y c l i c Shear S t r a i n and Number Cycles  of of  Computed Damping Ratios Shear Modulus as a Function of Shear S t r a i n and T o t a l Compaction Damping Ratios f o r Monterey Sand One Dimensional  Unloading  Curves  Hyperbolic Stress-Strain Loops (Undrained Test) Showing Softening due to Pore Pressure Rise S t r a i n and Pore Pressure Undrained C y c l i c Shear Test  during  X Figure  Page  J  4-1  Idealization  of  the  Response  Problem  97  4-2  Approximation by Lumped Mass System  100  4-3  Consolidation Model  107  4- 4  Flow Chart  114  5- 1  Properties of a S o i l Deposit  116  5-2  Input Base Acceleration  120  5-3  Surface Acceleration Response  122  5-4  S t r e s s - S t r a i n Response  5-5  Surface Acceleration Response Curve  5-6  Surface Acceleration  123  from  'SHAKE  1  Program  126  5-7  Simulated Earthquake Accelerogram  5-8  Surface  Acceleration  5-10  129  Response (No  Pore Pressure) 5-9  124  S t r a i n and Stress Responses f o r Layer 11 (No Pore Pressure) Surface Acceleration and Displacement Responses (Pore Pressure Generated)  131 133 136  xi  Figure 5-11  Page  Stress,  Strain  and  Pore Pressure  Responses f o r Layer 11 5-12  Stress-Strain Response f o r Layer  135 11 136  (Pore Pressure Generated) 5-13  Pore  Pressure  Distribution  at 139  Various Times (No Dissipation) 5-14  Pore Pressure Distribution Different k Values  for 141  xii  NOTATION  A  Amplitude of v i b r a t i o n  A^, A2. A-j  Constants f o r the shear modulus function  a^  Constants  t  &2  f o r the increase  in  shear  modulus due to compaction B^, B2, B 3  Constants f o r the shear modulus function  bj_, b2  Constants  f o r the increase i n l i m i t i n g  shear stress due to compaction Cit C2, C3, C4 Constants f o r the volume change function [C ]  Damping matrix  Cq  Equivalent  c±  Damping c o e f f i c i e n t f o r the i - t h layer  c  S p e c i f i c heat of a s o l i d  D  Damping r a t i o  e  D D  Maximum damping r a t i o  m a x  Relative  r  D50  viscous damping f a c t o r  T  n  e  density  maximum p a r t i c l e - s i z e of the smallest  f i f t y percent of a s o i l e  Void r a t i o  e"  C r i t i c a l void r a t i o l i n e  fi  Stress-strain-rate function  F  S t r e s s - s t r a i n function Shear modulus Maximum shear modulus I n i t i a l maximum shear modulus Maximum shear modulus at the n-th cycle Thickness of a s o i l  deposit  Thickness of the I-th layer Coefficient of the earth pressure at rest Bulk modulus of water One-dimensional rebound soil  modulus  of the  skeleton  S t i f f n e s s matrix S t i f f n e s s value of the i - t h layer Value of permeability Thermal conductivity Mass matrix Mass lumped at the top of the i - t h layer Coefficient of volume compressibility C r i t i c a l standard penetration  value  Number of cycles to l i q u e f a c t i o n Equivalent  number  of constant amplitude  shear stress cycles Porosity I n e r t i a l force vector In s i t u v e r t i c a l stress  Position vector Shear  strength  Time Time Volume change function Uniformity c o e f f i c i e n t Velocity Energy dissipated per cycle of loading F i n i t e difference approximation to cr  w  Displacement of the i - t h layer -  relative  to the base Displacement of the i - t h layer Base displacement Height measured from the bedrock Height of the top of the i - t h layer C r i t i c a l acceleration Constants f o r the numerical integration Damping c o e f f i c i e n t s Shear s t r a i n Hyperbolic shear s t r a i n Unit weight of water Amplitude of a shear s t r a i n cycle Increment  i n shear s t r a i n  Volume change increment Recoverable  component  of  the  volume  change increment Compaction compnent of the volume change Increment Strain Volumetric s t r a i n Recoverable component of  the  volumetric  strain Compaction  component  of  the volumetric  strain Temperature Mass densityStress Mean normal stress Octahedral normal stress Vertical  Stress  I n i t i a l v e r t i c a l stress Pore water pressure Minor p r i n c i p a l stress at f a i l u r e E f f e c t i v e stress Shear stress Amplitude of a shear stress cycle C y c l i c horizontal shear stress Equivalent  10  cycle  shear  stress  amplitude that i s developed C y c l i c shear stress amplitude required to  xvl  cause l i q u e f a c t i o n i n 10 cycles Tj  Maximum  n a x  shear stress that can be applied  to a s o i l element Maximum shear stress at the n-th cycle <j>' +1» +2*  E f f e c t i v e angle of shearing resistance ^3*  4*4 Constants f o r the volumetric function  compaction  xvil  ACKNOWLEDGEMENTS The author wishes to express his gratitude to his principal advice  adviser  Dr. W.D. Liam Finn  and guidance  research.  Thanks  during are  the entire also  due  Dr. R.G. Campanella, Dr. G.R. Martin their  f o r his invaluable course  of the  to Dr. P.M. Byrne,  and Dr. Y. Vaid f o r  encouragement and constructive discussion at various  stages of the research. The National Research Council of Canada the  financial  possible.  assistance  which  made  provided  this investigation  This assistance i s g r a t e f u l l y acknowledged. The author deeply  appreciated  the support  and  consideration given t o him by his wife during the period of his studies.  1  CHAPTER 1 INTRODUCTION  1.1  GENERAL CHARACTERISTICS OF SOIL LIQUEFACTION. It  has  relatively  been  loose  observed  sandy  soil  that  when  deposits  saturated  are  subjected  to  Imposed  by  r e p e t i t i v e or c y c l i c shear l o a d i n g , such as those an  earthquake,  turned  into  extreme  they  a  case  often  'quick' where  exhibit  or  great  the  liquefied  phenomena of b e i n g  condition.  earthquake  and  shocks  In  occur  the  In  an  a l l u v i a l r e g i o n , the pore water i n the l i q u e f i e d s o i l which i s under v e r y l a r g e p r e s s u r e may temporary pressure  geysers.  In  generated  is  be f o r c e d t o the s u r f a c e  the  less  severe  cases  excess  pore  original  occur  when  pressure i s d i s s i p a t e d . It  masses  the  o n l y a s m a l l f r a c t i o n of the  e f f e c t i v e normal p r e s s u r e and minor s e t t l e m e n t may  forming  is  generally  compact  vibration.  or  If  accepted  contract  the s o i l  In  that  volume  dry granular when  soil  subjected  i s s a t u r a t e d and drainage  prevented  then i t i s not d i f f i c u l t  to see t h a t the pore  will  a g r a n u l a r s o i l , the maximum s h e a r i n g  build  resistance  up.  For  that  can  be  developed  on  a  water  to  plane,  pressure  or  shear  s t r e n g t h , s, i s g i v e n by the f o l l o w i n g e x p r e s s i o n , s  =  (<r -<r )tan<t>' v  w  (i.-D  2  where  s  I s t h e shear s t r e n g t h ,  0* i s t h e t o t a l normal p r e s s u r e , v  cr $'  w  i s t h e pore water p r e s s u r e , and, i s the e f f e c t i v e angle o f i n t e r n a l r e s i s t a n c e .  As t h e pore p r e s s u r e r i s e s t h e e f f e c t i v e normal s t r e s s , 0^.-0^, decreases  causing  a  decrease  i n the shear s t r e n g t h , s .  In  the extreme case where the v i b r a t i o n a l d i s t u r b a n c e i s o f l a r g e magnitude and r e l a t i v e l y l o n g d u r a t i o n the pore water p r e s s u r e may r i s e t o such an extent t h a t t h e may  become  zero.  The  effective  condition  of  normal  zero  or  stress  near  c o n f i n i n g p r e s s u r e makes t h e s o i l behave l i k e a l i q u i d , little equation describe  or  no  shear r e s i s t a n c e can be developed  (1-1) and hence the the phenomenon.  tanks and o t h e r b u r i e d surface.  term  the next  1.2  i s used  to  floating  up  to  the  s l o p e and dam f a i l u r e s caused soil  ground by t h i s  have  been  i n numerous i n v e s t i g a t i o n s and they a r e d e s c r i b e d i n section.  DESCRIPTION OF SOME LIQUEFACTION CASE HISTORIES. In  liquefaction areas. of  according to  LIQUEFACTION  type o f s t r e n g t h r e d u c t i o n i n the s u p p o r t i n g reported  since  There have been i n s t a n c e s o f sewer  objects  Foundation,  zero  the  San  Francisco  phenomena  were  Earthquake  observed  The quake o c c u r r e d on A p r i l  about 8.2.  The f o l l o w i n g  r e p o r t o f the earthquake:  of  1906  soil  i n numerous a l l u v i a l  18, 1906 w i t h a magnitude  i s quoted  from  the  official  "One o f the most common phenomena i n  3 such  a  situation  was  the expulsion  of water i n jets from  apertures which suddenly appeared i n the f l a t  lying  grounds.  The water was usually thrown into the a i r f o r several feet; i n some  cases  i t was reported to be as much as 20 feet, and the  e j e c t i o n continued f o r several minutes a f t e r  the  earthquake.  The continuance of the e j e c t i o n a f t e r the shock indicates that an  e l a s t i c stress has been generated i n the saturated ground,  which thus found r e l i e f In the expulsion ...  The  of  contained  water  water i n I t s passage to the surface brought  considerable quantities  of  fine  sand,  which  from  up its  p r e v a i l i n g l i g h t bluish gray colour was evidently derived from considerable depth. the  On the flood p l a i n of the Salinas r i v e r ,  sand was recognized by the people of the neighbourhood to  be the same as that of a stratum of sand pierced by wells at a depth of 80 f e e t . "  Usually the sand brought up i n t h i s way  formed funnel-shaped sand craters on the ground surface. abundance of saturated a l l u v i a l f l a t s with such an  craters  gave  Indication of the size of the area which l i q u e f i e d as well  as of the s e v e r i t y of the l i q u e f a c t i o n problem. not many structures were b u i l t i n such areas at the  The  earthquake  and  which  may  the time  of  l i t t l e or no damage to structures due to  s o i l l i q u e f a c t i o n was reported. flows,  Fortunately,  have  been  Numerous  cases  of  earth-  due to l i q u e f a c t i o n , were also  recorded. Liquefaction of a saturated damage  to  buildings  and  sandy  deposit  causing  structures occurred i n the Nllgata  4 The quake took place at around 13  Earthquake of 1964. on  June  16,  1964  i n the Niigata and Yamagata area i n Japan. less than 7.5  The magnitude of the quake was Scale  and  effects  i n the  Hi enter  Intensity i n the Niigata City was around 8 i n  the  the Modified M e r c a l l i Scale. the  hours  of  the  A d e s c r i p t i o n and  analysis  quake have been published i n a d e t a i l e d  report compiled by the E d i t o r i a l Committee of "General on the  Niigata  Earthquake  observation was  that  structures(up to  Osaki(1968)  1964,"(1968).  of  modern  1964)  the Shinano r i v e r  buildings,  bridges  and  b u i l t on the a l l u v i a l formation around  were  damaged  more  readily  than  others.  a  whole  damage i n the superstructure, i n d i c a t i n g  severe  foundation f a i l u r e . upheaval  Moreover, f o r t h i s type of  of  soli  overturning  predicted by  theory,  shown  as  other  reckoned that two t h i r d s of the damaged buildings  much  usual  Report  An i n t e r e s t i n g  In these areas Just s e t t l e d , t i l t e d or overturned as without  of  in  on  the  the opposite side of t i l t i n g or  conventional  figure  1-la,  foundation was  deformation  pattern s i m i l a r to figure 1-lb  contrast  to  the  failure  settlement  not observed but a was  and  failure  obtained.  tilting  of  In heavy  superstructures as described, l i g h t buried structures, such as concrete s e p t i c tanks,  'floated' during the quake and remained  protruding above ground surface a f t e r the disturbance In  the  light  of these data, Osaki reasoned that the s o i l at  such s i t e s behaved l i k e a l i q u i d and foundation Also  an  soil area  ceased.  liquefied near  the  under same  he the  concluded  that  the  earthquake v i b r a t i o n .  location,  but  on  which  no  (a)  (b)  Fig. 1-1  Modes of (after  Foundation  Osaki, 1968 )  Failure  6 superstructures  were  built,  geysers erupted during the which  continued  for  a  mud  volcanoes  earthquake  and  and  these  temporary outbursts,  while ' a f t e r the tremor passed away,  f i n a l l y died down giving r i s e to sand c r a t e r s . Not only was liquefaction  under  loose  earthquake  could also l i q u e f y when the quake  were  sandy  both large.  Intensity  sand.  large  to  and  duration  of  the  Seed(1972) reported that i n the San 1971,  medium dense sand foundation beneath induced  susceptible  v i b r a t i o n s , medium dense sand  Fernando Earthquake of February 9,  Plant  deposit  the l i q u e f a c t i o n of a the  Jensen  Filtration  l a t e r a l movements i n the f i l l above the  As a r e s u l t of t h i s , extensive damage was done to  the  s i t e as well as to structures which were under construction. Failures  caused  r e s t r i c t e d to foundations flows  by of  soil  liquefaction  structures.  are  Numerous  not  earth-  and slope f a i l u r e s have occurred during earthquakes due  to l i q u e f a c t i o n of granular s o i l masses or seams or lenses similar  material.  Seed(1968)  of  i n his Terzaghi Lecture gave a  large number of examples of slope f a i l u r e s or s l i d e s r e s u l t i n g from s o i l  liquefaction  failure  of  each case.  dynamic  analysis,  Seed  and  he  examined  the  mechanism  of  Using a recently developed method of et  al  (1975)  reconstructed  the  mechanics of s l i d i n g of the upstream part of the embankment of the  Lower Fernando Dam,  which moved some 70 feet or more into  the  reservoir.  found  liquefaction  They  developed  near  that  an  the base  extensive  zone  of  of the embankment and  7 t h e i r conclusion was confirmed by the sand  bolls,  cracks  filled  observed  phenomena  with l i q u e f i e d sand i n the s l i d e  mass and spreading of the embankment s o i l 250 feet beyond toe  of  the  of the embankment.  1.3  SPONTANEOUS  LIQUEFACTION,  CYCLIC  MOBILITY  AND  LIQUEFACTION. With regards to the term controversy  on  i t s usage.  liquefaction  there  is a  There has been no d e f i n i t i o n f o r  i t i n terms of stress or s t r a i n nor  has  there  been  any i n  terms of other phenomenological parameters.  Castro(l969) Hazen(1920) was  Allen  "liquefies" the  i n h i s Ph.D. t h e s i s pointed out that the  first  one  to  use  the  word  to describe flow f a i l u r e s and he continued to use  term i n association with flow s l i d e s i n h i s t h e s i s .  wrote  :  "Liquefaction  or  flow f a i l u r e o f sand Is caused by  substantial reduction of i t s shear sand  actually  acting  within  compatible  flows the  with  spreading sand  the  However, Terzaghi and  He  become  reduced  strength.  mass  of  out u n t i l the shear stresses so  shear  Peck(l968)  The  used  small  that  strength the  term  they  are  of the sand". spontaneous  l l q u e f a c t i o n to describe t h i s phenomenon o f flow f a i l u r e which occurred shock.  mostly  i n loose  fine  sand when subjected to mild  They suggested that the metastable  structure  formed  by the sand grains was responsible f o r t h i s kind of behaviour. It  was  because  of  such  a metastable structure of the sand  8  aggregate that the f a i l u r e process could take place i n a short  time.  Indeed,  the term spontaneous  very  liquefaction i s  very appropriate i n describing t h i s type of flow f a i l u r e s . Several researchers i n the f i e l d of e.g..  Seed  soil  et a l at the University of C a l i f o r n i a and Finn et  a l at the University of B r i t i s h Columbia, c a r r i e d loading  undrained  tests  out  cyclic  on t r i a x i a l as well as simple shear  specimens of granular material and observed the  dynamics,  In general  that  pore pressure i n the specimens continued to r i s e , leading  to a decrease i n the e f f e c t i v e confining (or normal) pressure, and f i n a l l y reached a stage where large deformations under  constant  cyclic  "liquefaction"  was  phenomenon.  This  stress  also  used  is  in  quite  resistance of the s o i l decreased cycles  of  little  or no  However,  stress  until  on  so  that  suitable with  dilative  samples  term  with  this  since  the  number of  f a i l e d when i t offered large of  deformation.  the sand  response  at  specimen the peak  an  additional  cycle  f o r very of  shear  loose and stress  l i q u e f a c t i o n can increase the maximum shear s t r a i n by fifteen  per cent.  response suggestion  of  some  of  shear  f o r example very large s t r a i n s cannot develop  abruptly i n f a i r l y dense sand, although loose  The  increasing  at a l l to  the denseness  there are varying degree of shear  connection  i t finally  resistance  depending  amplitude.  occurred  near  ten to  In view of the above described d i l a t i v e samples,  Castro(1969)•  Dr. A. Casagrande,  following  the  referred to the process of  9  the gradual development of pore straining  as  cyclic  pressure  mobility  and  with  cyclic  reserved  shear  the  term  l i q u e f a c t i o n f o r the f a i l u r e of s o i l under shear when the loss of strength due to r i s e i n pore pressure was not regained (due to contractive response) even at large s t r a i n s . It i s well-known that v i b r a t i o n causes compaction i n a granular s o i l mass. cyclic  In other words, I t can  stresses, e s p e c i a l l y  granular s o i l mass always strain  (contraction).  shear  cause Each  an  though  f o r constant amplitude  stresses, increase  additional  s t r a i n always causes an a d d i t i o n a l  be  said  acting  -that on  a  i n volumetric  cycle of stress or  Increment  i n contraction  stress or s t r a i n the increment  i s a decreasing function of the number of  cycles.  I f the  s o i l mass i s saturated and drainage (volume change) prevented, a. portion of the confining (or normal) pressure a c t i n g on the s o i l skeleton would an  incremental  be transferred to the pore water, causing  rise  i n pore pressure.  Since the volume of  s o i l always decreases when drained, as described pore  pressure  above,  w i l l always r i s e i n the undrained case.  the As a  r e s u l t there i s always the p o s s i b i l i t y of making the e f f e c t i v e confining pressure zero by continuing the process the  stresses  or  strains.  word  mobility  does  cycling  Though the term c y c l i c mobility  describes the c y c l i c nature of the the  of  not  softening of the material l i k e  process  quite  adequately  give the notion of the gradual the word  liquefaction  does.  It may be more suitable to use the term c y c l i c l i q u e f a c t i o n J  10 Terzaghl  and  Peck(1968)  In  the recent e d i t i o n of  t h e i r book, " S o i l Mechanics In Engineering Practice", added new  section  on  'Liquefaction under  Reversals of Stress or  Strain* to describe the behavior of the r i s e i n pore in  a  soil  failure. Include  mass  due  to  cyclic  loading  pressure  which may  lead to  Thus the broad usage of the word ' l i q u e f a c t i o n ' flow  failures,  a  to  with sharp r i s e i n pore pressure and  c a t a s t r o p h i c a l l y large strains (spontaneous  liquefaction),  as  well as limited f a i l u r e , with incremental or continual r i s e i n pore under  pressure and comparatively smaller strains ( l i q u e f a c t i o n reversals  of  generalization.  As  w i l l be used i n  this  More  discussion  stress  or  strain),  is  not  an  over-  a matter of f a c t , the word l i q u e f a c t i o n broad  sense  throughout  this  thesis.  of the s i m i l a r i t i e s of these behaviors  will  follow i n the coming chapters.  1.4  PURPOSE AND In  subjected  SCOPE OF RESEARCH  the  to  liquefaction  c y c l i c loading  The volume  causes  the  described above, and the  pore  compaction water  dissipation  (consolidation)  decreases  comprehensive material model has  been  t h i s coupled problem of l i q u e f a c t i o n . and  of  a  soil  deposit  v i b r a t i o n there are two processes which work i n  opposite d i r e c t i o n s .  pressure  process  Seed(1975)  studied  tendency  under  pressure to r i s e , as  of  excess it.  presented  pore As  water  yet, to  no  explain  Recently, Martin, Finn  the mechanics of pore water pressure  11  generation relation  of a s o i l sample subjected between  pore water research,  shear  pressure a  model  was for  simple shear conditions has been relationship  that  includes  r e l a t i o n , volume compaction, dissipation.  Using  this  established.  liquefaction  constructed a  thick  horizontal  Liquefaction  may  deposit then  water  proposed  this  on  rise  comprehensive constructed  subjected  investigated  to  and  material  to  simulate  behavior  as  this  stress-strain  pressure  consolidation)  when  be  based  hysteretlc  pore  (including  In  granular s o i l under c y c l i c  model, a global mechanical model i s the  a  s t r a i n cycles, volume compaction and  increment  material  to c y c l i c loading and  base  of  a  motion.  the cummulatlve  consequences of the following i n t e r a c t i n g processes: 1. Dynamic response to base motion, 2. Increase i n volumetric  compaction s t r a i n  due  to  cyclic  straining, 3. Pore water pressure generation  as a r e s u l t of ( 2 ) , and,  4. Dissipation of excess pore water pressure (consolidation). Numerical techniques  for  the  solution  of the  liquefaction  problem are developed and an example problem i s solved showing that the model i s capable of simulating the of  a  soil  deposit  as  well  as  giving  permeability or c o e f f i c i e n t of consolidation potential.  actual the on  behaviour  influence  of  liquefaction  12 CHAPTER 2 . REVIEW OF PAST INVESTIGATIONS ON LIQUEFACTION OF SATURATED SANDS  2.1  LIQUEFACTION  POTENTIAL  OF  A  DEPOSIT  BASED  ON  FIELD  EXPERIENCE. As Earthquake,  early  1906,  as  i t was  after  observed  the  the  saturated  alluvial  sand  that  flats  aggregates  the main constituent of the s o i l deposits. at  Francisco  that the places that were most  susceptible to s o i l l i q u e f a c t i o n were where loose to very loose  San  I t was  were  recognised  time that these loose sands were responsible f o r the  liquefaction  of  the  soil  earthflows.  However,  descriptions of  the  soils  deposits apart  as  from  that  were  well  as  these most  f o r many qualitative  susceptible  to  l i q u e f a c t i o n , the values of the s o i l parameters, i f ever known at that time, were not reported. Florin field They  and  investigations studied  the  Ivanov(196l) into  the  liquefaction  published  soil  some  of t h e i r  liquefaction  problem.  behavior of r e l a t i v e l y loose  sandy s t r a t a which were around eight to ten metres (around feet)  i n depth  vibration.  using  explosives  as  an  30  energy source f o r  Throughout the i n v e s t i g a t i o n they used  a  charge  consisting of 5 kilograms of ammonite buried at a depth of  4.5  13 meters  (about  15  feet).  They observed that a f t e r a b l a s t ,  the loose saturated sandy s o i l deposit surface  settlements  after  liquefied  with  large  re-consolidation whereas i n dense  sands no evidence of l i q u e f a c t i o n could be found and l i t t l e or no surface settlement was observed. loose  Also, i n the case  of a  sandy deposit l i q u e f a c t i o n always occurred within a few  seconds of f i r i n g the explosives. experiments  they  As  a  result  of  these  were able to e s t a b l i s h a c r i t e r i o n f o r s o i l  l i q u e f a c t i o n using the magnitude of settlement due to b l a s t i n g as a key v a r i a b l e , namely:1. I f the average settlements  (due to the blast  created  by  the explosion of the 5 kg of ammonite described above) i n a  radius  of 5 meters (16  (around 3.5  inches), no  feet) was less than 8 to 10 cm  provision  against  liquefaction  should be made, and, 2. I f the r a t i o of settlements between successive blasts was greater than  1:0.6  l i q u e f a c t i o n of the stratum would be  very l i k e l y to occur. It was also reported that t h i s blasting  provided  quite  suffers the following  reliable  results.  of  standard  However, i t  disadvantages:-  1. Blasting has to be done investigation  method  which  may  on  the p a r t i c u l a r  site  under  be quite objectionable i n some  cases, 2. Uncertainty involved when applying the c r i t e r i o n to s o i l s with  properties  differ  from  those  described  i n the  14 literature, 3. Magnitude  and, and duration of v i b r a t i o n a l disturbance cannot  be taken into account. Koizumi(I966) examined the boring data f o r the conditions  in  the  Niigata  area  and a f t e r the 1964  before  Niigata earthquake i n an e f f o r t to study the changes density  due to the earthquake.  some  in  soil  From standard penetration N-  values f o r about 20 locations i n the N i i g a t a area that  soil  he  noticed  loose sand s t r a t a were compacted by the earthquake  while other not-so-loose sand s t r a t a were loosened.  Based on  these observations he reasoned that f o r the sandy s o i l at Niigata area a c r i t i c a l void r a t i o as a function of pressure  should e x i s t .  A  curve  was  proposed  c r i t i c a l standard penetration values, N , cr  overburden to  relate  which were used as  a measure f o r the c r i t i c a l void r a t i o at various depths, depth  of overburden  as shown i n Figure 2-1.  N__  deposits  with  N-values  above  Niigata.  given  or  below  However, the  loosening or compaction would occur i f the sand subjected  to  the  same quake again.  volume contraction occurred i n an  earthquake  the  be  the  result  N  c r  for  values,  deposit  were  Whenever compaction or saturated  sandy  deposit  pore water pressure would r i s e causing  the e f f e c t i v e stress i n the s o i l to would  as  would not s u f f e r any change i n density when subjected  to an earthquake s i m i l a r to that of  during  with  Thus, according  to Koizumi, a sand deposit having the same N values by  the  decrease.  Liquefaction  when the e f f e c t i v e stress approached a  N-Value 0  10  20  30  UO  50  •  -N  ig.2-1 Critical  N-Value  c r  versus Depth  (after Koizumi, 1966)  16 zero o r near zero v a l u e . also  be  this  the N  Consequently,  curve  cr  could  viewed as a c r i t e r i o n f o r s o i l l i q u e f a c t i o n .  criterion,  liquefaction  Osaki(1968) ' determined  the  Using  zone  of  of" sand d e p o s i t s i n h i s study i n o r d e r t o r e l a t e  the depth and t h i c k n e s s o f l i q u e f i e d s o i l s t r a t a t o the amount o f damage t o s u p e r - s t r u c t u r e s . a deposit  and  the  N  cr  The a c t u a l N-value curve  curve were p l o t t e d  i n the same graph,  and t h e zone o f l i q u e f a c t i o n was o b t a i n e d from the intersection,  the c h a r a c t e r i s t i c s at  a  general  potential  of  of  p u b l i s h e d s e v e r a l papers on  o f l i q u e f i e d sands i n an attempt t o a r r i v e  criterion a  points  2-2.  a and b, as shown i n F i g u r e  Kishida(1964,1969,l970)  for  level  f o r accessing deposit  f o l l o w i n g items determined  and  the  he  the occurrence  liquefaction  suggested  t h a t the  of l i q u e f a c t i o n of a  d e p o s i t :1. Earthquake i n t e n s i t y around V  to  VI  according  t o the  Japanese M e t e o r o l o g i c a l Agency I n t e n s i t y S c a l e . 2.  E f f e c t i v e overburden p r e s s u r e l e s s than 2 kg/sq.cm.  3. R e l a t i v e d e n s i t y l e s s than 15%, 4.  Saturated.  5. Coarse-grained s o i l w i t h 0.074<D^<2.0 mm. 6. S o i l g r a i n s a r e f a i r l y uniform w i t h U.<10. Q  From was o b t a i n e d .  item number 3 a curve o f N-values versus The curve, which has a s i m i l a r  N  curve,  Nrty.  curve which was o b t a i n e d by Koizumi i s a l s o  cr  is  shown i n F i g u r e 2-3.  usage  as  depth the  In the same f i g u r e the plotted  for  17 N-Value  Fig.2-2  Zone  of  Liquefaction  (after Osaki, 1968)  18 comparison. It i s quite obvious that these curves are applicable to  the  Niigata  Use of these curves to predict the l i q u e f a c t i o n deposit  located  elsewhere  would  potential  the  Nevertheless,  preliminary  these c r i t e r i a are  stage  of  site  LABORATORY INVESTIGATION  OF  quite  useful  i n v e s t i g a t i o n so that a  d i r e c t i o n f o r d e t a i l e d s i t e i n v e s t i g a t i o n may  2.2  relating  the s o i l properties, as well as those r e l a t i n g to the base  rock motion. in  of  be quite dubious since the  dynamic stresses developed depend on many parameters to  1964.  deposit subjected to the earthquake of  be  decided.  LIQUEFACTION  USING  SCALED  SAND MODELS. As dam Two  early as  1936  Casagrande(l936)  sections to demonstrate the model dam  sections, one  stability  Water  side  sandy  slopes.  tanks  mounted  l e v e l on the upstream side of the dam  was kept at a small distance below downstream  of  loose sand and one of dense, were  b u i l t of ordinary beach sand and placed i n casters.  constructed model  the  base was  the  crest  while  model  on  covered with a pervious  The  tanks  were shaken and i t was  the  filter  blanket to keep the l i n e of seepage away from the face of slope.  on  the  observed that the  loose sand model l i q u e f i e d ( f a i l u r e of flow s l i d e  type)  readily  However, no  while  the dense model remained i n t a c t .  d i r e c t l y useful numerical data were obtained testing.  from  the  very  model  N-Value 0  10  2.0  \/  ^  N  /  /  U0  — - Low Possib ility  /\  Y  30  \  \ y / / \  High Possibil ty  \>  /  s \  \ N. / / . «—s\ " \\ X / v // ^ /  ^\  N  Relative Density  'X / X  ^ s  dium \ P Dssibility'  Fig. 2-3  Probability (after  of  Liquefaction  Kishida, 1969)  20  Maslov(1957)  and  F l o r i n and  Ivanov(l96l)  model t e s t i n g , as well as f i e l d t e s t i n g , f o r the of  saturated  sandy  s o i l masses.  performed  liquefaction  Complete d e t a i l s of t h e i r  models and experimental procedures have not been published that accessment of boundary e f f e c t s cannot be made. their  already  generated  by, seismic  Maslov suggested the existence °crit»  w  n  o  s  e  of a  had been  disturbances.  critical  acceleration,  value depended on the s o i l density and type, and  i f the maximum  acceleration  of  of c *  would  or  vibrational  be  (re-consolidation)  concerned  process  established  a  criterion  water  a  comparison  pressure  between  obtained  section,  for l i q u e f a c t i o n based on the  amount of settlement due to b l a s t i n g at the s i t e . made  with the  of a l i q u e f i e d  s o i l stratum, though, as discussed i n the previous they  generated.  for various s o i l s i s shown i n Figure 2-4.  c r i t  The work of F l o r i n and Ivanov was mainly consolidation  the  f o r the p a r t i c u l a r  s o i l mass no excess pore water pressure graph  that  (vibrational)  disturbance was less than the required  A  However,  e f f o r t s seemed to have been focused on the d i s t r i b u t i o n  and d i s s i p a t i o n of excess pore water pressure  that  so  They  also  the d i s t r i b u t i o n of excess pore  theoretically  and that  obtained  through model t e s t i n g and seemingly good agreement was arrived at.  However,  neither  of these investigators attempted to  access the amount of excess  water  pressure  generated  function of the s o i l parameters and the acceleration duration of the v i b r a t i o n a l disturbance.  as a  l e v e l and  21  POROSITY (Vo)  Fig. 2-4  Critical A c c e l e r a t i o n Sands.  (After  for  Various  Masloy. 1957)  22 Yoshlml(1967) placed loose saturated sand i n a r i g i d  50cm  box,  25cm x 25cm  x  i n dimension, and applied h o r i z o n t a l  v i b r a t i o n to the box i n an e f f o r t e f f e c t on sand deposits. put over  the  sand  to  simulate  the  An impervious f l e x i b l e membrane was  i n the box so that a surcharge  applied by means of a i r pressure on the membrane.  two  stages,  pressure rose  an  initial  in a  liquefaction  stage  compaction  near-hyperbolic when  failure  consisted  stage i n which pore  fashion  occur  could be  Using such  a set up, he observed that the l i q u e f a c t i o n process of  seismic  and  in a  a  sudden  short  time.  Although quantitative r e s u l t s were not d i r e c t l y applicable to field  problems  conditions, describe  because  qualitative  the behavior  of  the difference  extrapolations  i n boundary  could  of a sand stratum  be  made  to  located at depth and  sandwiched between impervious l a y e r s . Recently,  Finn  saturated  et  al(1970,1971)  sands  performed  on  samples.  The size of sample used was approximately  by 18 Inches by 7 inches and input  was  along  the  using  reported  relatively  the d i r e c t i o n  longest  side  of  minimize, i f not eliminate, the boundary the  sample  by  the container.  of  tests  large sand 72 inches  base  motion  the sample so as to effects  Imposed  An electro-mechanical  on  servo-  controlled shake table was used as a source f o r imparting base motion to the sample. table  i s that  The advantage of using  such  a  shake  various c o n t r o l l e d input formats can be used.  The input can be e i t h e r a c c e l e r a t i o n , v e l o c i t y or displacement  23 controlled and the wave  form  sinusoidal,  or  records. that  triangular  of  the  square  input  or  can  be  irregular  Test equipments and procedures  were  sample uniformity was within reasonable  earthquake  developed  sinusoidal  base  input  the  ratio  acceleration amplitude to that of the base points  so  tolerance of +2%  Lee(1973).  v a r i a t i o n i n void r a t i o , see e.g., Emery, Finn and For  either  of  the  input  response  at  various  within the sample was also measured i n several samples  and was shown to be close to unity when the frequency of input was l e s s than 10 hz. constant  amplitude  Numerical  results  were  obtained  for  sinusoidal base motion input r e l a t i n g the  amplitude of acceleration to the number of cycles of v i b r a t i o n to cause l i q u e f a c t i o n , N , as shown i n Figure 2-5.  E f f e c t of  T  i n i t i a l e f f e c t i v e normal (or surcharge) pressure,  o-', on  N  V  Jj  was also examined and the r e s u l t i s as shown i n Figure The  general conclusion  that  can  2-6.  be obtained from  these model tests can be summarized as follows :1. High degree of saturation of the sand model i s required. 2. The looser the sand the  shorter  the  time  (or  smaller  number of cycles) to l i q u e f a c t i o n . 3. The  higher  the  i n i t i a l normal  e f f e c t i v e confining  stress)  the  pressure longer  (or the  Initial time  (or  larger number of cycles) to l i q u e f a c t i o n . 4. The (or  larger  the  acceleration input the shorter the time  smaller number of cycles) to l i q u e f a c t i o n . It should be pointed out that these conclusions agree  0.5  0.U  h  o  I  z  0.3  o  I—  <  DC LU  _l  0.2  LU O O <  0.1  0.0  10  Fig. 2-5.  1000  100  N, Effect to  CYCLES  of  TO  LIQUEFACTION  Acceleration  Liquefaction,  . I 10000  Level  (after Finn  on  Resistance  et al, 1971)  8.0  ~  6.0  Q.  LU O  <  U.0  o or 3  CO  2D  10000  10  N Fig. 2-6.  Effect to  CYCLES  L  of  TO  Surcharge  Liquefaction.  LIQUEFACTION  Pressure  (After  Finn  on  Resistance  et al, 1971)  ro  26 extremely  well with f i e l d observation.  In addition to t h i s ,  model t e s t i n g provides a basis f o r evaluating the e f f e c t of single variable on the l i q u e f a c t i o n p o t e n t i a l while from observations  the  effect  of  a  single  variable  is  a  field seldom  presented i n a clear-cut manner. In passing, mention should be made that i n most, not  a l l , model  t e s t i n g the r e s u l t s obtained are useful as a  stepping stone f o r p r e d i c t i n g f i e l d behavior. of material (or s o i l ) behavior predict  the  if  behavior  is  I f any theory  proposed  to  describe  or  of an a c t u a l structure (including s o i l  deposit), the theory should be able to predict  the  behaviour  of the model with better accuracy since, 1. the loading on the model can be controlled and measured, 2. the boundary conditions of the model are known, and, 3. the material  in  the  model can be made more homogeneous  and i s o t r o p i c than any f i e l d deposit i s l i k e l y to be. In  order  liquefaction  to  behavior  arrive of  at  some  sandy m a t e r i a l ,  samples under presumably uniform stress or should not  be  'theory*  overlooked.  testing strain  for  the  of small conditions  These are described i n the next  section.  2.3  LABORATORY INVESTIGATION OF LIQUEFACTION USING SMALL SAND SAMPLES During an earthquake  a typical  soil  element  in  a  27 deposit  Is  deformed  by  a  system  of dynamic stresses as a  r e s u l t of earthquake waves passing through The  most  significant  the  soil  medium.  set of dynamic stresses are considered  generally to be'caused by the propagation of shear waves the  bedrock  to  the  ground surface.  For a h o r i z o n t a l s o i l  deposit on a horizontal bedrock formation system  can  t h i s dynamic  be i d e a l i z e d as shown i n Figure 2-7.  the s o i l element i s subjected  to  stresses, cr^ and  K <r^ , where  normal stress  K  0  0  and  the  'at r e s t *  principal  cr^ i s the i n i t i a l e f f e c t i v e  0  Q  i s the c o e f f i c i e n t of earth pressure at  Q  During the earthquake the dynamic simple shear  system  as  developed  stress  i n Figure 2-7b i s imposed r e s u l t i n g i n the  stress system as shown i n Figure 2-7c. been  stress  Initially,  rest.  shown  from  to  apply  Test procedures  have  a uniform stress system s i m i l a r or  equivalent to that of Figure 2-7c and three o f the most useful and commonly employed ones are :1. C y c l i c Loading T r i a x i a l Test - This i s quite the  s i m i l a r to  conventional t r i a x i a l t e s t though a c y c l i n g deviator  load i s used load.  instead  of  the  unl-directional  The c y l i n d r i c a l saturated s o i l specimen i s f i r s t  consolidated under an ambient pressure of further  deviator  drainage  from  the  the  between  P  0  ^ D a:  a  n  d  °o~ D A0  then, with  %  sample prevented, a c y c l i n g  a x i a l load i s applied so that <r +/  <T o  e  r  i  axial o  d  i  c  a  l  stress  l • y  F  o  r  varies f u l  ly  saturated s o i l specimens the e f f e c t i v e stress applied can be represented should  as shown  i n Figure  2-8a  and  b.  It  be noted that the actual f i e l d conditions are not  a.  Fig. 2-7.  Idealised  b.  Field  Loading  c.  Conditions. CO  29  Fig. 2 - 8  Idealised  Triaxial  Conditions.  Test  30 exactly duplicated  since  the  actual  principal  stress  system rotates approximately +40° p e r i o d i c a l l y while that of  the  laboratory  test  conditions fluctuates by +90°.  The i n i t i a l stress systems are also d i f f e r e n t ; element  the  soil  i n the f i e l d i s K -consolidated while the sample o  i n the t r i a x i a l apparatus i s i s o t r o p i c a l l y  consolidated.  2. C y c l i c Simple Shear Test - In t h i s test a s o i l sample rectangular  cross-section  with  sides.  rigid  dimensionally  is  It  then  one-  t e s t i n g drainage i s prevented and a c y c l i n g shear  stress  to  the  vertical  consolidated  During  applied  a  confined i n a 'shear box'  stress of <T^.  is  under  is  of  top  and  0  bottom  of  the  sample.  Moreover, two sides of the 'shear box' rotates i n such way  so  that  the  uniform stress or s t r a i n condition as  shown i n Figure 2-7c stress  i s duplicated c l o s e l y .  conditions  in  the  simple  resembles those acting on a s o i l test  data  a  obtained  can  be  shear  element  applied  Since test  in  the  closely  the  field  d i r e c t l y to f i e l d  conditions without any c o r r e c t i o n . 3. C y c l i c Torsion Test - The set up test  is  similar  of  equipment  in  this  to that of c y c l i c t r i a x i a l t e s t .  But  the c y c l i c shear stresses are applied through the use torsional  load  applied  to both ends of the c y l i n d r i c a l  sample while the a x i a l load i s kept constant. and  Li(1972)  used  this  method  of  shear  test  apparatus  that  Ishihara  applying  deformation to the samples and they designed torsion  of  a  shear  triaxial  can be used to run  31  tests with constant confining pressure as well with  no  lateral strain.  the * second  test.  type  the  cyclic  triaxial  two  test  corresponds to the simple shear  However, more c o r r e l a t i o n of t h i s test  other  tests  It was claimed that the f i r s t  type of tests corresponds t c while  as  with  the  i s needed before the test data can be applied  to f i e l d conditions since the usage of average stress and s t r a i n to represent the varying stress and s t r a i n Inside  the c y l i n d r i c a l specimen under t o r s i o n needs more  substantiation. solved  the  an  of  1  field  Ishibashi  and  1974)  Sheriff  varying stress and s t r a i n f i e l d  cylindrical  ingenious  'donut-like strain  Recently  problem  Inside a s o l i d using  fields  specimen  under  sample c r o s s - s e c t i o n .  sample i n order to achieve a  torsion  by  They used a  uniform  shear  and consequently they were able to compare  t h e i r test r e s u l t s with those obtained from the  triaxial  and the simple shear apparatus based on the actual stress and s t r a i n values. Investigations of  sand l i q u e f a c t i o n using the above  mentioned types of laboratory test have been made by Seed  and  Lee(1966), Peacock and Seed(1968), Finn et al(1971), Ishlhara and Ll(1972) and Ishibashi and Sherif(1974). Most of the tests  performed  involved  c y c l i c shear stress, x^  t  during  the  test  the use  as shown  of in  a  constant  Figure  2-9.  amplitude Usually  the pore water pressure, cr , c y c l i c s t r a i n ,  7 t and the c y c l i c stress, T , were continuously monitored  by  32 using  electrical  transducers  and chart recroders.  records of simple shear and t r i a x i a l t e s t s  are  Figures 2-10  and 2-11.  portion  the test where the pore pressure  of  as  Typical shown  It can be seen that during the  the amplitude of c y c l i c s t r a i n developed i s  in  initial  r i s e s gradually, almost  zero  (of  the order decimal of a percent) but becomes increasingly large (+15  magnitude  to 20^)  towards l i q u e f a c t i o n f a i l u r e when the  pore pressure r i s e s more r a p i d l y .  It  sands  that  behaved  in  this  way  and  l i q u e f a c t i o n could be described  in  was  a  found  their  that  most  resistance to  simplified  manner  by  using the following variables :N:  number of stress cycles to cause l i q u e f a c t i o n .  Dt  r e l a t i v e density, defined as  L  r  D  =  where max  e  ** * " max a  r  :-  x 100%  6  (2-1)  min  e i s the current void r a t i o of the sand sample and  a n d  e  min  a r e  t n e  v  o  i  d  ratios  densest state r e s p e c t i v e l y . used  to  indicate  the  at  the  loosest  This variable i s  state  of  and  generally  denseness of the sand  specimen. td  !  c  y°ll  c  s  n  e  a  stress amplitude,  r  and,  C"y : the i n i t i a l e f f e c t i v e v e r t i c a l (or confining) s t r e s s . Q  It was and  the  combined  found that the c y c l i c stress amplitude,  initial and  effective  (T /crJ ) A  rs  can  vertical be  (j-' . vo as a  stress,  considered  T » d  can  be  single  33  Fig. 2 - 9  Cyclic Shear Wave  Forms  Stress  TEST 58 e = 0.617 2.0n  cc  «f nr f. r n f. r=Ml Hi h n I I i i 11 II I;II :.| i l l M i l l : i i l i I i.l.uhll • i.i'i t.111y11U11' 1  UJ  <  rfflllhl  —  UJ  nflWi'r¥rjr»r r rr i r  e w  ! PORE  i.oH  *  cc ~"  PRESSURE  i ri M i i i ii i i i i i i i i i ii i • i i i i i i i ii i i i i i i ii I i i i i i 11 I M i lr  iff  1  2 UJ=>  1  ce  ni  3  Z (O < UJ X  o-  SHEAR  1  cc  1 5  -10  STRAIN  x I z  <  — cc sS cc •10 <  SHEAR cc e  0.4-  <o  £ 2>  ii  0.2-  »  n fi ii  05 —X  0^ o p z < P I 0.2  £ S 0.4 b£  < c/)  . 2-10  II II II  II  II  M  II II II  ii  III II III II III II  STRESS  M  A II  II l| n M II M I ! ii ii ii Ii ii ii ii ii iii ii ii H II II II II II II II n in II it II  i  II  ii i II  il  II II II  v y  II II II II  H  II II  vhi  ii ii H i H n n II i II II n ii ,i M ii  II II II  H  II  vHy  A  A A II II II II  n n n  II II II  L-15 ni A.. Ml  III  ii :i ii i! H II n II II II  ii  H H  II  n  II II  II  III  H i n n ii  II  II II II II  i i i II n  A  ft  I!  II II  u uuuv uy nnJ U U  SHEAR  II  ii  FRICTION  Typical Cyclic  II  H  II  Ii II II  it  I 11 .,  II  H  II  H H  yv u  w  II II A:  fl II I  UJ X  V  i ii  v  y  =0.03 K g / c m '  Simple Shear Test Data  (after Finn et al, 1971)  36 Independent to present  variable.  Finn et al(197l) used a summary plot  the t e s t r e s u l t s r e l a t i n g N_ , D and L r  ) in  (T./(J'  d  a  vo  more e f f i c i e n t way as shown in' Figure 2-12, where each curve represents a constant (T,/o" ) value. It can be seen from o. vo f  this  figure  by N ,  that the resistance to l i q u e f a c t i o n , as measured  decreases as  T  :-  1. the r e l a t i v e density of the sample becomes smaller, 2. the c y c l i c shear amplitude becomes larger, and, 3 . the i n i t i a l  effective  vertical  (or  confining)  stress  becomes smaller. Although  the  c y c l i c t r i a x i a l t e s t imposes a stress  condition on the sand sample d i f f e r e n t from that i n the Finn et  al(l971)  found that i f due consideration i s  field,  given  to  i n i t i a l e f f e c t i v e normal pressure acting on the 45°-plane  the  in the t r i a x i a l t e s t specimen c o r r e l a t i o n between simple  shear  test  and  triaxial  test  is  the  cyclic  possible.  In  addition, the recent development of the novel Torsional Simple Shear Device by Ishibashi and Sheriff 1974)  makes  to  degree  run  liquefaction  confinement. AT^/<7^ ,  tests with varying  Moreover,  they  found  that  the  the octahedral normal stress, versus the number of liquefaction  al(1971).  since  it  of  possible lateral plot  of  the r a t i o of maximum change i n the shear stress to  ct  et  it  cycles  to  bears close resemblance to that obtained by Finn This c o r r e l a t i o n has a n o n - t r i v i a l  implication  indicates that the l i q u e f a c t i o n of sand produced i n  the laboratory i s not an a r t i f i c i a l ,  or  apparatus-dependent.  37  (07') = 2.0 0  X-  w  Kg/cm  Marked  1  for Edch  Line  .75 .21 I5  a> I  V\ V  .70  \.34  V  I0  \  .65  o <  y.50  \  .60  o > .55  .50  10 CYCLES  Fig. 2 - 1 2 Simple  100  TO  Resistance Shear  1000  LIQUEFACTION  to  10000  -T1,  Liquefaction  in  (After Finn et al, 1971)  38 phenomenon  and the data obtained i s good f o r a p p l i c a t i o n to  the analysis  of actual  soil  structures.  The method of  adapting l i q u e f a c t i o n data obtained In the laboratory to f i e l d problem i s discussed i n the following s e c t i o n .  2.4  SOME  METHODS  OF DETERMINING LIQUEFACTION POTENTIAL OF A  SOIL DEPOSIT. In sections 2.1 and 2.2 of t h i s chapter, some of the methods of determining l i q u e f a c t i o n p o t e n t i a l were  discussed,  for example :-  1.  Florin  and Ivanov(l96l) -  using  settlement  due to  blasting at the s i t e as a key v a r i a b l e ,  2. Osaki(l968)  - using  the c r i t i c a l  standard  penetration  values, N , as a c r i t e r i o n , c r  3. Kishlda(1969)  -  using  the  75%  relative  density as a  c r i t e r i o n f o r the occurrence of l i q u e f a c t i o n , and,  4. Maslov(l957) level,  c^Q-pif  - using the concept of c r i t i c a l acceleration to determine the p r o b a b i l i t y  of having  a  liquefaction, and  they  are not repeated  concerned with the a p p l i c a t i o n  here.  The following i s mainly  of laboratory  test  data to  f i e l d problems.  Castro(I969)  following  the  suggestion  of  Prof. A. Casagrande used the concept of c r i t i c a l void r a t i o to determine the l i q u e f a c t i o n p o t e n t i a l series  of stress-controlled  of a  sandy  soil.  A  t r i a x i a l undrained t e s t s , which  .3  were termed R t e s t s , were performed on representative from  the  field  for  various  39  samples  void r a t i o or r e l a t i v e density  values and at d i f f e r e n t i n i t i a l consolidation pressures, c r ^ . For  samples  that  developed  continuous pore pressure r i s e at  constant a x i a l load, (SPONTANEOUS LIQUEFACTION), the values of the  effective  against  minor  principal  2-13.  Figure a  ^ f *  the void r a t i o of the sample, e^,  through the points and a e~p l i n e  a^  stress,  was  w  e  r  A curve was drawn  obtained  as shown  For a s o i l element with consolidation pressure 2-13,  have a higher l i q u e f a c t i o n p o t e n t i a l than an element of  the same s o i l with <r^ and e^, represented by point the h o r i z o n t a l distance  AC  reduction  minor  pressure f  In  and void r a t i o e^, represented by point A i n Figure  would  e^  plotted  e  in  effective  and  BC  indicate  principal  r i s e , at l i q u e f a c t i o n .  B,  since  the amount of  stress,  or  pore  A s o i l element with (r*^ and  represented by point D, would not develop continuous  pressure  pore  r i s e at constant loading and i s said to be d i l a t i v e .  However, i n t h i s method of determining l i q u e f a c t i o n  potential  there are a number of factors not being taken into account  and  two of the important ones are :1.  Magnitude  and  duration  of the  disturbance  that cause  l i q u e f a c t i o n , and, 2.  E f f e c t of a further stress reversal on  the  sample  that  exhibit d i l a t i v e behavior. Nevertheless,  the  method i s a good guideline f o r recognizing  deposits with metastable spontaneous l i q u e f a c t i o n .  sand aggregates which  would  develop  0*7  41  Seed determining  Idriss(1967.1971)  and  the  liquefaction  s u b j e c t e d t o earthquake horizontal  shear  proposed  potential  forces.  of  a method f o r soil  deposit  The time h i s t o r y response o f  stresses,  T  ,  at  various  depths  was  xy e v a l u a t e d u s i n g dynamic a n a l y s i s .  Based  on  the  irregular 10  time h i s t o r y s t r e s s response the amplitude o f an e q u i v a l e n t cycles  o f constant c y c l i c shear s t r e s s was  each depth so t h a t a p l o t o f e q u i v a l e n t  then determined a t  10-cycle shear  stress  amplitude, T ^ Q ,  v e r s u s depth was  2-14,  From the l a b o r a t o r y c y c l i c simple shear t e s t  curve AB.  o b t a i n e d as shown i n F i g u r e  data, o r c y c l i c t r i a x i a l t e s t d a t a m u l t i p l i e d by a c o r r e l a t i o n f a c t o r , the shear s t r e s s amplitude t h a t caused l i q u e f a c t i o n i n 10  cycles, T^Q»  was  p a r t i c u l a r values of D curve  CD  determined  as shown i n F i g u r e 2-14  To  w i t h i n which  simplify  potential  the  Seed  T  d|  >T 0  was  is  no» »S«» e  procedure and  each  depth  using  of  Idriss(197l)  obtained.  determined z a  The zone o f  by  the  t o z^ i n F i g u r e  determining proposed  that  cyclic  shear  2-14.  instead of 10  s t r e s s , T,, , dlO approximation was used t o determine the e q u i v a l e n t N , o f constant c y c l i c shear o f amplitude, x . eq av the  constant amplitude  depth  liquefaction  approximating the s t r e s s response by an e q u i v a l e n t of  the  and cr^ a t the p a r t i c u l a r depth and the  r  l i q u e f a c t i o n , i f t h e r e i s one, interval  at  rt  cycles a direct cycles , However,  b a s i c method of d e t e r m i n i n g l i q u e f a c t i o n p o t e n t i a l was  not  changed. At  f i r s t glance,  it  would  seem  that  conclusions  42  Stress  Depth  FIG.Z-14  Method of Evaluating Liquefaction Potential (After Seed et a l , 1971)  43 about  l i q u e f a c t i o n p o t e n t i a l based on the c r i t i c a l void r a t i o  concept contradicts that using c y c l i c test to  data.  According  the c r i t i c a l void r a t i o concept an increase i n the i n i t i a l  consolidation  pressure  would  according  increase  while  reduced.  However, closer examination reveals that  in  the  two  cases.  the  latter  liquefaction  potential  'liquefaction p o t e n t i a l ' was  to  the  the p o t e n t i a l i s the  term  used to convey d i f f e r e n t meanings  In  actual  fact when  the  initial  consolidation pressure i s increased, a. the  critical  void  ratio  concept  Indicates  a  higher  p o t e n t i a l of f a i l u r e i n the spontaneous l i q u e f a c t i o n mode but  higher  failure, b. the  stresses  are  needed  to  initiate  a  while,  cyclic  magnitude  test of  data  indicates  disturbance  more  that  for  cycles  the  (smaller  liquefaction  same  of stresses are  needed to cause the pore pressure to r i s e to the level  such  failure  potential) though that the  d i l a t i v e response which has a s t a b i l i z i n g e f f e c t  may  be  reduced i s not mentioned. It  can  be  seen  that  while  the c r i t i c a l void r a t i o method  emphasizes the c l a s s i f i c a t i o n of f a i l u r e contractive  dilative  or  behaviour at f a i l u r e , the c y c l i c stress method i s  mainly concerned disturbing  modes,  with  stresses  the  magnitude  to  cause  and  duration  failure.  In  of  the  fact,  the  combination of the two methods would give the complete picture of l i q u e f a c t i o n f a i l u r e of s o i l structures when earthquake loading.  subjected  to  Even  though  the method of determining l i q u e f a c t i o n  p o t e n t i a l using c y c l i c test  data  i s able  to  provide  some  approximate solutions which are useful i n the design of actual of  actual  structures,  there  are s t i l l some questions l e f t  unanswered by t h i s method, e.g.:1. pvnamlo Response - In a r r i v i n g at values f o r N <' a  dynamic  analysis  using  equivalent  damping c o e f f i c i e n t s are used. The but  non-linearity the change  Seed  and  T  shear moduli and  and  Idrlss(l969)•  i n s o i l s t r e s s - s t r a i n i s approximated  of  shear  e f f e c t i v e normal stress and  modulus  due  to  changes  in  compaction (hardening) i s not  taken into account, 2. Time to Liquefaction  -  Since the time history  of  pore  pressure response cannot be determined by the method the time to l i q u e f a c t i o n cannot be determined, 3. Maximum  Strain  at  Liquefaction  -  As  a  result  of  approximating the actual s t r e s s - s t r a i n by a straight l i n e the one  value of maximum s t r a i n obtained i s not a and  to  realistic  obtain such s t r a i n value the actual s t r e s s -  s t r a i n should be used, and, 4. The Effect of Consolidation during pore pressure r i s e the  liquefaction  potential  of  a  soil  has not  on been  considered. It i s the purpose of t h i s thesis to t r y to a r r i v e at a in which account.  the above  mentioned  factors  can be  taken  method into  45  CHAPTER 3 MODELLING OF SATURATED GRANULAR SOIL UNDER CYCLIC SIMPLE SHEAR CONDITIONS  3.1  VOLUME  CHANGE  CHARACTERISTICS  DURING  CYCLIC  DRAINED  SHEAR. It  is  a  commonly  observed  phenomenon that a dry  granular s o i l compacts or reduces i t s volume when subjected cyclic  loading,  earthquake  e.g., machine  forces.  The  foundation general  vibration  observation  the  amount  of  testing  are quite  scattered  and  experiments  to  amount  of  enclose  the  acceleration  s o i l specimen during of  stress  the specimen was not known.  compaction level,  was or  Barkan(l962)  'Rigid' containers were used i n  course of v i b r a t i o n so that the state by  laboratory  uncorrelated and sometimes  Mogaml and Kubo(1953)t  and Bazant and Dvorak(1965).  experienced  that  Apart from t h i s , the quantitative  by many early investigators using  contradictory, see e.g.  their  is  the  volume reduction Increases with the magnitude  and duration of v i b r a t i o n . results obtained  and  from  performance of model as well as e x i s t i n g foundations  to  related  to  the  and  the  strain  Invariably, the frequency  and  v e l o c i t y l e v e l i n some cases, of the  input v i b r a t i o n .  It i s not d i f f i c u l t to r e a l i s e that r e s u l t s  obtained would be  very  much  apparatus-dependent  since  the  46  actual  stress  and  strain  experienced  by the s o i l specimen  depends not only on the magnitude of input v i b r a t i o n but on  the  also  dynamic response properties of the container-and-soil  system.  For example, i f one examines the  Schaffner(l962),  as shown i n Figure 3-1.  data  obtained  by  one may observe that  each of the f i n a l void r a t i o curves obtained  for  a  constant  acceleration l e v e l strongly suggests the presence of resonance of  the  container-and-soil  system  at  around  80 to 100  Hz.  Moreover, the e f f e c t of the superimposed v e r t i c a l component of o s c i l l a t i o n , which might be applying data.  vibration,  present  complicated  due  the  to  the  interpretation  Consequently, i t i s more r e a l i s t i c , i n  applying  laboratory  method  the  of  of the  sense  of  test data to the analysis of actual s o i l  structures, to consider the volume change of a s o i l element as a function of the stress and s t r a i n experienced by i t . In the dynamic response of a horizontal s o i l deposit to earthquake base system  motion  the  cyclic  simple  shear  caused by the propagation of shear waves from the base  to the surface of the deposit, as shown i n Figure 2-7. most  significant  investigators,  one.  e.g.,  device  compaction  to and  Silver  cyclic  and  Seed(197D.  Youd(1971#1972), used  determine  recently, Martin, Finn and mechanism  is  the  With t h i s i n mind many of the recent  Dmevich(1972,1972a) and shear  stress  the  relation  stress/strain  Hardin the  between conditions.  Seed(1975) studied the  and  simple volume More  fundamental  of s o i l l i q u e f a c t i o n and a r e l a t i o n s h i p between the  FIG. 3-1  Compaction  of Dry  (After Schaffner,  Coarse  1962)  Sand  48  volume change i n increment  in  proposed.  the  the  drained  undrained  case  and  case  the  during  The construction of a material  pore  pressure  cyclic model  shear to  was  include  both the volume change and pore pressure build-up phenomena i s made  possible  by  such  a  relationship.  Before going into  constructing the equations of the material model, the important  observations  concerning  the  various  volume changes under  drained c y c l i c shear are discussed f i r s t . In t h e i r investigations into the compaction subjected  to  S i l v e r and  Seed(197l)  the  amount  c y c l i c simple shear  of  and Seed and  volume  change  loading,  of  sand  Youd(1971,197), 2  Silver(1972)  (compaction)  observed  that  caused by c y c l i c  simple shear loading on dry granular s o i l depends on, 1. the r e l a t i v e density of the s o i l , 2. the magnitude of the c y c l i c shear s t r a i n , and, 3. the t o t a l number of s t r a i n cycles the s o i l  specimen  has  been subjected t o . Examples of experimental data obtained are as shown i n Figures 3-2  and  vertical particular  3-3.  It  effective test,  should  be  stress, does  noted  which  not  is  that the value of the constant  influence  the  during  volume  c h a r a c t e r i s t i c s s i g n i f i c a n t l y , as shown i n Figure 3-2. i s growing evidence that the compaction  is  caused  by  grains into a more 'stable  volume the  1  change There  change associated  slight  a  with  re-arrangement of the  configuration.  The  amount  of  .550  I  10  IO2  IO3  |0  4  10  Number of Cycles  FIG.3-2  Effect  of  Vertical Stress  on  Compaction (After Youd, 1972)  ^  i — T  "r  i i  i n  r—i"'  i  i i  111  • - 1  1  1  1 1 11 1  Cycle 10 D =45%  /  r  -  - •  /  -  '  8 0 % /  /  i  0.01  i  i  i  i i i t i  0.1  Cyclic  3-3  -  6 0 % /  i  t  i  i i i  1.0  Shear  Strain - 7^  y  (%)  Effect of Relative Density and Shear Strain on Compaction  (After Seed et al, 1972)  i  10.0  51 grain  re-arrangement  i s governed by how much adjacent grains  can move c y c l i c a l l y r e l a t i v e to one another,  and, since the  c y c l i c shear s t r a i n amplitude i s d i r e c t l y related to t h i s type of  relative  volume  displacement the strong dependence of compaction  change  surprising. the  on  cyclic  shear  strain  amplitude  i s not  Furthermore, mention should be made that most of  non-recoverable  volumetric s t r a i n (grain  caused by one-dimensional compression occurs  re-arrangement) i n the  initial  loading stage and that, 1. during  drained  cyclic  shear,  the v e r t i c a l  effective  the v e r t i c a l  effective  stress i s held constant, and, 2. during undrained c y c l i c  shear,  stress decreases, so  that  the effect  of v e r t i c a l  effective  stress  on the  •compaction curves* shown i n Figures 3-2 and 3-3 which are f o r constant c y c l i c s t r a i n amplitudes i s of minor importance. the other hand, the v e r t i c a l shear  modulus  and  so  effective  stress a f f e c t s  the compaction  which  results  This topic i s discussed i n section Seed and  the  has a s i g n i f i c a n t e f f e c t on the shear  strains which result from a given ground motion and on  On  from  therefore  such an e x c i t a t i o n .  3.2.  Silver(19?2) demonstrated the a p p l i c a b i l i t y  of data obtained from drained c y c l i c simple shear tests to the determination of t o t a l settlement of a model dry sand subjected  to harmonic base motion.  deposit  However, f o r non-uniform  52 c y c l i c loading, such as that Imposed by earthquakes, in t h e i r present form as shown i n Figures 3-2 directly that  applicable.  for  a  (relative  Martin,  soil  element  density)  the  compaction s t r a i n , Gyp» parameter  so  that  A£ p, caused  by  y  determined  -  total can  this  is  quite  on  a  vp  =  as  U ( 6  v  not  void  a  ratio  volumetric  strain  history  i n volume compaction,  strain similar  the  are  cycle, f to  theory  ,  &  can  be  the use of s t r a i n  of  plasticity.  An  c u r v e - f i t t i n g method was proposed to  describe the experimental  A e  used  additional  initial  accummulated  a further increment  an  based  or  be  hardening parameter used i n expression  given  and 3-3  Seed(1975) proposed  Finn and  with  test data  data:' V">  O" * 1  a  P  2 "' where  C^,  C  C^.  C^  experimentally.  A  l * a - 2 vt> (  C  €  )  T- "^  +  and  2  are  method  of  (3  constants  the  0.79, a  values  0.45  relative  density,  D, r  of 45$.  the  Figure  nature  process  3-4 0.80,  (in  of  at  It should be pointed out in  addition  the  d i s c r e t e steps) of the implied t o t a l  volume compaction s t r a i n function, e since  in  the  r e s p e c t i v e l y f o r a Crystal S i l i c a Sand  that the expression i s an empirical one and incremental  determined  of the four constants were found to be  and 0.73  la)  g r a p h i c a l l y presenting  function given by equation (3-la) i s as shown and  "  counting  v p  the  , has  to  be  accepted  number of shear s t r a i n  FIG. 3 - 4  Incremental Volumetric Strain (after Martin  et  al, 1974)  Curves  V*)  5^ cycles completed i s also of For  drained  amplitude  cyclic  simple  histories  characteristics  a  are  can  be  discrete shear  the when  volume the  compaction  four  associated with the p a r t i c u l a r s o i l are known. if  constants  For  example,  i s the t o t a l volumetric compaction at the end of the  1 i-th strain  cycle  i  and  and  are  vp volumetric  compaction  and  6  vp  =  €  ^  1  the  incremental  s s t r a i n amplitude  the i - t h s t r a i n cycle then, f o r equation  respectively for  (3-1),  *4p  +  (3-3)  Accordingly the t o t a l volumetric compaction s t r a i n completion such a  of  the  computation  Seed(l975)  and  a  N-th was  STRESS-STRAIN  after  cycle can be c a l c u l a t e d . carried  out  by  Martin,  the  In f a c t , Finn  and  comparison between calculated and measured  values was made as shown i n Figure  3.2  nature.  tests i n which the s t r a i n  monitored computed  incremental  RELATIONSHIP  3-5.  UNDER  CYCLIC  SIMPLE  SHEAR  CONDITIONS. From  a series of s t a t i c t r i a x i a l drained tests on a  sand Kondner and Zelasko(1963) concluded that the r e l a t i o n s h i p between the granular  axial  soil  stress,  under  axial  cr ,  and  axial  compression  represented by a hyperbolic curve,  strain, could  be  £ , of closely  55  to  I oI  £  0  1.2  I  I  5  I  10  Number of  Cycles  1  1  15  I  20  1  1 c  '5  0.8  GO  o 'xl •f  Measured  0.4  E  Computed  O  >  0  i  1  i  5  10  15  Number of  FIG. 3 - 5  20  Cycles  Comparison of Measured and Computed  Volumetric Strain  (after Martin et al, 1974)  56 e  c -  (3-4)  a +b  where a and b are constants determined The  from experimental data.  curve given by equation (3-4) i s shown i n Figure 3-6a and  i t can be seen that the reciprocals  of a  and b  give the  i n i t i a l tangent modulus and maximum stress l e v e l r e s p e c t i v e l y . In terms  of shear  stress, T  , and  shear  following equation can be written, see Figure  T  =  1  max  G  where G  m a x  T  .  +  s t r a i n , T , the 3-6b,  (3-4a)  •*  ^max  i s the i n i t i a l tangent modulus, and, i s the maximum shear stress that can be applied to  the sample, Hardin and Drnevich(1972a) showed that f o r sands  G „ and r max max  can be obtained from the following expressions,  max  =  .(2-973 1 +e  1+K.  -sin  max where G and T max max m o v  e  v  (3^5)  1230  n  a  v  4  - f f i T '.cr;  (3-6)  are i n p s i , '  i s the v o i d  r  ratio,  (T  i s the mean p r i n c i p a l e f f e c t i v e s t r e s s i n p s i ,  0-^.  i s the v e r t i c a l e f f e c t i v e s t r e s s i n p s i ,  K  i s the c o e f f i c i e n t o f e a r t h p r e s s u r e a t r e s t , and,  0  Q  4>'  i s the e f f e c t i v e angle o f s h e a r i n g r e s i s t a n c e .  FIG. 3-6  Hyperbolic  Stress-Strain  Curve.  58 Also, equation (3-^) can be rewritten into,  <fc  Where  G  m  i  •  = T/ /" , and  f  =  =  tabulated f o r various s o i l s .  °-  v  T  m Q  /G  m o  . and values of f  o f the sample  so  were  In addition, they used a hollow  c y l i n d r i c a l sand specimen and applied c y c l i c t o r s i o n ends  7)  that  to  both  a t y p i c a l s o i l element of the  specimen was subjected to a c y c l i c simple shear stress system. From the test data on a large number of s o i l s they to  were  able  show that 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 can be  used f o r s o i l loading ratio  and as  a  Furthermore,  elements subjected expressions function  to c y c l i c  simple  shear  f o r the shear modulus and damping  of  strain  level  were  obtained.  modifications to the i d e a l i z e d hyperbolic s t r e s s -  s t r a i n proposed were made and are b r i e f l y l i s t e d as follows»I , A hyperbolic  shear  strain,  ^  ,  i s defined  instead of equation ( 3 - 7 ) the following  so that  expression i s  used,  i +  G  max with  f.  2. Assumptions  h 'r  h  were  (3-8)  r.  . [ l + ae  r  J  (3-9)  made concerning the damping r a t i o D so  that at any cycle and s t r a i n l e v e l T ,  59 D  (3-10)  max with  defined s i m i l a r l y to (3-9)  a and b involved are d i f f e r e n t equation Numerical by  Hardin  from  those  appeared  in  (3-9).  values and  modifications  though the constants  f o r a's and b*s and D Drnevich(1972a).  the  expressions  are  m Q V  and  However, only  f even  are given with  the  applicable  to the  determination of average shear modulus, G, and damping  ratio,  D, f o r a complete s t r e s s / s t r a i n cycle, as shown i n Figure and 3-8.  Neither the actual s t r e s s - s t r a i n  nor the change i n shear modulus due to volume  curve  is  3-7  traced  compaction  and  change i n e f f e c t i v e stress are taken into account. Martin,  Finn  and  Seed(1975) suggested that during  the l i q u e f a c t i o n process the change i n shear volume  modulus  due  to  compaction and change i n e f f e c t i v e normal stress could  be taken into account by using the following expression,  hv where 1*  hv  f  a +  (3-U)  hf  i s the horizontal shear amplitude, i s the shear s t r a i n amplitude,  (Ty i s the current e f f e c t i v e normal stress, and, a and b are functions of volume compaction €vp*  as follows t —  strain,  Fig. 3-7  Equivalent Shear  Modulus  for  Initial Loading ON O  FIG. 3 ~ 8  Definition of Equivalent Shear Damping  Ratio  Modulus  for Symmetrical Loop  and ON  62 1  A  b = B, A  A , A^, B  1 #  2  experimental data..  (3-13)  l t  +A  €  ^2  (3-13)  W v p  1  with  2 3 vp  B  and  2  are constants determined  Comparing equations  from  (3-11), (3-12)  and  expressions  (3-4), i t can be seen that i n e f f e c t the i n (3-12) and (3-13) give the "hardening' effect  of  compaction  with  volume  -  i . e . , increase i n shear modulus and  maximum shear stress due to volume compaction. section,  a  general  stress-strain  In  relationship  the  next  i s proposed  based on the above described s o i l behavior.  3.3  MODELLING OF SATURATED GRANULAR SOIL UNDER CYCLIC  SIMPLE  SHEAR. 3.3.1  Hysteretic Hyperbolic S t r e s s - s t r a i n Relationship. In  order  to  s t r e s s - s t r a i n behavior cyclic  simple  a  of an  at expressions describing the idealized  shear conditions  Masing(1926) type that  arrive  soil  element  under  i t i s assumed here that the  of hysteretic r e l a t i o n s h i p i s followed  and  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 s i m i l a r to that  proposed by Martin, Finn and  Seed(1975)  i s used.  These are  discussed i n the following paragraphs. For the i n i t i a l loading phase of the  hyperbolic  stress-strain  the granular s o i l  relationship  formulated  by  63 Kondner and Zelasko(1963)  and  Drnevich(1972) i s used.  also  adopted  by  Hardin  and  Therefore, the i n i t i a l response of  the granular s o i l to simple shear  i s given  by the following  equation,  T  mo^  G  = 1  +  (3-14)  mo' mo  where T  i s the shear stress,  "t  i s the shear s t r a i n ,  Gmo^  i s the i n i t i a l maximum tangent modulus and  the maximum shear stress  that  can be  TL« mo i s  applied without  failure. Following  the suggestion  of Hardin and Drnevich x  i t i s assumed that these quantities G and TL are determined ^ mo mo n  by the i n - s i t u condition of the granular s o i l and are given by the following expressions,  G  mo  mo  mo where  =  ,(2.973 - e )  llt76o  1 + e  2  J 1  u  1+K„ 2 1-K ( -y^in*') - (  (3-15)  o  n  2  > 1 *  )  (3-16)  " are i n r psf, mo and lmo * e  i s the void r a t i o ,  <TQ i s the mean p r i n c i p a l e f f e c t i v e stress i n psf, i s the v e r t i c a l e f f e c t i v e stress i n psf, K  Q  i s the c o e f f i c i e n t of earth pressure at r e s t , and,  64  4>' i s the e f f e c t i v e angle of shearing resistance. Since  we  are considering the response  only,  i t w i l l prove more convenient  conditions  in  Then equation  terms  of  (3-15) i s  of to  horizontal layers express the  stress  the v e r t i c a l e f f e c t i v e stress,  a~y»  expressed i n the form,  Therefore, f o r a given i n - s i t u condition,  G  mo  ''•mo  = c*/^  (3-18)  =  (3-19)  pjffi  For unloading and reloading the Masing(1926) type of hysteretic behavior  is  i n i t i a l loading curve,  used. or  This  assumes  that  if  the  skeleton curve, which i s described  by equation (3-1*0 can be represented by,  X  s  =  f (D  "(3-20)  2  Then the unloading or reloading curve can be obtained from,  where (f ,T_„) i s the l a s t reversal point i n the s t r e s s - s t r a i n plot. Consider a s o i l undeformed  state  to  element a  shear  being strain  strained level  from  its  f,  with  a  81  65 .corresponding shear stress l e v e l unloaded  and reloaded i n the  I _.  The  element  i s now  opposite d i r e c t i o n so that the  shear s t r a i n or stress undergoes a complete reversal when the shear  f-^, with corresponding  s t r a i n reaches a f i n a l value of  shear stress, T ^ . Figure  3-9a.  The s t r e s s - s t r a i n  diagram  i s shown i n  (3-1^) while  Curve OA i s given by  curve  AB i s  given by, Tg  ~  = f (~ ) f  (3-22)  T&  2  2  The reversal curve AB i n t e r s e c t s the negative skeleton beyond followed  curve f  a  n  i t i s assumed  a  would  the skeleton  f^'-fa.  at  branch  when  d  that  of the  i s increased  the  stress-strain  path  be the continuation of the negative branch of  curve,  stress-strain  BC,  with  history.  (/, T ) Q  In  other  erased  from  the  when  the  last  words,  reversal point i s on the skeleton curve the absolute value the  stress  or  strain  can be used as a  l i m i t i n g value f o r  .switching back to the skeleton curve, thus 'erasing' (/"_, cl  from our c a l c u l a t i o n . Figure  3-9b.  T ) S a  A general loading history i s shown i n  Assuming that the current loading would bring  the s t r e s s - s t r a i n point, P to P', i t can be s t r e s s - s t r a i n branches, reversal points  of  (f .T a  AB s a  and  BP'A'  ) . (^ .^ ) D  sb  seen  that  are required.  and  ( T .X ) c  s c  two The  have  to  be recorded f o r the determination of the current branch of the stress-strain  curve, CPBP'A'.  (/„,!) C S C  branch CPB,  i s required f o r the  A* (a) first unloading  FIG. 3-9  (b) general reloading  Hysteretic  Characteristics  67 (3-23)  2 2 V  (fo.X^o) i s r e q u i r e d f o r the branch BP'A', a sa  2 While  (^b'tgt))  i  23) t o (3-24). P  passes  (3-24)  =  s  u  s  e  l  d  —  2'  a  s  c r i t e r i o n f o r s w i t c h i n g from  a  Note t h a t as soon as the s t r e s s s t r a i n  the  point  B  In  time.  this  the s t r e s s o r s t r a i n v a l u e o f the l a s t but one r e v e r s a l  p o i n t can be used the  point  the r e v e r s a l p o i n t s B and C would be  e l i m i n a t e d from our c a l c u l a t i o n a t the same case,  (3-  last  as a l i m i t i n g value f o r s w i t c h i n g  back  to  but one u n l o a d i n g or r e l o a d i n g branch and a l s o as a  c r i t e r i o n f o r ' e r a s i n g ' the l a s t and the l a s t but one r e v e r s a l p o i n t s from f u t u r e c a l c u l a t i o n . that  the  strain  above d e s c r i p t i o n  relationship  under  It  should  constitutes general  a  that  pointed  complete  loading  Herrera(1964) has shown  sequence.  be  and  stress-  unloading  combinations  Coulomb f r i c t i o n and l i n e a r e l a s t i c elements e x h i b i t t h i s of  behavior.  Newmark  and  Rosenblueth(1971)  out  of type  suggest  the  Masing model t o be a good r e p r e s e n t a t i o n o f s o i l b e h a v i o r  and  some e x p e r i m e n t a l support  f o r the  s u g g e s t i o n i s p r o v i d e d by  Marsal(1963). 3.3.2  Volume Compaction and Pore Pressure  Increments.  A t y p i c a l one d i m e n s i o n a l l o a d - d e f o r m a t i o n curve f o r a  laterally  confined  sand  specimen i s shown i n F i g u r e  where ab r e p r e s e n t s l o a d i n g , be u n l o a d i n g  and  ce  3-10a  reloading.  Log P  (a)  FIG. 3-10  One-Dimensional  Soil Model ON  00  69 I t i s generally observed that while ab i s quite d i f f e r e n t from be  and  cd  the l a t t e r  two are reasonably close together to  allow t h e i r representation by one curve, model  as  shown  i n Figure 3-10b  ce.  A  can be used to describe, to  some degree of approximation, the behavior of sand dimensional compression. be  assumed  conceptual  under  one  The spring element i n the model can  to be non-linear so that the i d e a l i s e d unloading  and reloading curve ce can be simulated. the model should  be  deformation curve  as  of a  The s l i p element i n  'hardening'  type  i n Figure 3-11&,  shown  with  Using  load-  i . e . , the force  required to cause s l i p increases with the magnitude displacement.  a  of  slip  such a conceptual model i t can be seen  that the volumetric s t r a i n , €  v >  i s given  by  the  following  equations;€  v  A£ where K  l r  v  = =  < lr K  +  K  s>'P  K .A lr  <3-25>  P*Pmax p<p  P  (3-26)  m a x  i s the modulus (non-linear) of the spring  element,  K_ i s the pseudo-modulus of the s l i p element, and, s p „ i s the past maximum v e r t i c a l e f f e c t i v e s t r e s s , 'max When  horizontal  cyclic  shear i s applied to a sand sample at  equilibrium with an v e r t i c a l load, e.g., c y c l i c test,  volumetric  strain  Since  shear  occurs at constant v e r t i c a l load so  that a load-deformation curve obtained.  simple  as  shown  i n Figure  3-12  is  the a d d i t i o n a l volume reduction caused by  c y c l i c loading i s non-recoverable an a d d i t i o n a l feature has to  FIG. 3-11  Force-Displacement  Relationship  of  Slip  Elements  -s3 O  71  Fig. 3-12 Volume Compaction Caused by Cyclic Loading  72 be incorporated into the s l i p element of the conceptual -  namely,  the  addition  of  s l i p component caused by c y c l i c  shear s t r a i n as described i n previous s e c t i o n .  25)  (3-26) have  and  *v A€  =  K  =  y  where €  < lr  +  K  s>'P  <Srp  +  P*Pmax  the  (3-  ^-27)  p-cpJJ^  K .AP lr  (3-1).  i s defined as i n equation  changes  Equations  to be replaced by the following,  £J  V  model  properties  (3-28)  Since  compaction  of the s o i l element, e.g. void r a t i o  and shear modulus, more or l e s s i n the  same  fashion  as  the  maximum past pressure, P ax' l ^ i " maximum past pressure, p!L , i s used i n equation (3-27) and (3-28) instead max a n  e c  u  v a  e n  t  m  V  °f  P™<, max .  The determination of p*~ max i s shown conceptually i n  v  Figure  QV  3-12. To account for the undrained  element,  is filled  with  these circumstances,  the same time  when  an  (grain s l i p ) s t r a i n , AG  water  and  drainage  a  soil  causing  pore water.  a  increment  in  volumetric  portion  pore  has  to  be  of i t s stress transferred to the this  phenomenon  l e v e l can be obtained from Martin, Finn and the  compaction  , i s induced by a shear s t r a i n c y c l e .  Description of  expression for  prevented.  volume reduction cannot occur at  As a r e s u l t , e l a s t i c s t r a i n i n the spring element relaxed  in  i t may be assumed that the conceptual model shown i n  Figure 3-10b Under  condition  water  at  the  grain  Seed(1975)»  and an  pressure increment, A(T » W  was  73 found to be,  . AC7  where n K  "  W  (  =  ' nK K  x  ^-29)  A€  w  i s the porosity of the s o i l element, w  i s the bulk modulus of water.  When water can be assumed to be equation  (3-29) can be ACT  It  j » lr' vp K  l r  should  be  incompressible  or  K  > > K w  =  w  pointed  stress  needed at a l l . dimensional  K  l r  .A  (3-30)  € v p  out that  during  the settlement or  This eliminates  maximum  (3-27) i s  i s never exceeded and equation  loading  r  s i m p l i f i e d to,  l i q u e f a c t i o n process of a hroizontal deposit the past vertical  i »  the determination  not  of one-  c h a r a c t e r i s t i c s of the granular material  under consideration, nevertheless, the one-dimensional rebound c h a r a c t e r i s t i c s , K^ , has to be q u a n t i t a t i v e l y known i n order r  that equation  (3-22) can be  used.  In the determination of the magnitude of volumetric compaction s t r a i n due to horizontal c y c l i c s t r a i n  (3-1)  and  (3-la) are used.  experimental data pressures  and  for a relative  Monterey  confirmed.  Sand  densities.  volumetric compaction s t r a i n conditions described  Recently,  on  Pyke(1973)  obtained  at various  normal  The independence  vertical  i n the previous  A comparison of the 6 ^  expressions  stress  under the  paragraphs was  function  of  again  f o r sands of  74 different  relative  obtained  densities  Pyke(1973)  by  can  into  be  one  made by p l o t t i n g data  graph.  A  'reduced'  incremental volumetric compaction s t r a i n , ^C^p, i s obtained by the  product  of  relative  density  i n per  cent,  D, r  and  incremental volumetric compaction s t r a i n , A C ^ , i . e . ,  (3-31)  D .A€ r vp  vp  The graphs are plotted i n Figure 3-13. the  curves  I t may be  seen  that  f a l l more or less to an unique pattern suggesting  the p o s s i b i l i t y of using incremental  one  shape  function,  U,  f o r the  s  compaction s t r a i n f o r sands at d i f f e r e n t r e l a t i v e  densities, namely, A €  vp  ^i-W^  =  (3-32) (3-33)  vp  (3-34) where ^ , ^ , ^ 2  and ^  are constants s i m i l a r to those used i n  equation (3-la).  In t h i s usage  i s inversely  to  density  r  the  relative  proportional  value, D , while i(/ , ^ 2  and ij;^ are  assumed to be constant f o r the same type of sand at densities. change  Martin, Finn and  increments  completion of However,  a  by  strain  Seed(1975)  assuming cycle  i t i s generally  that as  observed  different  calculated the volume these  shown  occur  i n Figure  at the  3-l4a.  that i n a drained c y c l i c  0.1 Shear  FIG. 3-13  0.2 Strain Amplitude —  'Reduced' Incremental  *f  %  Volumetric Strain ->3  77 shear t e s t on granular s o i l gradual volume compaction pressure build up i n the undrained stage  of unloading  from  occurs  during the  maximum shear s t r a i n .  In view of  t h i s , i t i s not unreasonable volume compaction from  test)  to assume a l i n e a r  value as  shown  assumed  i n order  that  change increment i s added every h a l f c y c l e . involved i n t h i s Martin, Finn and  3.3.3  3  thesis  the numerical  (3-3*0  the volume  For computations values  obtained by  and fy.  Hardening and E f f e c t of Pore Pressure Rise  cyclic  suggested  loading  that equation  the fundamental  behavior  of sands  i n simple shear, Martin, Finn and Seed  (3-11)  s t r e s s - s t r a i n behavior. and  to be  Seed(1975) were used to determine the values  In a study on under  of  The  by an expression s i m i l a r to that of equation  though the value of i f ^ i s halved  2  3-l4b.  i n Figure  magnitude of volume change increment i s s t i l l  of q/ , y  variation  increment with the reduction of shear s t r a i n  i t s maximum  determined  (or pore  might be used to represent the  In t h i s equation  the parameters  a  b are constants for a given load cycle but i n general are  functions of the volumetric compaction given by equations  (3-12)  strain,  (3-13).  and  the shear modulus, G , i n any  cycle  mn  ^yp. n d are a  The maximum value of of loading  n  can be  obtained i n the following manner, G  mn  =  < hv / dT  a  t  *=  0  (3-35)  78 From equation (3-H)t therefore, we obtain  mn = JK'  G  For  the case  3-36)  a  (  of i n i t i a l loading when no hardening  has yet  occurred we obtain mo  G  =  ' J o T /  A  (3-37)  l  A f t e r n cycles, substitution f o r a from equation  mn  G  i n which € cycle  "  /  A, - — 1  2  A  €  +  A  v  i n which a^ and a  3 vpJ  2  G  f o r fcr^ / A  m Q  1  to the nth of  equation  according to  are constants.  We can obtain an expression stress from equation  (3-11) by  The maximum shear stress, T t—oo  (3-38)  p  With algebraic manipulation  (3-12) and substitution of equation (3-37) we obtain  gives  €  i s the volumetric s t r a i n accumulated  of loading.  (3-12)  i s from equation  f o r the maximum shear  considering very large s t r a i n s .  , i n the nth cycle of loading as  (3-11),  = ! ^ / o During i n i t i a l loading when  (3-40) no  hardening has taken place  we  79 obtain  (3-19) gives  and a comparison with equation B  l  =  /PlK  1  Following a procedure find that a f t e r  n  tan "  (3  s i m i l a r to that used i n d e r i v i n g G  and b  vp  T  l  1  mo  b. + 1  IZt]  2  P^3) V B -I  i s the accumulated volumetric  strain  to  define  al(1975) used  of the constants A. and B  1  were used to determine  s t r a i n hardening. used because G  mQ  Herein only four hardening and  T  are determined  m 0  the Hardin-Drnevich equations.  and  1  a , 2  b^  (3-^3)  to  constants  are  independently using  The four hardening constants,  2  the  G ^ and  define the  and b , are obtained by f i t t i n g equations to  six  mo  and the remaining four constants were used  1  b,  the stress s t r a i n r e l a t i o n s h i p of sand.  1 1  a ,  and  are constants.  2  constants  X  we  6.  +  I t should be noted that Martin et  Two  m n  42)  cycles of loading  L  where again  "  (3-39)  r e s u l t s of constant s t r a i n c y c l i c loading  tests. In the case of undrained shear the increment i n pore water pressure, A<7", may w  and  (3-34).  be  computed using equations  In one load cycle the t o t a l change  (3-30)  i n effective  80 s t r e s s w i l l be Ac^, and the e f f e c t i v e s t r e s s changes from to  cr^ =  y0  0  The new l e v e l o f e f f e c t i v e s t r e s s  - Acr .  <f  <r^.  w  will  a l s o a f f e c t the maximum shear modulus, G , and the maximum mn shear s t r e s s , '^ * a p p l i c a b l e t o the next c y c l e o f l o a d i n g . xm  It i s difficult  to  apportion  g r a i n s l i p and volume compaction. that  the major c o n t r i b u t i o n  the hardening between  F i n n e t al(1970)  to  hardening  comes  c r e a t i o n o f more s t a b l e g r a i n t o g r a i n c o n t a c t s s l i p a t contact  points rather  than  from  Although the matter cannot be considered a t t r i b u t e most o f  the  hardening to  conclude t h a t hardening occurs a l s o  from  due t o  increased settled,  grain  slip  the  cyclic density.  herein,  we  and  we  so  i n undrained s h e a r .  the case o f undrained c y c l i c shear response represented by the p o t e n t i a l v o l u m e t r i c  In  the g r a i n s l i p  is  s t r a i n , ^-yp*  For undrained sands, t h e r e f o r e , moduli and maximum allowable  suggests  the maximum  shear s t r e s s e s  f o r the  shear  nth cycle  are r e l a t e d t o the i n i t i a l v a l u e s by the f o l l o w i n g e q u a t i o n s :  mn  =  "mn  VP  G. mo  'mo  (3-44)  2 FvpJ [^voj aT fc  1  L  +  K—Z" b  i  +  vp  u  Vvp  v  (3-45)  'vo  i n which cr^ 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 , and Q  is  the  current  vertical effective stress.  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 then g i v e n  The  ^  resulting  by the f o l l o w i n g :  81  r  G m n  1  3.4  COMPUTATION  OF  +  SOIL  G  m  .  f  (3-46)  n  BEHAVIOR  USING  THE PROPOSED SOIL  MODEL.  ~  .  In t h i s section the values of stress and s t r a i n and other  related variables, e.g., shear modulus, damping, volume  compaction and/or pore pressure build-up, etc., based  on  are computed  the equations proposed i n the previous section and  comparison i s made, whenever possible, with experimental data. A hypothetical sandy s o i l i s assumed  to have  the  following  properties:e  * K  Q  (T  y0  With  0.83 ; = 36.0° ; = 0.50 ; = 1600 psf. =  these and with the help of equations (3-5) and (3-6) the  i n i t i a l maximum shear modulus, G , and i n i t i a l maximum mo stress,  T , can be calculated.  To calculate the compaction  m Q  c h a r a c t e r i s t i c s values of 0.400,  0.790, 0.563  used  respectively.  f o r i j ^ , q> • ^3  a  2  n  d  and  O.730  are  These values  were obtained based on the discussion given i n section 3.4.1  shear  3.3.  Hypothetical Drained C y c l i c Simple Shear Test. In  the drained  cyclic  simple  shear  test  the  82 effective  stress  i s kept constant.  (3-44) and (3-45) can be  The values of a , a 1  2 >  Consequently,  equation  reduced to the following!-  b^ and b  2  take  the  values  of  0.754,  0,406, 0.550 and 0.500 respectively as discussed i n section 3.3. By using a s t r a i n increment of 0.005$ a drained c y c l i c simple  shear  test  calculated. section  The  3.3  are  with  constant  procedures followed.  shear s t r a i n amplitude i s  and  equations  The  first  is  well  in  four s t r e s s - s t r a i n  hysteresis loops were computed and the r e s u l t s Figure 3-15.  described  are  shown i n  I t can be observed that the e f f e c t of hardening  represented  by the increasing l e v e l of shear stress  required to bring the sample to the same peak s t r a i n  at each  cycle. The  average shear modulus or secant modulus i s also  computed f o r various shear s t r a i n l e v e l s and compaction values as shown i n Figure terms  of  3-l6.  Similar  larger 0.5% that  relationships  loading cycles are shown i n Figure 3-17.  shown that at larger strains hardening has e f f e c t on the average  1  cycle  while  in  The data  a proportionately  shear modulus.  For example, at  s t r a i n the average shear modulus a f t e r 5 cycles i s after  strain  twice  at 0.05% i t i s only 10% greater.  600.  Fig. 3-15  Hyperbolic Stress-Strain Loops Showing Effect of Hardening  (Drained Test)  .02  Fig. 3-16  .05 0.1 0.2 Cyclic Shear Strain Amplitude (%)  Shear Modulus as a function of Cyclic Shear Strain and Volumetric Strain  0.5  Fig. 3-17 Shear Modulus as a function of Cyclic Shear Strain and Number of Cycles  CO  •Vn  86 These observations indicate the d i f f i c u l t i e s and approximation inherent i n s e l e c t i n g a single shear modulus-strain curve using  in  dynamic  analysis  of  for  earth structures composed of  granular materials. The damping r a t i o defined by (see Figure 3-B). D  =  1  (i LQ)  area of loop APQ ' area of AOAB  i s also calculated.  v  For convenience  of evaluating  the  of the s t r e s s - s t r a i n hysteresis loop the volumetric strain  i s assumed  to  be  constant so  that  area  compaction  the  following  expression f o r the damping r a t i o can be obtained  D = | [ l + i ] [ l - log (l+Y)A] e  G where  Y  =  - -f-  (3-50)  f  -22-  (3-51)  Sin The  variation  amplitude,  of  damping  ratio,  D,  with  shear  f , at d i f f e r e n t values of compaction  strain  strain,  €  ,  are calculated and the r e s u l t s shown i n Figure 3-18. Pyke(l973) average shear Monterey  #0  apparatus. and 3-20, soil  modulus Sand  presented  and  obtained  damping from  data on the v a r i a t i o n of ratio  with  the Geonor  strain  simple  for shear  These experimental data are shown i n Figures 3-19 I t can be seen that the behavior of  model,  remarkable  has  namely  shown  in  Figures 3-16  s i m i l a r i t y to the experimental data.  the  and  proposed  3-18  bear  FIG.  3-18  Computed Damping Ratios  FIG.  3*19 Strain  Shear and  Modulus  as  a Function of  Total Compaction (after  Shear  Pyke, 1973)  Damping  68  Ratio  90 3.4,2  Hypothetical Undrained C y c l i c Simple Shear Test. In t h i s section the assumed s o i l data  the  previous  section  is  used.  In addition,  Seed(1975)  obtained by Martin, Finn and  mentioned  in  the power law  i s used f o r  the  one-  dimensional rebound s t r e s s - s t r a i n r e l a t i o n s h i p , namely, €  where  vr  =  ^Ko^K^  3-52)  (  i s the e l a s t i c (rebound) volumetric s t r a i n , <f  VQ  i s the past maximum  effective  v e r t i c a l stress  from  which unloading i s taking place, <7-y  i s the current e f f e c t i v e normal stress,  m, n and k  0.0025  respectively.  The set of unloading curves given by 3-21.  The  0.430, 0.620 and  are constants taken to be  2  and,  one-dimensional  (3-52) i s  rebound  obtained by d i f f e r e n t i a t i n g equation  modulus  shown i n Figure can  then  be  (3-52),  K vr (cr;) "" 1  The procedures i n similar  to  those  1  computing used  for  mk (<r; ) n  /  2  an  (3-53)  m  0  undrained  drained  cyclic  test.  There  additional step involved i n the c a l c u l a t i o n and calculation  (3-30).  of  test  that  are  i s one is  the  pore water pressure increments using equation  Small s t r a i n increments of  0.005%  a  re  again  used.  91 However, model  to  illustrate  during  amplitude  the  undrained  cyclic  test  softening  cyclic  shear  i s chosen.  behavior of the s o i l a  constant  stress  The f i r s t f i v e s t r e s s -  s t r a i n loops are shown i n Figure 3-22.  The decrease i n shear  modulus due to r i s e i n pore water pressure i s indicated by the increases i n s t r a i n amplitude generated by of  constant  amplitude  shear  representation of the development of  successive  stress.  pore  water pressure during c y c l i c loading are shown i n Figure  3-23.  progressive  increase  i n the  amplitude of c y c l i c shear  s t r a i n and the double curvature i n the development  curve,  which  are  strain  schematic and  The  shear  A  cycles  pore  water  very c h a r a c t e r i s t i c of actual  granular s o i l under undrained c y c l i c shear, are also from the model.  pressure  obtained  Recoverable  FIG. 3"2I  Volumetric  Strain -  e  One Dimensional Unloading (after Martin et al, 1974)  vr  Curves  93  FIG. 3 - 2 2  Hyperbolic Stress-Strain  Loops (Undrained  Test) Showing Softening due to Pore Pressure Rise  Time  FIG. 3~23  (sec)  Strain and Pore Pressure during Undrained Cyclic Shear Test  95 CHAPTER 4 MECHANICAL MODEL FOR THE ANALYSIS OF LIQUEFACTION  4.1  INTRODUCTION From the discussion of s o i l l i q u e f a c t i o n i n  chapters  i t i s seen  that  the pore  water  earlier  pressure,  <r ,  generated during shaking controls the occurrence and extent of soil failure. analysis  of  r i s e of pore normal  In view of t h i s , i t i s necessary to set up liquefaction pressure,  stress,  subjected  to  approximation,  in a an  p o t e n t i a l to calculate the maximum  or maximum  reduction  in  effective  saturated horizontal sand deposit when  oscillatory  liquefaction  base  motion.  As  an  f a i l u r e i s assumed to take place  when the e f f e c t i v e normal stress, o~ = <r^. -0^, y  zero,  an  o  i s reduced  to  or, f o r p r a c t i c a l purposes, to l e s s than f i v e per cent  of the i n i t i a l e f f e c t i v e normal stress, 0~ . yQ  A complete analysis of l i q u e f a c t i o n  should  include  the following processes, namely, 1. Dynamic response to base motion, 2. Increase i n volumetric compaction s t r a i n due to c y c l i c straining, 3. Pore water pressure generation as a r e s u l t of (2), and, 4. D i s s i p a t i o n of excess pore water pressure(consolidation). There  are  two main sources of n o n - l i n e a r i t y i n the  96  proposed analysis of l i q u e f a c t i o n : the n o n - l i n e a r i t y  of the  s o i 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 and that which a r i s e s from the coupling  of the above-mentioned four processes.  of the complex evident  coupling  of the processes,  As a r e s u l t  which  will  be  i n the following sections when the whole problem i s  formulated,  numerical  procedures  solution  of  the equations  equations  f o r numerical  have  to be  obtained.  computations  used  f o r the  For convenience, will  be  derived  immediately a f t e r each 'exact' equation i s obtained. 4.2  FORMULATION OF THE PROBLEM. A  level  sand  extent so that propagation can  be  deposit i s assumed to be o f i n f i n i t e of shear waves through the  treated as a one-dimensional problem.  without any confined  loss  to  a  of generalization  Consequently,  considerations  can be  s o i l column of unit cross-section through the  deposit, as shown i n figure 4-1. i s assumed  deposit  to be  completely  S o i l under the water  table  saturated so that when flow of  water occurs as a r e s u l t of s o i l compaction the d i r e c t i o n of flow  i s upward  towards  the ground  surface.  Governing  equations and relevant s o i l parameters w i l l be introduced  and  discussed i n the following sub-sections. 4.2.1  Dynamic Response to Base Motion Consider  a  s o i l element of thickness dz located at  height z from the base of the deposit as shown i n Figure Let  4-1.  x be the horizontal displacement, 1 the horizontal shear  FIG. 4-1  Idealization of the Response Problem  98  s t r e s s and y  the t o t a l mass  density  Dynamic e q u i l i b r i u m r e q u i r e s 2 p . d z . 3t  of  the  soil  element.  that,  ||.dz.l  0  ,  or,  ||  . For  visco-elastic  .o  -  materials  in  (4-x)  which the s t r a i n and s t r a i n  r a t e e f f e c t s are a d d i t i v e and/ non-coupled  the  stress-strain  r e l a t i o n s h i p can be r e p r e s e n t e d by, T  where  In  =  f^iT) +  (4-2)  f (f) 2  f  i s the shear  jr  i s the shear s t r a i n r a t e , and  f^  and f  strain,  are f u n c t i o n s  2  of f  and  the case where the s o i l p r o p e r t i e s  e q u a t i o n (4-2)  Since  if =  can be s u b s t i t u t e d  i f  into  f  respectively.  are c o n s t a n t w i t h depth (4-1)  d i r e c t l y to give,  ,  dZ  1£  =  2 *  .sir  _  A^ dz  9z  t[Ct)  f '(r) 2  2  =  (4-4)  and ,  and w r i t i n g  -$r> lCf) f  , and  (4-5) =  ^.f (/) 2  ,  99  Equation (4-3) can be rewritten i n the following form by using the r e l a t i o n s given by (4-4) and (4-5) »-  px  2 2 - f' l y i ) - f ' 2_(x) = 0 dz  J  t  t  f^=0 and ^2~ *  I t can be seen that when the shear modulus,  (4-6)  3z  G  the f a m i l i a r  w  equation  n  e  r  e  G  represents  f o r the dynamic  response of an undamped l i n e a r shear beam r e s u l t s , namely, j>x  -  G.^-J(X)  =  0  (4-7)  In general, the s o i l properties are functions of depth, (H-z), so that 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 given by equation (4-2) has  to be replaced by a more general equation which contains  s o i l parameters as functions of depth also. equation  similar  to  In t h i s case  (4-6) cannot be obtained i n a straight  forward manner since the p a r t i a l d i f f e r e n t i a t i o n in  equation  (4-3) would  involved  advantageous to apply numerical (4-1), finite  an  carried  more terms.  approximations  out  I t i s more to  equations  I t i s shown i n Appendix I that through the use of the difference method equation (4-1) can be d i s c r e t i z e d i n  the z-domain to give, m.  r  [M](X}  +  [C]{X)  +  [K]{X>  =  -  .'X  (4-8)  b  n so that a lumped mass system as shown i n Figure L m  used  f o r approximate  analysis.  matrices given i n the Appendix  I  4-2  can be  The components  of the  are rewritten  here f o r  100  m,  n  "'I "•^j'a&5S.'^5'iS^~ m; — — — — — — — —  X-  hi  »i k  c  i  wmmmmmm. m-.  h  l-c ,k 3  3  3  2  h  2  ^2*  K  2  Layer I  FIG. 4 - 2  Approximation Mass  System  by  Lumped  101 convenience :[M]  =  r  m  0  1  0 m 0 0  0 0  2  0 " 0 0  * •  •  (4-9a)  •  •  o  nj  [c] = ' l 2 C  0  -c  + C  c  0  2 3 +c  0 0  -3 3 4 c  c  +c  0  •  •  •  •  •  .  -c n  •  •  •  «  •  •  0 ' 0 0  •  •  •  -k n  k  <*  [K] =  V 2 "2 0  0  "2 2 3 k  k  k  k  0 0  +k  •  §  0 where  fyi-* i i  =  n c  i  b  = fA'ti 1  =  foU* C  I t should  be  1-2  h  +  b  h  A ) .b,/(h. 1 1 1  1  pointed out that  and width of the i - t h layer  C-  j,) are 1-2  1-T2  values  2  n  -  1  , and .  1—2  b^  p. i , J  fo(fi  »  (4-9c)  1-2  h^ and  while  (4-9b)  nj  i = 1  .f, i )  1-2  1—2  c  %+i i+l i+l  i . z . i).b./(h..f. i)  1—2  0 0 0  •  1-2  are the thickness f , ( f - i,z- i ) 1  1-2  and  1-2  referred to the centre of the 1-th  102 layer.  4.2,2  Hyperbolic  Stress-Strain Relationship, f^, and Material  Damping Values, f . 1  I t was suggested element  under  3  that  for a  soil  simple shear condition, such as i n the present  case of a horizontal s o i l from  i n chapter  deposit  transmitting  shear  waves  the base to the surface, a hyperbolic hysteretic s t r e s s -  s t r a i n representation i s a close approximation to the actual stress-strain.  Consequently,  the, f g function described i n  section 4.2.1 takes the form of equation (3-46) together the  Masing behavior given i n section  with  3.3.1.1.  When v i b r a t i o n a l energy i s being transmitted  through  a material medium a portion of i t i s dissipated i n t e r n a l l y due to  a  number  of mechanisms at the grain l e v e l , e.g. Coulomb  f r i c t i o n , p l a s t i c deformation, This  d i s s i p a t i o n ^ of energy  viscosity, which  dislocation, etc.  causes a decrease i n the  amplitude of v i b r a t i o n i s generally c a l l e d material damping or i n t e r n a l damping. into  two main  Generally, material damping i s  categories - a frequency dependent type and a  frequency independent type. equation  (4-2)  over  one  I n c i d e n t a l l y , i f one strain  s t r e s s - s t r a i n loop i s s i m i l a r to Figure  3-3a  classified  and  cycle, the loop  assuming AB  integrates that the  as shown  in  that the i n t e g r a l s involved are integrable,  one should be able to obtain the energy dissipated during s t r a i n cycle ABA, W  d  , as  the  103 JX .df  or,  =  yqch.df  w (ir,f) = d  Equation  (4-10)  jf (f).df  +  2  w'(ir) + vq(0  (4-2)  can  be  (4-n)  used  to describe a simple  l i n e a r v i s c o - e l a s t i c material, e.g., a Kelvin s o l i d where,  and  fAf) 1  =  c.'f  ,  f (jf)  =  K.r  .  2  (4-12)  where C and K are the v i s c o - e l a s t i c constants.  For  harmonic  v i b r a t i o n , f=Asin«/t, the energy dissipated per cycle, W,  can  d  (4-10).  be obtained from equation  Since,  M  j f^'f)  .df  =  Jc.f.r. dt o  2 j  and  f ( T ) .df 2  giving It  W  d  can  be  =  0 ,  =  C7TWA  2  of  f^  ( 1  ^_  1 3 )  seen that the energy dissipated per cycle, or the  damping depends on the frequency case  ^  where  the  w .  For  this  particular  damping stress i s l i n e a r with s t r a i n  rate, the term viscous damping applies. For f r i c t i o n a l  materials,  constitutive equation (4-2)  t^'f)  =  0,  the  components  of  can be written as follows :-  the  104 fM) ^ f (f) 2  where  f  =  +f  for  df>0  =  y -f  for  df<0  i s a p o s i t i v e quantity.  cycle can "be obtained  w  d  The energy dissipated per  s i m i l a r l y by (4-10) to give,  =  f  =  4.f .A  This type of damping  f (n.df 2  which  (4-i4)  is  called  Coulomb  Damping,  is  frequency independent. It the  can be observed from -6he above two examples that  constitutive  frequency  relation  dependent  dependent  stress  independent  (4-2)  damping,  component  damping,  or  is  can for  be  used  which  the  responsible,  structural  to  include  strain-rate  and  frequency  damping,  which  associated with the non-linear hysteretic s t r e s s - s t r a i n s t a t i c loading.  •  >-  is  under  •  In the equivalent l i n e a r analysis where a non-linear hysteretic  material  i s represented  by a l i n e a r v i s c o - e l a s t i c  model, damping has to be introduced a r t i f i c i a l l y through W W" vanishes. factor,  C eq ,  In t h i s case has  to  be  the  equivalent  obtained  by  viscous  equating  as  damping  the energy  dissipated per cycle i n the non-linear material and i t s l i n e a r equivalent, so that,  105 In a t r u l y non-linear analysis, as i s i n the present analysis, where the hysteretic s t r e s s - s t r a i n  relationship i s  used, a r t i f i c i a l viscous damping c o e f f i c i e n t given is  no  (4-15)  "by  longer needed, instead, the function f^ can be used to  describe the actual viscous material.  Or,  in  a  damping  effect  present  in  the  fashion s i m i l a r to that of the method  using equivalent viscous damping, i t can be used to compensate for the  the i n s u f f i c i e n t damping caused by actual  stress-strain W (f,T)  the uncoupling of one  though  d  this  approximation  by the function f^(f). into W  Note that  and W" i s not a  necessary  uncoupled  determination  form  the  of  first  damping  type  values.  of damping can be  obtained through v i s c o s i t y measurements while the l a t t e r can  be  of  i t brings about convenience i n discussion as well  as i n the experimental In  the  obtained  by  testing  the  type  material i n question i n a  quasi-static manner. For the present purposes,  i s assumed to be  of  the l i n e a r viscous type, i . e . ,  f (f) x  to  = Cf  (4-16)  take into account of the damping e f f e c t caused by the flow  of water occurring as a r e s u l t of the deformation skeleton,  see  e.g.  Biot(1956).  measure the v i s c o s i t y of saturated  Yen(1967) sands  near  of the attempted  to  liquefaction.  However, i t i s not clear whether the v i s c o s i t y values referred  soil  obtained  to the actual strain-rate dependent stress component  106  or to the ' a r t i f i c i a l ' equivalent viscous damping. 4.2.3  Volume Compaction and Pore Pressure Figure  4-3b  shows  in  chapter  3.  proposed  represents the s o i l container  the material The  model  that  spring-and-slip  while  the  fluid  was  element  inside the  represents the pore water of the granular material.  I t i s assumed composed  skeleton  Generation.  that  of two  component,  the t o t a l components,  volumetric namely,  and the compaction  strain,  the rebound  (slip)  £  v  , is  (elastic)  component,  ^vp»  i.e., €  =  v  Incremental  volumetric  € + vr  € vp  (4-1?) . "  compaction s t r a i n , &€.y^>  i s given by  the following expression,  e vp  =  Ae  +1-  ( f  a "*2V  +  * r 3  a  2  +^  (4vi8)  .  where  f  shear  s t r a i n l e v e l reached, i n the l a s t s t r a i n cycle (or h a l f  Q  i s the amplitude of shear s t r a i n ,  s t r a i n cycle) and the vj/s are s o i l constants cyclic  simple  shear  or t r i a x i a l t e s t s .  or the maximum  determined  In the case of dry  granular material an a p p l i c a t i o n of a shear s t r a i n cycle amplitude  with  f w i l l cause the t o t a l volumetric s t r a i n to change  by an amount s t r a i n , AE^,  from  equal  to  that  of the volumetric  compaction  since there i s n_o change i n the e l a s t i c component  when the e f f e c t i v e normal stays constant, i . e . ,  Grd. Surface  MR  02 0-  w +  ^dZ| oz  dz  / / / /  l i t  Base Rock  (a)  FIG. 4 - 3  ,  (b)  Consolidation Model  108  vr  A €  *  '  •  0  (4-19) and  A£y  I f the s o i l  is  =  vp  saturated  instantaneously,  and  relaxation  w i l l take place. vertical  &.  The  stress  is  drainage  of  result  cannot  take  place  the e l a s t i c s t r a i n component is  that  a  portion  of  the  transferred to the pore water so that an  increment of pore water pressure,  A(T , r  occurs,  as explained  i n chapter 3 . A<T  W  where  K< lr  ir'  A €  vp  ( 4  -  2 0 )  soil  D e t a i l discussion of the volumetric s t r a i n caused  b,y c y c l i c loading  4.2.4  K  i s the one-dimensional'rebound modulus of the  skeleton.  repeated  -  is  given  in  section  3.3.2  and  is  not  here.  Reconsolidation. The  pressure compaction  equation  generation can  as  for a  reconsolidation with pore water result  of  internal  volumetric  be derived i n the same manner as that f o r the  conventional consolidation equation.  The assumptions used i n  the d e r i v a t i o n are i 1.  Complete saturation.  2.  Negligible compressibility of pore water.  3.  One  4.  One dimensional flow of water.  dimensional compression.  109 5.  V a l i d i t y of Darcy Law. Consider a  section  as  shown  typical  soil  i n figure 4-3a,  element  the  unit  cross-  where v i s the v e l o c i t y of  flow of water into the s o i l element. flow,  of  For the  continuity  of  difference i n volume of flow of water into and out  of the element should equal the volume change of  the  element  i t s e l f , so that we have, v + & z  - vldt  £  =  f^.dz.dt  -  U  <*-> 2 L  I f the flow s a t i s f i e s Darcy's Law, v  where  k  =  -k  i s the c o e f f i c i e n t of permeability,  f  i s the unit weight of water , and w tj i s the excess pore water pressure. W  If  the  expression  f o r v i s substituted into equation  the following equation i s obtained,  W  From equation  Since,  (4-17) i t can  Ae^ =  jp-.A<r  lr  v  be seen that,  (4-21)  110 = V - (A(T - ACT)  (4-25)  lr  During the l i q u e f a c t i o n of a horizontal s o i l deposit the t o t a l overburden pressure i s not changed, i . e . ,  (4-23)  so that equations  (4-26)  , (4-24) ,  (4-27)  and  can  be combined to give, +K  When  the  coefficient  of permeability,  depth and i f ( i « / ^ ) K  k  r  -^'  w  i s denoted by  (4  -  28)  k, i s constant with  c , y  equation  (4-28)  t  can be s i m p l i f i e d to, •Wvt  ^vr.  (4-28)  Equation  or  (4-29) i s the desired equation governing  the d i s s i p a t i o n of excess pore water pressure with pore  water  compaction. compaction €^ =0,  pressure  generation  due  I t can be seen that by strain  function  continuous  to i n t e r n a l volumetric  setting  the  volumetric  i d e n t i c a l l y equal to zero, i . e . ,  the well known Terzaghi's Equation of Consolidation i s  obtained. Inspite of the extra term i n equation  (4-29),  i t is  s t i l l a t y p i c a l parabolic p a r t i a l d i f f e r e n t i a l equation i n one  Ill  dimension and i s i d e n t i c a l i n form to that f o r the of  conduction  heat i n a bar with heat sources d i s t r i b u t e d along the bar.  For example, Sneddon(1957) gave the following equation f o r the conduction of heat i n s o l i d s ,  =  3t where Q  (cv  2  i s the temperature  k  (4-30)  H(r,e,t)  +  0  j>c  at a point with p o s i t i o n vector r,  i s the thermal conductivity,  j). i s the mass density of the material, c  i s the s p e c i f i c heat of the material, and,  K  =  k/(pc).  The function  H(f,0,t)  gives  the  generated per unit volume per unit time. gives  us  the  conduction.  rate This  of is  amount  Physically  r i s e of temperature quite  consistent  of  heat  (H/j)c)  when there i s no  with  the  physical  meaning implied i n the l a s t term of equation (4-29), namely, the  rate  of r i s e of pore pressure when there i s no drainage!  Note that with the help of equation (4-18) the equation that  (4-28)  the  expressed  or  internal as  a  Tf^  volumetric  in  compaction  strain  can  be  function of the t o t a l number of shear s t r a i n  .  method  soil  element  has  been  The f i n i t e difference scheme f o r equation  (4-29) i s derived i n Appendix integration  term  (4-29) can be t h e o r e t i c a l l y evaluated so  cycles (or h a l f cycles) to which the subjected,  last  III.  When  i s used to determine  a  step  by  step  the dynamic response  of the s o i l deposit the shear strains can be computed f o r a l l  112  the  layers  f o r each  time  increment  generation and d i s s i p a t i o n can he  and the pore pressure  incorporated.  Discussion  of the solution procedures i s g i v e n i n the next section, v  4.3  SOLUTION SCHEME. A  numerical  scheme  using a step by step method i s  used to solve the governing equations. since  [K] i s a  In  (4-8)  equation  function of the displacement vector, {X}, a  d i r e c t integration procedure i s used so that the displacement, v e l o c i t y and acceleration values at any time t can be to  values  at  (t-At),  where At i s the time step used i n the  d i r e c t integration procedure. numerical  integration  There are several methods f o r  of equations s i m i l a r to equation (4-8)  and a comparison of the methods based on t h e i r accuracy present  has  Newmark's  used  i n the  discussed i n Appendix I I . volumetric  method(1959)  of  step  compaction  The actual Also,  function,  particular  For the by  step  equations  value every  time  the  i n equation €  shear s t r a i n ,  ,  f ,  In  the as  increase  a in  associated with the n  3-14.  (4-29),  i s handled  same layer decreases from a maximum value, f„max „, J  and  computation are derived and  discontinuous function, which has an incremental  Figure  and  was chosen f o r i t s unconditional s t a b i l i t y i n the  analysis of l i n e a r systems. procedures  stability  been given by Bathe and Wilson(1973).  prupose  integration  related  as  shown  in  t h i s case the f i n i t e difference method of  solution becomes quite s u i t a b l e .  The relevant equations  and  113 expressions  arising  from  the use  of the f i n i t e difference  method of solution f o r equation (4-29)t or equation (4-28) f o r the more general case, are given i n Appendix I I I . A b r i e f outline of the step by step solution  scheme  can be given as follows :1.  on the current values of f,  Based  and <r at time t , y  the shear modulus i s calculated as described  i n section  4.3.2. 2.  The matrix [K] i n equation (4-8) i s then updated.  3.  With  the base acceleration value at (t+At) new values f o r  {X} , {X} , {X} and {f}  can be calculated using a  direct  integration method. 4.  When  a  maximum  layer  has i t s shear s t r a i n decreasing from the  value  incremented  i t s volumetric  according  compaction  to the equations  strain  is  described i n  section 4.2.3. 5.  Current value of pore water pressure, j * , w  using  the f i n i t e  difference method  are  as  calculated  discussed i n  Appendix I I I . 6.  Current e f f e c t i v e normal stresses, g-^, are then calculated "by 0~ - (7"-(7^, f o r use i n the next time step. v  v  A computer program has been written solution 4-5.  based  on  this  scheme and the o v e r a l l flow chart i s shown i n figure  A l i s t i n g of the computer program i s given i n Appendix  114  READ IN DATA  NUMERICAL INTEGRATION FOR X, X, X & f ETC.  PORE PRESSURE, VOLUME COMPACTION, IP REQD. CHANGE SHEAR MODULUS & OTHER PARAMETERS. RECONSCNIDATION CALCUL/LTION  PRINT OUTPUT &/0R TAPE OUTPUT FOR PLOTTING.  MORE [NCREMENTS ?.  YES  NO  FIG. 4 - 4  Flow  Chart  115  CHAPTER 5 APPLICATION OF THE MECHANICAL MODEL - AN ANALYSIS OF LIQUEFACTION OF A SOIL DEPOSIT  5.1  A  HYPOTHETICAL  SOIL  DEPOSIT  The method of analysis developed i n e a r l i e r chapters was applied to a hypothetical s o i l deposit described section.  The  equivalent Assumed  same  deposit  was also analysed by using the  l i n e a r analysis so that  values  f o r the  in this  comparison  can  be  deposit are as shown i n Figure  The v a r i a t i o n of r e l a t i v e density with depth was made to that of the Niigata deposit analysed by Seed et that  the  problem  this  thesis  5-1.  similar  al(l96l) so  i s as close to r e a l s i t u a t i o n as possible.  However, mention should be made that i t i s not of  made.  the  intention  to project the r e s u l t of the analysis to the  actual f i e l d condition at N i i g a t a . The maximum and minimum void r a t i o s were assumed be  1.00  and 0.50  so that the void r a t i o s at various r e l a t i v e  density can be calculated.  For convenience, the bulk density  of the s o i l above water table, ground surface  which  was  5*0  feet  beneath  i n t h i s case, was assumed to be 122.0 pcf while  the buoyant unit weight was assumed to be 60.0  pcf.  deposit  i n order  was  to  sub-divided  into  a p p l i c a t i o n of the previously  14  layers  described  method  of  The s o i l that  analysis  Ground Surface  o  Water Table =====  "T^5J0 ft  or 50  SAND 200.0ft  T  T  =  1  2  buoy.  2  s  PCf  0  6  0  P  0  100 H a.  c f  UJ  a  150  ' ^  Bedrock ^  ^  ^  ' ^  ^  200  ^  1  (a)  FIG. 5-1  Properties  of  a  Soil  Deposit  could  be made.  Table 5-1. up.  The properties of each layer are as shown i n  Notice the numbering of the layers i s from bottom  E f f e c t i v e shearing a n g l e w a s assumed to  the c o e f f i c i e n t of earth pressure at r e s t , K ,  be was  Q  30  and  calculated  from an equation given by Blshop(1958). K  1.  =  Q  The i n i t i a l mean e f f e c t i v e  stress  e f f e c t i v e v e r t i c a l stress, <r  0  the  use  of  obtained. system  ,  as  (5-2)  and  shown  Column 9 gives  calculated  based  the  and (3-6) initial  the  initial  maximum  shear  i n columns 7 and 8 of Table 5-1 the  masses  of  the  were  lumped  mass  on a unit cross-section s o i l column  while column 10 gives the i n i t i a l spring s t i f f n e s s e s . of damping values i s discussed i n section  5.2  from  , by,  v  equations (3-5) mo  m Q  calculated  0  maximum shear modulus, G , stress,t  <r£ was  ^.(l.+2.K ) .<4  <T = Through  (5-1)  - s i n d)'  Choice  5*2.  RESPONSE CHARACTERISTICS OF THE SOIL DEPOSIT. In order to investigate the steady-state response of  the s o i l deposit to harmonic base motion through the ij> this  proposed  obvious  step-by-step  use  of  non-linear analysis the volume change parameter  used i n equation is  the  (4-29)  was  set to zero.  The reason  for  since transient response i s included i n the  integration  technique,  steady-state  vibration  TABLE 5-1 • VALUES FOR LATER APPROXIMATION TO SOIL DEPOSIT LATER NO.  THICKNESS ft.  DEPTH TO CTR ft.  EFF VERT STRESS psf  D r  14  5.0  2.5  50.  7-5 15.0  305.0 760.0  13 12  5.0 10.0  1210.0  50.  11  25.0  1810.0  10  10.0 10.0  35.0  2410.0  9 8  10.0 10.0  45.0 55.0  7 6  10.0 20.0  5  G mo kip/ft  T  _  m.  K  mo psf  i slug  1 kip/ft 0  594. 938.  85. 212.  9.5 19.0  118.8  338. 506.  28.4  50.  1184. 1448. 1733. 2010.  674.  3010.0  55. 60.  37.9 37.9  118.4 144.8  841.  37.9  173.3 201.0  65. 70.  2283.  1009.  37.9  228.3  65.0  3610.0 4210.0  2557.  37.9  80.0  5110.0  2923.  20.0  100.0  6310.0  75. 80.  1177. 1428.  3369.  1764.  4  20.0  120.0  7510.0  85.  2099.  3 2  20.0 20.0 30.0  HO.O 160.0 185.0  8710.0 9910.0 11410.0  85.  3813. 4106.  56.9 75.8 75.8  255.7 146.1 190.6  2434.  75.8  85.  4380.  2770.  85.  4700.  3189.  75.8 94.8  1  % 50.  187.6  190.6 205.3 219.0 156.7  119 occurs  only  a f t e r a few cycles of v i b r a t i o n and i f hardening  or softening, i . e . change i n maximum shear modulus, maximum  shear  stress,  T^,  f a r as  m n  Moreover, to minimize  possible the duration of the transient v i b r a t i o n ,  the amplitude of base acceleration input was made to gradually  increase  i n the f i r s t few cycles to the required l e v e l and  kept constant, thereafter. that  and  i s allowed f o r every cycle the  steady state would never be achieved. as  G »  After a few t r i a l s  i t was  found  i f the amplitude of base acceleration i s made to r i s e as  a sine function within the f i r s t two and h a l f cycles transient v i b r a t i o n has i t s least linear  system  influence.  The  response  of non-  depends not only on the frequency of the input  loading but also on the magnitude.  Consequently,  complete  response c h a r a c t e r i s t i c s requires a l o t of computer runs. view  of t h i s , an a r b i t r a r i l y fixed base acceleration l e v e l of  0.065g was  used.  A time acceleration plot f o r the base input  i s as shown i n Figure  5-2.  Even though the hysteretic quite  cent,  soil  model  (3-30), at  addition,  due  inside  the s o i l  to the stepwise-llnear  structure.  nature  integration procedure a small amount of viscous  points.  0.10 per  s t r a i n l e v e l s around  viscous damping i s s t i l l needed to take into account of  the e f f e c t of the flow of water  needed  contributes  large damping, e.g. 20 to 30% i n terms of damping r a t i o  defined by equation  In  In  of the  damping  is  to eliminate high frequency v i b r a t i o n at the reversal Small amount of viscous damping was added through f  n  Time  FIG. 5 - 2  Input Base  (sees)  Acceleration  121 defined i n equation viscous  (4-2).  s t i f f n e s s values, k c  a  i  ft  are  i t  proportional  a  value  of  vibrations  given  0.0008  was  used surface acceleration  was  obtained  0.002  by equation  gives  while  p =0.002  It should be noted that changed  but  the  high  at r e v e r s a l points are eliminated.  2 %.  for  the  By  viscous  (4-24) i t can be shown that a /$  value of 0.0008 gives approximately while  initial  (5-3)  c a l c u l a t i n g the energy dissipated per cycle as  the  k  the amplitude of response i s s l i g h t l y  force  to  jS- i,t=0  88  response as i n Figure 5-3b.  frequency  linear  t=0* namely,  ±  response as shown i n Figure 5-3& gave  that  damping can be used and that the value of the viscous  damping c o e f f i c i e n t s , c  When  Also i t was assumed  a damping  ratio  of  0.6%  For a l l analyses described i n t h i s  chapter a ^5 value of 0.002 was  used.  Step-by-step i n t e g r a t i o n procedures were c a r r i e d out for at least 10 cycles of v i b r a t i o n and i t was  found  most cases steady-state response could be reached. 3b  shows  a  typical  acceleration.  time  of  the  acceleration  Figure 5-5.  Figure  5-  of the surface  Twelve frequencies were used and  amplitude to  that  of  computed for each frequency. in  response  in  A t y p i c a l s t r e s s - s t r a i n response of a layer i s  as shown i n Figure 5-4. ratio  history  that  of  the  steady-state  the  base  acceleration  The r e s u l t obtained i s  the  surface Input  was  plotted  For comparison purposes, the "SHAKE" program  122  cu u o>  in  c o  Q> U U <  Time  u  c o  ^  o  2.  f\ /I  /I  /I  /I  /I  0  v. 0) O u  (sees)  -Z  VJ VI  -4. _L  2. Time  FIG. 5 - 3  VI VI VI VJ 3.  4.  (sees)  Surface Acceleration Response  123  .RESPONSE OF HORIZONTAL SOIL DEPOSIT TO HARMONIC BASE ACC. STRESS-STRAIN FOR LAYER 12  0.04  0  0.02  Shear  FIG.  5-4  0.02  0.04  Strain (%)  Stress - Strain  Response  Steody State Surface Acceleration Base Acceleration  125 written by Schnabel et al(1972) was also used and i s plotted i n the same f i g u r e .  the  result  It can be seen that while the  equivalent l i n e a r analysis used i n the "SHAKE" program gives a pronounced  peak  at  the  resonant  frequency,  the non-linear  analysis gives a much less pronounced peak i n the response  curve.  A  frequency-  plausible explanation i s that when both  the equivalent l i n e a r and the non-linear models are to  the  same  s t r a i n l e v e l , the former has a larger amount of  s t r a i n energy which i s responsible curve.  subjected  Apart  from  this,  f o r the peaked  the two methods give  close a c c e l e r a t i o n r a t i o within the frequency  response reasonably  range of 2 to 10  hertz. A t y p i c a l surface acceleration by  the  "SHAKE"  program  response  calculated  i s shown i n Figure 5-6.  Comparing  with Figure 5-3, i t can be observed that the equivalent  linear  analysis gives the usual s i n u s o i d a l surface a c c e l e r a t i o n while in the non-linear case the acceleration response i s no sinusoidal.  Jennings(1964)  method to obtain the response system  with  non-linear  used of  a a  numerical  a  obtained  sinusoidal one. by  interesting  integration  single-degree-of-freedom  hysteretic  force-deflection  r e l a t i o n s h i p and observed that the a c c e l e r a t i o n not  longer  response  was  This i s i n agreement with the r e s u l t  the proposed  non-linear  analysis.  It  is  to note from Figure 5-3 and Figure 5-6 that both  the response computed by the 'SHAKE' program and that computed by the proposed non-linear analysis have more or less the same  127 phase s h i f t with respect to the input acceleration wave form.  5.3  MAGNITUDE OF PORE PRESSURE RISE AND  ITS  EFFECT  ON  THE  DYNAMIC RESPONSE OF THE DEPOSIT. In  this  used as the previous  section  base  input  section.  a simulated earthquake record  to  The  the  deposit  described  in  was the  acceleration values were calculated  from a random acceleration function  suggested  by  Bogdanoff,  Goldberg and Bernard(I96I), namely, J x (t)  = • .1. jjt.exp(-5jt).cos(wjt+ej)  g  where j^y ^ l  w  <  S^i  2  w  t>0  (5-4)  are r e a l p o s i t i v e numbers, <  3  w  <  4  w  ' • • » J » W  62t Ojt • • » t 0j  »  a n d  are  real  random  variables  uniformly d i s t r i b u t e d over the i n t e r v a l 0 to 27T. A large number of terms can be used when of  an  accelerogram  is  needed.  best  representation  In a study of the seismic  response of structure-foundation systems Parmelee et followed  the  and  al(l968)  suggestion of Bogdanoff et al(196l) and assumed to  accelerogram. function used was  be  constant  Ten  terms  when were  generating  used  and  a  simulated  the acceleration  given by,  10 x ( t ) = O.50.t.exp(-0.333t).£ cos(wjt+elj) g  1  in which the frequencies, Wj, and phase angles, in  Table 5-2.  t>0 8-j, are  (5-5) shown  In t h i s section, the same accelerogram scaled  128 to various values of maximum acceleration Figure 5-7  analyses.  was  used  shows a plot of the accelerogram.  can be seen that a f t e r 15 seconds the amplitude of relatively  for - a l l It  motion  small so that the analyses can be stopped at  is this  time without losing any s i g n i f i c a n t data. For the I n i t i a l parameter,  was  comparison could be 5-8a  Figure  still made  shows  computer  the  kept  run,  Identically  with  results  using  the  computer  volume  change  zero so that a  obtained  otherwise.  surface acceleration obtained by the  non-linear analysis while Figure 5-8b by  the  shows the same  program "SHAKE".  obtained  In both cases the  acceleration values of the base motion were scaled so that the maximum acceleration was point  0.065g.  Though  direct  both  i s no  responses.  material  softening  the  be  present  Thus, so f a r , i t can be seen when there  hardening  due  to  compaction  nor  material  due to r i s e i n pore water pressure, both the method  of d i r e c t integration of the and  to  comparison of the surface acceleration responses i s not  possible the same general c h a r a c t e r i s t i c s seem to in  point  method  of  wave  non-linear  equations  propagation with strain-compatible  material constants give more or  less  the  harmonic  as well as to random e x c i t a t i o n .  show the  typical  strain  dynamical  and  stress  same  response  to  Figure 5-9a and b  response  of  a  layer  obtained by the non-linear a n a l y s i s . In  order  that  the  pore pressure generated during  shaking might be accounted f o r the respective «J/.  values  for  FIG. 5 ~ 7  Simulated  Earthquake Accelerogram  TABLE 5-2 FREQUENCIES AND PHASE ANGLES FOR THE SIMULATED EARTHQUAKE ACCELEROGRAM 3  UK (rad/sec)  8. (radians)  1  6.00  3.7663  2  8.00  1.3422  «J  3  10.00  4.8253  4  0.2528  5  11.15 12.30  6>  13.25  1.8834  7  14.15  1.3320  8  16.20  1.7852  9  17.35  0.1517  10  19.15  2.4881  4.5204  TABLE 5-3 %-VALUES FOR DIFFERENT LAYERS Layer No.  Relative Density it)  +1  14  50.  0.0000  13  50.  0.5200  12  50.  0.5200  11  50.  0.5200  10  55.  0.4727  9  60.  0.4333  8  65.  0,4000  7  70.  0.3714  6  75.  0.3467  5  80.  0.3250  4  85.  0.3059  3  85.  0.3059  2  85.  0.3059  1  85.  0.3059  131  Time  FIG. 5-8  (sec.)  Surface Acceleration Response (no Pore Pressure )  132 each  layer  were calculated based on the Inverse proportional  r e l a t i o n s h i p between 4^ i n section 3.3.2 soil  deposit  Figure 3-11, Since  a  n  d  r e l a t i v e density, D , as discussed r  of chapter 3. has  Assuming  that  the present  the volume change c h a r a c t e r i s t i c s given by  the vj/^ values l i s t e d i n Table 5-3  are obtained.  the top layer (layer 14) i s above the water table i t s  volume change parameter \\>^ was set to zero so that excess pore water pressure due to compaction before,  the acceleration  maximum of were  0.065g was  could not be generated.  values  obtained.  As  were scaled down so that a In t h i s analysis the  layers  assumed to be separated by impervious boundaries and r e -  d i s t r i b u t i o n of pore pressure was not allowed.  As previously  discussed the integration procedures were carried out up to 15 seconds. and  Figure 5-10a  and b show "the surface  acceleration  displacement responses and Figure 5-Ha, b and c show the  s t r a i n , stress and pore water pressure of the 11th layer which i s situated between depths of 20.0 to 30.0  feet.  shows a s t r e s s - s t r a i n plot f o r the 11th l a y e r . which  liquefaction  occurred  marked i n a l l the response  pressure  The  time  5-12 at  was around 8.5 second and i t i s  plots.  From Figure 5-8a and 5-10a the pore  Figure  r i s e s during  i t can be  the  seen  liquefaction  surface acceleration tends to be attenuated.  The  that as process  relatively  large acceleration response that occurs at a time of 9 seconds when  pore  pressure i s not generated i s not seen i n Figure 5-  10a where pore water pressure  i s taken  into  account.  As  133  FIG. 5 - 9  Strain and Stress Responses for Layer II  (no Pore  Pressure)  FIG. 5-10  Surface Acceleration and Displacement Response  (Pore Pressure Generated)  0.6 — .  — [  T  ...,  ,  I  ,  ,  ,  1  — i  1  0.2  r  i  i  J\\ *A A ~ — — - \ A/ \/ v~  •A  0.4  1  0  c  to  —  -0.2  (a)  -0.4 -0.6 400.h  1  T  1  1  I  I  1  1  1  1  1  1  1  1  r  1  1  1  r  I  1  1  1  : 1  1  1  T — | — i — i — i — r  T  1  1  1  I  I  L  1  T  r  ^2000. a.  CO CO  £ 1000. Q_ <P i— O  (c)  0.  5.  Time (sec)  FIG. 5-11 Stress, Strain and Pore Response  15.  10.  for Layer 11  Pressure  FIG. 5-12  Stress - Strain Response for Layer II (Pore  Pressure  Generated )  137 expected  the acceleration  response  within  the f i r s t three  seconds was not much affected by the Increase i n pore pressure which as shown i n Figure 5-13 did not r i s e  to a  high  value  within t h i s time i n t e r v a l . 5-9a and 5-Ha shows the dramatic e f f e c t of  Figures the  pore pressure r i s e on the  layer.  s t r a i n response  the 11th  In the case where no pore pressure was generated the  peak s t r a i n response was +0.05$ same  of  peak  at around  2  seconds.  was also obtained i n the case where pore pressure  was generated.  In the l a t t e r  case,  however,  the +0.05$  s t r a i n response at around 2 seconds was dwarfed by the s t r a i n at 8 seconds, a shown  short  i n Figure 5 - l l a .  the input  base  maximum  strain  some doubt on the analysis  before  increasing  acceleration  response are decreasing. the  time  +0.25$  liquefaction,  as  The softening behaviour of the s o i l  l s well i l l u s t r a t e d by the when  The  strain as well  response  even  as the stress  Moreover, the large difference i n  levels  developed  i n the two cases casts  validity  of using  an equivalent  f o r the study  linear  of s o i l l i q u e f a c t i o n as was done by  Seed and I d r l s s ( I 9 6 7 ) . The offset of the s t r a i n response of the 11th layer, and hence the surface displacement response, the  to one side of  zero axis a f t e r l i q u e f a c t i o n i s probably due to the s l i g h t  non-symmetry of the input base motion. stress-strain  response  Figure 5-12 shows the  of the 11th layer and i t can be seen  that even a f t e r l i q u e f a c t i o n the stress and s t r a i n were  still  138 describing  the  describe. reveals  kind  Close  of  loops  which  examination  of  they were supposed to  the  stress-strain  plot  the s l i g h t one-sidedness of the stress response a f t e r  9 seconds.  The reason i s that at such low l e v e l of confining  stress the f a i l u r e stress i s very small so that a small stress increment near t h e . f a i l u r e stress  can  induce  a  very  large  strain. 5-13  Figure  shows the d i s t r i b u t i o n of pore pressure  in the 200-foot deep s o i l pore  water  pressure  was  deposit.  Initially,  an increasing  the  excess  function of  depth.  However, as time went on the s o i l layers situated and. .50  feet  pressure. strains  depth  interval  Consequently, at  later  water pressures. liquefaction  between  10  began to develop higher excess  these  layers  developed  larger  times which led to even higher excess pore The  interaction  occurred  in  the  seemed  11th  to  go  on  until  layer which i s situated  between 20 and 30 feet.  5.4  EFFECT OF PORE PRESSURE DISSIPATION. To  dissipation  investigate on  the  effect  of  pore  pressure  the l i q u e f a c t i o n process the complete coupled  analysis, i . e . including dyanmic response, pore water pressure generation and d i s s i p a t i o n , was performed on the s o i l described  earlier.  deposit  The simulated earthquake input described  i n section 5*3 of t h i s chapter was again used.  Analyses were  carried out using three d i f f e r e n t values of the  permeability,  Pore Pressure (psf) (Note •• values for I. sec. curve should be divided by I0<)  FIG. 5-13  Pore  Pressure Distribution at  Various Times  (no  Dissipation)  140 k.  In  each  analysis,  the  value  of the permeability was  assumed to be constant f o r the entire depth Even  though  the  permeability capability  computer  values was  not  program  that  vary  used  since  from  a  matter  of  the  set  layer  deposit.  up  to  accept  to  layer  the  the purpose was to obtain a  general idea on the influence of As  was  of  pore  pressure  dissipation.  f a c t , such investigation i s very l i k e l y the  f i r s t of i t s kind carried out to date(1975). Terzaghi permeability  and  values  encountered i n analyse three  Peck(l96?)  for  engineering different  approximately  types  practice.  cases  0.000003, 0.0003 and 0.03 these  various  listed  feet  with per  correspond  a  of  soil  It  was  to  of  that  are  decided  permeability second  range  values of  respectively  and  impervious  fine  nearly  sands, medium permeable sands and very permeable clean From  here  onwards  these  these  cases  are  sands.  three cases w i l l be referred to as  case D l , D2 and D3 r e s p e c t i v e l y . for  to  shown  Pore pressure in  Figure  distribution  5-l4a,  b  and  c  respectively. Without pore pressure liquefied  at  around  t h i s chapter.  When  was  soil  layer  dissipation  of  pore  water  pressure  allowed but with a very low value of permeability, as  of  liquefaction  agrees with the  11  8 second as described i n section 5«3 of  case Dl, l i q u e f a c t i o n occurred at zone  dissipation,  was  boundary  the  12th  s h i f t e d upwards. conditions  of  the  layer,  in  i . e . the  This behaviour problem  since  Pore 0  Pressure IOOO.  (psf) 2000.  Pore  Pressure  0  1000.  \  \• I  50.  •  / / •  100  100.  I50.h  I50.h  (psf) 2000.  •  I /  •  Pore Pressure .0  1000.  (psf) 2000.  •  • .  / / / /  1/ •  •  /  //  •  • •  3  CD  01 CO CD O CO  200  200.  6  (b)  (a)  FIG. 5-14  200.  Pore Pressure  Distribution  (c)  for Different  k Values  142 drainage  can only  take place at the top surface.  When the  permeability value of the s o i l was made larger, as i n case D2, liquefaction  d i d not occur  distributions second  respectively.  similar  pore  pressure  occurred at much l a t e r time, as shown by the 815-second  and  though  curves  in  Figures  5-l4a  and  When the permeability value was made 0.03 feet  per second the pore pressure f o r s o i l s t r a t a situated depths  0  deposit  between  and 100 feet a c t u a l l y decreased during the i n t e r v a l  of time 8 to 15 seconds. statements  b  could  be  In general  terms  the following  made i n r e l a t i o n to the p a r t i c u l a r s o i l  analysed.  1. Relatively low permeability values per  second)  tends  to s h i f t  (around  0.000003  feet  the zone of l i q u e f a c t i o n  upwards and d i s s i p a t i o n of pore pressure does not improve the  resistance  to l i q u e f a c t i o n  failure  of the  soil  deposit to any s i g n i f i c a n t extent. 2. Medium to large permeability value (0.0003 to 0.03 feet per second)  Improves  the resistance to l i q u e f a c t i o n f a i l u r e  by d i s s i p a t i n g the high pore pressure pressure  may  decrease  instead  generated.  Pore  increase a f t e r the peak  disturbance i s passed even though there are t r a i l i n g weak disturbances fed to the deposit. that  I t may  be  visualized  i f the earthquake motion consists of strong pulses  separated by weak ones the maximum pore pressure near the surface would be even smaller. 3. Within the time i n t e r v a l of 15 seconds  the magnitude and  143  d i s t r i b u t i o n of pore pressure below the 100-foot depth i s l i t t l e affected, i f any, by the d i s s i p a t i o n process. the  earthquake  and  uniform  duration,  waves consist of a r e l a t i v e l y low-valued  amplitude  say  If  60  acceleration  waves  with  long  to 80 seconds, l i q u e f a c t i o n of deeper  zones i s possible. Quantitative strictly  applicable  section 5*1. terms  of  the  r e s u l t s obtained i n t h i s to  the  chapter are  p a r t i c u l a r deposit described i n  Qualitative extrapolation to other deposits i n effects  of some of the s o i l parameters on the  l i q u e f a c t i o n process can be made.  In addition, the numerical  solution of the coupled problem Is shown to be practicable well  as  feasible.  Computer time f o r numerical integration  for any one case described i n t h i s for  a  as  section,  i . e . integrating  time history of 15 seconds at 0.01 second i n t e r v a l , i s  approximately 15 seconds CPU time.  144  CHAPTER 6 SUMMARY AND CONCLUSIONS /  6.1  SUMMARY The problem of s o i l l i q u e f a c t i o n has  concern  been  a  main  i n the design of foundations f o r earthquake-resistant  structures  on  saturated  granular  soils.  Previous  investigations indicated that f o r horizontal s o i l deposits the important  factor  accumulated  i s the, amount  of pore  during the period of shaking.  water  pressure  In t h i s thesis i t  i s assumed that f o r the saturated granular s o i l comprising the deposit,  l i q u e f a c t i o n occurs when the e f f e c t i v e normal stress  reaches a value of 0.05 of i t s model  initial  value.  A material  i s proposed whereby such changes i n pore water pressure  can be calculated by following the s t r a i n history element  subjected to shear deformation.  assumed to be of i n f i n i t e on  of a  soil  The s o i l deposit i s  extent and f i n i t e thichness  resting  a firm bedrock i n order that s i m p l i f i c a t i o n of the problem  to one of the shear-beam type derived  i s possible.  Equations are  to describe the coupled processes involved during the  progress of l i q u e f a c t i o n , namely, 1. Dynamic response to base motion, 2. Increase i n volumetric straining,  compaction  strain  due to  cyclic  145  3. Pore water pressure generation as a r e s u l t of (2),  and,  4. Dissipation of excess pore water pressure(consolidation). Solution  of  the  numerical i n t e g r a t i o n .  problem  was  achieved  a  lumped-mass-system  integration method dynamical  using  The f i n i t e difference method was used  to a r r i v e at an approximation of the continuous by  by  was  equation.  soil  i n a r a t i o n a l way. then  used  to  deposit  A step-by-step  solve  the  resulting  A separate f i n i t e difference scheme was  used to solve the consolidation equation. Several example problems were solved f o r an  assumed  deposit by varying some of the s o i l parameters so that factors a f f e c t i n g the l i q u e f a c t i o n process can be studied.  6.2  CONCLUSIONS  Conclusions  reached  in  this  research  are  as  follows:1. The hyperbolic h y s t e r e t i c s t r e s s - s t r a i n volumetric  compaction  parameter,  can  be  'strain  used  serving  effectively  relationship  with  as a 'hardening' to  describe  the  l i q u e f a c t i o n behavior of s o i l under dynamic loading. 2. Using  the  proposed  s o i l model f o r a deposit a r e l a t i v e l y  • f l a t ' frequency response i s obtained. that  the  resonant  capacity  This  indicates  of the non-linear model i s  lower than i t s l i n e a r equivalent. 3. When pore water pressure accumulated during shaking reaches  146 A value of about 90.% of stress,  very  large  the  strains  initial  account  doubt  on  of  equivalent  but without  pore pressure r i s e .  the v a l i d i t y linear  of  analysis  concerning l i q u e f a c t i o n  normal  develop In comparison with  those computed f o r the same deposit into  effective  using  This throws some strain  to  taking  compatible  obtain  information  potential.  For r e l a t i v e l y low permeability values (0.000003 to feet  per second)  i n the deposit  a l o c a l i z e d type of  l i q u e f a c t i o n i s produced so that s t i f f n e s s deterioration  and  strength  due to pore pressure r i s e f o r a p a r t i c u l a r  layer occur at a much faster rate than others. usual  case  strata zone  when  In the  the r e l a t i v e densities of the surface  are lower than of  0.0003  those  l i q u e f a c t i o n occurs  of near  the deeper strata the the surface,  i.e.  at depth between 20.0 to 50.0 f e e t . R e l a t i v e l y large permeability values (around 0.03 feet second)  reduce  the tendency  strength deterioration near  the  surface  per  to l o c a l i z e s t i f f n e s s and  to a single layer. the maximum  pore  For water  strata pressure  accumulated during the time of earthquake distrubance, i s reduced and larger magnitude base motion  i s require  to  In general, the r a t i o of the pore pressure accumulated  to  cause  liquefaction.  the i n i t i a l e f f e c t i v e normal stress decreases with after  the  50.0  feet  mark.  depth  For s o l i deposits without  e s p e c i a l l y weak layers l i q u e f a c t i o n at depth greater than  14? 100.0 feet i s very u n l i k e l y . For s o i l s t r a t r a located  50.0  feet  or more  beneath the  ground surface, the pore pressure accumulated i s not much affected  by  the  dissipation  process even  permeability value i s as high as 0.03 Therefore,  when  a  permeable  feet  when  the  per second.  s o i l deposit i s shaken by  earthquake motions consisting of r e l a t i v e l y low, uniform amplitude  and  long  duration  acceleration  waves,  l i q u e f a c t i o n at larger depths ( i . e . 50.0 and 100.0 is  3 In  feet)  possible.  SUGGESTIONS FOR FURTHER RESEARCH view  of  the p a r t i a l s t a b i l i z a t i o n e f f e c t produced by  sandy material confining  due  to  i t s dllatant  pressures, incorporation  behavior  at low  of dilatancy  into 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 s needed so that the analysis may be continued near and at l i q u e f a c t i o n (zero  or  near  zero e f f e c t i v e s t r e s s ) . A  thorough  investigation  of  the volume  c h a r a c t e r i s t i c s i s necessary i n order compaction function 1. a  continuous  can be generalised  that  compaction the volume  to Include,  v a r i a t i o n of volumetric compaction as a  function of s t r a i n or stress, and, 2. e f f e c t of a general state of s t r e s s . A correlation  study with f i e l d data should be c a r r i e d out.  148  LIST OF REFERENCES  Ambraseys, N., and Sarma, S., I969. "Liquefaction of S o i l s Induced by Earthquakes." B u l l , of the Seism. Soc. of Amer., v o l . 59» no. 2, A p r i l 1969,  PP.  651-664.  Barkan, D.D., I 9 6 2 , "Dynamics of Bases and Foundation." McGraw H i l l Book Co., Inc., I962. Bathe, K.J., and Wilson, E.L., I973t " S t a b i l i t y and Accuracy Analysis of Direct Integration Methods." I n t e r n l . J . Earthquake Engg. and Struct. Dyn., v o l . 1, 1973.  PP.  283-291.  Bazant, Z., and Dvorak, A., 1965. "Effects of Vibrations on Sand and the Measurement of Dynamic Properties." Proc. 6th I n t e r n l . Conf. S o i l Mech. & Found. En*rg., v o l 1, 1965. PP. 161-164. Blot, M.A., 1956, "Theory of Propagation of E l a s t i c Waves i n Fluid-Saturated Porous S o l i d . I - Low-Frequency Range." J . Acoustical Soc. of Amer., v o l . 28, no. 2, March 1956, pp. 168-178.  -  Bogdanoff, J.L., Goldberg, J.E., and Bernard, M.L., I 9 6 I , "Response of a Simple Structure to a Random Earthquake Type Disturbance." B u l l . Seism. Soc. Amer., v o l . 511 no. 2, A p r i l 1961,  PP. 293-310.  Casagrande, A., 1936, "Characteristics of Cohesionless S o i l s A f f e c t i n g the S t a b i l i t y of Slopes and Earthf i l l s . J . Boston Soc. C i v i l Engrs., Jan. 1936. 1 1  149  Castro, G., 19^9. "Liquefaction of Sands." Ph.D. Thesis, Harvard Univ., Cam., Mass., I969. E d i t o r i a l Committee of the "General Report.", I968, "General Report on the Niigata Earthquake of 1964." Finn, W.D.L., 1972, " S o i l Dynamics - Liquefaction of Sands." Proc. I n t e r n l . Conf. Microzonation f o r Safer Construction Research arid Application, Seattle, Washington, USA, Oct.,  1972, pp. 87-111.  Finn, W.D.L., Emery, J . J . , and Gupta, Y.P., 1970, "A Shake Table Study of the Liquefaction of Saturated Sands during Earthquakes.", Proc. 3rd European Sym. Earthquake Engg., Sept. 1970, PP. 253-262. Finn, W.D.L., Emery, J . J . , and Gupta, Y.P., 1971, "Liquefaction of Large Samples of Saturated Sand on a Shake Table.", Proc. 1st Can. Conf. Earthquake Engg., Van., UBC, May  1971, PP. 97-110.  Finn, W.D.L., Emery, J . J . , and Gupta, Y.P., 1971a, " S o i l Liquefaction Studies Using a Shake Table." Closed Loop, System Corporation, Fall/Winter, 1971. Finn, W.D.L., Pickering, D.J., and Bransby, P.L., 1971. "Sand Liquefaction i n T r i a x i a l and Simple Shear Tests." J . S o i l Mech. & Found. Div. ASCE, v o l . 97, SM 4, A p r i l  1971. PP. 639-659.  t  F l o r i n , V.A., and Ivanov, P.L., 1961, "Liquefaction of Saturated Sandy S o i l s . " Proc. 5th I n t e r n l . Conf. on S o i l Mech. and Found. Engg.,  v o l . 1, 1961, pp. 107-111.  Forssblad, L., I967, "New Method f o r Laboratory S o i l Compaction by Vibration." Highway Research Board, no. 177, 1967, pp. 219-224.  150  Hardin,. B.O., and Dmevich, V.P., 1972, "Shear Modulus and Damping i n S o i l s : Measurement and Parameter E f f e c t s . " J . S o i l Mech. & Found. Div., ASCE, v o l . 98, SM 6, June 1972, pp. 603-624. Hardin, B.O., and Drnevich, V.P., 1972a, "Shear Modulus and Damping i n S o i l s : Design Equations and Curve s." J . S o i l Mech. & Founds. Div., ASCE, v o l . 98, SM 7, July  1972, pp. 667-692.  Hazen, A., 1920, "Hydraulic F i l l e d Dams." ASCE Transaction, v o l . 83, 1920, pp. I9I7-I945. Herrera, I., I 9 6 5 ,  "Modelo Dinamicos para Materiales y Estructuras d e l Tipo Masing." 'Boletin Sociedad Mexicana de Ingenieria Sismica, 3(1), 1965, PP. 1-8.  Ishibashi, I., and Sherif, M.A., 1974, " S o i l Liquefaction by Torshional Simple Shear Device." J . Geot. Engg. Div., ASCE, v o l . 100, no. GT 7, July, 1974,  pp. 871-904.  Ishihara, K., and L i , S.L., 1972, "Liquefaction of Saturated Sand i n T r i a x i a l Torsion Shear Test." S o i l s and Foundation, v o l . 12, no. 3. June 1972, pp. I939. Iwan, W.D., I 9 6 6 , "A Distributed-Element Model f o r Hysteresis and i t s Steady-State Dynamic Response." J . Appl. Mech., v o l . 33, no. 4, I966, pp. 893-900. Jennings, P.L., 1964, "Periodic Response of a General Yielding Structure." J . Engg. Mech. Div., ASCE, v o l . 90, EM~ 2, A p r i l 1964  pp. 131-195.  151 Kishlda, H., 1964, "Damage to Reinforced Concrete Buildings i n Niigata City with Special Reference to Foundation Engineering." S o i l and Foundation, v o l . 7, no. 1, March 1964, pp.71-88. Kishida, H., 1969, "Characteristics of Liquefied Sands during Mino-Owari, Tohnankal and Fukui Earthquakes." S o i l s and Foundations, v o l . 9, no. 1, March 1969t  PP. 75-92.  Kishlda, H., 1970, "Characteristics of Liquefaction of Level Sandy Ground during the.Tokachioki Earthquake." S o i l s and Foundations, v o l . 10, no. 2, June 1970,  pp. 103-111.  Koizumi, Y., 1966, "Changes i n Density of Sand Subsoil Caused by the Niigata Earthquake." S o i l and Foundation, v o l . 6, no. 2, June 1966, pp. 38-44. Kondner, R.L., and Zelasko, J.S., I963. "A Hyperbolic Stress-Strain Formulation f o r Sands." Proc. 2nd Pan Amer. Conf. S o i l Mech. & Found. Engg., 1963, pp. 289-324. Marsal, R.J., 1963. Internal Report, Institute of Engineering, National University of Mexico, Mexico, 1963. Martin, G.R., Finn, W.D.L., and Seed, H.B., 1975, "Some Fundamental Aspects of Liquefaction under C y c l i c Loading." J . Geot. Engg. Div., ASCE, v o l . 101, no. GT 5,  May 1975, PP. 423-438.  Maslov, N.N., 1957, "Questions of Seismic S t a b i l i t y of Submerged Sandy Foundations and Structures." Proc. 4th I n t e r n l . Conf. S o i l Mech. & Found. Engg.,  v o l . 1, 1957, PP. 367-372.  152  Maslng, G., 1926, "Eigenspannlngen und Verfestigung beim Messing." Proc. 2nd Internl Cong. Applied Mech., 1926,  PP. 332-335.  Mogami, R., and Kubo, 1953, "The Behaviour of S o i l during Vibration." Proc. 5th I n t e r n l . Conf. S o i l Mech. and Found. Engg., v o l . 1, Zurich, 1953, PP. 182-185. Newmark, N.M., 1959 "A Method of Computation for S t r u c t u r a l Dynamics." J . Engg. Mech. Div., ASCE, v o l . 85, no. EM3, July 1959. Newmark,.N.M., and Rosenblueth, E., 1971, "Fundamentals of Earthquake Engineering." Prentic H a l l , Inc., Eaglewood C l i f f s , N.J., USA,  1971. PP. 162-163.  Osaki, Y., 1968, "Building Damage and S o i l Condition." General Report on the Niigata Earthquake of 1964. Committe of "General Report", I968, pp. 355-383. Parmelee, R.A., Perelman, D.S., Lee, S.L., and Keer, L.M.,  1968,  "Seismic Response of Structure Foundation System." J . Engg. Mech. Div., ASCE, v o l . 94, EM 6,  Dec. 1968, pp. 1295-1315.  Peacock, W.H., and Seed, H.B., I968, "Sand Liquefaction under C y c l i c Loading Shear Conditions." J . S o i l Mech. & Found. Div., ASCE, v o l . 94, SM 3, May 1968, pp. 689-708.  Pyke, R.M., 1973. "Settlement and Liquefaction of Sands under M u l t i D i r e c t i o n a l Loading." Ph.D. Thesis, Univ. of C a l i f . , Berkeley, C a l i f . , 1973, P. 280.  153 Schnabel, P.B., Lysmer, J . , and Seed, H.B., 1972, "A Computer Program for Earthquake Response Analysis of Horizontally Layered S i t e s . " Earthquake Engineering Research Centre, Report Number EERC 72-12, Dec. 1972. Seed, H.B., I968, "Landslides during Earthquakes due to Liquefaction." J . S o i l Mech. & Found. Div., ASCE, v o l . 94, SM 5, Sept. 1968, pp. 1053-1122. Seed, H.B., 1972, "Dams and S o i l s . " The San Fernando Earthquake of February 9, 1971, and Public P o l i c y . Special Subcommittee of the Joint Committee on Seismic Safety, C a l i f o r n i a Legislature, July 1972. Seed, H.B., and I d r i s s , I.M., I 9 6 7 , "Analysis of S o i l Liquefaction : Niigata Earthquake." J . S o i l Mech. & Found. Div., ASCE, v o l . 93, SM 3, May  1967, PP. 83-108.  Seed, H.B., and I d r i s s , I.M., 19^9, "Influence of S o i l Conditions on Ground Motions during Earthquakes." J . S o i l Mech. & Found. Div., ASCE, v o l . 95, SM 1, Jan. I969, PP. 99-137. Seed, H.B., and I d r i s s , I.M., 1971, "Simplified Procedures f o r Evaluating S o i l Liquefaction Potential." J . S o i l Mech. & Found. Div., ASCE, v o l . 97, SM 9, Sept. 1971, pp. 1249-1273. Seed, H.B., and Lee, K.L., I 9 6 6 , "Liquefaction of Saturated Sands during Cyclic Loading." J . S o i l Mech. & Found. Div., ASCE, v o l . 92, SM 6, Nov. 1966, pp. 105-134. Seed, H.B., and S i l v e r , M.L., 1972, "Settlement of Dry Sands during Earthquakes." J . S o i l Mech. & Found. Div., ASCE, v o l . 98, SM 4, A p r i l 1972, pp. 381-397.  154  S i l v e r , M.L., and "Volume Changes J . S o i l Mech. & Sept. 1971, PP.  Seed, H.B., 1971. i n Sands during C y c l i c Loading." Found. Div., ASCE, v o l . 97. SM 9. 1171-1182.  Sneddon, I.N., 1957. "Elements of P a r t i a l D i f f e r e n t i a l Equations." McGraw H i l l Book Company, Inc., 1957. State Earthquake Investigation Commission, 1908, "The C a l i f o r n i a Earthquake of A p r i l 18, 1906." Report of the State Earthquake Investigation Commission, Published by the Carnegie Institute of Washington, 1908. Terzaghi, K., and Peck, R.B., I968, " S o i l Mechanics i n Engineering P r a c t i c e s . " John Wiley Publishing Company, F i r s t Corrected P r i n t i n g , 1968. Yen,  B.C.,  I967,  "Viscosity of Saturated Sand near Liquefaction." Proc. I n t e r n l . Sym. Wave Propag. & Dyn. Properties of Earth Materials, New Mexico, I967. PP. 677-888. Yoshimi, Y., 1967. "An Experimental Study of Liquefaction of Saturated Sands." S o l i and Foundation, v o l . 7. no. 2, June 1967, PP. 20-32.  155  APPENDIX. I LUMPED MASS SYSTEM APPROXIMATION FOR NON-LINEAR SHEAR BEAMS The equilibrium equation, equation i n chapter 4 i s re-written here f o r  -  s$  (4-1),  obtained  convenience,  8  =0  When the s o i l parameters are functions of depth, (H-z), where H  i s the t o t a l  depth  of  the s o i l deposit and  z  i s the  height from the surface of the bedrock, equation (4-2)  has to  be replaced by a more general equation, namely,  and also, Equations  T  =  f^JT.z)  f  =  j>(z)  ( 1 - 2 ) and (1-3)  +  (1-2)  f (f,z) 2  (1-3) can be substituted  into  (1-1)  to  give,  Consider, intervals,  f o r example, so  that  the  z-domain  O^^z^z^  i s sub-divided into  . , . <z =H , n  and  length of each i n t e r v a l be denoted by h^, i . e . ,  l e t the  156 Also, j> f  l e t x., 1  at  f. and /  p. represent  1  the values of x,  i , e  t and  c  ,. and denote the p a r t i c u l a r  2,i'  function  at  z^  by  ' f  The  1  2,i  central  =  f  2 i' i> (f  ^- )  z  difference  6  quotient  for  Tz^2^' ^  a  z  ^  z  i^  s  given by,  A2 H f  (ir  (f  2.1 i  "  +  ,i-i>/<  f 2  i h  i l * i> +  h  +  and a s i m i l a r expression can be written f o r  f^(f,z)  when  into  these  expressions  are  following equation can be i(h  Now,  i+1  +h ).^ .x i  1  i  substituted  so  that  (1-3)• the  obtained,  +  (f^..^ - f  +  <2,i-i " 2,i+*> f  1 § i + i  )  f  =  0  tt- ) 8  write,  t (f,z) z  tM,z) -L-  =  .t  (1-9)  3x dz  Since f  i-Jr  (  Consequently, equations  x  i * i-l x  )  /  h  i  (1-9) and (I-10) l)  fo(f? 1  1-2  I f a new symbol i s defined as follows :-  ( I  give us,  -  1 0 )  157  i  (1-12)  I  h. . f.  ( I - 1 1 ) can be w r i t t e n i n a simple  equation  f  2,i-i  *'  i '  k  (  i  x  "  i - i  x  form, (1-13)  }  S i m i l a r l y , i t can be shown t h a t ,  f  2.1+4  F  '  k  (1-14)  i+l*( i+l " i ^ x  x  (1-15)  ^i'( i ~ i-l^ x  f  i.i+i  c  '  ,=  x  i+l*( i+l x  i)  A.Z.  where,  c  Equations equation  (1-17)  i  (1-13)  (1-16)  i ^  x  1  (I-l6)  through  ( 1 - 8 ) and  1-2  the  can  following  be  substituted  'discretized•  into  dynamic  equation i s obtained, X  V  x  i  +  c +c  (-^  i  -c  i + 1  i + 1  ).^  i-1  X.. X  >  i 1 +  *i-l +  (" i i k  k  +  i l  k  - i l) k  +  X  +  x  Note  that, i(h  i + 1  +h ).j> i  i  i i+l  *  =  0  (1-18)  158 and since the average density within the i n t e r v a l i s usually known  instead  of  the density  value  at  z^, m^ i s "better  calculated by, V  *\fi-i- i  =  h  ^/i+i^i+l  +  I t should be pointed out that the f i n i t e difference (1-18)  holds  i n general f o r i=2,3,4,... , n - l .  observing that when i = l the variables x  For i = l and  and x  Q  i n equation (1-18) are i n fact the v e l o c i t y and of  the base  Q  be  used  appearing  displacement  respectively so that the symbols x^ and x^ are  used instead. subscript  1 9 )  equation  i=n, i . e . at both end points, the same equation can by  "  ( I  And,  of n+l  when  i=n  the variables  are set i d e n t i c a l l y  equal  with  the  to zero.  system of n equations thus obtained can be represented  by  A a  simple matrix equation![M].{'x} where x  i'  x  {x}, i  a  n  d  +  [C].{x}  +  [K].{x>  =  {P > b  {xj- and {x} are column vectors whose components x  i Sive  us  the  acceleration,  velocity  displacement respectively of the point at z^, {x}  =  (x  x  £  x^  .  .  . x )  T  x  {x}  =  (x  x  2  x^  .  .  . x )  T  x  {x}  =  (x  x  2  x  .  .  . x )  T  x  {P } = b  (1-20)  (k x +c x 1  b  1  b  3  0  0  n  n  n  .  .  . 0)  T  and  159 where  the s u p e r s c r i p t  T  r e f e r s t o the transpose o f the row  v e c t o r i n t o a column v e c t o r and x^ and x^ r e f e r t o v e l o c i t y and displacement The  matrices  the  base  respectively. [M],  [C]  and  [K] as o b t a i n e d from  e q u a t i o n (1-20) are as f o l l o w s  / [M] =  m  0  1  0  m  0  •  0 '  •  0  0  2  0  0  m^  •  •  •  0  0  0  ,  0  •  (I -21)  •  *  •  m  n j  [C] =  c  c  1 +  " 2 0 C  [K] =  ~2 C  0  2 3  " 3  2  C  + C  - 3  c  c  c  0  •  •  0  0  •  •  0  •  •  0  •  •  •  •  - n  3 4 + c  •  *  •  •  0  0  0  0  0  0  •  •  0  -k  0  t  •  0  •  •  0  •  •  •  " l  + k  k  " 2 0 k  2 k  2 +  k  -k  3  3  3  k +k 3  4  •  •  •  •  0  0  0  0  C  c  (I -22)  n -  k  vector  displacement  vector,  i s preferred. {X)  and  {X}  L e t {x} be  the  -23)  n  F o r base a c c e l e r a t i o n type o f l o a d i n g the displacement  (I  relative  be the r e l a t i v e corresponding  160 velocity  and acceleration vectors and l e t x^ be the variable  describing the displacement of the base  as  a  function of  time, i t can be seen that  {x}  Similar  expressions can be  acceleration. equation  (1-24)  When  written  equation  f o r the v e l o c i t y  and  i s substituted  into  (1-24)  (1-20), i t can be shown that the following equation  i s obtained,  [Ml{X)  +  [ C ] . { X )  kx  [K].{X}  +  m  I t i s quite apparent obtained  that  the f i n i t e  (1-25)  b  n  difference  equation  i s the same as that f o r the dynamic response  system of n dashpots.  discrete Thus,  a  masses lumped  connected mass  by  model  springs  can be  used  of a and to  approximate the continuous shear beam model f o r the dynamic response  of the s o i l  column discussed i n chapter 4,  The  value of the masses are obtained by lumping h a l f the mass of each  layer  to i t s corresponding nodes, z^ ^ and z^.  s t i f f n e s s e s of the springs and the damping the dashpots are obtained from f relationship, and  (I-l?),  equation  1  and f  2  The  coefficients for  of the s t r e s s - s t r a i n  (4-2), and through the use of (1-12)  Also, the width of the s o i l column can be made  161 d i f f e r e n t from u n i t y and i t s e f f e c t can be taken i n t o account by m u l t i p l y i n g the c o r r e s p o n d i n g width o f l a y e r , b^, i n t o the appropriate c^  and  terras.  So t h a t , i n g e n e r a l ,  the components  m^,.  k^ o f the m a t r i c e s [M], [C] and [K] are g i v e n by the  following i  -  m  c  k  i  *«/>i-*'V i  =  h  £-j>i+£' i+l- i+l  +  b  i = tilh-i'H-tf'*! i  =  f  2 ^ i - i '  Z  i - ^ -  h  / ih-h-tf b  i  '/  ( I  -  2 6 )  <Vi-*>  where i r e f e r s t o the l a y e r number o f the lumped mass system.  162  APPENDIX II NUMERICAL INTEGRATION OF THE DYNAMIC RESPONSE EQUATION The dynamic response equation for. the  lumped  mass  system i s re-written here f o r convenience,  [M].{X} + [C].{X} + [KJ.{X}  =  {P}  (II-l)  where {p} i s the i n e r t i a force vector, •m.  {p}  II.1  -  (II-2)  Incremental S t r e s s - s t r a i n . Let  where  t  T  t  + At,  ( H - 3 )  i s an instant of time at which a l l the quantities  i n equation ( I I - l ) increment.  are known  and  equation  (II-l)  Since  At  is a  small  time  has to he s a t i s f i e d at  every instant, we have,  [M] {X} + [CJ {X} + [K] {X} T  T  T  T  T  T  [M] {x} + [C] {x} + [K] { } t  t  t  t  t  X  t  = (P} =  {P}  where the subscript r e f e r s to the instant of the  (II-5)  T  p a r t i c u l a r quantity takes on i t s value.  (II-6)  t  time  at  which  Since the mass  163 matrix i s a constant matrix, [M] so that  =  T  [M] {X} T  T  [M]  t  =  [ M],  [MJ {x}  -  t  =  t  (X}  [M]  - {X}  T  t  However, the [Cj and [Kj matrices depend on the shear of the s o i l , as a relation,  result  of  ways  non-linear  strain  stress-strain  and thus on the displacement of the masses so that  (II-7) cannot  equations s i m i l a r to that of these  the  (II-7)  be  obtained  matrices i n a s t r a i g h t forward manner. to  arrive  at  expressions  similar  There are two  (II-7) using  to  A crude method i s to replace [K]^  approximation methods.  for  by  [K]^, so that, [K] {X} T  -  T  [K] {X} t  = [K]  t  t  {x}  T  - {x}  1  t  (II-8)  A better method i s as follows »Define,  {F}  so that,  [K] {X}  [K].{X}  T  T  -  (H-9)  [K] {X} t  t  Wt 3F; 5X  = where  [K],  -  {X},  W  t  jj  [ K] .|{X} t  3F^  ax. J-  T  -  ^ >t} X  (11-10) (11-11)  164 and 9).  i s the component of the column vector defined i n ( I I Therefore, i n general, [K] {X} T  -  T  [K] {X> t  =  t  LK] .|{X} t  (11-12)  - {X} |  T  t  where [K]^ can take the value of [K]^ f o r crude approximation or i t can be  evaluated  approximation.  according  Similar  to  procedures  (11-11)  f o r better  can be used to obtain  IF equation (II-6) i s subtracted from equation ( I I 5) and taking note of (II-7) and (11-12), an incremental form for the dynamic response equation can be obtained as, [M]|{X} +  T  -  +  {X} j t  [Kj [{x> t  T  [C] |{x} t  -{x>  t  J  T  -  {X} J t  = {P} - {P} T  T  (11-13)  I t can be seen that the only unknown i n t h i s equation are the acceleration, v e l o c i t y and displacement vectors As  such,  at  time  T.  the numerical procedures proposed by Newmark(1959)  can be applied to the equation so that {X}^,  {x}  T  and { x }  T  can be expressed i n term of {*X} , {X}^ and { x } . t  II.2  t  Step by Step Integration In Newmark's method of step by step integration two  parameters,  o< and  displacement at time  |3 are used T  so  that  can be expressed  the v e l o c i t y and i n terms  of the  acceleration, v e l o c i t y and displacement at time t, and of the  165 unknown  acceleration  at time  T.  The r e l a t i o n s  can be  written as follows :{X}  T  =  {X} + (l-cx).At.{X} + *.At.{X}  {X}  T  = {X} + At.{X}  t  t  t  + (i-/5).(At)?{X}  t  '+ Newmark( 1959)  j8.(At)?{x} and /3=i  proposed that  (11-14)  T  t  (11-15)  T  f o r unconditionally-  stable integration procedure, which i n c i d e n t a l l y to  a  constant  average  acceleration method of i n t e g r a t i o n .  I f «=(l/2) and /6=(l/6) method  as proposed  are used  interval  and /3=^  °<=i  the acceleration  At=T-t {Q)  acceleration  are used.  increment  during  the time  i s represented by {Q}^,. i . e . ,  =  T  the l i n e a r  by Wilson and Clough(1952) i s obtained.  For the present purpose If  corresponds  (x}  T  - {x}  (11-16)  t  and i f the following simplifying symbols are used,  then  {a}  t  =  At.{x}  {b}  t  =  At.{X}  equations  (II-14)  t  (H-17)  t  + 2-.(At)?{X} and  (11-15)  (11-18)  t  can be written into a  simple form, namely, {X} -{X} T  {X}  T  t  {X>  t  =  { a ) + c*.At.{Q}  =  {b} + y5.(At)?{Q}  t  (II-19)  T  t  T  (11-20)  166 Equations (II-16), (11-13)  into  (11-19) and  to give  (11-20)  can be  substituted  a matrix equation involving only one  unknown, {Q} , namely, [D] {Q} t  where  [D]  and  {P}  t  T  =  T  {P}  (11-21)  T  + p. (At) .[K]  =  [M] +  ot.At. [C]t  =  { P } - {P} - [ C ] { a } T  Consequently,  t  t  {Q}  T  t  =  t  - [K] {b) t  t  calculated  T  (11-23) (11-24)  [Dj^P}^,  and {x}^, { x } and { X ) can be T  (11-22)  2  using  equations  (11-16), (11-19) and (11-20). In  the numerical  step  by  step  integration,  a  t y p i c a l sequence of calculations that have to be performed i n one time step, At, i s as follows »1.  Based on { x }  t >  {X} and { X } at time t , (11-17) and t  t  (II-  18) are evaluated,  2.  {a}  t  =  At.{x}  {b}  t  =  At,{x}  t  Based on {x}^ and { x } , {f} t  m  t  {i}  t  where  +  t  t  i.(At) .{x} 2  t  and { f } are calculated by, t  =  CaJ{x}  (11-25)  =  [^]{x}  (11-26)  t  t  167  l/h  =  [a]  -l/h  0  1  l/h  2  2  0  0  0  0  - l / h ^ l/h^  0  0  -lA 3.  With  the {i\  {t}  and  t  evaluated according to  values,  t  the  [C]  appropriate  n  [ K ]  and  t  l/h  n  t  are  expressions i n  appendix I . 4.  [D]^. i s then calculated using equation (11-22).  5.  Using  the  base acceleration value at time T, x  b T  , i t is  possible to evaluate, r -\  T  bT m n  and 6.  {Q}  T  (P}  =  T  {P} - {P} - [CJ {a} T  i s solved  velocity  and  by  t  (11-24)  displacement  t  so  that  t  - [K] {b} t  the  vectors, {X) , T  t  acceleration, {x)  T  and fX}  T  can be found by, {X}  n  {X)  =  {X}  + { a } +o<.At.{Q}  (X}  =  {X}  + {b} + |3.(At) .{Q} .  T  T  7.  t  t  t  The steps 1 through 6 are increment.  T  2  t  repeated  T  f o r the next  time  168 An option i s b u i l t into the analysis to bring the s t r e s s strain  point  I f the  strains  integration rates,  closer  to the actual s t r e s s - s t r a i n curve.  and  strain  determined  are accepted as true s t r a i n s and  the  stress-strain  relationships  can  be T  t C ] . {X} .  T  However,  T  and  used  r e s t o r i n g force, [ K ] . { X ) , T  rates  to and  true  each strain  stress-strain-rate  determined the  in  the  'true'  'true' damping force,  these r e s t o r i n g and damping forces  do not necessarily s a t i s f y the equilibrium equation 5).  An  artificial  'external' force, { P . ™  by the following can be applied to  the  },  system  (II-  defined so  that  equilibrium i s restored, {P  e r r  The  {P_„_  for  the  }  =  {P}  - [M].{X}  T  T  - [C] .{X}  } can be added to {P}m next time increment.  T  i n step  T  - [K] .{X} T  number  CX  .  5  for  I t should be pointed  out  that t h i s method of correction forces i s only good the magnitude of {P-„_  T  where  } i s small as i n the present case  X  of step-by-step integration with small time increment.  169  APPENDIX I I I FINITE DIFFERENCE EQUATIONS FOR PORE PRESSURE DISSIPATION In  the  general  case  where  the permeability and  compressibility of the s o i l are functions of  depth  equation  (4-19) i s used,  *cr  w  3  IF  =  K  ae^  d<r„  k  l r - ^ ^ ' 1 7  )  +  K  lr-"iT  Following the same approach as i n appendix I, i t i s assumed  that  the  z-domain  i s divided into n i n t e r v a l s so  that, 0<z,<z <z l c. } o  .  Q  .  ,<z =H, n  where H i s the t o t a l depth of the deposit.  For convenience  the following notations w i l l be used :hj  to denote the length of the i-th i n t e r v a l , i . e . , h  W  1 ,T  i  =  Z  i" i-1 Z  (IH-2)  to denote the f i n i t e difference approximation to W  i,T  at z=z^ and t=T, or, =  °w< f > z  T  (IH-3)  170 Moreover, attention should be given to the following m  v  k^  1 #  and  ^  which are used to denote the  values at the centre of their  values  €vp^ i^  at  are  z  the  used.  the  i-th  interval.  symbols  respective  To  indicate  'nodes* the symbols m (z ), k(z,) and v 1 i m i s the coefficient of volume v  compressibility so that, 1/K  m  (III-4)  lr  For the f i n i t e difference approximation of equation (III-l)  consider f i r s t the term on the l e f t hand side of the  equation  and  use  a  forward  difference  quotient  for  approximation so that, w St t=T z=z,  acr  ( W  A  central  i.T AT-  W  +  J  difference  i.T  quotient  (III-5)  ) / A T  will  approximation of the f i r s t term on the  be  used  for  the  right  hand  side  of  equation ( I I I - l ) , namely, r  .A( * l r 3z 2f  2£) t=T dz  w  W L  i+1*'  z=z  - W  1+1,T l  J  i+l  V * l+1 (  And  for  the  1,T  W  111  h  +  *  h 1  " 1-1,T h. W  (III-6)  )^ (z ) v  1  second term on the right hand side of equation  ( I I I - l ) a backward f i n i t e difference quotient w i l l  be  used,  17.1  l  3t  r  m ( v  -  T  Z l  ) .AT T  v *  (III-7)  It can be shown that by combining equations ( I I I - l ) ,  (III-5),  ( I I I - 6 ) and (III-7)  of the  a  finite  difference  equation  following form can be obtained,  W  " MVl+l.I  i.T AT +  +  where,  R,  W  1.T  =  1  N  a  n  It  d  W  can  +  i+l/  h  f,T=  be  1.T  seen  U  that  1  I  -  »  8  (III-9)  2  (h  ^  I  h  *w' i+l  -  (  / — i  — —  =  1  A M  "l-l.T * <Vl>-Wi.T]  +  l  +  AT  r . ~  hi)»n (z ) h v  1  (111-10)  1  f  (  1  1  1  -  U  )  equation (III-8) resembles the usual  f i n i t e difference s o l u t i o n of a one dimensional consolidation problem  except  •generation  f o r the  of  1  pore  values of m  and £  instead  the  v  of  v p  last  water  term  which  pressure.  represents the In order that the  at the centre of the layer can  respective values at the nodal  be  points  used the  following average values are used,  ffi  / . v< l> = z  m  v,i+l' H-l * h + h  1 + 1  h  l  m  v,i' l h  (111-12)  172 vp 1T m (z, )_ v IT c  A£  vp,l+1,T mv.i+l  A G  vp,i,' mv . i  (111-13)  Expressions (III-10) and ( I I I - l l ) can be modified  into the  following, 2. AT r  AW g  l.T  The  w- v,i+l- i+l  = 4.  final  form  (m  h  vp.i+l.T mv,i+l  *i  +  m  (111-14)  v,i'V' i h  Srp,i,T  A  +  (111-15)  m  of the equation  v.l  is still  (III-8)  expressions (III-9), (111-14) and (111-15) w i l l the analysis.  be  but  used i n  LISTING  OF LIQUEFACTION  PROGRAM  -  A U G , 75  LIQPG  1  173 APPENDIX IV C**************************************** C*  INPUT  CARDS  FOR LIQUEFACTION  PROGRAM  *  C**************************************** C C  (20AU)  (1  1-  TITLE  CARD)  2.  NMA T , N P T Y P E , N V O L , N L A Y E R , N L I N E L , N D A M P , N B T , N B B , N C R V WHAT = NUMBER OF MATERIAL TYPES.  C C C  NPTYPE  C C C C C C C C C C C C C C C C C  NVOL NLAYER NLINEL NDAMP NBT NBB NCRV  .  . (20IU)  (1  C  =  PRO B K E M T Y P E NUMBER, 1 FOR DYNAMIC RESPONSE ONLY, 2 FOR DYN. RSSP. PLUS PORE PRESS. R I S E OR V O L . CHANG 3' F O R D Y N . R E S P . PLUS PORE PRESS. PLUS DISSIPATION. = 0 FOR PORE PRESSURE RISE. = 1 FOR VOLUME CHANGE (SETTLEMENT). = NUMBER OF L A Y E R S . = 0 FOR NON-LINEAR MATERIAL (CHANGE MODULI) = 1 FOR LINEAR MATERIAL (CONSTANT MODULI) = 1 FOR DAMPING C O E F F I C I E N T S INPUT. = 2 FOR DAMPING PROPPORTIONAL TO MASS AND S T I F F N E S S . = 0 FOR ZERO PORE PRESSURE AT T O P BOUNDARY. =1 F O R ZERO G R A D I E N T AT T O P BOUNDARY. = 0 FOR ZERO PORE PRESSURE AT BOTTOM BOUNDARY, ='1 F O R ZERO G R A D I E N T AT BOTTOM BOUNDARY. = 0 F O R N O T C A L L I N G S U B R O U T I N E TO C H E C K S T R A I N REVERSAL = 1 F O R C A L L I N G S U B R O U T I N E TO CHECK FOR S T R A I N REVERSAL  C C  3.  C  MATYP(I),  NSUBD(I),  H(I) ,  ( 2 1 4 , 2 F 1 0 . U) ,  WIDTH (I)  (NLAYER  C C C C C  U.  C  MATYP(I) NSUBD(I) H(I) WIDTH(I) A1 ( I ) , A 2 ( I ) D1 ( I ) , D 2 ( I )  = M A T E R I A L T Y P E NUMBER O F I = NUMBER O F S U B D I V I S I O N OF = THICKNESS OF I - T H LAYER. = WIDTH OF I - T H L A Y E R . , A 3 (I) , B 1 (I) , B 2(I) , B 3(I) , C 1 ,DEN-H (I) , DEN-V (I) ,PERM (I) ,  - T H LAYER. I - T H LAYER  FOR PP  CARDS) DISSIPAT  (I) , C 2 (I) ,C3 (I) ,C4 (I) , ALFA (I) ,BETA (I) .  C (6F12.4)  C C C  5 .  NEQ,INTYP,NC,NCPR,NCPRM,NPCON,NCONT  R  NPLD,DT  —  ( 3 * N M AT (5I4,F10.3)  CARDS (1 C A  NEQ  = 0 = 1 = 2  FOR CONSTANT AMPLITUDE HARMONIC INPUT. F O R EARTHQUAKE RECORD INPUT. F O R G R A D U A L L Y I N C R E A S I N G A M P L I T U D E TO T Y P E  C C  INTYP  = =  0 1  FOR WILSON AND CLOUGH'S METHOD OF I N T E G R A T I O N . FOR NEWMARK'S METHOD O F I N T E G R A T I O N .  C C  NC NCPP  = =  NUMBER NUMBER  C C  NCPRM NPCON  = =  NUMBER OF INTEGRATIONS F O R P R I N T I N G MAXIMUM V AL. 1 F O R P R I N T I N G OUT F I N A L V A L U E S FOR C O N T I N U I N G NEXT  C C C  "  = 0  C C C C C  C  OF INTEGRATIONS OF INTEGRATIONS  FOR EACH  PRINT OUT  FOR NO P R I N T I N G OUT INITIAL  1.  NCONT  =  0  F O R ZERO  NPLD  = = =  1 0 1  FOR NON-ZERO I N I T I A L CONDITIONS. F O R NO O U T P U T ON U N I T 7 F O R P L O T . F O R OUTPUT ON U N I T 7 FOR P L O T .  CONDITIONS.  T  L I S T I N G OF LIQUEFACTION  PROGRAM - AUG, 75  LIQPG  2  174 C 6. FOR NEQ=1, NCARD,NREC,NFTS,RDFR C (2014) (1 CARD) C NCARD = NUMBER OF EQ REC CARDS. C NREC = NUMBER OF EQ REC PER CARD. C NFTS = 0 GOR NUMBER OF GEE'S INPUT. C . RDFR = FACTOR TO BE MULTIPLIED TO ACC RECORDS. C ' =1 FOR FPS UNIT INPUT. C FMT C (20A4) C FORMAT FOR EQ REC INPUT. C NEQ=0 C NEQ=0, NC, NFTS,AMP,OMEGA C (2014) C NC NUMBER OF SINE CYCLES. C NCFTS C NFTS = 0 FOR NUMBER OF GEE'S EQ INPUT C =1 FOR FPS UNIT INPUT C AMP = AMPLITUDE OF SINUSOIDAL INPUT C OMEGA - CRICULAR FREQ. OF S I N . INPUT. C 7. T I T L E C (20A4) (1 CARD) C DESCRIPTION FOR BASE MOTION INPUT. C C 8. NEQ=0 S NCONT=0 THAT IS ALL NO MORE CARDS REQD. C NEQ=0 & NCONT=1 I N I T I A L VALUES INPUT ACCORDING TO TM1VAL. C NEQ=1 & NCONT=0 INPUT EARTHQUAKE RECORDS. C NEQ= 1 5 NCONT=1 INPUT INITIAL. VALUES AND EQ. REC. C** DYNAMIC ANALYSIS FOR LIQUEFACTION OF A HORIZONTAL SOIL DEPOSIT C SIX TYPES OF PROBLEMS CAN BE SOLVED, NAMELY: C 1. DYNAMIC RESPONSE OF SOIL DEPOSIT ONLY C 2. DYNAMIC RESPONSE WITH PORE PRESSURE GENERATION C 3 . DYNAMIC RESPONSE WITH PORE PRESSURE GENERATION AND DISSIPATION C AND FOR EACH OF THE ABOVE EITHER C I) LINEAR CASE (NLINEL= 1) , I . E . STIFFNESS ETC I S CONSTANT C I I ) NON-LINEAR CASE (NLINEL=0) C SO NPTYPE=1 CORRESPONDS TO TYPE ONE PROBLEM, ETC. C LCPPG IS USED TO CONTROL CALLING OF PORE PRESS GEN SUB-PROG. C LPPDIS I S USED TO CONTROL CALLING OF DISSIPATION SUB-PROG. C LPLOT I S USED TO CONTROL OUTPUT FOR PLOTTING PROGRAM. C INTYP IS USED TO INDICATE TYPE OF INTEGRATION SCHEME USED: C INTYP=0 WILSON AND CLOITGH'S METHOD IS USED C INTYP=1 NEWMARK'S METHOD OF INTEGRATION IS USED. C LDOF1 IS USED TO INDICATE ONE-DEGREE OF FREEDOM C LOGICAL* 1 LVDIR ( 2 0 ) ,LREV ( 2 0 ) ,LRG ( 2 0) ,DODEVP ( 2 0 ) , L L I Q ( 2 0 ) , LREAD,LCPPG,LPPDIS,LPLOT,LDOF1,LLIN,LSTOP,LMCLG, LJRV,LCKF3,LIT,LPMAX COMMON /NCONTL/N,NDAMP,INCR,INTYP,NVOL,NBT,NBB . /LCONTL/LDOF1,LLIN,LSTOP,LMCLG,LJRV,LCKFB,LIT,LPMAX . /LOGCOM/LVDIR,LR EV,LRG,DODEVP,LLIQ . / P R O P T Y / T I T L E ( 2 0 ) , P R O P ( 2 0 , 1 7 ) , H ( 2 0 ) ,DEPTH ( 2 0 ) ,WIDTH ( 2 0 ) ,WT ( 2 0 ) , M A T Y P ( 2 0 ) ,SIGZ ( 2 0 ) ,CC ( 2 0 ) ,GG ( 2 0 ) , NSUBD ( 2 0 ) , GAMR ( 2 0 ) ,  L I S T I N G OF LIQUEFACTION PROGRAM  - AUG,  75  LIQPG  3  175 TAUR(20) ,CM(40) ,SK(40) ,GMN(20) ,TMN(20) /DYNEQT/DYM (40) /PEP (20) ,FORER(20) /CALV AL/DIRN (20) ,XD2 (20) , ABXD2 (20) ,XD1 (20) ,XD0 (20) ,BD(20) , BA (20) ,BB (20) , GD2 (20) , GD 1 (2 0) , G A MM A (20) , G A AO (20) , GAMMAX (2 0) ,TAXY(20) ,EVP (20) ,SK1R (20) , DPP (20) ,PP(20) /CPVALU/BAP(20) ,BBP(20) ,XD2P (20),XD1P (20) ,XD0P (20) ,GD2P (20) , GD1P (20) ,TAXYP (20) /REV PA R/NRPT (20) ,GGAMR(20,8) ,TTAUR(20,8) ,DIRNR(20,8) DIMENSION FMT (20) , AE (10) , DTJMIN (17) DATA LREAD/.TRUE./,LCPPG/.FALSE./, LPPDIS/. FALSE./,LPLOT/. F A L S E . / DATA TNOT/0.0/,UG1/0.0/,NCPRC/0/,NRECC/0/,NPRMC/0/  MM  c 1  2  9001  1003 1004 7  C** C C c** 1001 8  READ(5,1) T I T L E FORMAT(20 AU) WRITE (6,2) T I T L E FORMAT(1H1,10X,20A4) READ (5,3) NMAT,NPTYPE,NVOL,NLAYER,NLINEL,NDAMP,NBT,NBB,ITERRV, ICKFB FORMAT (2014) I F (NLINEL.EQ. 1) LLIN=.TRUE. I F (NPTYPE.GE.3) LPPDIS=.TRUE. I F (NLAYER. EQ. 1) LDOF1=. TRUE. LJRV = . F A L S E . LCKFB = .FALSE. I F (ITERRV. EQ. 0) LJRV=.TRUE, I F (ICKFB. GE.1) LCKFB=. TRUE. I F (LLIN) WRITE (6, 9001) FORMAT ('0** LINEAR MATERIAL IS USED **•) N = NLAYER N2=N*2 N21=N2-1 IF(NLAYER.GE.1) GO TO 1003 WRITE (6,6) FORMAT(1H0,•** ERROR 2, ERROR IN INPUT **•) STOP DO 1004 1=1,N READ(5,7) MATYP (I) ,NSUBD (I) ,H (I) , WIDTH (I) FORMAT(2I4,2F10.4) I F (NMAT.GE.1) GO TO 1001 WRITE (6,4) FORMAT (1 HO,'** ERROR 1,ERROR IN INPUT * * « ) STOP ORDER OF INPUT CONSTANTS ARE: A 1,A2,A3,B1,B2,B3,C1,C2,C3,C4,D1,D2, DEN-H,DEN-V,PERM.,ALFA,BETA WRITE INPUT DATA WRITE (6,8) NMAT,NPTYPE NVOL NLAYER,NLINEL,NDAMP,NBT, NBB, ITERRV, ICKFB FORMAT ('0 NMAT =',13,' ; NPTYPE =',13,' ; NVOL =',13, ' ; NLAYER = « , I 3 , ' ; NLINEL =',I3,' ;'/' NDAMP =',13, « ; NBT =',I3,« ; NBB =',I3, ; ITERRV = ' , 1 3 , ' ; ICKFB =',13) WRITE(6,9) f  r  1  L I S T I N G OF LIQUEFACTION PROGRAM - AUG,  75  LIQPG  4  176 9  FORMAT(1H0,'***********************'/* i  C C  * * * * * * * * * * * * * * * * * * * * * * *  * MATERIAL  PROPERTIES * ' /  1  I F NDAMP=1, CC INPUT AS FORCE/UNIT VELOCITY 2, ALFA=PROP(1, 16) AND BETA=PROP (1,17) DO 1002 I=1,NMAT READ(5,5) (DUMIN (IJ) , I J = 1 , 17) FORMAT (6F12.U) WRITE(6,10) I , (DUMIN (I J) , IJ= 1,17) FORMAT (1X, 'MATERIAL TYPE', 13,' : A " S AND B«'S : ' , 1 P6 E1 4. 3/20 X , • C ^ S AND D " S :',6E14. 3/20 X, • DEN. - H = , 0 P F 5 . 1 , « P C F ; ' , 3 X , •DEN.-V = ' , F 5 . 1 , « P C F ; PERM. =',1PE10.3,' ; ', * ALFA =',E11.3,3X,'BETA = ',E11.3) DO 1012 IJ=1,N IF (MATYP (IJ) . NE. I) GO TO 1013 DO 1014 IK=1,17 PROP ( I J , I K ) =DUMIN (IK) CONTINUE CONTINUE CONTINUE WRITE(6,11) FORMAT(1H0,'*********************/1X,'* LAYER PROPERTIES *'/1X,  5 10  ,  1014 1013 1012 1002 11 m  12 13 C C C C C  1009  1016  1007 C C  t********************i/j  WRITE (6, 12) FORMAT (2X 'LAY.NO. MATYP NSUBD WRITE (6,13) ( I , MATYP (I) , NSUBD (I) ,H(I) FORMAT(1X,I4,3X,2I8,2F10.2) f  THICK , WIDTH (I)  ASSIGN LOWER MOST BOUNDARY AS PROBLEM' BOUNDARY FORM DIAGONAL MASS MATRIX FIND DEPTH AND SIGZ AT CENTRE OF EACH LAYER DO 1007 1=1,N I J = N-I+1 IK = IJ+1 DT2 = 0. 5*H (IJ) UGI = DT2*PROP (I J , 13) CA = DT2*PROP ( I J , 14) IF (I.GT. 1) GO TO 1009 WT(N) = UGI DEPTH (N) = DT2 SIGZ(N) = CA GO TO 1016 WT(IJ) = UGI+UG1 DEPTH(IJ) = DEPTH (IK) +DT1 + DT2 S I G Z ( I J ) = SIGZ (IJ+1) +CA+BK DT1 = DT2 UG1 = UGI BK = CA CONTINUE CHANGE MASS UNIT TO LB-SEC.SQ./FT. UNIT CA = 1./32. 172 DO 1033 1=1,N  WIDTH') ,I=1,N)  L I S T I N G OF LIQUEFACTION  PROGRAM - AUG, 75  LIQPG  177 1033 C C C** C  WT(I) = WT(I) *CA*WTDTH(I) FIND DAMPING VALUE AND STIFFNESS VALUE FOP. EACH LAYER FORM STIFFNESS AND DAMPING MATRIX. NOTE HALF BAND WIDTH = 2  CALL SYSCKM (NCONT) C C 15  16 24 25 26 27  WRITE I N I T I A L PROPERTIES OF THE LUMPED MASS SYSTEM WRITE(6,15) FORMAT(1 H1,• * * I N I T I A L PROPERTIES OF THE LUMPED MASS SYSTEM** /1H0,' NO. ALFA BETA»,10X,»G•,11X,•C»,11X,•M' 7X,'DEPTH',8X,•SIGZ'/) WRITE (6, 16) (I,PROP (1,16) ,PROP (1,17) ,-GG (I) ,CC(I) , WT (I) ,DEPTH (I) S I G Z ( I ) ,1=1 ,N) FORMAT(1X,14,1P7E12.3) WRITE (6,24) FORMAT('0***INITIAL STIFFNESS MATRIX***') J = 2*N-1 WRITE (6, 25) (SK (I) ,1=1, J , 2) FORM AT(' DIAGONAL TERMS :•/(10X,1P10E11.3)) WRITE(6,26) (SK (I) ,1=2, J , 2) FORMAT(' OFF-DIAGONAL TERMS :•/(10X,1P10E11.3) ) WRITE(6,27) FORMAT ('0***INITIAL DAMPING MATRIX***') WRITE(6,25) (CM (I) ,1=1, J , 2) WRITE(6,26) (CM (I) ,1=2, J,2)  C C** . START DYNAMIC ANALYSIS C READ (5,18) NEQ,INTYP,NC,NCPR,NCPRM,NPCON,NCONT,NPLD,DT 18 FORMAT(8I4,3F10.0) WRITE (6,28) NEQ,INTYP,NC,NCPR,NCPRM,NPCON,NCONT,NPLD 28 FORMAT ('6*** CONTROL NUMBERS FOR DYNAMIC RESPONSE ANALYSIS ***• 10X,'NEQ =',I4,» ; INTYP =',I4,« ; NC =',14 ' ; NCPR =',I4,» ;•/10 X,'NCPRM =',14,' ; NPCON = 14,' ; NCONT =',14,' ; NPLD =',I4,' ;•) LPLOT = .FALSE. RDFR = 1. I F (NPLD. EQ. 1) LPLOT=. TRUE. I F (NCPRM.EQ.0) NCPRM=NCPR I F (NEQ. NE. 1) GO TO 1056 READ(5,30) NCARD,NREC,NFTS,RDFR,FMT 30 FORMAT(3I4,F10. 0/20A4) WRITE (6,29) NCARD,NREC,NFTS,RDFR,FMT 29 FORMAT('0*** EARTHQUAKE RECORD USED AS INPUT ***'/10X,'NCARD =• 14,' ; NREC =',I4,» ; NFTS =',I2,' ; MULTIPLYING 'FACTOR =',F5.2/10X,'FORMAT OF INPUT : ',20A4) NCINP = NCARD*NREC GO TO 1057 1056 READ(5,22) NCINP,NFTS,AMP,OMEGA 22 FORMAT (214,2F10.4)  L I S T I N G OF LIQUEFACTION PROGRAM - AUG,  LIQPG  75  1?8  31 , 1057 3333  19 1050 23 20  WRITE(6,31) NCINP,AMP,OMEGA,NFTS FORMAT(* 0*** SINUSOIDAL BASE INPUT I S USED * * * » / 1 0 X , . * NUMBER OF BASE ACC. INPUT =',15,' ; AMPLITUDE = », 1PE11.3,' ; OMEGA =\,E11.3,' ; NFTS = ' , 0 P I 2 , « .») NREC = 1 0 READ (5,3) NS.UBDT WRITE (6,3333) NSUBDT I F (NSUBDT.EQ.0) NSUBDT=1 F O R M A T ( » 0 * * * * * SUB-DIVISION OF TIME INTERVAL FOR », •CONSOLIDATION =',I3,» ****••) ACSC = RDFR I F (NFTS.EQ.0) ACSC=32.17*RDFR I F (DT.LE.0.1) GO TO 1050 WRITE(6,19) DT F O R M A T ( » 0 * * * TIME INTERVAL TOO LARGE, D T = « , F 6 . 3 , » ! » * * * ' ) STOP READ (5,1) T I T L E WRITE(6,23) T I T L E FORMAT(«0*****»,20A4,•*****»/) WRITE (6,20) DT FORMAT (* 0** TIME INTERVAL USED, D T = , F 6 . 3 , » SEC. ***») IF(NCONT.EQ.0) GO TO 1052 CALL TM1VAL(UG1,DT,TNOT) ,  C  C  I N I T I A L CYCLE  C  1052 21  WRITE(6,21) FORMAT(1H1,* OUTPUT VALUES ARE IN THE FOLLOWING ORDER : * * ' / » ** LAY.NO.,ACC., VEL., DISPL., GAMMA., TAU, EVP, PP. . • **«//)  C  C  START R E P E T I T I V E INTEGRATION,  (PP GEN. AND  PP DIS. I F  C  100  1051 1034  UG1 = 0.0 I F (. NOT. LREAD) GO TO 1034 CALL GETACC(AE,FMT,DT,AMP,OMEGA,NEQ,NREC,NCINP,TNOT) DO 1051 1=1,NREC A E ( I ) = ACSC*AE(I) LREAD = .FALSE. NRECC = NRECC+1 NCPRC = NCPRC+1 NPRMC = NPRMC+1 INCR = INCR+1 ITERNO = 0 CA = INCR TIME = CA*DT+TNOT DDT = DT I F (INCR.GT.NC) GO TO 9999 DGI = AE (NRECC) -UG 1 UGII = UGI UG1 = AE (NRECC) I F (NRECC. NE. NREC) GO TO 1023 NRECC = 0  REQUIRED)  LISTING OF LIQUEFACTION PROGRAM -• AUG, 75  LIQPG 179  LREAD = .TRUE. 1023 CONTINUE C C************** 1054 CALL SYSCKM(NCONT) r j * * * * *********** C 1029 DO 1063 1=1,N BAP (I) = BA(I) BBP(I) = BB(I) XD2P(I) = XD2(I) XD1P (I) = XD1 (I) XDOP (I) = XDO (I) GD2P(I) = GD2(I) GD1P (I) = GD1 (I) GAMMAO (I) = GAMMA (I) 1063 TAXYP(I) = TAXY(I) C CALL NUMERICAL INTEGRATION SCHEME C c****************  LIT = .FALSE. 1064 CALL INTSCH(DDT,UG1,UGII) C**************** C I F (LJRV) GO TO 1062 C ( 3 * * * * * * * * * * * * * * * *  1065 CALL CHEKRV (DT,DDT,UG1,UGT,UGII,ITERNO,NRETN) C**************** C GO TO (1062,1054,1064,1065),NRETN WRITE(6,93) NRETN 93 FORMAT (' 0 * * SOMETHING WRONG, NRETN = ,I4,» ***•) STOP C 1062 CONTINUE IF (NPTYPE. LT. 2) GO TO 1035 LSTOP = .FALSE. C ,  c***************  CALL PPGEN(TIME,LCPPG,LSTOP) IF(LSTOP) GOTO 9 9 9 9 C*************** c  1061 C  I F (.NOT.LPPDIS) GO TO 1058  c****************  CALL PPDISN(DT,NSUBDT,LCPPG) C**************** C GO TO 1035 1058 I F (.NOT.LCPPG) GO TO 1035 DO 1059 1=1,N  7  LISTING OF LIQUEFACTION PROGRAM - AUG, 75  LIQPG  8  180  1059 C 1035 1666  C C  IF (.NOT.DODEVP(I)) GO TO 1059 PP (I) = PP (I) +DPP (I) CONTINUE , DO 1666 1=1,N IF (PP(I) .GT,SIGZ(I)) PP (I) =SIGZ (I) CONTINUE  •  .  LPMAX = .FALSE. CALL MAXVAL (TIME,UG1)  IF(LPLOT) WRITE (7) TIME, UG1 , (I, ABXD2 (I) , XD1 (I) , XDO (I) , GAMMA (I) , TAXY (I) , EVP (I) ,PP(I) ,I=1,N) 1055 I F (NCPRC.LT.NCPR) GO TO 1667 NCPRC = 0 C DISPLACEMENT IS REFERRED TO THE TOP OF THE SOIL LAYER WRITE(6,90) TIME,UG1 , GG (1) 90 FORMAT (1H0, •TIME = ,F6.3,' S E C ; ACC. = ,F8.t»,' FPSS', 1P2E13.4) IF (NPTYPE.LE. 1) WRITE (6,91) (I, ABXD2 (I) ,XD1 (I) ,XD0 (I) , GAMMA (I) , TAXY (I) ,1=1, N) 91 FORMAT (1X,IU, 1P5E1 2. 3) IF (NPTYPE. GT. 1) WRITE(6,92) (I,ABXD2(I) ,XD1 (I) ,XD0 (I) ,GAMMA (I) , TAXY (I) ,EVP (I) ,PP (I) ,1=1 ,N) 92 F0RMAT(1X,IU,1P7E12.3) 1667 I F (NPRMC.LT.NCPRM) GO TO 100. NPRMC = 0 LPMAX = .TRUE. CALL MAXVAL (TIME,UG1) . GO TO 100 9999 IF (NPCON.NE.1) GO TO 9998 WRITE(8) TIME,UG 1 , DT DO 1060 1=1,N 1060 WRITE (8) I,XD2 (I) ,XD1 (I) ,XD0 (I) ,TAXY (I) , EVP (I) , PP (I) ,GAMMA (I) , GAMMAO (I) , GAMR (I) ,TAUR (I) , DIRN (I) 9998 I F (LPLOT) END FILE 7 C LPMAX = .TRUE. CALL MAXVAL (TIME,UG1) C WRITE (6,150) INCR 150 FORMAT(* 0******************************************* */ • * NORMAL TERMINATION OF PROGRAM *'/ » * NUMBER OF INCREMENTAL CALCULATION =',14,' *•/ ,  •  ,  * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * i )  STOP C *********************** C * END OF MAIN PROGRAM * C *********************** END BLOCK DATA LOGICAL*1 LVDIR (20) ,LREV (20) ,LRG (20) , DODEVP (20) ,LLIQ (20) , LDOF1,LLIN,LSTOP,LMCLG  L I S T I N G OF LIQUEFACTION PROGRAM - AUG, 75  LIQPG  9  181 COMMON /NCONTL/N,NDAM P,INCR,INTYP,NVOL,NBT,NBB /LCONTL/LDOF1,LLIN,LSTOP,LMCLG /LOGCOM/LVDIR,LREV,LRG DODEVP,LLIQ /PROPTY/TITLE (20) ,PROP(20, 17) ,H (20) ,DEPTH(20) ,WIDTH(20) ,WT (20) , MATYP(20),SIGZ(20) ,CC(20) ,GG (20) ,NSUBD(20) ,GAMR(20) , TAUR (20) ,CM (40) , SK (40) ,GMN (20) ,TMN (20) . /DYNEQT/DYM (40) , PEP (20) ,FOREP (20) . /CALVAL/DIRN(20) ,XD2(20) ,ABXD2(20) ,XD1 (20) ,XD0(20) , B D ( 2 0 ) , BA (20) ,BB(20) ,GD2 (20) ,GD1 (20) ,GAMMA (20) , GAM MAO (20) , . GAMMAX (20) ,TAXY (20) , EVP (20) ,SK1R (20) ,DPP (20) ,PP (20) . /CPVALU/BAP (20) ,BBP(20) ,XD2P (20) ,XD1P(20) ,XD0P(20) ,GD2P(20) , GD1P(20) ,TAXYP(20) . , /REVPAR/NRPT (20) ,GGAMR (20, 8) ,TTAUR (2 0,8) , DIRNR (20,8) DIMENSION BLOC(360) ,BLOC2 (160) EQUIVALENCE (BLOC (1) , DIRN (1) ) , (BLOC2 (1) ,BAP (1) ) DATA INCP/0/, LDOF1/.FALSE./, L L I N / . F A L S E . / , LSTOP/.FALSE./, LMCLG/.FALSE./, L V D I R / 2 0 * . F A L S E . / , LEEV/2 0*.FALSE./, LRG/20*.FALSE./, DODEVP/20*.FALSE./, LLIQ/20*.FALSE./, DEPTH/20*0.0/, WT/20*0.0/, S I G Z / 2 0 * 0 . 0 / , CC/20*0.0/, GG/20*0.0/, BLOC/360*0.0/, BLOC2/16O*0.0/, FORER/20*0.0/, NRPT/20*0/, GGAMR/20*0.0/, TTAUR/20*0.0/ END SUBROUTINE TM1VAL(TNOT,UG1,DT) COMMON " . /NCONTL/N,NDAMP,INCR,INTYP,NVOL,NBT,NBB . /PROPTY/TITLE(20) ,PROP(20,17) ,H(20) ,DEPTH (20) ,WIDTH (20) ,WT(20) , HATYP(20) ,SIGZ (20) ,CC (2 0) ,GG (20) ,NSUBD (20) ,GAMP. (20) , .. TAUR (20) ,CM (40) , SK (40) , G M N (20) ,TMN (20) . /CALVAL/DIRN (20) ,XD2 (20) ,ABXD2 (20) ,XD1 (20) ,XD0 (20) ,BD(20) , BA(20) ,BB(20) ,GD2(20) ,GD1 (20) ,GAMMA(20) ,GAMMAO(20) , GAMMAX (20) , TAXY (20) ,EVP (20) , SK1R (20) , DPP (2 0) ,PP (20) . /REVPAR/NRPT (20) ,GGAMR (20,8) ,TTAUR (20,8) , DIRNR (20,8) WRITE (6,21) FORMAT ('0*** SUBROUTINE TM1VAL I S CALLED ***•/) READ (8) TN0T,UG1,DT1 WRTTE(6,22) TN0T,UG1,DT1 FORMAT(* TIME =',F7.3,' SEC.; A C C = « , F 7 . 3 , ' FPSS; TIME', » INTERVAL = , F 5 . 3 , ' SEC.'/) TEM = ABS (DT-DT1) I F (TEM. GT. 0. 0005) GO TO 888 WRITE (6,23) FORMAT('ONO. ACC',9X,'VEL•,8X,'DISP',7X,'TAUXY•,9X,'EVP•, 10X,*PP',7X,'GAMMA',6X,'GAMMAO',8X,'GAMR »,8X,'TAUR', 7X,'DIRN'/) . . . „ .  21  22  f  ,  23  C DO 101 1=1,N READ (8) a,XD2 (I) , XD1 (I) ,XD0 (I) , TAXY (I) ,EVP(I) ,PP(I) ,GAMMA (I) , GAMMAO (I) , GAMR (I) , TAUR (I) , DIRN (I) I F ( J . NE.I) GO TO 887 BA (I) = XD2 (I) *DT BB(I) = (XD1 ( I ) + 0 . 5*BA (I) ) *DT  LISTING OF LIQUEFACTION PROGRAM - AUG, 75  LIQPG 10  182 101 24 888 25 887 26 999  14 15  11  17 12 18 102  WRITE(6,24) J,XD2(I) ,XD1 (I) ,XD0(I) ,TAXY(I) ,EVP(I) ,PP(I) , GAMMA (I) , GAMMAO(T) , GAMR (I) ,TAUR(I) ,DIRN(I) FORMAT (1X,12,1X,1P10E12.4,0PF6.1) RETURN' WRITE (6, 25) FORMAT (»0**** SOMETHING WRONG, D T " S ARE NOT EQUAL !!! ****«) GO TO 999 WRITE(6,26) FORMAT(* 0**** SOMETHING WRONG, I NOT EQUAL TO J I ! ! ****») STOP END SUBROUTINE GETACC{AE,FMT,DT,AMP,OMEGA,NEQ,NREC,NCINP,TNOT) LOGICAL*1 LDOF1,LLIN,LSTOP,LMCLG,LO,L1,L2,L3 COMMON . /NCONTL/N,NDAMP,INCR,INTYP,NVOL,NBT,NBB . /LCONTL/LDOF1,LLIN,LSTOP,LMCLG DIMENSION AE (10) ,FMT(20) DATA LO/.FALSE./,L1/.FALSE./,L2/.FALSE./,L3/.FALSE./ , IF (L3) GO TO 999 IF (INCR. GT. NCINP) GO TO 99 > I F (LO) GO TO 11 IF (L1) GO TO 12 IF (L2) GO TO 13 IF (NEQ. NE. 0) GO TO 14 LO = .TRUE. GO TO 11 • I F (NEQ.NE.1) GO TO 15 L1 •= .TRUE. GO TO 12 . FNEQ = FLOAT (NEQ) RNPI = FNEQ*3.1415926 RNPI = 1./RNPI L2 = .TRUE. GO TO 18 FI = INCR T = FI*DT DO 17 1=1,10 T = T+DT TEM = T+TNOT TEM = TEM*OKEGA AE(I) = AMP*SIN (TEM) CONTINUE RETURN READ (5, FMT, END=99) (AE (IJ) ,IJ= 1 , NREC) RETURN I F (TNOT.LT.DT) GO TO 13 WRITE(6,102) TNOT FORMAT(•0**** GRADUALLY INCREASING AMPLITUDE REQUIRED, BUT TNOT *, • IS NON-ZERO, TNOT = ,F6.2,« SECS. ****') GO TO 99 F I = INCR T = FI*DT DO 19 1=1,10 ,  13  L I S T I N G OF LIQUEFACTION PROGRAM - AUG,  75  LIQPG  11  183  C  21 20 19 99 22 999  T = T+DT TEM = T*OMEGA ACCT =. AMP*SIN (TEM) TEM = TEM*RNPI TEM IS DIVIDED BY N PI I F (TEM. GE. 1.0) GO TO 21 TEM = TEM*1.570796 TEM = SIN (TEM) ACCT = ACCT*TEM GO TO 20 LO = .TRUE. L2 = .FALSE. A E ( I ) = ACCT CONTINUE RETURN DO 22 1=1,10 A E ( I ) = 0.0 L3 = .TRUE. RETURN END SUBROUTINE SYSCKM(NCONT) LOGICAL*1 LVDIR (20) ,LREV (20) ,LRG (20) ,DODEVP (20) ,LLIQ (20) , LDOF1,LLIN,LSTOP,LMCLG,LJUMP1,LJUMP2,L3,L4 COMMON . /NCONTL/N,NDA MP,INCR,INTYP,NVOL,NBT,NBB . /LCONTL/LDOF1,LLIN,LSTOP,LMCLG . /LOGCOM/LVDIR ,LR EV,LRG,DODEVP, LLTQ. . /PROPTY/TITLE(20) ,PROP(20,17) ,H(20) ,DEPTH(20),WIDTH(20) ,WT(20) , MATYP{20) ,SIGZ(20) ,CC(20) ,GG(20) ,NSUBD(20) ,GAMR(20) , TAUR (20) ,CM (40) ,SK (40) ,GMN (20) ,TMN (20) . /CALVAL/DIRN(20) ,XD2(20) , ABXD2(20) , XD1 (20) ,XDO (20) , BD (20) , BA(20) , BB ( 20) ,GD2(20) ,GD1 (20) ,GAMMA (20) ,GAM M AO (20) , GAMMAX (2 0)",TAXY (20) ,EVP (20) , SK1R (20) , DPP (20) ,PP (20) . /REVPAR/NRPT (20) ,GGAMR (20,8) , TT AUR (20,8) , DIR NR (20,8) DIMENSION FK (20) DATA LJUMP1/.FALSE./,LJUMP2/.FALSE./,KPRT/0/ DATA RLIM1/0.0500/, RLIM2/0.05/  C C  THIS PORTION  FINDS SHEAR  MODULUS FOR  A LAYER  c IF C C 27 15  C C  (LJUMP1) GO TO  14  HERE IS FOR LINEAR CASE OR FIRST INCREMENT DO 15 1=1,N GMN (I) = PROP (1,1) TMN(I) = PROP(I,4) GG (I) = GMN (I) I F (.NOT. LLIN) LJUMP1 = .TRUE. GO TO 50 HERE IS NON-LINEAR CASE, AND  INCR  > 1  LISTING OF LIQUEFACTION PROGRAM - AUG, 75  LIQPG 12  184 14  17  C C 19 2005 29 28  31  30  DO 16 1=1,N I F (.NOT. LVDIR (I) ) GO TO 16 IF (.NOT.DODEVP(I)) GO TO 17 EVPI = EVP (I) CEP 1 = 1. +EVPI/ (PROP (1,2) + PROP (1,3) *EVPI) CEP2 = 1.+EVPI/(PROP (1,5)+PROP (1,6) *EVPI) RATS = 1.-PP (I)/SIGZ (I) SATT = RATS IF (RATT.LT. RLIM1) RATT=RLIM 1 TMN(I) = P R O P (1,4) *CEP2*RATT RATT = SQRT (RATS) IF (RATT. LT.RLIM2) RATT=RLIM2 GMN (I) = PROP (I, 1) *CEP1*RATT CONTINUE J = NRPT(I) GMNT = GMN (I) TMNT = TMN (I) DIRNI = DIRN (I) GAMI = GAMMA (I) TAUI = TAXY (I) GAMAB = ABS (GAMI) TAUAB = ABS (TAUI) IF (. NOT. LRG (I) ) GO TO 20 LRG(I) = -FALSE. HERE REVERSAL OCCURS, INCREASE NRPT BY 1 I F (J.LT.R) GO TO 29 . IF (J.EQ.8) GO TO 29 WRITE (6,2005) INCR,I,J FOR MAT (' **** AT INCR',14," FOR LAY.»,I3,« » ****!)  STOP J = 7 NRPT (I) = J GO TO 40 TEM1 = GAMR(I)*GMNT/TMNT TEM2 = TAUR (I) /TMNT IF (J.EQ.O) GO TO 30 TDIF = TEM2-TTAUR (I, J) TPRD = DIRNR (I, J) *TDIF IF (TPRD.GT.0.0) GO TO 31 J = J+1 NRPT (I) = J DIRNR (I, J) = SIGN (1. ,TDIF) GGAMR (I, J) = TEM 1 TTAUR (I, J) = TEM2 GO TO 40 J = 1 NRPT (I) = 1 DIRNR (1,1) = SIGN (1. ,TEM2) GGAMR (1,1) = TEM 1 TTAUR(I,1) = TEM 2 GO TO 40  J > 8 AND =',19,  LISTING OF LIQUEFACTION PROGRAM - AUG, 75  LIQPG 13 185  C 20 C C C  C 21 C  I F (J.EQ.O) GO TO 18 CHECK FOR REDUCTION IN NRPT, I E . CHECK FOR ENVELOP CURVES :NOTE THAT THE LAST BUT ONE REVDRSAL POINT IS ALSO THE LIMIT FOR SWITCHING TO ENVELOP CURVE I F (J.GT. 1) GO TO 21 TEMP = ABS (TTAUR (I, 1) ) *TMNT IF (TAUAB.LT.TEMP) GO TO 40 J = 0 NRPT (I) = 0 GO TO 18  TEMP = TTAUR (I,J-1) *TMNT CHECK FOR CONSISTENCY BETWEEN DIRN AND TTAUR'S TDIF = TEMP-TTAUR (I, J) *TMNT TPRD = DIRNI*TDIF IF (TPRD) 24,25,25 24 KPRT = KPRT+1 J = J-1 NRPT (I) = J IF (KPRT. GT. 50) GO TO 20 IF (J.NE.O) GO TO 26 WRITE(6,2004) INCR,I,J GO TO 20 26 J1 = J+1 WRITE (6,2004) TNCR,I,J, (TTAUR(I,K) ,K=1,J1) ,TAUI 2004 FORMAT (1X,314,1P9E13.3) GO TO 20 25 I F (TDIF) 22,23, 23 23 . IF (TAUI.LT.TEMP) GO TO 40 J = J-2 NRPT (I) = J GO TO 20 22 IF (TAUI.GT.TEMP) GO TO 40 J = J-2 NRPT (I) = J GO TO 20 40 TEMP = GGAMR(I,J)*TMNT/GMNT GGT = 1. 0/(1.0 + 0.005*GMNT*ABS (GAMI-TEMP)/TMNT) GO TO 116 18 GGT = 1.0/(1.0+0.01*GMNT*GAMAB/TMNT) 116 GGT = 0.98*GGT*GGT*GMNT RLG = 1./PROP (1,1) IF (GGT*RLG.LT.0.005) GGT=0.005/RLG GG(I) = GGT C IF (TAUAB. LT. TMNT) GO TO 16 KPRT = KPRT+1 IF (KPRT.GT.50) GO TO 16 IF (KPRT. EQ. 50) GO TO 45 IF (TAUI) 42,43,43 42 I F (DIRNI.GT.0.0) GO TO 16 WRITE(6,2001) INCR,I  LISTING OF LIQUEFACTION PROGRAM - AUG, 75  LIQPG  186 2001 43 2002 45 2003 C 16 C  FORMAT(' ** AT I N C R ' , 1 4 , ' LAY.• , 1 3 , • TAXY < -TMN *) GO TO 16 I F (DIRNI.LT.0.0) GO TO 16 WRITE(6,2002) INCR,I . , FORMAT (' ** AT I N C R ' , 1 4 , ' L A Y . ' , 1 3 , • TAXY > TMN') GO TO 16 WRITE (6,2003) FORMAT(* *** NO MORE WRITING FOR |TAXY| > |TMN| ***') CONTINUE  50  CONTINUE  C  THIS PORTION FORMS SYSTEM DAMPING  C C** C C  CC AND G VECTORS CONTAIN LAYER DAMPING AND STIFFNESS VALUES C AND K ARE BANDED DAMPING AND STIFFNESS MATRICES  C  57 56  59 C 58 52 C  AND STIFFNESS MATRICES  IF (LJUMP2) GO TO 58 LJUMP2 = .TRUE. AA1 = PROP(1,16) BB1 = PROP (1,17) DO 56 1=1,N IF (NDAMP.EQ.1) GO TO 57 CC(I) = BB1*GG(I) GO TO 56 CC(I) = PROP (I, 16) CONTINUE DO 59 1=1,N 12 = 2*1 I2M1 = 12-1 CCI = CC (I) *WIDTH (I)/H (I) CCI1 = CC (1 + 1) *WIDTH (1+1)/H (1+1) CM(I2M1) = CCI IF (NDAMP.EQ. 2) CM (I2M1) =CM (I2M1) +AA1*WT (I) IF (I.GE.N) GO TO 59 CM(I2M1) = CM (I2M1)+CCI1 CM (12) = -ecu CONTINUE DO 52 1=1,N FK(I) = GG (I) *WIDTH (I)/H (I) DO 51 1=1,N 12 = 2*1 I2M1 = 12-1 FKI = FK(I) FKI1 = FK (1+1) SK (I2H1) = FKI IF (I.GE.N) GO TO 51 SK(I2M1) = SK (12M1) + FK11  .  LISTING OF LIQUEFACTION PROGRAM - AUG, 75  LIQPG 15  187 SK(I2) = -FKI1 CONTINUE IF (INCR. LT. 113) GO TO 9999 IF (INCR.GT. 120) GO TO 9999 ; WRITE(6,2999) INCR, (GAMMA (K) , K= 1 , N) WRITE (6,2999) INCR, (TAXY(K) ,K=1,N) FORMAT (1X,I4,1P10E12.3/(5X,10E12. 3) ) RETURN END SUBROUTINE SOLB2 COMMON /NCONTL/N,NDAMP,INCR,INTYP,N70L,NBT,NBB . /DYNEQT/DYM (40) ,PEP(20) ,FORER(20) DIMENSION T(20) C DYM BECOMES THE LOWER TRIANGULAR MATRIX C T IS THE UPPER TRIANGULAR MATRIX C NOTE L33 IS DYM (5) , ETC, T12 IS T ( 1 ) , T23 IS T (2) ETC. IF (ABS(DYM(1)).GT.1.0E-75) GO TO 103 WRITE (6, 10) DYM (1) 10 FORMAT («0*** ZERO ENCOUNTERED IN ROW 1 !!! *** DYM (1)=•,E12.4) STOP 103 CONTINUE PEP(1) = PEP (1)/DYM (1) IF(N.EQ.1) RETURN DO 101 1=2,N IM1 = 1-1 MN = 2*IM1 NN = MN+1 MM = MN-1 AMN = DYM (MN) • T (IM1) = AMN/DYM (MM) DYM (NN) = DYM (NN)-T (IM1) *AMN IF (ABS (DYM (NN) ) .GT. 1. OE-75) GO TO 104 WRITE(6,11) I , NN , DYM (NN) 11 FORMAT (* 0*** ZERO ENCOUNTERED IN ROW',13,' !!! *** DYM(',I2, •)=',E12.4) STOP 104 CONTINUE PEP(I) = (PEP(I) -AMN*PEP(IM1) )/DYM(NN) 101 CONTINUE C BACK SUBSTITUTION DO 102 1=2,N MN = N-I + 1 NN = MN+1 102 PEP,(MN) = PEP(MN) -T(MN) *PEP(NN) RETURN END SUBROUTINE CHEKSV(DT,DDT,UG1,UGI,UGII,ITER NO,NRETN) LOGICAL*1 LDOF1,LLIN,LSTOP,LMCLG,LJUMP,LITE,LCPPG,LPRT, LVDIR (20) , LREV (20) , LRG (20) ,DODSVP (20) , LLIQ (20) , LRGL(20),LJRV,LCKF3,LIT,LPMAX COMMON /NCONTL/N,NDAMP,INCR,INTYP,NVOL,NBT,NBB /LCONTL/LDOF1,LLIN,LSTOP,LMCLG,LJRV,LCKFB,LIT,LPMAX 51 C C C C C2999 9999  LISTING OF LIQUEFACTION PROGRAM - AUG, 75  LIQPG 16  188  C C  102 101  103  104 100  111 9 124 C  /LOGCOM/LVDIR,LREV,LRG,DODEVP,LLIQ /PROPTT/TITLS (20) , PROP (20, 17) , H (20) , DEPTH (20) , WIDTH (20) ,WT (20) MATYP (2 0) ,SIGZ (20) ,CC (20) ,GG (20) , NSUBD (20) , G AMR (20) , TAUR (20) ,CM (40) ,SK{40) ,GMN (20) ,TMN (20) /DYNEQT/DYM (40) ,PEP (20) ,FORER (20) /CALVAL/DIRN (20) ,XD2(20) ,ABXD2(20) ,XD1 (20) , XDO (20) ,BD(20) , BA (20) , BB (20) ,GD2 (20) ,GD1 (20) ,GAMMA (20) ,GAMMAO (20) , GAMWAX(20) ,TAXY(20) ,EVP(20) ,SK1R(20) ,DPP(20) ,PP(?0) . /CPVALU/BAP (20) ,BBP (20) ,XD2P (20) ,XD1P (20) ,XD0P (20) ,GD2P (20) , GD1P (20) ,TAXYP (20) DIMENSION NREV (20) ,GAMRO (20) ,TAURO (20) DATA JCALL/0/,LJUMP/.FALSE./,LPRT/.TRUE./,JPRINT/0/, LRGL/20*. FALSE. / IF (INCR. EQ. 50) LPRT=. FALSE. I F (LJUMP) GO TO 100 IF (JCALL.GE.1) GO TO 101 DO 102 1=1,N LVDIR (I) = . FALSE. DO 103 1=1,N IF (GAMMA (I) .EQ.0.0) GO TO 103 IF (LVDIR (I) ) GO TO 103 DIRN(I) = SIGN (1. , GAMMA (I) ) LVDIR (I) = .TRUE. CONTINUE JCALL = JCALL+1 I F (JCALL. EQ. 1) GO TO 131 DO 104 1=1,N • I F (. NOT. LVDIR (I) ) GO TO 100 CONTINUE LJUMP = .TRUE. I F (ITERNO.NE.O) GO TO 111 NDT = 0 NDDT = • 10 NTDT = 0 ITERNO = ITERNO+1 IF (ITERNO.LT.11) GO TO 124 WRITE (6,9) ITERNO FORMAT(' ***ITERNO =',14,' ***») STOP CONTINUE CHECK GAMMA FOR REVERSAL LITE = .FALSE. DO 105 1=1,N LREV (I) = .FALSE. IF (LRGL (I) ) GO TO 105 IF (.NOT.LVDIR(I) ) GO TO 105 TEMY = GAMMA (I)-GAMMAO (I) IF (TEMY. EQ. 0.0) GO TO 105 TEMY = DIRN (I) *TEMY IF (TEMY. GT. 0.0) GO TO 105 LREV (I) = .TRUE.  LISTING OF LIQUEFACTION PROGRAM - AUG, 75  LIQPG 17  189 LRG (I) = .TRUE. GAMRO (I) = GAMR (I) TAURO (I) = TAUR (I) GAMR (I) = GAMMAO (I) TAUR (I) = TAXYP(I) GAMMAX (I) = GAMMAO (I) DIRN (I) = -DIRN (I) IF (ABS(GAMMAO(I)) .LT.1.E-8) GO TO 105 LITE = .TRUE. 105 CONTINUE IF ( (.NOT.LITE) .OR.NTDT. EQ. 10) GO TO 131 C SINCE RE-CALCULATION IS REQUIRED, CHANGE GAMR'S AND DIRN'S BACK C TO ITS VALUES BEFORE CHECKING FOR REVERSAL. DO 123 1=1,N IF (. NOT. LREV (I) ) GO TO 123 GAMR (I) = GAMRO (I) TAUR (I) = TAURO (I) DIRN(I) =-DIRN(I) 123 CONTINUE C NREV IS A NO. USED FOR INDICATING NEAREST TENTH OF (DT) FOR C (GAMMA-DOT)=ZERO. DO 106 1=1,N I F (. NOT. LREV (I) ) GO TO 106 IF (ABS (GAMMAO (I) ) .GT. 1. E-8) GO TO 120 NREVI = 1 0 GO TO 121 120 ' TEMY = GD1P (I) TEMP = TEMY-GD1 (I) I F (ABS (TEMP) .GT. 1. E-20) GO TO 133 IF (LPRT) WRITE (6,3) I NCR , ITERNO, I 3 FORMAT(12X,'INCR =»,I4,», ITERNO =',13,', LAYER',13, ', (GD1-GD1P)=0«) NREVI = 1 0 GO TO 121 133 TEMY = 10.0*TEMY/TEMP NREVI = IFIX (TEMY) IF (NREVI.LE.O) NREVI=0 IF (NREVI.GE.10) NREVI=10 121 CONTINUE NREV (I) = NREVI 106 CONTINUE C DO 107 KNT=1,11 K = KNT-1 DO 108 1=1,N IF (.NOT.LREV(I)) GO TO 103 IF (NREV (I) .EQ*. K) GO TO 109 108 CONTINUE 107 CONTINUE 109 NDT = K IF (NDT.GT.NDDT) NDT=NDDT C IF (NDT. NE. 0) GO TO 112  LISTING OF LIQUEFACTION PROGRAM - AUG, 75  LIQPG 18 190  4  113  112  C 134  117  5  116  DO 113 1=1,N LRG (I) = - FALSE. GAMMA (I) = GAMMAO(I) IF (LRGL (I) ) GO TO 113 IF (. NOT. LREV (I) ) GO TO 11 3 I F (NREV (I) . NE. 0) GO TO 113 LRG (I) = .TRUE. LRGL (I) = .TRUE. DIRN(I) = -DIRN(I) IF (LPRT) WRITE(6,4) INCR,ITERNO,I,NDT FORMAT(12X,«INCR =«,I4,', ITERNO =',13,»: ,» CHANGES MOD. , NDT =«,I2) GAMR (I) = GAMMAO (I) TAUR (I) = TAXYP (I) CONTINUE CALL SYSCKM(O) NRETN = 3 LIT = .TRUE. RETURN IF(ITERNO.EQ.1.AND.NDT.EQ.10) GO TO 134 FN = FLOAT (NDT) UG1 = UG1-UGII DDT = 0.1*FN*DT UGII = 0.1*FN*UGI UG1 = UG1+UGII  ',13,  LIT = .TRUE. CALL INTSCH (DDT,UG1,UGII) DO 116 1=1,N LRG (I) = . FALSE. I F (LRGL (I) ) GO TO 1 16 IF (. NOT. LREV (I) ) GO TO 116 I F (NPEV(I) .EQ.NDT) GO TO 117 TEMY = GAMMA (T)-GAMMAO (I) IF (TEMY. EQ. 0.0) GO TO 116 TEMY = DIRN (I) *TEMY IF (TEMY. GT. 0.0) GO TO 116 LRG (I) = .TRUE. LRGL (I) = . TRUE. I F (LPRT) WRITE (6,4) INCR,ITERNO,I,NDT IF (LPRT) WRITE (6,5) I, ABXD2 (I) , XD1 (I) ,XDO (I) , GAMMA (I) , TAXY (I) ,EVP (I) ,PP (I) FORMAT(20X,I4,1P7E12.3) DIRN (I) = -DIRN (I) GAMR(I) = GAMMA (I) TAUR (I) = TAXY (I) CONTINUE NTDT = NTDT+NDT NDDT = 10-NTDT FN = FLOAT(NDDT) DDT = 0.1*FN*DT UGII = 0.1*FN*UGI  LISTING OF LIQUEFACTION PROGRAM - AUG, 75  LIQPG 19 191  131  DG1 = UG1+UGII NRETN = 2 I F (NTDT.GE.10) NRETN=4 RETURN NRETN = 1 DO  132  100 C  IK  103 102  121  132  1 = 1 , N .  LRGL (I) = .FALSE. RETURN END SUBROUTINE INTSCH(DT,UG1,UGI) LOGICAL*1 LJUMP,L2,L3,L4,LDOF1,LLIN,LSTOP,LMCLG, LJRV,LCKFB,LIT,LPMAX COMMON . /NCONTL/N,NDAMP,INCR,INTYP,NVOL,NBT,NBB . /LCONTL/LDOF1,LLIN,LSTOP,LMCLG,LJRV,LCKFB,LIT,LPMAX . /PROPTY/TITLE(20) , P R O P ( 2 0 , 1 7 ) ,H(20) ,DEPTH(20) ,WIDTH(20) , W T ( 2 0 ) , MATYP (20) ,SIGZ (20) ,CC (20) ,GG (20) ,NSUBD (20) ,GAMR ( 2 0 ) , TAUR(2 0) ,CM (40) , SK (40) , GMN (2 0) , TMN (20) . /DYNEQT/DYM (40) ,PEP (20) ,FORER(20) . /CALVAL/DIRN(20) , X D 2 ( 2 0 ) , A B X D 2 ( 2 0 ) , X D 1 (20) , X D 0 (20) , B D ( 2 0 ) , BA (20) , B B (20) , G D 2 (20) , G D 1 (20) , G AMMA (20) , GAMMAO ( 2 0) , GAMMAX(20),TAXY(20) , EVP ( 2 0 ) , S K 1 R ( 2 0 ) , D P P (20) , P P ( 2 0) . /CPVALU/BAP (20) ,BBP(20) , X D 2 P ( 2 0 ) , X D 1 P ( 2 0 ) , X D 0 P ( 2 0 ) , G D 2 P ( 2 0 ) , G D 1 P (20) ,TAXYP (20) . /REVPAR/NRPT(2 0 ) ,GGAMR (20,3) , T T A U R ( 2 0,8) , DIRNR (20,8) DIMENSION FCC E M (20) , FOCEC (20) , FOCEK (20) ,FORERO(20) DATA LJUMP/.FALSE./ IF (LJUMP) GO TO 100 LJUMP = .TRUE. DFINT = 0 . 5 AFINT = 0.16666667 IF (INTYP. EQ. 1) AFINT=0.25 CONST 1 = DFINT*DT C O N S T 2 = AFINT*DT*DT FORM RIGHT HAND SIDE OF DYNAMICAL EQUATION. DO 1 0 1 1 = 1 , N TEM1 = 0.0 TEM2 = 0.0 DO 102 IJ=1,3 I I = I+IJ-2 =  2* (.1-1)  +IJ-1  IF (IK.LE.O) GO TO 103 IF (IK.GE.2*N) GO TO 103 TEM 1 = T E M 1 + C M (IK) *BAP (II) T E M 2 = T E M 2 + SK(IK) *BBP(II) CONTINUE CONTINUE PEP(I) = -WT (I) * U G I - T E M 1 - T E M 2 IF (.NOT. LCKFB) GO TO 1 0 1 IF (.NOT.LIT) GO TO 1 2 1 PEP (I) = PEP (I)-FORERO (I) GO TO 1 0 1 PEP (I) = PEP(I)-FORER (I)  LISTING OF LIQUEFACTION PROGRAM - AUG, 75  LIQPG 20 192  101 C  105 C C C C  C  107  106  115 116 114 C  CONTINUE FORM LEFT HANDSIDE OF DYNAMICAL EQUATION DO 10S 1=1,N I J = 2*1-1 DYH (IJ) = WT (I) +CONST1 *CM (IJ) +CONST2*SK (IJ) I J = IJ+1 DYM(IJ) = CONST1*CM(IJ)+CONST2*SK (IJ) CONTINUE CALL SOLB2 CALCULATE ACC., VEL.,DISPL., GAMMA AND TAXY ETC. DO 106 1=1,N TEH 1 = PEP (I) XD2 (I) = XD2P (I) + TEM1 ABXD2 (I) = XD2 (I) +UG1 XD1 (I) = XD1P (I) + BAP (I) +TEM1*C0NST1 XDO (I) = XDOP (I) + BBP (I) +TEM1 *CONST2 BA (I) = XD2 (I) *DT BB(I) = XD1 (I) *DT+BA (I) *CONST1 IF (I.NE. 1) GO TO 107 TEM 1 = 1./H(1) GD2 (1) = XD2 (1) *TEM 1 GD1 (1) = XD1 (1) *TEM 1 GAMMA (1) = 100. 0*XD0 (1) *TEM1 GO TO 106 TEM 1 = 1./H(I) GD2(I) = (XD2 (I) -XD2 (1-1) ) *TEM1 GD1 (I) = (XD1 (I)-KD1 (1-1) ) *TEM1 GAMMA (I) = 100. 0* (XDO (I)-XDO (1-1) ) *TEM1 CONTINUE DO 114 1=1,N GMNT = GMN (I) TMNT = TMN (I) IF (LLIN) GO TO 116 J = NRPT (I) I F (J.EQ.O) GO TO 115 TEM 1 = 0. 01* (GAMMA (I) -GGAMR (I, J) *TMNT/GMNT) TAXY (I) = TTAUR (I, J) *TMNT +GMNT*TEM1/(1.0+0.5*GMNT*ABS(TEM1)/TMNT) GO TO 114 TEM 1 = 0.01*GAMMA(I) TAXY (I) = GMNT*TEM1/(1.0 + GMNT*ABS (TEM 1) /TMNT) GO TO 114 TAXY (I) = 0.01 *GMNT*GAMMA (I) CONTINUE DO 111 1=1,N FOCEM (I) = WT (I) *ABXD2 (I) FOCEC(I) = CC (I) *WIDTH (I) *GD1 (I) FOCEK(I) = TAXY (I) *WIDTH (I)  LISTING OF LIQUEFACTION PROGRAM - AUG, 75  LIQPG 21  193 IF (I.EQ.N) GO TO 113 11 = 1+1 FOCEC(I) = FOCEC (I) -CC (11) *WIDTH (11) *GD1 (11) FOCEK (I) = FOCEK (I) -TAXY (11) *WIDTH (11) 113 CONTINUE TEM 1 = FOCEH(I) TEM2 = FOCEC(I) TEM3 = FOCEK (I) TEM4 = TEM1+TEM2+TEM3 FORERO(T) = FORER (I) FORER (I) = TEM4 TEM4 = ABS(TEM4) TEM5 = ABS (TEM1) +ABS (TEM2) +ABS (TEM3) IF (TEM4.LE.TEM5) GO TO 111 LSTOP = .TRUE. WRITE (6,2001) INCR,.I,GDI (I) ,GD2 (I) ,TEM 1, TEM2, TEM 3 , FORER (I) 2001 F0RMAT(1X,2I4,2X,1P6E14.4) 111 CONTINUE IF (LSTOP) STOP RETURN END SUBROUTINE PPGEN(TIME,LCPPG,LSTOP) LOGICAL* 1 LVDIR (20) ,LREV (20) , LRG (2 0) , DODEVP (20) , LLIQ (20) , LLPR(2 0) , LUPMP (20) ,LUPMC (20) ,LNUPK (20) ,LXZY (20) , LPRINT,LCALL,LCPPG,LSTOP COMMON . /NCONTL/N,NDAMP,INCR,INTYP,NVOL,NBT,NBB /LOGCOM/LVOIR,LREV,LRG,DODEVP,LLIQ . /PROPTY/TITLE(20) ,PROP (20,17),H (20),DEPTH (20),WIDTH (20) ,WT (20) , MATYP (20) ,SIGZ (20) ,CC (20) ,GG (20) ,NSUBD (20) ,GAMR (20) , TAUR (20) ,CM (40) , SK (40) ,GMN (20) ,TMN (20) . /CALVAL/DIRN (20) ,XD2 (20) ,ABXD2 (20) ,XD1 (20) ,XD0 (20) ,BD(20) , BA (20) ,BB (20) ,GD2 (20) ,GD1 (20) ,GAMMA (20) ,GAMMAO (20) , GAMMAX (20) ,TAXY (20) ,EVP(20) ,SK1R(20) ,DPP (20) ,PP (20) . /CPVALU/BAP (20) ,BBP(20) ,XD2P(20) ,XD1P(20) ,XD0P(2O) ,GD2P(20) , GD1P (20) ,TAXYP (20) DIMENSION YMXP (20) ,DEVP (20) ,DEVPP (20) ,DPPP (20) ,TFYMP (20) , DEVPT(20) ,YMXT(20) ,VCC1 (20) ,RDENK (20) DATA LCALL/.FALSE./ C ,  c  **********  C C C  INITIALIZATION FOR FIRST CALL NOTE PP(I) AND EVP (I) HAS TO BE INITIALIZED IN THE MAIN PROGRAM ********** IF (LCALL) GO TO 105 LCALL = .TRUE. NPRINT = 0 LPRINT = .TRUE. NWLIQ = 0 VCC2 = 0.790 VCC3 = 0.5625 VCC4 = 0.730 POWTM = 0.43  LISTING OF LIQUEFACTION PROGRAM - AUG, 75  LIQPG 22  194  106 C  Q  C C C C C C 105  C C 101  131 C 102 C C 103 C  POWTN = 0.620 SK2T = 0.0025 POWTMC = 1.-POWTM DO 106 1=1,N YMXT(I) = GAMMA (I) BD(I) = GAMMA (I) YMXP (I) = O.'O TFYMP(I) = 0 . 0 DPPP(I) = 0 . 0 LUPMP(I) = .TRUE. VCC1 (I) = 0. 5*PROP (1,7) TEMP = POWTN-POWTM TEMP = POWTM*SK2T*SIGZ (I)**TEMP RDENK(I) = 100.0/TEMP SK1R(I) = RDENK (I) *SIGZ (I) **POWTMC LLPR(I) = • FALSE. RETURN **********  CHECK GAMMA FOR CROSSING ZERO VALUE AND DETERMINE GAMMAX LXZY = .TRUE. FOR GAMMA CROSSING ZERO STRAIN AXIS LUPMP = .TRUE. FOR GAMMA INCREASING IN MAG. DURING LAST INTERVAL LUPMC = .TRUE. FOR GAMMA INCREASING IN MAG. DURING CUR. INTERVAL LNUPK = .TRUE. FOR NEW PEAK OCCURING DURING CURRENT INTERVAL ********** CONTINUE LCPPG = .FALSE. DO 122 1=1,N DPP (I) = 0 . 0 LUPMC (I) = .FALSE. DODEVP (I) = .FALSE. LXZY (I) = .FALSE. TEMP = BD (I) TEMC = GAMMA (I) TPRD = TEMP*TEMC IF (TPRD) 102,103, 101 GAMMAO AND GAMMA ARE ON THE SIDE, CHECK FOR GAMMAX TEMPI = ABS(YMXT(I)) IF (ABS (TEMC) .LE.TEMPI) GO TO 131 YMXT(I) = TEMC IF (ABS (GAMMAX (I)) . LT. ABS (YMXT (I) ) ) GAMMAX (I) =YMXT (I) I F (ABS (TEMC) .GT. ABS (TEMP) ) LUPMC (I) =. TRUE. GO TO 104 CONTINUE GAMMA CHANGES SIGN HERE - I.E. CROSSES ZERO, AXIS, STOP PPGEN LXZY (I) = .TRUE. GO TO 104 CONTINUE EITHER GAMMA OR GAMMAO OR BOTH EQUAL 0 IF (TEMP. NE. 0.0) GO TO 102  LISTING OF LIQUEFACTION PROGRAM - AUG, 75  LIQPG 23  195  104  121 122 C  c C C  LOPHC(I) = .TRUE. GO TO 121 LNUPK(I) = .FALSE. I F (LUPMP (I) . AND. . NOT. LUPMC (I) ) LNUPK (I) = . TRUE. IF (LUPMC (I)) GO TO 121 LCPPG = .TRUE. DODEVP (I) = '.TRUE. LUPMP(I) = LUPMC (I) CONTINUE  **********  CHECK FOR PRINTING GAMMAX ********** IF {.NOT. LCPPG) GO TO 888 DO 107 1=1,N I F (LLPR (I) . AND.LXZY (I) ) GO TO 109 107 CONTINUE GO TO 108 109 I F (. NOT- LPRINT) GO TO 108 NPRINT = NPRINT*1 IF (NPRINT.GT.20) LPRINT=.FALSE. WRITE (6,5) INCR 5 FORMAT('0 INCR =',I4,» GAM-MAX VOL CHANG " TIME •) DO 110 1=1, N IF (. NOT. LLPR (I) ) GO TO 110 WRITE (6,6) I,YHXP(I) ,DEVPP(I) ,DPPP(I) ,TFYMP(I) 6 FORMAT (1 OX, 14 ,1P3E 12. 3, 0PF8. 2) LLPR (I) = -FALSE. 110 . CONTINUE 108 CONTINUE C  c  **********  C  CALCULATE VOLUME CHANGE OR PORE PRESSURE INCREMENT  Q  114  115 116  PRES CHAN',  **********  DO 130 1=1,N DEVP(I) = 0 . 0 DPP (I) = 0 . 0 IF (.NOT.DODEVP(I)) GO TO 130 IF (. NOT. LNUPK (I) ) GO TO 116 TEMC = EVP (I) TEMP = ABS (YMXT (I) ) IF (TEMP.GT. 1.E-10) GO TO 115 DODEVP(I) = .FALSE. YMXT (I) =0.0 DEVPT(I) = 0 . 0 GO TO 130 TPRD = VCC1 (I)*(TEMP-0.79*TEMC+0.5625*TEMC*TEMC/(TEMP+0.73*TEMC) ) IF (TPRD.LT.0.0) GO TO 114 DEVPT (I) = TPRD I F (LUPMC(I)) GO TO 130 TEMP = ABS (YMXT (I) ) IF (TEMP.LT. 1.E-10) GO TO 130  LISTING OF LIQUEFACTION PROGRAM - AUG, 75  LIQPG 196  C  C 117  TEMC = ABS (BD (I) -GAMMA (I) ) IF (LXZY (I)) TEMC=ABS (BD (I) ) DEVP(I) = TEMC*DEVPT (I)/TEMP EVP (I) = EVP (I) +DEVP (I) IF (NVOL.EQ.O) GO TO 117 DPP (I) = DEVP(I) PP (I) = EVP (I) GO TO 118  '. •  ,  CALCULATE PORE PRESSURE INCREMENT SGVP = SIGZ (I)-PP (I) I F (SGVP.GT.0.05*SIGZ (I) ) GO TO 112 SGVP = 0.05*SIGZ (I) I F (LLIQ (I) ) GO TO 112 LLIQ(I) = -TRUE. WRITE ( 6 , 3 ) INCR, I 3 FORMAT(' **** AT INCR =',I4,» ; LAYER',13,' LIQUEFIED ! ****') NWLIQ = NWLIQ+1 N02 = N/2 IF (NWLIQ.LE.N02) GO TO 112 WRITE(6,200U) 2004 FORMAT('0** HALF OF DEPOSIT HAS LIQUEFIED ! **') LSTOP = .TRUE. GO TO 888 112 DEVPS = 0.01*DEVP(I) TEMC = RDENK(I)*SGVP**POWTHC DPP(I) = TEMC*DEVPS SK1R(I) = TEMC 118 CONTINUE IF (. NOT. LXZY (I) ) GO TO 130 LLPR(I) = .TRUE. YMXP(I) = GAMMAX (I) GAMMAX (I) = 0.0 YMXT(I) = GAMMA(I) DEVPP (I) = DEVPT (I) DPPP(I) = DPP (I) 130 CONTINUE 888 DO 889 1=1,N 889 BD(I) = GAMMA (I) RETURN END SUBROUTINE PPDISN(DT,NSUBDT,LCPPG) C SOLUTION OF ONE-DIMENSIONAL HEAT OR CONSOLIDATION PROBLEMS. C MAXIMUM 20 LAYERS C FINITE DIFFERENCE SOLUTION C PROP (I,J) = PROPERTIES OF SOIL LAYER FROM MAIN PROGRAM C SK1R(I) = ONE-D REBOUND MODULUS FOR THE LAYER C H(I) = THICKNESS OF SOIL LAYER C NSUBD (I) = NUMBER OF SUBDIVISIONS IN THE LAYER C PPO (I) = PORE PRESSURES AT THE START OF TIME INTERVAL C DPP (I) = INCREMENTS OF PORE PRESSURE FOR THE TIME INTERVAL C REFERRED TO MID-LAYER C DT = TIME INTERVAL FOR ITERATION  LISTING OF LIQUEFACTION PROGRAM - AUG, 75  LIQPG  197 C C C C C C C  C 101  130 131 132 133 C  15 102 103 49 C C  N NTB,NBB DZ(I) CV(I) AM (I)  NUMBER OF SOIL LAYERS 0 FOR W=0 AT ALL TIMES AT TOP OR BOTTOM BOUNDARY, = 1 FOR (DW/DN) =0 AT ALL TIMES. = DEPTH INTERVAL = COEFFICIENTS OF CONSOLIDATION C.V*DT/(DZ (I) **2) —  LOGICAL*1 LVDIR (20) ,LREV(20) ,LRG (20) ,DODEVP(20) ,LLIQ(20) , LDOF1,LLIN,LSTOP,LMCLG,LCPPG,LJUMP,L3,L4 COMMON /LCONTL/LDOF1,LLIN,LSTOP,LMCLG /NCONTL/N,NDAMP,INCR,INTYP,NVOL,NBT,NBB /LOGCOM/LVDIR,LREV,LRG,DODEVP,LLIQ /PROPTY/TITLE(20) ,PROP (20,17) ,H (20) ,DEPTH (20) ,WIDTH (20) ,WT (20) HATYP (20) ,SIGZ (20) ,CC (20) ,GG (20) ,NSUBD (20) ,GAMR (20) , TAUR (2 0) ,CM (40) ,SK (40) ,GMN (2 0) ,TMN (20) /CALVAL/DIRN (20) ,XD2(20) , ABXD2(2 0) ,XD1 (20) ,XD0 (20) ,BD(20) , BA(20) ,BB(20) ,GD2 (20) ,GD1 (20) , GAMMA (20) , GAMMAO (20) , GAMMAX (2 0) ,TAXY (20) ,EVP (20) , SK1R (20) , DPP (20) ,PP (20) DIMENSION W (100) ,TW (100) ,CV (20) , AM (20) ,DZ (20) DATA LJUMP/.FALSE./, W/100*0.0/, TW/100*0.0/ IF (LJUMP) GO TO 25 LJUMP = .TRUE. WRITE(6,101) FORMAT (•0*** DESCRIPTION OF CONSOLIDATION PROBLEM ***«/) I F (NBT. EQ. 0) WRITE(6,130) IF (NBT. EQ. 1) WRITE(6,131) IF (NBB. EQ. 0) WRITE(6,132) I F (NBB. EQ. 1) WRITE(6,133) FORMAT (' AT THE TOP BOUNDARY THE CONDITION IS ( W=0 )') FORMAT (' AT THE TOP BOUNDARY THE CONDITION IS ( DW/DN = 0 )') FORMAT (' AT THE BOTTOM BOUNDARY THE CONDITION IS ( W=0 )•) FORMAT (' AT THE BOTTOM 30UNDARY THE CONDITION IS ( DW/DN = 0 ) ' NTOT IS THE TOTAL NUMBER OF SOIL POINTS FOR PP CALCULATION NTOT = 0 DO 15 1=1,N IF (NSUBD (I) . EQ. 0) NSUBD(I)=2 FN = FLOAT(NSUBD(I)) DZ (I) = H (I) /FN NTOT = NTOT+NSUBD(I) CONTINUE NTOT = NTOT+1 WRITE (6, 102) NTOT FORMAT (» TOTAL NUMBER OF SOIL POINTS FOR PP CALCULATION =',I4) IF (NTOT.LE.100) GO TO 49 WRITE(6,103) FORMAT ('0*** NTOT IS GREATER THAN 100, PROCESSING STOPS.***') STOP CONTINUE INITIALIZING THE VARIABLES DO 16 IJ=1,NTOT  LISTING OF LIQUEFACTION PROGRAM - AUG, 75  LIQPG 26  198  16 105 C C  135 201 50  109 14 107 108 C C 25  17  C  19  W (IJ) = 0. 0 TW(IJ) = 0.0 CONTIN.UE WRITE (6,105) FORMAT('0*** INITIAL PROPERTIES OF SOIL LAYERS ***•/) CHECK THAT ALL AM'S ARE LESS 0.5 LSTOP = .FALSE. IF (NSUBDT.LE.10) GO TO 201 WRITE (6,1 35) NSUBDT FORMAT(•0***** NSUBDT IS TOO LARGE !! =•,14,' ***•*•) STOP RFNDT = FLOAT (NSUBDT) DDT = 0.01605*DT/RFNDT DO 14 1=1,N TEMP = DZ (I) TEMP = TEMP*TEMP AH (I) = DDT*PROP ( I , 15) *SK1R (I)/TEMP IF (AM (I).LT.0.5) GO TO 14 WRITE(6,109) I FORMAT (•0*** AM VALUE IS TOO HIGH FOR LAYER',13,' ***•) LSTOP = .TRUE. CONTINUE WRITE(6,107) FORMAT ('0 LAY.NO. K-S K1R-S THICK. •AM-S NSUBD-S'/) WRITE (6, 108) (I, PROP (I, 15) ,SK1R(I) ,H(I) ,AM(I) , NSUBD (I) ,I=1,N) FORMAT (1X,I6,4X,1P2E12.3,0PF12.2,1PE12.3,3X,0PI6) I F (LSTOP) STOP RETURN PORE PRESSURE REDISTRIBUTION STARTS HERE CONTINUE DO 37 K=1,NSUBDT NP = 0 DO 26 1=1,N NN = NSUBD(I) TEMP = DZ (I) TEMP = TEMP*TEMP AMI = DDT*PROP (1,15) *SK1R (I) /TEMP I F (AMI.LT.0.5) GO TO 17 LSTOP = .TRUE. GO TO 50 CONTINUE DO 26 IJ=1,NN NP = NP+1 IF (I.NE. 1) GO TO 19 IF (IJ.NE. 1) GO TO 21 CALCULATION FOR BOTTOM BOUNDARY IF (NBB. EQ. 0) GO TO 18 TW(1) = 2. *AMI* (W (2)-W (1) )+W (1) GO TO 18 IF ( I J . NE. 1) GO TO 21  LISTING OF LIQUEFACTION PROGRAM - AUG, 75  LIQPG 27  199 C  C 21 18 26 C  27 C C  CALCULATION FOR INTERFACE OF 2 LAYERS TEM 1' = DZ (I) /SK1F. (I) TEMP = DZ (1-1)/SK1R (1-1) TEMP = TEMP+TEM1 CONST = 2.*DDT/TEMP A = PROP(I-1.,15)/DZ(I-1) B = PROP (I, 15) /DZ (I) TW(NP) = CONST* (A*W (NP-1) + B*W (NP+1) - (A + B) *W (NP) )+W (NP) GO TO 18 CALCULATION FOR POINTS IN A LAYER TW(NP) = AMI*(W(NP-1) + W (NP+1)-2. *W (NP) ) + W (NP) CONTINUE CONTINUE CALCULATION FOR TOP BOUNDARY IF (NBT. EQ. 0) GO TO 27 NP = NP+1 TW(NP) = 2. *AMI* (W (NP-1)-W (NP) )+W (NP) CONTINUE  ADD DPP TO HT(I) TO GIVE W (I) NP = 0 DO 30 1=1,N DPI = DPP (I) *RFNDT NN = NSUBD (I) .IF (T.EQ. N) NN=NN+1 DO 30 IJ=1,NN NP = NP+1 IF (NP.EQ.LAND.NBB.EQ.O) GO TO 33 IF (NP.EQ.NTOT.AND. NBT. EQ.O) GO TO 33 IF (I.NE. 1. AND.IJ. EQ. 1) GO TO 32 C NORMAL SOIL POINTS, BOTTOM BOUND. PT. AND TOP BOUND. PT. W(NP) = TW(NP)+DPI GO TO 33 C INTERFACE 32 W (NP) = TW (NP) +0 . 5 * (DPI + DPP (1-1) ) 33 CONTINUE 30 CONTINUE 37 CONTINUE C FIND PP(I) WHICH IS AVERAGE OF W (I) FOR A LAYER NP = 1 DO 35 1=1,N NP = NP-1 NN = NSUBD (I) +1 FN = FLOAT (NN) A = 0.0 DO 36 IJ=1,NN NP = NP+1 36 A = A+TW (NP) 35 PP(I) = A/FN C WRITE (6,9999) (TW (K) , K= 1 , NP) C999 FORMAT (' OWT (HP) : , 1P5E14. 3, (/9X, 1P5E1 4. 3) ) C WRITE (6,9998) (DPP (K) , K= 1 , N) C998 FORMAT (' DPP (I) :»,1P5E14.3, (/9X,1P5E14.3)) 1  LISTING OF LIQUEFACTION PROGRAM - AUG, 75  LIQPG 200  C C997 C C996  120 C  101  102  103  104  150 C C 200 2001  WRITE (6, 9997) (W (K) , K=1 , NP) FOR MAT(' W(NP) :»,1P5E14.3, (/9X,1P5E14.3)) WRITE (.6,9996) (PP (I) ,T= 1, N) FORMAT (' PP (I) :», 1P5E14.3, (/9X,1P5E14.3)) RETURN END SUBROUTINE MAXVAL(TIME,UG1) LOGICAL*1 LDOF1,LLIN,LSTOP,LMCLG,LJRV,LCKFB,LIT,LPMAX COMMON . /LCONTL/LDOF1,LLIN,LSTOP,LMCLG,LJRV,LCKFB,LIT,LPMAX . /NCONTL/N,NDAMP,INCR,INTYP,NVOL,NBT,NBB . /CALVAL/DIRN (20) ,XD2 (20) ,ABXD2 (20) ,XD1 (2 0) ,XD0 (20) ,BD (20) , . BA(20) ,BB(20) ,GD2(20) ,GD1 (20) , GAMMA (20) ,GAMMAO(20) , GAMMAX (20) ,TAXY (20) , EVP (20) ,SK1R (20) ,DPP (20) ,PP (20) DIMENSION ACCX(20) ,GAMX(20) ,TMAX(20) ,EVPX(20) ,PPMX(20) , TFM (20,5) DATA ACCX/20*0.0/, GAMX/20*0.0/, TMAX/20*0.0/, EVPX/20*0.0/, PPMX/20*0.0/, TFM/100*0.0/, ACCINX/0.0/, T1/0.0/, TP/0.0/ IF (LPMAX) GO TO 200 TEM = ABS (UG1) IF (ACCINX.GE.TEM) GO TO 120 ACCINX = TEM T1 = TIME CONTINUE DO 150 1=1,N TEM = ABS (ABXD2(I)) IF (ACCX (I) .GE.TEM) ACCX(I) = TEM TFM (1,1) = TIME TEM = ABS (GAMMA (I) ) IF (GAMX (I) .GE.TEM) GAMX (I) = TEM TFM (1,2) = TIME TEM = ABS (TAXY (I) ) IF (TMAX (I) .GE.TEM) TMAX(I) = TEM TFM (1,3) = TIME TEM = ABS (EVP (I) ) IF (EVPX (I) .GE.TEM) EVPX(I) = TEM TFM (1,4) = TIME TEM = ABS (PP (I) ) I F (PPMX(I) .GE.TEM) PPMX(T) = TEM TFM (1,5) = TIME CONTINUE RETURN  GO TO 101  GO TO 102  GO TO 103  GO TO 104  GO TO 150  PRINT OUT MAXIMUM VALUES BEFORE CONTINUE WRITE (6,2001) TP,TIME FORMAT (  TERMINATION  LISTING OF LIQUEFACTION PROGRAM - AUG, 75  LIQPG 29 201  .  »0********************************************* • * MAXIMUM VALUES OCCURRED -',F8.3,' TO',F8.3,'  / SEC *'/  • ******************************************* ' LAY. TIME ABS ACC -TIME STRAIN TIME ', •STRESS TIME VOL STRAIN TIME PORE PRESS'/) DO 2000 1=1, N WRITE (6,2002) I, TFM ( 1 , 1) , ACCX (I) ,TFM (I , 2) , GAMX (I) ,TFM (I, 3) , TM AX (I) , TFM (1,4) ,EVPX(I) ,TFM (1,5) ,PPMX (I) FORMAT (14,1X,5 (0PF7. 2, 1PE11. 3) ) WRITE (6,2003) ACCINX,T1 FORMAT('0** MAX BASE ACC =',1PE11.3,' : TIME =',0PF5.2, ' SECS **') ;  2000 2002 2003 C C  RESET VARIABLES TO ZERO ACCINX = 0,0 DO 160 1=1,N ACCX(I) =0.0 GAMX(I) = 0.0 TMAX(I) = 0.0 EVPX(I) =0.0 PPMX(I) =0.0 160 CONTINUE TP = TIME RETURN END EXECUTION TERMINATED  $SIG  

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0063017/manifest

Comment

Related Items