UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

A non-linear dynamic finite element analysis Quong, Wayne 1983

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

Item Metadata

Download

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

Full Text

C.  A NON-LINEAR DYNAMIC FINITE ELEMENT ANALYSIS by WAYNE QUONG B.A.Sc.  U n i v e r i t y Of B r i t i s h Columbia, 1977  A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in THE FACULTY OF GRADUATE STUDIES Department Of C i v i l E n g i n e e r i n g  We a c c e p t t h i s t h e s i s as c o n f o r m i n g t o £heyrequired s t a n d a r d  THE UNIVERSITY OF BRITISH COLUMBIA J u l y 1983  ©  Wayne Quong, 1983  i  In p r e s e n t i n g requirements  this thesis f o r an  of  British  it  freely available  Columbia,  understood for  by  that  Library  shall  for reference  and  study.  I  for extensive copying of be  her  copying or shall  g r a n t e d by  not  be  of  The U n i v e r s i t y o f B r i t i s h 1956 Main Mall Vancouver, Canada V6T 1Y3  Date  DE-6  (3/81)  of  Columbia  make  further this  thesis  head o f  this  my  It is thesis  a l l o w e d w i t h o u t my  permission.  Department  the  representatives. publication  the  University  the  h i s or  f i n a n c i a l gain  the  I agree that  s c h o l a r l y p u r p o s e s may  department o r  f u l f i l m e n t of  advanced degree a t  agree t h a t p e r m i s s i o n for  in partial  written  i i ABSTRACT  A two-dimensional predicting  the  structures to behavior  stress  and  seismic  of  approach  f i n i t e element  the  permanent  loading  soil  method  of  displacements  i s presented.  i s modelled  analysis for  by  The  level  of both t h e shear s t r a i n and t h e mean normal  shear  modulus  relationships  dependency  governing i n i t i a l  is  based  with  stress.  on  the The2  hyperbolic  l o a d i n g and u n l o a d i n g b e h a v i o r ,  l e a d i n g t o a h y s t e r e t i c type energy bulk  inelastic  an i n c r e m e n t a l l i n e a r  i n which t h e tangent shear modulus i s v a r i e d  strain  of e a r t h  dissipation.  The  tangent  modulus i s v a r i e d w i t h t h e l e v e l of the mean normal  stress  o n l y , and h y s t e r e t i c e f f e c t s a r e not c o n s i d e r e d . The  i n c r e m e n t a l l i n e a r e q u a t i o n s of motion of t h e s t r u c t u r e  are s o l v e d u s i n g the Newmark s t e p - b y - s t e p i n t e g r a t i o n in  procedure  t h e time domain a l l o w i n g t h e s t r e s s e s and d i s p l a c e m e n t t o be  computed.  A f t e r each time  modulus a r e r e - e v a l u a t e d . hyperbolic  shear  step  the  tangent  shear  and  bulk  H y s t e r e t i c damping as a r e s u l t of the  stress-strain  law i s i n h e r e n t i n the model.  V i s c o u s damping may a l s o be i n c l u d e d . The a n a l y s i s i s a p p l i e d t o a number of dams and s l o p e s the  earthquake  induced  displacements  a r e compared w i t h those  p r e d i c t e d by a s i m p l e r Newmark s i n g l e degree plastic  analysis.  As w e l l , a comparison  p r e d i c t i o n of d e f o r m a t i o n of structure, a clay  dam  of  freedom  rigid  i s made w i t h M a k d i s i ' s  embankments.  For  a  clay  slope  t h e o v e r a l l d i s p l a c e m e n t s a r e of s i m i l a r o r d e r . struture  the non-linear  and  finite  element  For  results  i n d i c a t e t h a t t h e Newmark type methods a r e o v e r l y c o n s e r v a t i v e . The  more  rigorous  multi-degree  d i s t r i b u t i o n of d i s p l a c e m e n t s embankment t o be o b t a i n e d .  of freedom a n a l y s i s a l l o w s the  w i t h i n and on the s u r f a c e  of the  iv TABLE OF CONTENTS Page ABSTRACT  i i  LIST OF TABLES  vi  LIST OF FIGURES  vii  ACKNOWLEDGEMENTS  ix  Chapter 1  1  INTRODUCTION  1.1  P r e l i m i n a r y Remarks  1.2  Scope of T h e s i s  Chapter 2  CRITICAL REVIEW OF CURRENT METHODS IN THE ASSESSMENT OF SEISMIC PERFORMANCE OF EARTH STRUCTURES  2.1  Pseudo-static  2.2  Newmark A n a l y s i s Procedure  2.3  Seed Lee I d r i s s Proceduce  Chapter 3  Analysis  THE FINITE ELEMENT METHOD  3. 1  Introduction  3.2  Constitutive relationships  3.3  Chapter 4  3.2.1  Bulk modulus s t r e s s dependency  3.2.2  H y s t e r e t i c h y p e r b o l i c Shear stress-strain Relationship  21  F o r m u l a t i o n of s t r u c t u r e s t i f f n e s s , mass and damping m a t r i c e s . 3.3.1  P a r t i a l S t i f f n e s s matrices  3.3.2  Mass m a t r i x  3.3.3  Damping  matrix  NUMERICAL ANALYSIS OF NON-LINEAR DYNAMIC RESPONSE  4.1  General  4.2  E q u a t i o n of motion  44  4.3  I n c r e m e n t a l e q u a t i o n of motion  4.4  Step by Step I n t e g r a t i o n  4.5  Correction  4.6  Procedure f o r s t r a i n r e v e r s a l  4.7  Summary of procedure  Chapter 5  factor occurence  APPLICATIONS OF THE NON-LINEAR FINITE ELEMENT METHOD  5. 1  Introduction  5.2  S i n g l e Element Comparison Analysis  5.3  Dynamic Response a n a l y s i s of C l a y S t r u c t u r e s 5.3.1  5.3.2  Chapter 6  w i t h Newmark  Comparison w i t h t h e Newmark and the M a k d i s i - S e e d A n a l y s i s 5.3.1.1  Clay Slope  Comparison  5.3.1.2  C l a y Dam Comparison  Comparison of H y p e r b o l i c and E l a s t i c - P l a s t i c Shear s t r e s s - s t r a i n Laws  SUMMARY AND CONCLUSIONS  6.1  Summary  6.2  Conclusions  6.3  S u g g e s t i o n s f o r F u r t h e r Research  REFERENCES APPENDIX I D e r i v a t i o n of element s t i f f n e s s m a t r i x  vi  LIST OF TABLES Table I  S t a t i c and Dynamic S o i l  Properties  vii LIST OF FIGURES F i g u r e 2-1  Pseudo-Static  F i g u r e 2-2  Forces on a S l i d i n g  F i g u r e 2-3  I n t e g r a t i o n of E f f e c t i v e A c c e l e r a t i o n Time H i s t o r y  F i g u r e 2-4  T y p i c a l Shear moduli and Damping R a t i o s f o r Sands  F i g u r e 3-1  H y p e r b o l i c Shear s t r e s s - S h e a r S t r a i n R e l a t i o n s h i p  F i g u r e 3-2  H y p e r b o l i c R e l a t i o n s h i p under G e n e r a l  F i g u r e 4-1  Program  F i g u r e 5-1  A c c e l e r a t i o n Response and Time h i s t o r y f o r San Fernando Earthquake  F i g u r e 5-2  S i n l g l e F i n i t e Element Model  F i g u r e 5-3  Shear D e f o r m a t i o n R e l a t i o n s h i p  F i g u r e 5-4  S i n g l e F i n i t e Element - Newmark ' Displacement Comparison  F i g u r e 5-5  S i n g l e F i n i t e Element Model Newmark A n a l y s i s Displacement H i s t o r y Comparison  F i g u r e 5-6  F i n i t e Element G r i d of C l a y Slope  F i g u r e 5-7  F i n i t e Element G r i d of C l a y Dam  F i g u r e 5-8  Slope S t a b i l i t y A n a l y s i s on C l a y Slope  F i g u r e 5-9  Slope S t a b i l i t y A n a l y s i s on C l a y Dam  F i g u r e 5-10,  C l a y Slope N o n - l i n e a r - M a k d i s i - Newmark Displacement Comparison  F i g u r e 5-11  C l a y Slope N o n - l i n e a r F i n i t e Element Newmark A n a l y s i s D i s p l a c e m e n t H i s t o r y Comparison  F i g u r e 5-12  D i s p l a c e d G r i d of C l a y Slope  F i g u r e 5-13  N o n - l i n e a r - M a k d i s i - Newmark D i s p l a c e m e n t Comparison ( C l a y Dam)  F i g u r e 5-14  Clay Dam D i s p l a c e m e n t Time H i s t o r y  method Block  Loading  Flowchart  Structure  Structure  Structure  vi i i Figure  5-15  D i s p l a c e d G r i d of C l a y  Figure  5-16  E l a s t i c - P l a s t i c A p p r o x i m a t i o n of H y p e r b o l i c Curve  Figure  5-17  Hyperbolic - E l a s t i c P l a s t i c Displacement Comparison ( C l a y S l o p e )  Figure  5-18  Hyperbolic - E l a s t i c - P l a s t i c Displacement H i s t o r y Comparison  F i g u r e I-1a  Q u a d r i l a t e r a l Element  Figure  U n i t Square  I-1b  Dam  ix  ACKNOWLEDGEMENTS  The  author  wishes t o thank h i s p r i m a r y  advisor,  Professor  P e t e r M. Byrne f o r h i s u n f a i l i n g guidance and encouragement, and for  making c r i t i c a l c o n t r i b u t i o n s towards t h e c o m p l e t i o n  thesis. Professor  D.L.  Anderson  was  of  valuable  of the  assistance  d u r i n g t h e p r e p a r a t i o n and review of t h e t h e s i s . The  author  wishes t o e x p r e s s h i s g r a t i t u d e t o h i s  for t h e i r continued The  support  parents  and p a t i e n c e d u r i n g t h i s p e r i o d .  f i n a n c i a l a s s i s t a n c e p r o v i d e d by a N a t u r a l S c i e n c e s and  Engineering  Research C o u n c i l of Canada p o s t g r a d u a t e  is gratefully  appreciated.  scholarship  1  CHAPTER 1 INTRODUCTION 1.1  P r e l i m i n a r y Remarks The f r e q u e n t o c c u r r e n c e of d e s t r u c t i v e  the  earthquakes  during  past few decades has caused b i l l i o n s of d o l l a r s i n p r o p e r t y  damage and  loss  important  that  of  thousands  earth  of  lives.  structures  It  is  therefore,  l o c a t e d i n an a c t i v e  seismic  r e g i o n be d e s i g n e d t o w i t h s t a n d s a f e l y t h e e x p e c t e d  earthquake  motion  earthquake  f o r that  region.  In  measuring 6.6 on t h e R i c h t e r Fifty-eight  people  February  scale  earthquake,  occurred  in  California.  were k i l l e d , over two thousand i n j u r e d , ' and  1500 b u i l d i n g s were damaged beyond This  1971, a  said  to  be  C a l i f o r n i a i n t h e past  80  years,  safe  occupational  the most  severe  caused  an  to  levels. occur i n  estimated  five  b i l l i o n d o l l a r s i n damage. Of  particular  geotechnical  importance  was  f a i l u r e and near t o t a l c o l l a p s e of t h e Lower San  the  Fernando  and downstream movement of t h e Upper San Fernando Dam. the  earthquake  the water  level  below the c r e s t of the dam. upstream  shell  of  the  San  which  material.  p a r t of t h e downstream.  same  Prior to  occurred  i n the  Fernando Dam f o l l o w i n g the  earthquake r e s u l t e d i n a f r e e board of cracked  Dam,  i n the r e s e r v o i r was 35 f e e t  The s l i d e  Lower  partial  only  4  to  5  feet  of  At the same t i m e , t h e Upper dam which forms reservior  Fortunately  complex,  displaced  some  6  feet  t h i s movement d i d not cause a r e l e a s e  of water from t h e r e s e r v o i r of t h e Upper Dam.  If  i t had, the  2  resulting  overtopping  of  what remained of the Lower Dam  have caused c o n s i d e r a b l e damage and l o s s of l i f e  (Seed 1979).  F i v e y e a r s p r e v i o u s , based on e x t e n s i v e s t a t e procedures,  a  consulting  board,  a  might  of  the  d e s i g n agency, and  review  board had deemed the d e s i g n of the Lower San Fernando Dam adequate a g a i n s t any earthquake t o which i t c o u l d be But  the  margin  by  Fernando Dam  was  this  other  and  which  to  be  subjected.  the t o t a l c o l l a p s e of the Lower  a v o i d e d was recent  consequences, (Osaki 1966,  art  uncomfortably events  small.  with  Because  associated  Seed e t . a l . 1966)  San of  disastrous  t h e r e was  a  need  f o r a r e t h i n k i n g of the earthquake r e s i s t a n t d e s i g n p h i l s o p h y of earth  and  analysis  rock  fill  dams  and  the  development of a l t e r n a t e  procedures.  Up t o t h i s t i m e , the s t a n d a r d method f o r 40 y e a r s for  more  e v a l u a t i n g the s a f e t y of e a r t h dams or embankment s l o p e s to  earthquake f o r c e s has been the analysis.  In  this  method  p o t e n t i a l s l i d e mass i s horizontal  force,  so-called the  represented  determined  by  earthquake  pseudo-static  method  e f f e c t of the earthquake on a by the  c o e f f i c i e n t and the t o t a l weight of the the  or  an  equivalent  product slide  of mass.  a  static seismic Assuming  f o r c e a c t s permanently on the s l o p e m a t e r i a l i n  one d i r e c t i o n o n l y , and t h r u the c e n t r o i d the mass, c o n v e n t i o n a l s l o p e s t a b i l i t y a n a l y s i s are a p p l i e d t o e v a l u a t e the safety  a g a i n s t movement.  nature  of  A f a c t o r of s a f e t y of l e s s than u n i t y  i m p l i e s movement, but because oscillating  factor  of  the  transient  and  randomly  of earthquake motions, a c o l l a p s e ( i n terms  3  of  an e x c e s s i v e d i s p l a c e m e n t c r i t e r i a ) may  sense,  a  factor  of  safety  not o c c u r .  In  cannot determine adequate  this  seismic  performance. I t i s now g e n e r a l l y a c c e p t e d t h a t the magnitude  of r e l a t i v e  d i s p l a c e m e n t s i s a more l o g i c a l and b e t t e r  criterion  factor  performance of e a r t h  of s a f e t y f o r a s s e s s i n g the dynamic  structures.  The a l l o w a b l e d i s p l a c e m e n t s may  vary  than  from  a  the  few  i n c h e s t o a few y a r d s depending on the f u n c t i o n a l a s p e c t s of the s t r u c t u r e s concerned. for  predicting  Newmark (1965) f i r s t proposed a procedure  permanent  d u r i n g earthquakes u s i n g degree  of  freedom  deformations a  simplifying  approach.  in  a  more  elaborate  p r e d i c t i n g displacements.  earth  rigid  embankments  plastic  More r e c e n t l y , dynamic  a n a l y s e s of s o i l s t r u c t u r e s (Seed included  in  et. and  a l . , 1973) encompassing  single response  have  been  procedure i n  These p r o c e d u r e s w i l l be d i s c u s s e d i n  Chapter 2. Dynamic response a n a l y s e s as p r a c t i s e d today had i t s o r i g i n i n the p i o n e e r i n g attempts of Seed and ifniversity  of  California  at  his  Berkeley  co-workers to  at  the  predict  the  a c c e l e r a t i o n s i n h o r i z o n t a l s o i l d e p o s i t s by bedrock motions.  The  Seed  Approach  (Dynamic  s t r e s s p a t h Method) was  used t o e v a l u a t e the s e i s m i c performance of on soil  soil  deposits.  deposit  subjected  to  to  earthquake  structures  founded  For such s t r u c t u r e s , the performance of the act  as  earthquake  an  adequate  motions  and  foundation  base  when  the m o d i f i c a t i o n of the  bedrock motion as i t propogates through the s o i l  are  important  4  considerations.  It  has  been i d e n t i f i e d t h a t a m p l i f i c a t i o n or  d e - a m p l i f i c a t i o n (depending on the c h a r a c t e r i s t i c s of deposit)  occurs  a d d i t i o n , there top  of  the  seismic  waves  soil  t r a v e l through s o i l .  i s a tendency f o r the r e s u l t a n t  motion  at  In the  d e p o s i t t o be a f f e c t e d by the n a t u r a l frequency of  the d e p o s i t . a deposit  when  the  The predominant p e r i o d of the motion a t the t o p of  i s g e n e r a l l y longer  than the predominant p e r i o d of the  bedrock m o t i o n . The p r e d i c t i o n of s e i s m i c deposits  has  been  well  response of h o r i z o n t a l l y  researched  (Schnabel  S t r e e t e r e t a l , 1974; F i n n e t a l , 1977). that  soil  p r o p e r t i e s vary  assumptions  that  are  loading,  equation  of m o t i o n .  Often s o i l horizontally consider dams,  and  structures  layered  the  These methods vary  loss  methods  cannot  deposits,  be and  two or even t h r e e d i m e n s i o n s .  where  the  cross-sectional  variability  of  dimension.  Where  other  in  the  soil the  dimension  have  analyses been  assume  in  of  soil  used  of the  during  the  t o i n t e g r a t e the  modelled  adequately  i t becomes For  the  as  important to  example,  in  zoned  shape as w e l l as the s p a t i a l  properites  necessitate  a  second  t h i r d dimension i s much l a r g e r than the  two and p r o p e r t i e s i n t h i s t h i r d  two d i m e n s i o n a l  methods  made , the m o d e l l i n g  n o n - l i n e a r i t y and p o s s i b l e s t r e n g t h dynamic  These  a l , 1972;  i n the v e r t i c a l d i r e c t i o n and remain  uniform i n the h o r i z o n t a l d i r e c t i o n . simplifying  et  layered  a r e adequate.  considered  dimension  are  uniform,  The e f f e c t s of the t h i r d  i n two d i m e n s i o n a l  analyses  by  5 Makdisi cross  (1976).  The predominant p e r i o d of the two  s e c t i o n was a l t e r e d a c c o r d i n g t o a l e n g t h of s t r u c t u r e t o  height r e l a t i o n s h i p . to  dimensional  Dynamic response a n a l y s e s were g e n e r a l i z e d  two dimensions w i t h t h e development of the computer  QUAD4,  LUSH and FLUSH.  T h e i r method of a n a l y s i s i s based on an  e q u i v a l e n t l i n e a r e l a s t i c approach.  Linear e l a s t i c  not a l l o w computation of d i s p l a c e m e n t s conjunction  with  programs  methods  do  d i r e c t l y , but a r e used i n  l a b o r a t o r y and a d d i t i o n a n a l y t i c a l procedures  to determine d i s p l a c e m e n t s .  These p r o c e d u r e s w i l l be  discussed  i n g r e a t e r d e t a i l i n Chapter 2. Streeter  proposed t h e f i r s t t r u e n o n - l i n e a r a n a l y s i s where  a Ramberg-Osgood r e p r e s e n t a t i o n of t h e s t r e s s - s t r a i n b e h a v i o r of s o i l was used.  The a n a l y s i s can o n l y be used i n one d i m e n s i o n a l  total stress applications.  The development of a two d i m e n s i o n a l  n o n - l i n e a r dynamic a n a l y s i s i s t h e next l o g i c a l s t e p . h e r e i n , i s presented  1.2  The  work  w i t h t h i s aim i n mind.  Scope of t h e T h e s i s The  thesis  response a n a l y s i s  presents  a  to predict  n o n - l i n e a r f i n i t e element dynamic deformations  during  earthquake  f o l l o w i n g b a s i c assumptions a r e made w h i l e  formulating  loadings. The the model: 1)  Theory of l i n e a r i n c r e m e n t a l e l a s t i c i t y i s applicable.  2)  S o i l behaves  isotropically.  6 3)  Plane s t r a i n c o n d i t i o n s p r e v a i l .  4)  The earthquake ground motion i s i d e n t i c a l a t a l l p o i n t s a l o n g the base of the s t r u c t u r e .  The s t r e s s - s t r a i n b e h a v i o r of s o i l under shear  loading  is  assumed t o be a h y p e r b o l i c 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 t o that  used  by  Martin,  Finn  and  Seed  b e h a v i o r as proposed by Duncan e t . The  dynamic  analysis  (1975).  Volume change  a l . , (1980)  is  followed.  i s c a r r i e d out u s i n g an i n c r e m e n t a l time  s t e p approach where f o r each element the s t r e s s - s t r a i n c u r v e followed  in  an  i n c r e m e n t a l manner.  With  this  is  approach,  permanent d e f o r m a t i o n s can be e v a l u a t e d d i r e c t l y . The  seismic  travelling  excitation  shear  waves  a c c e l e r a t i o n time h i s t o r y element  is  which at  assumed are  the  to  defined  rigid  be by  base  of  vertically a  specified the  finite  model.  In Chapter 2, a d e s c r i p t i o n of c u r r e n t methods i n a s s e s s i n g seismic  instability  and  given.  Their  assumptions,  main  predicting  permanent d e f o r m a t i o n s i s analysis  procedures  and  l i m i t a t i o n s are c r i t i c a l l y r e v i e w e d . Chapter linear  3  dynamic  constitutive  and  4  finite  behavior,  deals  w i t h the f o r m u l a t i o n of the non-  element  analysis.  Assumptions,  s t i f f n e s s and damping p r o p e r t i e s of the  s o i l m a t e r i a l s , and the method of s o l u t i o n of the  equations  of  motion a r e p r e s e n t e d . The method i s used t o e v a l u a t e the s e i s m i c b e h a v i o r of some  7  typical  earth  Chapter  5.  results  structures. Displacement  obtained  investigators.  from  A brief  These  results  are  predictions  are  procedures  developed  summary,  conclusion,  f u r t h e r r e s e a r c h a r e p r e s e n t e d i n Chapter 6.  presented  compared by  in  against earlier  suggestions f o r  8  CHAPTER 2 CRITICAL REVIEW OF CURRENT METHODS IN THE ASSESSMENT OF SEISMIC PERFORMANCE EARTH STRUCTURES  The s e i s m i c performance of e a r t h embankments and dams received  considerable  attention  i n recent  c u r r e n t l y s e v e r a l methods b e i n g u s i n g for  this  thesis,  By way of i n t r o d u c t i o n t o  these  methods  practice  deformations  later  chapters  w i l l be examined h e r e i n .  assumptions made i n the f o r m u l a t i o n , a n a l y s i s  procedures,  l i m i t a t i o n s of each method w i l l be c r i t i c a l l y  2.1  There a r e  i n engineering  a s s e s s i n g t h e s e i s m i c s t a b i l i t y and p r e d i c t i n g  of e a r t h s t r u c t u r e s . in  years.  have  The and  reviewed.  P s e u d o - s t a t i c Method In  this  method  the  randomly  oscillating inertia  caused by an earthquake i s r e p r e s e n t e d horizontal  by an  forces  equivalent  static  f o r c e and a c o n v e n t i o n a l s l o p e s t a b i l i t y a n a l y s i s i s  a p p l i e d t o e v a l u a t e t h e f a c t o r of s a f e t y a g a i n s t c o l l a p s e . equivalent the d e s i g n  The  s t a t i c h o r i z o n t a l f o r c e i s d e t e r m i n e d by m u l t i p l y i n g seismic  coefficient  p o t e n t i a l s l i d e mass.  by  the t o t a l  weight  of the  The method i s i l l u s t r a t e d i n F i g u r e 2-1.  The e f f e c t i v e n e s s of t h e method depends upon, amongst other things,  the  selection  of a d e s i g n s e i s m i c c o e f f i c i e n t  V a l u e s f o r t h e s e i s m i c c o e f f i c i e n t i n t h e d e s i g n of worldwide  have  varied  from  .05  t o .20 .  value.  earth  dams  The s t a n d a r d  North  American p r a c t i c e i s t o use s e i s m i c c o e f f i c i e n t v a l u e s of 0.05,  FIG. 2-1  PSEUDO-STATIC  METHOD  10  0.10,  0.15, i n areas of low, medium and  selection  Prior to  accordance  these  with  field  1971, v e r y  few  The  experience  i n which  p s e u d o - s t a t i c approach  to  in  been  no  base t h e adequacy of the  failed  cause of f a i l u r e was  stability.  to  predict  the  the  build-up  of  pore  i n t h e embankment and t h e l o s s of s t r e n g t h  from these pore p r e s s u r e s . occur  designed  slope  of t h e Lower and Upper San Fernando Dam (Seed, 1979).  primary  pressures  dams  F o r t h i s reason, t h e r e had  p s e u d o - s t a t i c method f o r a s s e s s i n g s e i s m i c  failures  The  p r i n c i p l e s , have been s u b j e c t e d t o very  s t r o n g earthquake s h a k i n g .  The  seismicity.  of a s e i s m i c c o e f f i c i e n t v a l u e seems l a r g e l y based on  past experience.  real  high  The s e i s m i c  instability  water  resulting which  can  i n l o o s e c o h e s i o n l e s s s o i l s depends upon; t h e peak ground  a c c e l e r a t i o n and frequency  content  of  t h e earthquake  i n i t i a l s t r e s s s t a t e of t h e s o i l , r e s i d u a l pore water These  factors  cannot  be  represented  equivalent s t a t i c h o r i z o n t a l force.  motion,  pressures.  i n any r a t i o n a l way by a  Seed (1979)  suggests  that  the p s e u d o - s t a t i c method can o n l y be used w i t h any a s s u r a n c e f o r soil  materials  suffer  significant  s t r e n g t h or  s t i f f n e s s l o s s d u r i n g dynamic l o a d i n g .  Under these  conditions,  factors  of  that  safety  do  not  a g a i n s t movement i n t h e range of 1.0 t o 1.2  are c o n s i d e r e d adequate f o r s e i s m i c  2.2  stability.  Newmark A n a l y s i s Procedure Newmark (1965), and Seed (1966) have  both  criticized  the  concept of f a c t o r of s a f e t y as a means of a s s e s s i n g t h e probable  11  performance  of  an e a r t h dam d u r i n g an e a r t h q u a k e .  A f a c t o r of  s a f e t y of l e s s than u n i t y i n a s t a t i c a n a l y s i s i s not a c c e p t a b l e as i t i m p l i e s c a s t a s t r o p h i c safety  of  less  displacements.  However  factors  than u n i t y a r e a c c e p t a b l e i n dynamic  analyses,  s i n c e t h e earthquake f o r c e s a c t s f o r a s h o r t time and in  direction,  prediction  alternate  o n l y s m a l l d i s p l a c e m e n t s may occur and these may  be q u i t e a c c e p t a b l e . a  of  on  The p s e u d o - s t a t i c  method does not  t h e magnitude of d i s p l a c e m e n t s .  cannot be viewed as an adequate method  provide  Therefore i t  for evaluating  seismic  performance. In  t h e Rankine L e c t u r e of 1965, Newmark f i r s t o u t l i n e d t h e  b a s i c elements of a deformations  of  procedure  for predicting  the  potiential  an embankment s l o p e due t o earthquake f o r c e s .  I t was assumed t h a t s l o p e  f a i l u r e would be i n i t i a t e d and outward  movement would b e g i n  occur  potential  slide  to  mass  were  i f the  slide  were mass  modelled  by  ( F i g u r e 2-2).  reversed. along  Newmark  i t s failure  stop  when  on  a  the  inertia  proposed t h a t t h e movement of surface  could  t h e movement of a r i g i d block  be  adequately  on a i n c l i n e d  plane  By computing an a c c e l e r a t i o n a t which t h e i n e r t i a  f o r c e s become s u f f i c i e n t l y h i g h t o cause s l i d i n g integrating  forces  l a r g e enough t o overcome the y i e l d  r e s i s t a n c e and t h a t the movement would forces  inertia  the e f f e c t i v e  acceleration  to  begin  and  on t h e r i g i d b l o c k i n  excess of t h i s y i e l d a c c e l e r a t i o n as a f u n c t i o n of time f o r the duration  of  t h e earthquake motion ( F i g u r e 2-3),  v e l o c i t i e s and  u l t i m a t e l y d i s p l a c e m e n t s of t h e r i g i d block c o u l d be determined. Newmark  (1965)  presented  a  chart  for  computing  such  12  FIG. 2-2  FORCES O N  SLIDING  BLOCK  13  4  Time  FIG. 2-3  INTEGRATION OF  EFFECTIVE ACCELERATION TIME HISTORY  14  displacements.  In  displacements  the  were  normalized to a  computed  maximum  v e l o c i t y of 30 i n / s e c . that  the  development  earthquake  other  The  proposed  by  Newmark  used  For  of  an  (and  a  earthquake .5g  and  the  motions  a  maximum  e s s e n t i a l l y the same  not be n e c e s s a r i l y  embankment  true  e s t i m a t e of the  slope,  the  equations  shown on the c h a r t t o g i v e an upper  acceleration  are  the  The  maximum  motion.  chart,  conservative  be used.  earthquake  have  T h i s may  bound f i t t o the data) may values  four  the  narrow s c a t t e r i n the data i n d i c a t e s  motions  earthquakes.  permanent d e f o r m a t i o n  using  a c c e l e r a t i o n of  number of s i g n i f i c a n t p u l s e s . for  of  The  maximum  velocity  and  v a l u e s a p p r o p r i a t e t o the d e s i g n resistance  coefficient  is  d e f i n e d as the v a l u e of the h o r i z o n t a l s e i s m i c c o e f f i c i e n t which will  give  a f a c t o r of s a f e t y e q u a l t o u n i t y i n a p s e u d o - s t a t i c  a n a l y s i s of the embankment.  The  resistance coefficient  can  be  c o n s i d e r e d the y i e l d a c c e l e r a t i o n of the s l i d e mass. Evidence  had been p r e s e n t e d t h a t s l i p i n dense sands, when  subjected to uniform a c c e l e r a t i o n s , occurs along a t h i n zone.  The  s l i d i n g mass may  r e s t i n g on a i n c l i n e d p l a n e .  be c o n s i d e r e d analogous t o a b l o c k Shaking  t e s t s performed by Goodman  and Seed (1966) on s m a l l s c a l e embankments of dry demonstrated  the  surface  dense  sands,  v a l i d i t y of the fundamental p r i n c i p l e s of the  Newmark approach. Refinements t o a l l o w f o r throughout  the  the  variations  in acceleration  embankment and s l i d e mass were proposed by Seed  and M a r t i n (1966), Ambraseys and  Sarma  (1967),  and  Seed  and  15  Makdisi  (1978).  analysis  on  On  the  basis  embankments  acceleration  time  of  two  subjected  history,  to  average  dimensional a  given  response  earthquake  induced a c c e l e r a t i o n  time  h i s t o r i e s f o r a number of p o t e n t i a l s l i d i n g masses were computed by Seed and M a k d i s i . average  induced  procedure.  Deformations  accelerations  are  time  determined  histories  Design c u r v e s p r e s e n t e d by  Seed  using  and  and  the  Newmark's  Makdisi  show  good agreement w i t h Ambraseys and Sarma r e s u l t s . The approach  determination  d i s p l a c e m e n t v a l u e s u s i n g a Newmark  i s s t r a i g h t f o r w a r d enough.  freedom system model.  of  A complex  multi-degree  of  i s p r e s e n t e d by a s i m p l e s i n g l e degree of freedom  But  in  doing  so,  i n t e r p r e t a t i o n of r e s u l t s .  there  arises  the  problem  Where does t h i s d i s p l a c e m e n t  of  occur?  I t can be seen e i t h e r a s , the h o r i z o n t a l movement or the s l i d i n g movement a l o n g a s l i p s u r f a c e of the f a i l u r e mass. the  Newmark  approach  displacement.  gives  Ideally,  some  maximum  idea  of  displacement  In any the  as  case, average  well  as  a  d i s t r i b u t i o n of d i s p l a c e m e n t s w i t h i n an embankment i s d e s i r e d .  2.3  Seed Lee I d r i s s A n a l y s i s Procedure Seed and h i s co-workers  Berkeley procedure method undergone account  have  developed  at the a  University  comprehensive  of  California,  dynamic  analysis  for p r e d i c t i n g deformations i n earth s t r u c t u r e s . initially  proposed  refinements for  dynamic  (Seed forces  by et.  Seed  (1966)  al.,1973),  and  later  endeavours  induced by the earthquake and  The has to the  16  e f f e c t of s t i f f n e s s l o s s due t o dynamic c y c l i n g . been used response  with of  a  reasonable  success  number  earth  of  in  back  The method has calculating  dams (Seed e t .  a l . , 1975).  Commonly r e f e r r e d t o as the Seed Dynamic s t r e s s p a t h the  p r o c e d u r e i s summarized by the f o l l o w i n g a)  the  Appproach,  steps:  S e l e c t a p p r o p r i a t e c r o s s s e c t i o n of e a r t h s t r u c t u r e t o be  used  in  a n a l y s i s and model w i t h a f i n i t e element  grid. b)  Determine a p p r o p r i a t e s t a t i c and dynamic p r o p e r t i e s of the  s o i l t o used i n s t a t i c and dynamic f i n i t e  element  a n a l y s i s. c)  Determine the s t a t i c s t r e s s e s which e x i s t e d b e f o r e the earthquake.  Insitu  stresses  are  evaluated using a  c o n v e n t i o n s t a t i c f i n i t e element a n a l y s i s . d)  S e l e c t a d e s i g n earthquake time h i s t o r y .  e)  Compute the dynamic shear s t r e s s selected  elements  within  the  time  histories  for  cross section using a  t w o - d i m e n s i o n a l computer response a n a l y s i s . f)  S u b j e c t r e p r e s e n t a t i v e samples of the material  to  the  combined  effects  s t r e s s e s and the superimposed dynamic  soil  structure  of i n s i t u shear  static  stresses  and d e t e r m i n e t h e i r e f f e c t s i n terms of development of potential g)  From  the  strains. knowledge  of  the  soil  deformation  and  18  shown i n F i g u r e 2-4. These c u r v e s a r e used i n an i t e r a t i v e the  e q u i v a l e n t l i n e a r method.  procedure  known  as  In t h i s method, s u c c e s s i v e l i n e a r  problems a r e s o l v e d u n t i l the modulus and damping c o r r e s p o n d s t o the the  average dynamic shear s t r a i n max  shear s t r a i n ) .  The  ( u s u a l l y taken as 65 p e r c e n t of  f i n a l linear solution with  c o m p a t i b l e s o i l p r o p e r t i e s i s taken as true  non-linear  response.  The  an  strain-  approximation  equivalent  linear  to  a  solution  c o r r e s p o n d s t o the assumption of an e l l i p t i c h y s t e r e s i s l o o p f o r c y c l i c loading.  C u r r e n t l y , the computer  programs  QUAD4,  LUSH  and FLUSH based on the e q u i v a l e n t l i n e a r method, are used i n two dimensional  dynamic  response a n a l y s i s .  QUAD4 and LUSH are the  most f r e q u e n t l y used programs f o r s o i l s t r u c t u r e s such as s l o p e s and e a r t h dams. dynamic  FLUSH i s commonly  soil  structures  used  interaction  response of embedded n u c l e a r r e a c t o r  for  the  problems  solution such  structures  to  as  of the  earthquake  loading. Since  the  final  iteration  with  strain  compatible s o i l  p r o p e r t i e s i s p u r e l y e l a s t i c , the permanent d e f o r m a t i o n s  caused  by earthquake s h a k i n g cannot be computed d i r e c t l y from t h i s type of  analysis.  As i n a l l e l a s t i c a n a l y s e s the f i n a l d e f o r m a t i o n s  return to zero. procedure from  has  computed  computed f i e l d they  To c i r c u m v e r t t h i s  a  laboratory  been developed ( s t e p s f t o h) t o p r e d i c t dynamic  elastic are  shortcoming,  shear  strains used  for  bear  stress no  histories.  strains  While  the  r e l a t i o n t o s t r a i n s i n the  d e r i v i n g the  s t r a i n compatible s o i l  FIG. 2-4  TYPICAL SHEAR MODULI AND DAMPING RATIOS FOR SANDS  20  properties.  Stresses obtained  from  the a n a l y s i s  using  these  s t r a i n c o m p a t i b l e p r o p e r t i e s a r e assumed t o be r e p r e s e n t a t i v e of stresses  i n the ground.  As i n d i c a t e d i n s t e p f , t h e computed  s t r e s s e s a r e used t o e s t i m a t e lies  a  serious  inconsistency  computed a r e c o n s i d e r e d obviously  there  There  earth  accurate  are  i n t h e Seed p r o c e d u r e ; s t r e s s e s while  the  with  the  I t cannot  not  are not,  the  s t r u c t u r e due t o a pseudo resonance e f f e c t . period  of  the  uncommon  response T h i s can  earthquake  motion  n a t u r a l p e r i o d of the e a r t h s t r u c t u r e .  occur  in materials,  properties are highly non-linear is  strains  l i n e a r methods may o v e r e s t i m a t e  pure l i n e a r e l a s t i c systems, resonance i s a event.  Herein  two other p o s s i b l e shortcoming t o t h i s method.  occur when t h e predominant coincides  deformations.  i s a one t o one c o r r e s p o n d e n c e .  F i r s t , equivalent of  permanent  that  computed  s t r e n g t h r e s i s t a n c e of t h e s o i l .  real where  and  possible  the s t i f f n e s s  s t r a i n dependant. stresses  In  exceed  Second, i t the dynamic  T h i s cannot a c t u a l l y happen as  at most the dynamic s t r e s s can e q u a l t h e dynamic r e s i s t a n c e . The able  to  u l t i m a t e aim of a dynamic response a n a l y s i s predict  permanent d e f o r m a t i o n .  d i r e c t l y by a n o n - l i n e a r  i s to  be  T h i s can o n l y be done  a n a l y s i s i n the time  domain.  A  non-  l i n e a r f i n i t e element a n a l y s i s i s p r e s e n t e d i n t h e f o l l o w i n g two two  chapters.  21  CHAPTER 3 THE FINITE ELEMENT METHOD 3.1  Introduction Several mathematical  evaluate  the seismic  techniques  performance  of  s t r u c t u r e s founded on s o i l d e p o s i t s . has to  been  earth  developed  to  embankments, and  The f i n i t e  element  method  been commonly used i n g e o t e c h n i c a l e n g i n e e r i n g problems due t h e ease and a c c u r a c y w i t h which geometry  properties  can be m o d e l l e d .  s o i l behavior, i n i t i a l the  finite  and  varying  nonlinearity  in  stress-strain  s t a t i c s t r e s s e s , and boundary  element  soil  More r e c e n t l y , t h e method has been  used t o e v a l u a t e t h e e f f e c t s of  As  have  conditions.  method i s t h e o n l y a v a i l a b l e t e c h n i q u e  which w i l l a l l o w a r i g o r o u s assessment of the n o n - l i n e a r i t y s o i l , t h i s method w i l l adopted i n t h i s The  finite  element  method  of  thesis.  has  consequence of t h e advent of h i g h - s p e e d  been  developed  digital  as  computers  a  and  has been e x t r e m e l y s u c c e s s f u l i n s o l v i n g many s t a t i c and dynamic problems  i n continuum  mechanics.  I t s application to s t a t i c  a n a l y s e s of e l a s t i c c o n t i n u a has been d e s c r i b e d Clough  (1962)  and  and P e n z i e n (1975). The f i n i t e discretization  Wilson  and  i t s e x t e n s i o n t o dynamic a n a l y s e s by Clough T h i s method w i l l be b r i e f l y o u t l i n e d h e r e .  element method may be d e s c r i b e d as proceduce  whereby  assemblage of d i s c r e t e e l e m e n t s . an i n f i n i t e  by  a  numerical  continuum i s i d e a l i z e d as an The n u m e r i c a l t e c h n i q u e a l l o w s  degree of freedom system  to  be  transformed  to a  22 finite  degree of freedom system.  i n terms i n element s i z e predictions  is  selection  displacements,  will  stresses,  permit  system  accurate  and  strains  of the  The d i s p l a c e m e n t f o r m u l a t i o n of the f i n i t e  element  method  actual  of  and  Proper modeling of the  continuum.  employed  herein.  distribution  within  displacements  This an  such  compatibility  and  displacement  field  assumes  element  that  in  been  internal  terms  certain  completeness has  an  of  required  are  displacement the  conditions  satified.  assumed  for  nodal  any  Once  the  element  of  p r e s c r i b e d geometry and once the c o n s t i t u t i v e p r o p e r t i e s of elements  are  determined,  v i r t u a l work theorem t o element.  The  it  derive  stiffness  is the  matrix  on  the  p o s s i b l e w i t h the a i d of the stiffness  matrix  represents  the  of  the  stiffness  p r o p e r t i e s a s s o c i a t e d w i t h the d i s p l a c e m e n t s at the nodes of the element.  The  obtained  by the proper a d d i t i o n of i n d i v i d u a l element s t i f f n e s s  m a t r i c e s by  stiffness  the  matrix  direct  of  stiffness  the  entire  method.  continuum  Solution  s t i f f n e s s equations s a t i s f i e s e q u i l i b r i u m ( i n a global The  advantage of t h i s d i s c r e t e m a t h e m a t i c a l  t h a t f o r the dynamic problem be  of  is  the  sense.)  formulation i s  the e q u i l i b r i u m of the  system  may  e x p r e s s e d by a s e t of o r d i n a r y d i f f e r e n t i a l e q u a t i o n s r a t h e r  than a set of p a r t i a l d i f f e r e n t i a l e q u a t i o n s , w h i l e i n a problem of points  static  the p a r t i a l d i f f e r e n t i a l e q u a t i o n s are reduced t o a set  algebraic equations. in  the complete  I f the d i s p l a c e m e n t s  of  a l l nodal  assemblage are d e s i g n a t e d by the v e c t o r  23  {r} , t h e c o r r e s p o n d i n g nodal f o r c e s by the v e c t o r {R}, and t h e stiffness  matrix  of  t h e e n t i r e system by the m a t r i x [K] ; the  s t a t i c e q u i l b r i u m e q u a t i o n may be expressed  i n the form  [K] {r}= {R}  The  finite  (3-1)  element  numerical techniques. prescribed  method has many advantages over o t h e r  Different  or  coupled  themselves.  e l a s t i c systems by u s i n g  displacement  functions,  elements  with  increases.  This  a c c u r a c y can be o b t a i n e d . the  Any  displacement,  I t can be shown t h a t properly  selected  the f i n i t e element method converges t o  the exact s o l u t i o n as t h e number of elements system  can be  boundary c o n d i t i o n can be handled, and t h e  boundaries can be v e r y i r r e g u l a r i n shape. for  properties  from one element t o t h e next and/or can have v a r y i n g  p r o p e r t i e s w i t h i n t h e element stress,  material  indicates  used  t o model  a  t h a t any d e s i r e d degree of  The g o v e r n i n g c o n s t r a i n t  here  being  c o m p u t a t i o n a l c o s t s a s s o c i a t e d w i t h the i n c r e a s e d n u m e r i c a l  o p e r a t i o n s performed  on l a r g e r s t i f f n e s s m a t r i c e s and v e c t o r s .  Two d i m e n s i o n a l systems may be r e p r e s e n t e d by various triangle  shapes. (CST),  The has  s i m p l e s t , t h e t h r e e node c o n s t a n t been  investigators.  Displacements  vary  through  linearly  function results element,  hence  elements  used  extensively  i t s name  strain earlier  w i t h i n t h i s element a r e assumed t o  t h e element.  i n constant  by  of  strain CST.  This linear and  Higher  stress order  displacement within  the  displacement  24  displacement functions  f o r t r i a n g l e s , and q u a d r i l a t e r a l elements  have been i n t r o d u c e d i n t o s o i l  dynamic  Miller  respectively.  (1971), and Seed (1969)  The  plane  quadrilateral  the p r e s e n t a n a l y s i s . eight  at  each  'isoparametric' interpolation  the  derived  functions  two four  from  t o map  translational corner the  use  The  of  the  the  same shape  element.  The  element need not be r e c t a n g u l a r but can be of any  shape.  structures  t o be m o d e l l e d w i t h o u t d i f f i c u l t y .  This  feature  allows  irregular  shaped  As w e l l , the four  element may be s p e c i f i e d such t h a t the l o c a t i o n of any two  a d j a c e n t nodes i s the same, making the element t r i a n g u l a r . stiffness as  of term  the q u a d r i l a t e r a l element  arbitrary  node  and  degrees  nodes.  as a r e used t o d e f i n e d i s p l a c e m e n t s w i t h i n quadrilateral  Finn  The p l a n e q u a d r i l a t e r a l element p o s s e s s e s  of  is  by  i s o p a r a m e t r i c element i s used i n  degree of freedom, namely  freedom  analyses  a  matrix  four  node  isoparametric  of a t r i a n g u l a r element, o b t a i n e d t r e a t i n g i t element  elements,  words, by s p e c i f y i n g  and  a  CST.  for  the  procedure  the  fourth  node,  Or i n o t h e r  the  From t h i s r e s u l t , q u a d r i l a t e r a l s  plane  for  the node number of the t h i r d node t o be the  can be used i n the a n a l y s i s . matrix  following  reduces t o t h a t of a CST.  same of the node number of actually  The  The f o r m u l a t i o n of  quadrilateral  element  and/or CST's  the  isoparametric  is  stiffness element  is  p r e s e n t e d i n Appendix I . The dynamic thesis consists  finite  element  analysis  of t h r e e major p a r t s  presented  is  this  25  a.  M o d e l l i n g the c o n s t i t u t i v e b e h a v i o r of  b.  soil materials.  F o r m u l a t i o n of s t i f f n e s s , mass and damping m a t r i c e s  c.  S o l u t i o n of the e q u a t i o n s of  Parts of  a  motion  and b w i l l be d i s c u s s e d i n the r e m a i n i n g s e c t i o n s  t h i s chapter.  P a r t c w i l l be d i s c u s s e d i n c h a p t e r  4.  26  3.2  Constitutive  relationships  A fundamental behaves  assumption made i n the a n a l y s i s i s t h a t  isotropically.  This  allows  two-dimensional  stress-  s t r a i n b e h a v i o r t o be d e s c r i b e d by two e l a s t i c parameters each increment. and for  In the present a n a y l s i s , the  shear  Dynamic propagating  loading shear  should  is  waves.  essentially  parameters  due  to  vertically  The shear waves induce dynamic shear  deformations  deformations  displacements.  G  reasons.  s t r e s s e s and shear shear  during  modulus  t h e b u l k modulus B were s e l e c t e d as t h e e l a s t i c the f o l l o w i n g  soil  make  up  on a  B e a r i n g i n mind  the d e p o s i t .  great that  part  the  of  'failure  Therefore, the  overall  condition'  be based on t h e magnitude of r e s i d u a l d i s p l a c e m e n t s t h a t  can be t o l e r a t e d .  A proper  stress-strain  law  should  control  shear  deformation  i n o r d e r t h a t a p r e d i c t i o n on f a i l u r e can be  made.  Based on t h i s argument, i t was d e c i d e d t h a t shear modulus  be taken as one of the e l a s t i c parameters, which can be if  near  failure  conditions  occurs  reduced  using a simple h y p e r b o l i c  model. D u r i n g dynamic l o a d i n g a  soil  element  often  experiences  near f a i l u r e c o n d i t o n , i e . , reaches near maximum shear s t r e n g t h . During  this  condition,  higher  shear  deformation  w i t h o u t an accompanying i n c r e a s e i n v o l u m e t r i c bulk  modulus  remains  constant.  l a b o r a t o r y t e s t s on c l a y and sand. and  v  can  be  can  strain  occur  as the  T h i s b e h a v i o r i s observed i n I n a s t r e s s - s t r a i n model,  used as t h e v a r i a b l e s .  E  And from e l a s t i c i t y , we  27  know, B=  Near  E 3(1-2*)  the  (3-2)  c o n d i t i o n , where E i s r e d u c i n g , i f v i s kept  failure  c o n s t a n t , from e q u a t i o n  (3-2) B w i l l  increase i n volumetric s t r a i n . to  keep  B  constant  relationships  been w e l l developed al.  resulting  i n an  T h i s can be a v o i d e d by v a r y i n g v  (Byrne e t a l 1982).  the bulk modulus d i r e c t l y as one well,  reduce,  I t i s s i m p l i e r t o use  the e l a s t i c  f o r determining  parameters.  As  bulk modulus v a l u e s have  from e x t e n s i v e i n v e s t i g a t i o n s by Duncan e t .  (1980). In  strain  order to a r r i v e at expressions behavior  of  an  idealized  describing  soil  the  element under g e n e r a l  l o a d i n g i t i s assumed t h a t h y p e r b o l i c shear s t r e s s - s h e a r relationship  similar  to  that  used  by  (1975), and t h a t volume change b e h a v i o r as et.  a l .  (1980) i s f o l l o w e d .  stress.  .strain  M a r t i n , F i n n and Seed proposed  by  Duncan  Tangent v a l u e s of shear modulus  i s v a r i e d w i t h t h e l e v e l of both t h e shear s t r a i n and normal  stress-  the mean  V a l u e s of bulk modulus i s v a r i e d w i t h t h e l e v e l  of mean normal s t r e s s o n l y . in the f o l l o w i n g s e c t i o n s .  These r e l a t i o n s h i p s  are  discussed  28  3.2.1  Bulk  Studies volume  with  by  change  accurately  with  pressure  K -P  B=  b  a l . , (1980)  o f most that  pressure,  mobilized.  dependency  et.  by a s s u m i n g  increasing  normal  stress  Duncan  behavior  confining  strength  Modulus  soils  shown  c a n be m o d e l l e d  the bulk  modulus  and i s independent  Values  o f B have  confining  pressure  c a n be a p p r o x i m a t e d  have  that the  reasonably  of the s o i l  varies  of the percentage  been  found  and  to  i n terms  increase  o f t h e mean  by t h e e q u a t i o n ,  ( ^ f  a  of  (3-3)  Pa'  where, B  = tangent  Kb  modulus  = bulk  modulus  parameter  m = bulk  modulus  exponent  o =  normal  mean  m  P  bulk  a  = atmospheric  effective  stress  pressure,  included  non-dimensional  equation,  consistent  o  The  tangent  with  bulk  2  expressed  a  in units  a n d B.  modulus  Ag, + A o + A p 3Ae  B=  m  t o have  i s defined  by  (3~4)  3  v  where  Aa,, A a , and A a 2  principal  stresses,  volumetric  strain,  Ae = A e + A e v  x  alternatively  a r e the changes  3  and  which  Ae  v  i s  the  f o r the plane  in  values  corresponding strain  condition  of  the  change i n i s ,  (3"5)  y  equation  the  ( 3 - 4 ) c a n be w r i t t e n ,  29  (3-6)  B=  where, A o  i s t h e change i n mean normal  m  stress.  As w i l l be d i s c u s s e d i n g r e a t e r d e t a i l i n c h a p t e r non-linear  response  of  s o l v e d u s i n g the Newmark This  earth  4, t h e  s t r u c t u r e t o s e i s m i c motions i s  step-by-step  integration  procedure.  procedure computes d i s p l a c e m e n t s and s t r a i n s a t the end of  consecutive  time  intervals  At  assuming  the  the  p r o p e r t i e s have been h e l d c o n s t a n t d u r i n g t h e increment.  material At t h e  end of a time i n t e r v a l t h e v a l u e s of bulk moduli must be updated for  t h e e f f e c t s of s t r e s s change.  The bulk modulus f o r an s o i l  element a t any time can be determined as f o l l o w s :  1.  At time T, t h e b u l k modulus B from e q u a t i o n 3-3 as i s known.  2.  D u r i n g time s t e p At u s i n g bulk modulus v a l u e s calculate e  3.  v  a t time T + A t .  At time T + At c a l c u l a t e the new b u l k modulus B.  30  3.2.2  H y s t e r e t i c h y p e r b o l i c Shear Relat ionship  stress-strain  A comprehensive survey of t h e f a c t o r s a f f e c t i n g  t h e shear *  moduli  of  soils  and e x p r e s s i o n s f o r d e t e r m i n i n g t h i s p r o p e r t y  have been p r e s e n t e d by H a r d i n and D r n e v i c h study,  an  empirical  equation  v a l u e s of maximum shear modulus  (1970).  In  their  was p r e s e n t e d t o determine t h e G  •  max  Their  equation  i s as  follows: G  =  m a x  320.8-Pn ( 2 . 9 7 3 - e ) - (OCR)° o ' / (1+e) Pc7' 2  f  (3-7)  ;  m j  (  where, = Maximum shear modulus  Gmax  e  = void  OCR  ratio  = overconsolidation ratio a = parameter  t h a t depends on t h e p l a s t i c i t y  index  of the s o i l a  = mean normal e f f e c t i v e  m  stress.  For c l a y s , Seed and I d r i s s used an e q u a t i o n of t h e form, ( S )• ( c o n s t a n t )  Gmax=  where S  u  (3-8)  0  i s t h e u n d r a i n e d s h e a r i n g s t r e n g t h of t h e c l a y  Laboratory  and  in-situ  test  data  performed  by  i n v e s t i g a t o r s have found the c o n s t a n t v a l u e t o v a r y from  several 1000 t o  3000. For t h e i n i t i a l strain  l o a d i n g phase the h y p e r b o l i c shear  stress-  r e l a t i o n s h i p f o r m u l a t e d by Kondner and Z e l a s k o (1963) t o  31  model  the response  The  same  triaxial shown  of g r a n u l a r  hyperbolic tests  this  behavior. For  dimensional  loading,  reversal The  starting  T  =  1  Idriss  assumed the  i s used et.  isotropic  initial  m o x Tmax  +  T  is  for clays,  a l .  (1978)  response  up  a n d i s shown curve  as  used. cyclic  on c l a y s  has  stress  i s  i n general  two  shear  material,  on t h e s t r e s s - s t r a i n  + (T~  shear  of the s t a t i c  by t h e e q u a t i o n ,  point  G  by  The e f f e c t  the  to given  i n simple  relationship  performed  included.  soil  to  the  in Figure  i s (0, r  5 t  first 3-1.  ).  (3~9)  ST  -y  "max /mox T  ult1  'st  Rf  where, T Tmax r| 0  =  m  a  x  shear  strain  shear  = the f a i l u r e  ratio  The  - sign  the  negative  and + s i g n  Tmax  =1  )(T* )  of the stress  are applied  direction  i s the second  m a x  strength  = the s t a t i c  ST  Rf  7  stress  = the shear  t  T  = shear  2  +  of the  soil  to loading  "  U  of s t r a i n , e  y  )  as determined  c = x  the normal  strain  i n the  x-direction  e =  the normal  strain  i n the  y-direction  7  = the shear  Tmax  i s given  strain  and  i n t h e x-y  t h e same  sign  by, (3-10)  2  where,  y  i n the p o s i t i v e  respectively.  invariant  y  soil  as 7  plane. xy  .  32  Shear strain  FIG. 3-1  HYPERBOLIC SHEAR STRESS-SHEAR STRAIN  RELATIONSHIP  33  For u n l o a d i n g  and r e l o a d i n g ,  hysteretic  behavior  development  of a  analysis  t h e Masing  has been  used  one-dimensional  by  Lee  dynamic  f o r s a t u r a t e d sand d e p o s i t s .  skeleton  (1977),  stress  B r i e f l y described loading  of  i n the  effective  the Masing c r i t e r i o n assumes t h a t i f t h e i n i t i a l or  type  ( 1 9 2 6 )  here, curve,  ( 3 - 9 ) can be  c u r v e , which i s d e s c r i b e d by e q u a t i o n  r e p r e s e n t e d by, T  =  f(7max  > st r  Then t h e u n l o a d i n g  where ( 7 , r r  and  where r t s  r  equation  branches  of a  m a x  hysteretic  loop  a r e t h e same  c u r v e w i t h both t h e s t r e s s and s t r a i n s c a l e s  The  . Further  curve  stress-strain  means t h a t t h e u n l o a d i n g  ( 3 - 1 2 )  by a f a c t o r of two and t h e o r i g i n  G  from, ( 3 - 1 2 )  ) i s the l a s t r e v e r s a l point i n the  reloading  point.  ( 3 - 1 1 )  or r e l o a d i n g c u r v e can be o b t a i n e d  Geometrically  skeleton  i s s e t t o zero  - f(7.  2  plot.  )  tangent  translated  by  to the r e v e r s a l  modulus a f t e r s t r e s s r e v e r s a l i s equal t o  t o t h i s Lee assumes t h a t ,  described  equation  ( 3 - 1 2 )  i f the  stress-strain  i n t e r s e c t s an e x t e n s i o n of  the s k e l e t o n c u r v e , t h e s t r e s s - s t r a i n p a t h f o l l o w s t h e c u r v e u n t i l t h e r e i s a r e v e r s a l of l o a d i n g a g a i n . strain  curve  increased  intersects  I f the s t r e s s  t h e c u r v e of the p r e v i o u s  the s t r e s s - s t r a i n path then f o l l o w s  the l a t t e r  skeleton  load c y c l e , stress-strain  curve. The  stress-strain  model  requires  that  previous  r e v e r s a l p o i n t s f o r each s o i l l a y e r or element be  stress  recorded,  in  34  order  to  previous have  determine previous  unloading  occurred.  performed may  on  be  whether  Because  sands  or  or  reasonable  the  present  reloading  very  clays  more  intersection  to  to  of  the  stress  strain  little  experimental  verify  the  simplify  the  Lee's  skeleton  or  curves  work has  been  assumptions,  reloading  and  it  unloading  behavior.  In described  r  Then  by e q u a t i o n  =  it  f ( 7  is  obtained  *  m  ,±r,  w  assumed  =  a T  (3-9)  a  a  .  The  as  by,  and  reloading  curve  can  be  while curve  of  point  value  Tui  A  + r  asymptote  hyperbolic  the  A  is  unloading  .  and Upon  equal  with  A  in Figure  A  the  of  the  from t h a t  3-2.  reloading  of  +  C u r v e OA i s  the  is  r . B  shear  point.  (3-14),  reversal  u ) t  undeformed  corresponding  at  The given  where  (3-14)  the  has  the  at  the  soil  equal  to  point  B the  new  stress-strain  strength  T  its  Equation  a  another to  from  by e q u a t i o n  (7 , T ) .  origin  and  7 ,  given is  r  strained  now u n l o a d e d  shown  (7 ,r )  A  is  AB i s  r  being  level  is  (7 ,r ), t  (3"14)  element  diagram  resetting  nature.  unloading  element  effect  strength  the  strain  point  the  represented  loading  )  r  soil  reversal  reversal  initial  (3-13)  that  shear  stress-strain by  is  the  )  t  <7max " 7 r , ± T  f  to  stress  (3-9)  if  from,  Consider state  analysis,  curve  This  essentially  simplified  hysteretic  in  35  FIG. 3-2  HYPERBOLIC RELATIONSHIP UNDER GENERAL  LOADING  36  The of  shear modulus at any t i m e , i s the v a l u e of the  the  stress-strain  corresponding to calculated  by  d e s c r i b i n g the whichever  3.3  that  curve time.  evaluating initial  at The  the  or  the  stress-strain  tangent  shear  differential  after  tangent  reversal  of  point  modulus the  is  equation  hyperbolic  curve,  i s a p p r o p r i a t e at t h a t i n s t a n t i n t i m e .  F o r m u l a t i o n of s t r u c t u r e s t i f f n e s s , mass and damping matrices 3.3.1 An  linear  Partial Stiffness element  stiffness  Matrices matrix  i s , f o r a g i v e n geomtry, a  f u n c t i o n of the terms i n the  ( {6o}=  [D]{6e} ).  In  stress-strain  general  the D m a t r i x i s a f u l l  m a t r i x , w i t h 36 independent terms. stiffness  matrix  must  [D] m a t r i x t o be as w e l l . analysis  For s t a b i l i t y ,  be p o s i t i v e As  considers isotropic  terms(which  may  plane  strain  earlier,  the  s o i l b e h a v i o r under the  be  Thus  only  expressed  m a t e r i a l parameters) are r e l e v a n t i n isotropic  the  [D], 6 by 6  element  d e f i n i t e , t h i s r e q u i r e s the  mentioned  but p r a c t i c a l case of plane s t r a i n . independent  matrix  the  6  in  restricted of  terms  present  linear elastic relation  present  the of  study.  21 two The  between s t r e s s  and s t r a i n can be w r i t t e n as: {a}=  [D]{e}  For an increment {Aa}= where  (3-15) change t h i s can be w r i t t e n a s :  [D]{Ae}  '  (3-16)  37  For  {Ae}  the  {Aa}  i s the  the  written  incremental  plane  strain  incremental  strain  stress  condition  (A£  vector  the  ,A£  x  (Ao  [d]  ,A7 )  y  xy  ,Aa  x  ,Ar )  y  x y  matrix  is  commonly  as, ~\-v  fDl=  v  E (1+*)(1-2*)  v  0  \-v  0  0  where  vector  D=  For  Young's  the  modulus  simple  independent  ratio  programs  did by  is  written  commonly  a  do  multiplying  *=  et.  any  Poisson's  case  the  more  Young's  than  The  ratio.  where  a l .  constant.  f  [k] =  and  E,  (Duncan  not  1-2*  isotropic  modulus  Poisson's  0  (3-17)  there  modulus  1980),  is  only  with  a  constant  non-linear  analysis  merely  change  the  D  matrix  stiffness  matrix  of  an  is  as:  [B] [D][B]dA  (3-18)  T  element,  Likewise, load  step  the  r  displacement  Appendix  i t  can  be  i s obtained  slightly  more  Poisson*s  ratio  complex  obtained  directly  and  by  element  ~Area  [B]  one  for  the  isoparametric  I.  seen  that  merely case  Young's by  matrix  a  by is  the  new  stiffness  multiplying isotropic  modulus.  a  new of  the  at  constant.  conditions with  The  multiplication  by  matrix  [D] old  each A  varying  cannot [D],  be From  38  equation  (3-18),  regenerated  each  at  undertaken,  each  element load  stiffness  matrix  For  type  step.  the  must  be  of a n a l y s i s  t h i s would be a q u i t e e x p e n s i v e p r o c e s s .  E  v  and  are not independentterms of the [D] m a t r i x : we cannot w r i t e , E dp dE  [D]=  +  yfdp  Id?  We can, however, w r i t e , B 'dD dB  [D]=  ' + G 'dD dG  '  (3-19)  where G and B are the shear and bulk modulus,  dD dB  The  above  1  4/3  -2/3  1  -2/3  4/3  0  0  and,  (3-20),(3-21)  0  m a t r i c e s are commonly r e f e r e d to as the c o n s t i t u t i v e  p a t t e r n s f o r the G, B model. S i n c e the element s t i f f n e s s m a t r i x i s a l i n e a r f u n c t i o n the  terms  of  i n D, and s i n c e the terms i n D are a l i n e a r f u n c t i o n  of the independent moduli G & B, i t stiffness  matrix  follows  that  the  element  ( f o r a g i v e n geometry) i s a l i n e a r f u n c t i o n of  the independent moduli G & B.  Thus we can w r i t e ,  for  the  G-B  model, [k]= B [ S ] + R  Where G and increment developed  B  are  (3-22)  G[S ] G  tangent  values  at  the  beginning  and e v a l u a t e d a c c o r d i n g t o c o n s t i t u t i v e i n chap  3.2.  of  an  relationships  39  S  B  i s t h e element ' p a r t i a l s t i f f n e s s m a t r i x ' f o r B= 1 and G= 0,  and S  B  S  i s element p a r t i a l s t i f f n e s s m a t r i x f o r B= 0 and G= 1.  G  and S  G  a r e o b t a i n e d from:  f [B]  [S ] = R  •/Area  f [B]  [S Jr  •'Area  dD dB  [B]dA  dp dG  [B]dAr  For i s o p a r a m e t r i c q u a d r i l a t e r a l s ,  (3-23)  r  (3-24)  the i n t e g r a l s  24) have t o be e v a l u a t e d n u m e r i c a l l y .  (3-23)  and (3-  The procedure i s o u t l i n e d  i n Appendix I . The advantage of t h e p a r t i a l s t i f f n e s s approach i s t h a t t h e shear  and bulk  stiffness  e v a l u a t e d j u s t once. matrix  matrices  f o r each element need be  At each l o a d s t e p ,  an  element  stiffness  i s o b t a i n e d by f i r s t m u l t i p l y i n g t h e i r p a r t i a l  stiffness  m a t r i c e s by t h e i r r e s p e c t i v e modulus v a l u e s and adding  t h e two  m a t r i c e s , as i n d i c a t e d by e q u a t i o n It  should  be  noted  that  (3-22). f o r i s o p a r a m e t r i c elements t h e  s t r a i n s and t h e r e f o r e t h e s t r e s s e s v a r y consequently If  this  t h e moduli  was  taken  over  t h e element  and  a r e not c o n s t a n t throughout t h e element.  into  consideration  the p a r t i a l  stiffness  approach would not be v a l i d and t h e r e - g e n e r a t i o n of t h e element s t i f f n e s s m a t r i x a t each s t e p would be s t r e s s e s and s t r a i n s representative  of  required.  I n s t e a d the  a t t h e element c e n t r o i d were taken as being the values  i n the e n t i r e  element and the  40  moduli  a p p r o p r i a t e o n l y t o t h e c e n t r o i d was used t o d e s c r i b e the  s t r e s s - s t r a i n b e h a v i o r of t h e e n t i r e element. allows  partial  t h e i r present  3.3.2  assumption  element s t i f f n e s s m a t r i c e s t o be r e p r e s e n t e d i n  form, e q u a t i o n s  (3-23) and (3-24).  Mass M a t r i x  I t i s p o s s i b l e t o d e v e l o p a mass element  This  which  i s consistent  with  matrix the  'for each  adopted  displacement  i n t e r p o l a t i o n f u n c t i o n and mass d i s t r i b u t i o n w i t h i n the (Cooke  1975).  The  resulting  elements  to  l e a d s t o a d i a g o n a l mass diagonal the  and  A lumped mass a p p r o x i m a t i o n  i s assumed  mass  rotational  overestimated,  matrix  be  mass  coupling  , where t h e mass of the  and  no  coupling.  When  i s used, Lsymer (1979) has observed of  the  results  i n an  individual  matrix  elements  underestimation  h i g h e s t n a t u r a l f r e q u e n c i e s of t h e system. consistent  possesses  c o n c e n t r a t e d a t t h e nodal p o i n t s ,  matrix  inertia which  element  mass m a t r i x i s l i k e i t s element  s t i f f n e s s m a t r i x i n t h a t i t i s banded properties.  finite  On o t h e r  frequencies  m a t r i x may be s u i t a b l e .  of  For s i t e the  Good r e s u l t s have been shown by Penzien approximation.  I t has the  s a v i n g s i n computer s t o r a g e and computation  required f o r matrix  the  e i t h e r approach f o r d e r i v i n g a mass  (1969) when u s i n g the lumped mass advantages  is  l e a d s t o an o v e r e s t i m a t i o n of the same  response problems where t h e response i s governed g r e a t l y by natural  that  of the  hand,  order of magnitude of t h e h i g h e s t n a t u r a l f r e q u e n c i e s .  lower  the  calculations.  time  41  For the p r e s e n t a n a l y s i s , o n e - t h i r d of triangular  element  and  one-fourth  of  the  mass  of  each  the  mass  of  each  q u a d r i l a t e r a l element are lumped a t t h e i r r e s p e c t i v e nodes. mass of any one node i s the sum  of  the  contributions  s u r r o u n d i n g elements t o t h a t p a r t i c u l a r node.  The  of i t s  The m a t r i x [M] i s  as f o l l o w s :  [M] =  m.  0  0  0  0  0  —  0  0  m  2  0  0  0  0  -  0  0  0  m  3  0  0  0  -  0  0  0  0  m  0  0  -  0  0  0  0  0  m  5  0  -  0  0  0  0  0  0  m  -  0  0  6  (3-25)  0 0  0  0  0  0  0  -  m  where m i s the mass of the node a s s o c i a t e d w i t h the  i t h degree  of freedom, and n i s the t o t a l number of degrees of freedom.  3.3.3 When material internally  Damping M a t r i x vibrational medium, due  a to  energy portion  a  is of  being  t r a n s m i t t e d through a  i t s energy  number of mechanisms.  v i s c o u s type of damping.  T h i s d i s s i p a t i o n of  is  dissipated  One of t h e s e i s a energy  causes  a  decrease i n the a m p l i t u d e of v i b r a t i o n and can be b r o a d l y termed 'material  damping'.  In the e q u i v a l e n t l i n e a r a n a l y s e s by Seed  42  (1969) and Lsymer (1969) where n o n - l i n e a r h y s t e r e t i c m a t e r i a l i s represented has  by a l i n e a r v i s c o - e l a s t i c  t o be i n t r o d u c e d a r t i f i c i a l l y  v i s c o u s type damping. the  present  model,  material  through a frequency  where  a  hysteretic  r e l a t i o n s h i p i s used, a r t i f i c i a l v i s c o u s damping needed  to  model  dependent  In t r u l y n o n - l i n e a r a n a l y s i s ,  analysis,  the  different  types  of  damping  as  is  in  stress-strain is  no  longer  m a t e r i a l damping.  Damping i s i n t r o d u c e d i n t o the p r e s e n t a n a l y s i s t o d e s c r i b e actual  v i s c o u s e f f e c t s due  grains. degree  t o the presence of water i n the  S t u d i e s by F i n n e t . of  viscous  the  a l . , (1979) have shown  that  soil some  damping i s r e q u i r e d t o s t a b i l i z e systems a t  r e v e r s a l p o i n t s , where  there  are  abrupt  changes  in  modulus  Rayleigh  to  produce  values. Damping  expressions  orthogonality  in  introduced  modal  by  superposition  solutions  to  dynamic  s t r u c t u r e response, have been used i n e q u i v a l e n t l i n e a r by  Seed  (1960).  analyses  A R a y l e i g h type damping e x p r e s s i o n  the o r t h o g o n a l i t y  feature  important  present  i n the  present  analysis.  The  f o l l o w i n g r e l a t i o n s h i p i s used f o r each element; + b [k]  i n which [ c ] , [m]  [k]  matrices  respectively  for  e  matrix corrresponds and b are g i v e n  by:  (3-26)  e  and  e  effects  the  used  e  viscous  in  is  e  model  not  analysis)  [ c ] = a [m]  to  is  (although  e  are the damping, mass and  stiffness  element  stiffness  e.  The  t o i t s v a l u e at time= 0.  element The  parameters  a  43  The  a =  Xg'to,  (3-27)  b =  X/w,  (3-28)  value  X  of  e  expressed  i n percentage of c r i t i c a l  r e p r e s e n t s the damping r a t i o f o r element e. equal  to  evaluated  the by  fundamental  frequency  of  the a n a l y s i s at time t=0  damping  The parameter CJ, i s the  system  (or based on G  max  and  is  ).  The  complete damping m a t r i x [C] of the e n t i r e s t r u c t u r e i s o b t a i n e d from  the  individual  s t i f f n e s method. damping  matrices  element  Given the form  damping of  matrices  equation  by the d i r e c t  (3-26),  element  and t h e r e f o r e the s t r u c t u r e damping m a t r i x i s  assumed t o remain c o n s t a n t d u r i n g dynamic  loading.  44  CHAPTER 4 NUMERICAL ANALYSIS OF NON-LINEAR DYNAMIC RESPONSE 4.1  General In  t h i s c h a p t e r t h e n u m e r i c a l t e c h n i q u e used t o  solve  the  f i n i t e element m o d e l l i n g of the n o n - l i n e a r dynamic reponse of an earth  structure  t o earthquake  motions i s p r e s e n t e d .  p r a c t i c a l cases a s t a t e of p l a n e s t r a i n can be assumed for  I n many so  that  a n a l y s i s purposes t h r e e d i m e n s i o n a l s t r u c t u r e s can i d e a l i z e d  by a f i n i t e element system r e p r e s e n t i n g t h e c r o s s s e c t i o n of the structure.  Modelling  of t h e c r o s s s e c t i o n by a f i n i t e  element  system of t r i a n g l e s and q u a d r i l a t e r a l s r e q u i r e s , the f o r m u l a t i o n of the  c o n s t i t u t i v e s t r e s s s t r a i n b e h a v i o r of t h e s o i l m a t e r i a l , and d e r i v a t i o n of the s t i f f n e s s and  damping  properties.  This  has been o u t l i n e d i n c h a p t e r 3. In  earthquake response a n a l y s e s of l i n e a r s t r u c t u r e s , many  e a r l y i n v e s t i g a t o r s have (Wilson  (1962)  involves  the s o l u t i o n  represented  by  used  ,Clough  the  and of  free  the mode-superposition Penzien  (1975)).  the c h a r a c t e r i s t i c vibration  response  This value  of  method, method problem  t h e system,  f o l l o w e d by t h e t r a n s f o r m a t i o n of the d i s p l a c e m e n t s t o t h e mode shapes the  of t h e system.  T h i s procedure uncouples t h e response of  system, so t h a t t h e response of each mode may  independently  of  the others.  The  a n a l y s i s i s c a l l e d the step-by-step direct  numerical  their original  integration  form.  of  second  method,  be e v a l u a t e d  method of dynamic and  involves  the  the e q u i l i b r i u m equations i n  45  The main advantage of t h e mode s u p e r p o s i t i o n method i s t h a t the response  of a system may be o b t a i n e d w i t h good  considering  only  step-by-step retained,  by  a few of t h e lower normal modes, w h i l e i n t h e  method  Penzien  accuracy  a l l generalized  (1969).  coordinates  must  be  On t h e o t h e r hand, the e v a l u a t i o n of  the c h a r a c t e r i s t i c v a l u e problem and t r a n s f o r m a t i o n t o the mode shapes  a r e major  computational  s t e p - b y - s t e p method.  Recent  1969)  equivalent  have  used  an  method t o approximate response  solved  problems  not r e q u i r e d i n t h e  investigators  (Seed  and  linear viscoelastic  the non-linear  i n t h e frequency  solution. domain  Idriss  iterative  The  dynamic  i s based  on the  assumption of l i n e a r s t r u c t u r a l b e h a v i o r , cannot be used  for a  n o n - l i n e a r system. The be  s t e p - b y - s t e p i n t e g r a t i o n method, on t h e other hand can  applied  response  to non-linear  systems.  In  this  approach,  i s c a l c u l a t e d f o r a s h o r t time increment  w i t h known c o n d i t i o n s , t o e v a l u a t e t h e c o n d i t i o n s time. by  At, s t a r t i n g at  a  later  The i n c r e m e n t a l l i n e a r n a t u r e of t h e system i s c o n s i d e r e d  assuming  step,  the  and  properties  linear by  behavior  making  proper  p r i o r t o each s t e p .  throughout  each s u c c e s s i v e time  modifications  to  the  linear  I n t h e p r e s e n t s t u d y , where the  s t i f f n e s s p r o p e r t i e s of the s t r u c t u r e behave n o n - 1 i n e a r l y , being s t r a i n and s t r e s s dependent, t h e s t e p - b y - s t e p i n t e g r a t i o n method i s used. The complete displacement,  response  velocity  i s obtained  by  using  t h e known  and a c c e l e r a t i o n a t t h e end of one time  46  i n t e r v a l as the i n i t i a l c o n d i t i o n s f o r the next p r o c e s s i s c o n t i n u e d s t e p - b y - s t e p from i n i t i a l to  the  c o m p l e t i o n of s e i s m i c motions.  interval.  The  static conditions  The dynamic e q u i l i b r i u m  c o n d i t i o n i s s a t i s f i e d a t t h e b e g i n n i n g and  end  of  each  time  interval. The  equations  together with  their  of  motions  solution  by  for non-linear  structures,  a  integration  step-by-step  procedure w i l l be p r e s e n t e d i n the f o l l o w i n g  4.2  sections.  E q u a t i o n of motion As  mentioned  ground motion structure.  p r e v i o u s l y i t i s assumed t h a t the earthquake  i s i d e n t i c a l a t a l l p o i n t s a l o n g t h e base Spatial  variations  in  the  of the  ground motion a r e not  c o n s i d e r e d i n the p r e s e n t a n a l y s i s . The dynamic e q u a t i o n of motion of nodal the  rigid  earthquake  points  ground motion can be e x p r e s s e d i n t h e m a t r i x form:  of  the  derivation  Newmark and R o s e n b l u e t h  (4-1 )  i s g i v e n i n Z i e n k i e w i c z (1971) and  (1971).  The  matrices  are  follows: [M]  the  base f o r the f i n i t e element system when s u b j e c t e d t o  [M] {x} + [C] {x} + [K] {x.} = {R}  Details  above  = the mass m a t r i x ( s e c t i o n 3.3.2)  [C]  the damping m a t r i x ( s e c t i o n 3.3.2)  [K]  the s t i f f n e s s m a t r i x ( s e c t i o n 3.3.2)  defined  as  47  {x}  = v e c t o r of nodal d i s p l a c e m e n t s  {R}  = i n e r t i a force vector  r e l a t i v e t o base  For t h e lumped mass system the i n e r t i a f o r c e v e c t o r {R} i s (R)  = -{M^x  b  - {M\%  (4-2)  where, {M}  x  =(m,0m 0---m  0  - - m 0 )  (4-3)  = (0 m, 0 m  m  - - 0 m )  (4-4)  2  and {M}  Y  2  - - - 0  {M} and {M} a r e the v e c t o r s of the mass of t h e nodes  associated  w i t h the x and y degrees of freedom r e s p e c t i v e l y (m  i s the mass  x  of  Y  t h e node a s s o c i a t e d w i t h t h e i t h degree of freedom which i s  in the x d i r e c t i o n ) . mass  matrix  Therefore,  [M], x  D  and %  {M} + {M} = d i a g o n a l x  Y  of t h e  a r e r e s p e c t i v e l y t h e h o r i z o n t a l and  v e r t i c a l components of t h e ground a c c e l e r a t i o n . The s t i f f n e s s m a t r i x f o r each time  interval  of  the step  by  procedure d e s c r i b e d i n c h a p t e r complete  stiffness  matrix  finite step  3.3  and  [K] o f  element  during  method i s o b t a i n by t h e i n appendix  the e n t i r e  stiffness  method.  damping m a t r i x i s o b t a i n e d .  In a  I.  The  structure i s  o b t a i n e d from t h e i n d i v i d u a l element s t i f f n e s s m a t r i c e s direct  any  by t h e  s i m i l a r manner t h e s t r u c t u r e  These m a t r i c e s a r e of order N x  N,  48  where N i s the number of degrees of freedom. are  banded  and  symmetric,  o n l y m a t r i c e s of s i z e N x M need be  c o n s i d e r e d , where M i s the h a l f  bandwidth  T h i s g r e a t l y reduces the computer s t o r a g e  4.3  I n c r e m e n t a l e q u a t i o n of Equation  S i n c e the m a t r i c e s  of  the  structure.  requirements.  motion  (4-1) must be s a t i s f e d at every i n s t a n t i n t i m e .  L e t T = t + A t , where At i s a s m a l l time  interval  then, [ M ] { x } + [ C ] {x} + [ K ] {x} = {R} t  t  t  t  t  [M] {x} + [C] {x} + T  The  T  T  subscript  [K] {x} =  T  T  refers  to  t  (4-5)  {R}  (4-6)  t  T  the  T  instant  of  p a r t i c u l a r q u a n t i t y t a k e s on i t s v a l u e .  time  S i n c e the  a t which mass  the  matrix  i s constant matrix, [M]  = [M]  T  t  =  [M]  so t h a t , [M] {x} T  - [Mj {x}  T  t  t  = [M] ({x} -{x} T  T  )= [M]{Ax}  t  (4-7)  where, ({x} and  T  - {x} )=  {Ax}  t  (4-8)  T  the damping m a t r i x i s assumed t o remain c o n s t a n t  throughout  the a n a l y s i s , [C] {x} T  where,  T  - [C] {x} t  t  = [ C ] ( { x } - { x } )= [C] {Ax} T  T  t  T  (4-9)  49  ({x}  - {x} )= {Ax}  T  t  (4-10)  T  However, the [K] m a t r i x depends on the l e v e l of shear s t r a i n the  s o i l , as a r e s u l t of the n o n - l i n e a r s t r e s s - s t r a i n  and thus e q u a t i o n s obtained  for  similar  this  to  (4-7)  matrix in a straight  a r e two ways t o a r r i v e a t e x p r e s s i o n s approximation  and  methods.  (4-9)  relation,  can  not  forward manner.  similar  to  in  be  There  (4-9)  using  A crude method i s t o r e p l a c e [ K ]  T  with  [K] , so t h a t , [K] {x} T  - [K] {x}  T  t  = [K] ({x} -{x} )=  t  t  T  t  [K] {Ax}  (4-11)  f  where, {x} )=  Ux} -  t  T  The  tangent  {Ax}  T  s t i f f n e s s p r o p e r t i e s d e f i n e d at the b e g i n n i n g of the  time i n t e r v a l are used. i s changing [K]  D u r i n g any  time i n t e r v a l the [k] m a t r i x  as the s t i f f n e s s p r o p e r t i e s change w i t h s t r a i n .  matrix  should  r e p r e s e n t some average s t i f f n e s s p r o p e r t i e s  d u r i n g the time i n t e r v a l .  T h i s can be o n l y  done  by  iteration  because the d i s p l a c e m e n t  at the end of the time increment  on  As  these  properties.  iteration  is  not  always  will  be  desired.  (4-5) from e q u a t i o n  ( 4 - 7 ) , (4-9) and equation  of  which i s as  (4-11),  For  this  motion follows:  for  reason  the  incremental  4.3, the  Subtracting  (4-6) and t a k i n g note of  gives  depend  e x p l a i n e d i n chapter  approximate method i s used i n the p r e s e n t a n a l y s i s . equation  The  form  equations of  the  the time i n t e r v a l s t a r t i n g at time t ,  50  [M] {Ax} + [C] {Ax} + [ K ] { A x } =  {AR}  Matrix  is a  T  T  equation  t  (4-12)  T  which  (4-12)  T  s e t of second  d i f f e r e n t i a l e q u a t i o n s , may be reduced t o a r e c u r r e n c e if  an  assumption  i s made  regarding  the  equation  variation  of  a c c e l e r a t i o n of each node w i t h i n the time i n t e r v a l A t . end,  numerical  applied  to  {x} ,{x} ,{x} T  T  procedures  the  equation  at  T  developed  time  so T,  order  To  the this  by Newmark (1959) can be  that  unknown  displacements  can be expressed i n terms of known  d i s p l a c e m e n t s {x} ,{x} , { x } a t time t . t  4.4  t  t  Step by Step I n t e g r a t i o n In Newmark's 0 method  parameters,  a  and  displacement at acceleration,  0  time  of  are  T  step used  can  velocity  be  and  by so  step that  the  expressed  displacement  unknown a c c e l e r a t i o n a t time T.  integration  in at  velocity terms  time  of  two and the  t , and of  The e q u a t i o n s f o r v e l o c i t y  and  d i s p l a c e m e n t a t time T a r e as f o l l o w s :  (x}  T  = {x}  t  + (1-a)At{x}  {x}  T  = {x}  t  + At{x}  + aAt{x}  t  (4-13)  T  + (.5-/3) ( A t ) { x } 2  t  + 0(At) {x} 2  t  T  (4-14)  At i s r e f e r e d t o as the 'time s t e p ' of i n the i n t e g r a t i o n . Newmark integration  (1959) procedure  proposed that  a  for a  unconditionally  =  and  1/2  0  =  1/4.  stable This  51  corresponds  to  integration.  When v a l u e s of a = 1/2 and /3 =  linear  a  variation  increment.  constant  of  average  acceleration  The l i n e a r  variation  acceleration  is of  proposed by W i l s o n and Clough (1962).  1/6  assumed  method  are  used,  over  acceleration  of a  the time  method  was  Both approaches have been  i n c o r p o r a t e d i n t o the computer program and a r e o p t i o n s a v a i l a b l e t o the u s e r . The  l i n e a r v a r i a t i o n of a c c e l e r a t i o n method may  instability dependent  of on  the the  solution. size  of  the  This time  instability step  At  l e a d t o an is  used  i n t e g r a t i o n , m a t e r i a l p r o p e r t i e s , and s i z e of the f i n i t e grid.  The W i l s o n 8  solution  method  which  provides  ( W i l s o n et a l , 1973) has been w r i t t e n  stability  usually in  the  element in  the  i n t o the computer  program. I f the f o l l o w i n g s i m p l i f y i n g e x p r e s s i o n s a r e used, {a}  t  = At{x}  t  {b}»  = At{x}  t  then e q u a t i o n  T  Equations 12).  T  +.5(At)Mx}  + aAt {Ax }  t  = {b}  Transferring  . (4-17)  T  + 0(At) {Ax}  T  (4-18)  and  T  (4-18) are s u b s t i t u t e d  i n t o e q u a t i o n (4-  a l l forms a s s o c i a t e d w i t h known v a l u e s t o the  r i g h t hand s i d e l e a d s t o the f o l l o w i n g  T  can w r i t t e n a s :  2  t  (4-17)  [D] {Ax}  (4-16)  t  (4-13) and e q u a t i o n ( 4 - 1 4 )  {Ax} = {a} {Ax}  (4-15)  = [P]  T  equation, (4-19)  52  where, [D]  T  = [M] + aAt[C] + / 3 ( A t ) [ K ]  [P]  T  = ( R } ~ (R), - [C] ( a } - [ K j { b }  (4-20)  2  and  T  Equation  t  (4-19)  (4-21)  t  a  static  increment  e q u i l i b r i u m r e l a t i o n s h i p and may be s o l v e d by  matrix  inversion  and m u l t i p l i c a t i o n . {Ax}  with  incremental  (4-22)  T  displacements  o b t a i n e d from e q u a t i o n s This  {Ax}  T  time  step  velocities  {Ax}  to  vary  i n some  by the v a l u e s of t h e a and 0 ) ,  and (2) t h e s t i f f n e s s p r o p e r t i e s of t h e  s t r u c t u r e remain c o n s t a n t d u r i n g any time s t e p and a r e e q u a l its  values  at  t h e b e g i n n i n g of t h e time s t e p .  Because t h e constant  when  the  structure  time  damping  f o r the d u r a t i o n  of  as  good  s t e p At i s chosen t o be s m a l l . matrix  i s assumed  to  remain  the a n a l y s i s , inherent e r r o r s  a r i s i n g from assumption (2) a r e not p r e s e n t i n t h e modeling damping p r o p e r t i e s .  to  N e i t h e r of the  two assumptions a r e s t r i c t l y c o r r e c t but can be viewed approximations  T  i n c l u d e s two s i g n i f i c a n t  (1) t h e a c c e l e r a t i o n i s assumed  d e s c r i b e d manner, (as determined the  and  (4-17) and (4-18).  n u m e r i c a l a n a l y s i s procedure  assumptions:  during  to  Consequently,  = [D)fMPl  T  i s equivalent  of  53  As soil  d i s c u s s e d i n Chapter 3.3, t h e s t i f f n e s s p r o p e r t y of the  i s determined by t h e tangent  with  the  level  the tangent  moduli  s t r a i n and mean normal s t r e s s , and by  b u l k modulus which i s v a r i e d w i t h  stress only. calculated  of shear  shear modulus which i s v a r i e d  t h e mean  normal  Element shear s t r a i n and mean normal s t r e s s v a l u e s at  the end  according  to  of  their  t h e time s t e p a r e used t o compute constitutive  relations.  The  [K]  m a t r i x f o r t h e next time s t e p i s based on these newly c a l c u l a t e d moduli.  I n t h i s way t h e i n e l a s t i c b e h a v i o r of s o i l  by an i n c r e m e n t a l l i n e a r  approach.  In n o n - l i n e a r problems t r u e convergence seldom that  i s modelled  occurs,  in  t h e i n c r e m e n t a l f o r c e s a p p l i e d a r e not e q u i l i b r a t e d by the  incremental s t r e s s e s .  I n order t o r e s o l v e t h i s i n c o n s i s t e n c y an  error c o r r e c t i o n i s incorporated  into  the a n a l y s i s  This  is  d i s c u s s e d i n t h e next s e c t i o n .  4. 5  Error Correction A  tangent  s t i f f n e s s method has been used t o f o r m u l a t e the  incremental f i n i t e  element m a t r i x e q u a t i o n  method t h e r e a r e b a s i c a l l y  1.  (4-12).  With  this  t h r e e approaches t h a t can be t a k e n .  Accept t h e s t r e s s e s : Accept  {6c}=  [D ]{6e}, L  where [ D ] i s t h e tangent L  stress-  s t r a i n m a t r i x used f o r the l a s t i t e r a t i o n , see e q u a t i o n (320).  C a l c u l a t e the new [ D ]  attained.  L  In  this  approach  based  on  the  equilibrium  stress is  level  maintained  54  throughout the a n a l y s i s , c u r v e i s not f o l l o w e d .  however  the  true  stress-strain  The r e s u l t i n g i n c o n s i s t e n c y between  the computed s t r a i n and exact s t r a i n s i s a c c e p t e d . 2.  Accept the s t r a i n s : Accept  strains  computed from n o d a l d i s p l a c e m e n t s o b t a i n e d  by the s t e p by s t e p i n t e g r a t i o n procedure and c a l c u l a t e the new  tangent  Equilibrium is 3.  modulus is  based  the  the  strain  level.  v i o l a t e d , but the t r u e s t r e s s - s t r a i n curve  followed.  Accept the s t r a i n s and a p p l y c o r r e c t i o n  forces:  The t r u e c u r v e i s f o l l o w e d , and any e r r o r i n e q u i l i b r i u m at the end of the time s t e p i s c a l c u l a t e d and a p p l i e d b e g i n n i n g of the next time s t e p .  at  the  In t h i s way any e r r o r s i n  e q u i l i b r i u m do not accumulate. A l l the above methods can b e n e f i t from m u l t i p l e For  iterations.  example, from the i n i t i a l tangent moduli used f o r the  i n t e g r a t i o n of the time s t e p , the tangent moduli a t the the  time  s t e p can be o b t a i n e d .  end  be d e s c r i b e d  as  a  of  I t e r a t i o n i s then performed  the time s t e p u s i n g the average of the two tangent m o d u l i . approach may  first  step-secant  approach.  on  This While  i t e r a t i n g i s d e s i r a b l e i n terms of i m p r o v i n g the a c c u r a c y of the incremental  approach,  the  number  of a d d i t i o n a l i t e r a t i o n s at  each time s t e p would p r o p o r t i o n a l l y i n c r e a s e analysis. elements  T h i s may is  the  cost  used.)  the  prove t o be e x p e n s i v e when a l a r g e number of  involved.  ( A l t h o u g h t h e r e i s a t r a d e o f f here i n  t h a t one can g e n e r a l l y use l a r g e r time increments i f are  of  iterations  For t h i s r e a s o n , i t e r a t i o n i s not c o n s i d e r e d i n the  55  present a n a l y s i s .  Approach (3) has  the advantage  of  following  the a c t u a l s t r e s s - s t r a i n c u r v e w h i l e at the same time s a t i s f y i n g equilibrium.  In  this  sense, and  when m u l t i p l e i t e r a t i o n s are  not performed, approach (3) i s the most a c c u r a t e approach  has  in greater The the  been b u i l t  i n t o the a n a l y s i s and  method.  This  i s discussed  here  detail.  s t r a i n s determined i n each i n t e g r a t i o n are a c c e p t e d  true  strains  determine the restoring  the s t r e s s - s t r a i n r e l a t i o n s are used to  'true r e s t o r i n g ' f o r c e  forces  equation (4-3). defined  and  do  as  not  [K]  necessarily  An set of a r t i f i c i a l  {x} .  T  However  T  satisfy  the  these  equlibrium  ' e x t e r n a l ' f o r c e s , {P  by the f o l l o w i n g can be a p p l i e d t o the  system  },  so  that  equilibrium i s restored, {P  The  {P  }=  err  {R}  - [M]{x}  T  ) i s added t o {P}  err  analysis.  All  terms  s t r a i g h t forward manner. i s evaluated  T  in  T  - [C] { x } T  [K] {x} T  (4-23)  T  f o r the next time increment  of  the  the above e q u a t i o n are computed i n a The  deserves comment.  manner i n which the term [ K ] A l l other  terms  in  {X}  T  the  T  above  e q u a t i o n are computed i n a s t r a i g h t forward manner. A b a s i c assumption of the f i n i t e element method i s t h a t a l l forces  are t r a n s m i t t e d  at the n o d a l p o i n t s .  any  node i s the sum  of the c o n t r i b u t i o n s of i t s  one  elements t o t h a t p a r t i c u l a r node. then,  are  The  nodal f o r c e s representing  the e l e m e n t s .  The  s t a t e of s t r e s s  The  t o t a l f o r c e at surrounding  restoring forces the sum  within  an  [K]  y  of the s t r e s s e s element  can  {x}  T  of be  56  represented forces.  by  a  statically  equivalent  system of nodal p o i n t  From t h e s t r e s s e s {a}, nodal f o r c e s a r e o b t a i n e d by, (4-24)  As  was  the d e c i s i o n e a r l i e r i n e v a l u a t i n g t h e s t i f f n e s s  f o r an element, Chapter 3.3., t h e s t r e s s e s and centroid  of  a t the  the element a r e taken t o be r e p r e s e n t a t i v e of the  average s t r e s s e s and s t r a i n s i n t h e element. stress  strains  matrix  vector  I n t h i s case,  the  can be removed from w i t h i n t h e i n t e g r a l s i g n , as  shown above. The t o t a l contributions node.  force of  at  one  node  i t s surrounding  The r e s t o r i n g f o r c e s  representing  any  the  sum  of  T  then,  the. stresses  manner i n which element s t r e s s e s  sum  elements t o t h a t  [K] {x} T  i s the  of the  particular  a r e nodal  forces  of t h e elements.  are c a l c u l a t e d  is  The  described  below. From t h e a c c e p t e d s t r a i n s the c o r r e s p o n d i n g s t r e s s e s can be d e t e r m i n e d as f o l l o w s : 1.  Calculate incremental {Ae}  ;  strains:  = [B; ]{Au}i  (4-25)  where, {Au}| i s t h e v e c t o r of i n c r e m e n t a l  nodal  displacements  a s s o c i a t e d w i t h element i . [B i ]  i s t h e v a l u e of the s t r a i n d i s p l a c e m e n t  element i , a t i t s c e n t r o i d , see Appendix I .  matrix  57  {Ae}jthe  i n c r e m e n t a l s t r a i n v e c t o r (A£  for  ,&Z ,&y ) Y  Y  element i . 2.  C a l c u l a t e increment s t r e s s e s : {Aa},=  [D ]{Ae}j  (4-26)  L  where, [D ]  i s the s t r e s s - s t r a i n  L  shear  and • bulk  moduli,  matrix, of  based  element  i n t e r a t i o n time s t e p , see e q u a t i o n  on  tangent  i , d u r i n g the  (3-20).  {Ao}j i s the i n c r e m e n t a l s t r e s s v e c t o r (Aa  ,Aa  x  y  ,Ar  xy  )  f o r element i . 3.  C a l c u l a t e t o t a l element s t r e s s e s : For  any  element, s t r e s s e s at time T are equal t o the  summation of the i n c r e m e n t a l s t r e s s e s up t o They are o b t a i n e d as (a}  Ti  = {a}  ti  time  T.  follows:  + {Aa}  (4-27)  Ti  where, T = t + At {a}  i s the t o t a l s t r e s s (o , o , r ) , a t time t  f  x  Y  xy  { A a } i s the i n c r e m e n t a l s t r e s s f o r the l a s t time T  4.6  Procedure  for strain reversal  A complete d e s c r i p t i o n of under 3.2.2.  general  loading  the  occurrence stress-strain  strain  relationship  and u n l o a d i n g has been g i v e n i n Chapter  On u n l o a d i n g or r e l o a d i n g of a s o i l element,  stress-shear  step  behavior  the  shear  of the m a t e r i a l i s governed by a  58  newly d e f i n e d  curve,  equation  (3-14),  where  i t s origin  l o c a t e d a t the shear s t r e s s shear s t r a i n r e v e r s a l p o i n t , ( 7 , see  Figure  3-2.  For  two  reasons  a n a l y s i s can determine when s t r a i n reversal  a  new  it  shear  reversal  modulus  at  ,T ), R  i s important t h a t the occurs:  s t r e s s - s t r a i n branch i s f o l l o w e d , and  i s a d i s c o n t i n u i t y of  is  any  strain  (1)  upon  (2) t h e r e reversal  location. Shear for  strain,  as d e f i n e d i n Chap 3.2.2  i s re-written  here  convenience,  7  (4-28)  m a x  The  rate  of shear s t r a i n , d i f f e r e n t i a t i n g e q u a t i o n (4-29)  w i t h r e s p e c t t o t i m e , i s d e f i n e d by,  7 ,  (  max  7  X  y )  ( 7 x y )  +  U  x  ~  C y H e x  ~  (4-30)  f y )  m a x  Based on d i s p l a c e m e n t and v e l o c i t y v a l u e s {x} and element  {x},  s t r a i n s and shear s t r a i n r a t e s are c a l c u l a t e d  by, {e}; = [ B j ] {u}j  (4-31)  where, {u}| i s the v e c t o r of n o d a l d i s p l a c e m e n t s a s s o c i a t e d w i t h element i . i s the s t r a i n v e c t o r ( £ , i - y , 7 ) f o r element i . x  x y  and, (4-32) where,  59  {u}j i s t h e v e c t o r of n o d a l v e l o c i t i e s a s s o c i a t e d  with  element i . {e}j  =  i s the  rate  of s t r a i n v e c t o r ( £ , J-y ,7 ) f o r x  xy  element i . S t r a i n r e v e r s a l can be d e f i n e d as a decrease recent  |7  m a x  | v a l u e , see F i g u r e 3-2.  A l o c a l maximum or minimum  i n t h e v a l u e of shear s t r a i n o c c u r s when 7 Therefore,  from the most  i s equal t o zero.  m a x  s t r a i n r e v e r s a l o c c u r s i n an element when t h e r e i s a  change i n s i g n i n the r a t e of shear s t r a i n 7  from  m Q X  one  time  s t e p t o t h e next time s t e p . L e t T = t + At S t r a i n r e v e r s a l o c c u r s when, 7 'max,  -7 * mri  < 0  (4-33)  'maxj  where, t h e s u b s c r i p t r e f e r s t o t h e i n s t a n t particular Now  that  a  criterion  point  for strain  has  been  considered.  reduce  the  time  interval  Since  a  At  of the  w i t h subsequent i t e r a t i o n when they o c c u r , i n order t o 'points'.  time s t e p of t h e a n a l y s i s s h o u l d be reduced  order t o l a n d r i g h t on them. Salgado  the  reversal  o b t a i n a b e t t e r a c c u r a c y when t u r n i n g these the  which  i s c h a r a c t e r i z e d by a d i s c o n t i n u i t y i n modulus,  i t may be d e s i r a b l e t o analysis  at  q u a n t i t y t a k e s on i t s v a l u e .  e s t a b l i s h e d , one f u r t h e r p o i n t s h o u l d be reversal  time  (1980)  on  a  one  With degree  system, t h i s was e a s i l y a c h i e v e d .  earlier  Ideally  sufficiently in  investigations  by  of freedom, e l a s t i c - p l a s t i c With m u l t i - d e g r e e of  freedom  60  systems,  where  more  than  one  strain  elements or l a y e r s ) may occur w i t h i n a necessarily  at  precisely  reversal ( i n different time  t h e same  time,  r e v e r s a l i n one p a r t of t h e system and may time  interval  cause  strain  interval  reversal  and later  reasonable  reversal  points  approach would  which  be  where  strain  i n the same  i n another p a r t of the  system, c o s t c o n s i d e r a t i o n s make t h i s approach A  and not  does  impractical.  require  landing  on  t o make a p r e d i c t i o n r e g a r d i n g the  segments of t h e time i n t e r v a l f o r which s t r a i n r e v e r s a l has and has  not o c c u r r e d .  W i t h i n the time i n t e r v a l of s t r a i n  t h i s can be approximated acceleration.  by  assuming  linear  variation  in  The time of r e v e r s a l can be e x p r e s s e d a s ,  At,  •  (TW >  'At  T  ' Tmaxj  and,  a  reversal,  (4-34)  — ^maxt  A t = At - A t ,  (4-35)  2  where, At,  i s t h e time segment ' b e f o r e ' s t r a i n  reversal  At  i s t h e time segment ' a f t e r '  reversal  2  strain  T h i s a p p r o x i m a t i o n can be used on s t r a i n r e v e r s a l with  iteration  performed  f o r t h e time  interval  c o r r e s p o n d i n g t o the s t r a i n r e v e r s e d element w i t h At,  time i n t e r v a l .  the  element,  and  the  r e m a i n i n g p o r t i o n of t h e time check  procedure  segment shortest  A f t e r the i n t e g r a t i o n f o r t h i s time segment  has been completed, s t r a i n r e v e r s e d modulus i s a s s i g n e d appropriate  elements,  and  i n t e g r a t i o n i s performed on the  interval.  subsequent  t o the  The  strain  segmenting of A t  2  reversal  further  into  61  two s m a l l e r segments i f r e v e r s a l o c c u r s , i s performed This  type  of  approach  dimensional s o i l layer halfbandwidth=  2)  has  been  systems  cost  (max  efficient  degree  of  on A t . 2  when  one  freedom=  20,  have been a n a l y s e d , w i t h t h e l i m i t i n g of the  s i z e of A t , and A t t o At/10, 2  e x o r b i t a n t l y expensive  (Lee, 1977).  for t y p i c a l l y  T h i s approach can be  large  degree  of  freedom  systems t h a t r e s u l t from f i n i t e element m o d e l l i n g . A  crude  approximation  e s t i m a t i n g , on s t r a i n closer  to  determined  to  this  reversal(s),  the beginning from e q u a t i o n  approach  whether  can be made by  reversal  occurred  o r t h e end of the time i n t e r v a l .  As  (4-34),  i f , A t , < .5«At, assume s t r a i n r e v e r s a l a t t i f , A t , > .5«At, assume s t r a i n r e v e r s a l a t T  I f A t , < ,5-At f o r any s t r a i n r e v e r s e d elements, moduli  reversed  a r e a s s i g n e d t o these elements and the complete time  At i s r e p e a t e d . assumed  Any a d d i t i o n a l  t o have  reversed  closer  to  t h e end  of  strain  reversed  When s t r a i n  occcurs  the time i n t e r v a l , t h e s t r a i n  reversed  I t e r a t i o n i n t h i s case  crude approach m u l t i p l e at  time  any  time  time this  step are  and  procedure  t o handle the o c c u r r e n c e s of s t r a i n r e v e r s a l has analysis.  step  of  With  avoided,  adopted i n t h e p r e s e n t  a  t h e next  i s not n e c e s s a r y .  subiterations  most,  elements a r e  reversal  moduli a r e not a s s i g n e d t o these elements u n t i l step.  step  a t t h e end of t h e time s t e p and a r e  a s s i g n e d the a p p r o p r i a t e m o d u l i .  step  strain  i s r e p e a t e d once.  This been  62  4.7  Summary of A  brief  procedure outline  of  solution  scheme  increment  i s g i v e n as  1.  Based on the c u r r e n t v a l u e s of  f o r any g i v e n time  follows: 7  ,  T  , a t time t ,  the tangent shear modulus G and tangent bulk modulus B are c a l c u l a t e d as d e s c r i b e d i n c h a p t e r 3.2.1  and  3.2.2  respect i v e l y . 2.  The m a t r i x [K] i n e q u a t i o n (4-20) and updated,  3.  and {b}  substitute  i n t o e q u a t i o n (4-21) t o set up [P] and  t  then  3.3.  {a}  t  e q u a t i o n s (4-15) and  (4-16), set  [D]  S o l v e f o r { A x } u s i n g e q u a t i o n ( 4 - 2 2 ) , and  subsequently  {Ax} .and {Ax}  (4-18)  T  T  5.  is  Calculate  up 4.  see c h a p t e r  (4-21)  from e q u a t i o n s (4-17) and  T  Acceleration,  v e l o c i t y and d i s p l a c e m e n t s of the nodes  at the end of the increment are o b t a i n e d from:  6.  {x}  T  = {x}  t  +  {Ax}  T  {x}  T  = {x}  t  + {Ax}  T  {x}  T  = {x}, + {Ax}  T  Based on d i s p l a c e m e n t and v e l o c i t y v a l u e s at of  increment  {x}  T  and  the  {x} , element s t r a i n s and T  s t r a i n r a t e s are c a l c u l a t e d  end shear  by,  {e} = [B|]{u}  (4-31)  {e} = [B;]{u}  (4-32)  and,  7.  Check t o see  if  there  are  S t r a i n r e v e r a l i s determined  any by,  strain  reversals.  63  a change i n s i g n o f :  If  there  beginning  are of  s t r a i n r e v e r s a l s and they occur a t the the  time  step,  the  integration  is  r e p e a t e d f o r t h e time s t e p , by r e p e a t i n g s t e p s 1 t o 8. 8.  Calculate  stresses  according  to  the  appropriate  s t r e s s - s t r a i n law. 9.  In order t o b r i n g the s t r e s s - s t r a i n the  actual  stess-strain  ' e x t e r n a l ' f o r c e , {P {Perr  err  curve,  closer  an  T  (4-23),  - [C] {X } - [ K ] { x } T  to  artificial  }, d e f i n e d by e q u a t i o n  }= { R ) " [M]{x} T  point  T  T  i s c a l c u l a t e d and added t o {P} i n s t e p 3 f o r the  next  increment. When increment the  next  step  9 has been completed,  i s finished. time  the a n a l y s i s f o r one time  The e n t i r e p r o c e s s may be  interval.  The  process  c o n s e c u t i v e l y f o r any d e s i r e d number of the  complete  response  can time  be  repeated  for  carried  out  increments;  thus  h i s t o r y can e v a l u a t e d f o r the n o n - l i n e a r  system. A computer program has been w r i t t e n based on t h i s  solution  scheme and the o v e r a l l f l o w c h a r t i s shown i n f i g u r e 4-1.  64  Read i n Data  Initialize Values  S t a r t of do l o o p  Assign BULK and SHEAR Modulus f o r each Element at each time s t e p B u i l d Master STIFFNESS and DAMPING M a t r i c e s [M] & [C] at each time s t e p  redo time s t e p  Set up Dynamic eqns [D] & {R}  Solve f o r DISP.,VEL.,ACC.  C a l c u l a t e STRAINS at t h e c e n t r o i d of each Element  Are t h e r e any s t r a i n Reversals  Yes  65  No Calculate Stresses according to S t r e s s - S t r a i n Law  S t o r e S t r e s s and Strains i f required  More Increments No STOP  FIG 4-1  PROGRAM  FLOWCHART  Yes  66  CHAPTER 5 APPLICATIONS OF THE NON-LINEAR FINITE ELEMENT METHOD 5.1  General The  non-linear  earlier  chapters  systems.  dynamic  was  response  applied  to  a  analysis  developed  number of s o i l  The same s t r u c t u r e s were a l s o a n a l y s e d  in  structure  using  Newmark  type methods t o e s t i m a t e movements so t h a t a comparison  c o u l d be  made on the p r e d i c t i o n of permanent d i s p l a c e m e n t s . I f comparison should  be  based  characteristics. soil  of r e s u l t s are t o be m e a n i n g f u l , the a n a l y s e s on  similar  In the Newmark  stiffness  analysis,  (Chapter  i s assumed t o behave i n a r i g i d p l a s t i c manner.  match  damping 2.2) To  the  closely  the r i g i d p l a s t i c b e h a v i o r , the n o n - l i n e a r f i n i t e element  method employs a l i n e a r e l a s t i c - p l a s t i c strain  law.  This  is  done  by  the  form  approaches.  of  plastic  shear  setting  h y p e r b o l i c r e l a t i o n s h i p e q u a l t o 0.01 in  and  the  or l e s s .  deformation  is  stress R  f  -  shear  v a l u e i n the  M a t e r i a l damping  inherent  in  both  In the f i n i t e element a n a l y s i s v i s c o u s damping must  be i n c l u d e d f o r n u m e r i c a l r e a s o n s .  A nominal but s m a l l v a l u e of  2  p e r c e n t of c r i t i c a l damping i s used t o s t a b i l i z e the s o l u t i o n  at  t u r n i n g p o i n t s of  procedure  for  stress-strain  done  in  comparison  and  this  structures i s evaluated While  The  Makdisi-Seed  p r e d i c t i n g d i s p l a c e m e n t s i s a l s o performed.  matching of s t i f f n e s s strictly  curves.  may  damping  characteristics  cannot  The be  c a s e , as the dynamic response of e a r t h using not  be  an  equivalent  linear  method.  as v a l i d , M a k d i s i and Seed have  67  shown good agreement w i t h the r i g i d t h i s , t h e i r approach i s i n c l u d e d A  hypothetical  i n the s t u d y . 1.  The  a h e i g h t of 150  acceleration  value  s t r u c t u r e s , o b t a i n e d from p s e u d o - s t a t i c 5.3.1) the  are  same  the  same.  value  of  displacements  of  for  geometrical  c l a y dam  ft.  and  two  of  the  structures  The  The  be  order  different.  In  p o s s i b l e d i f f e r e n c e s , the n o n - l i n e a r different soil  Station  1971  The  spectra  are  a c c e l e r a t i o n of Hughes  .35g  record  and was  The  to  on  15 seconds  rock. of  The the  inherently  investigate  element  the  method  at  Lake  digitized  is  earthquake had  of  an  of  .18  i s an earth  Hughes  rigid  acceleration  o b t a i n e d on rock and  As a  large  base  response a maximum  sec.  The  appropriate structure  most s i g n i f i c a n t motions occur i n the record.  may  Fernando earthquake  a predominant p e r i o d  f r e e f i e l d motion t o use when the base rest  taken  remaining  not  recorded  used as the  5-1.  This  structures.  a c c e l e r a t i o n time h i s t o r y and shown on F i g u r e  predicted  displacement f i e l d s  finite  component  12 i n C a l i f o r n i a was  motion.  Lake  N21E  given  consideration  In a l l the methods of a n a l y s i s , the San 9,  Chapter  i n f l u e n c e of the  c o n s i d e r e d i n the Newmark a n a l y s i s .  February  (see  w i l l be the same. no  to  different  the  s o i l mass of the s t r u c t u r e on the s l i d e mass i s  a p p l i e d t o the two  were used  two  acceleration,  differences.  of  s l o p e s of 2  analyses  c o r r e c t as t h e r e i s  significantly  view  Newmark a n a l y s i s i m p l i e s t h a t  yield  the  cannot be s t r i c t l y  The  In  i n the c o m p a r a t i v e s t u d y .  c l a y s l o p e s t r u c t u r e and  Each has  yield  block analogy.  first  p e r c e n t a g e of  the  68  FIG. 5-1  ACCELERATION AND  TIME  SAN  FERNANDO  RESPONSE  HISTORY  SPECTRA  FOR  EARTHQUAKE  69  permanent  deformation  should  occur  within  this  time,  (for  s t r u c t u r e s s u b j e c t e d t o t h i s motion) a d u r a t i o n of 15 seconds i s used  in  the  analyses.  The  time  earthquake r e c o r d i s used as the  interval  time  step  of the d i g i t i z e d (At= .02)  of  the  dynamic f i n i t e element a n a l y s e s . As  a  preliminary  v e r f i c a t i o n of the n o n - l i n e a r method, a  s i n g l e q u a d r i l a t e r a l f i n i t e element  comparison  with  Newmark's  procedure i s p r e s e n t e d i n the next s e c t i o n .  5.2  S i n g l e element comparison w i t h Newmark a n a l y s i s A  single  q u a d r i l a t e r i a l f i n i t e element r e s t i n g on a r i g i d  base w i t h a s s i g n e d s t a t i c and dynamic s o i l p r o p e r t i e s are in  F i g u r e . 5-2.  Nodes 3 and 4 a r e a l l o w e d t o t r a n s l a t e i n the  horizontal d i r e c t i o n only. elastic-plastic of  this  a  linear  manner t o shear as noted i n the e a r l i e r  section  of  s t a t i c shear s t r e s s Q whereby  on  The  chapter.  characteristics  shown  the st  element  The  in  shear s t r e s s - shear s t r a i n  element can  behaves  be  a r e shown i n F i g u r e 5-3. considered  a  'static  bias',  s h a k i n g , the a c c u m u l a t i o n of d i s p l a c e m e n t w i l l  t o be i n the d i r e c t i o n of the s t a t i c f o r c e .  The  The f o r m u l a t i o n  tend of  the mass m a t r i x assumes the mass of the element i s lumped a t the nodes,  see  conditions  Chapter is  3.3.2.  essentially  The a  finite  lumped  element  mass,  under these  spring  dashpot  m e c h a n i c a l model, where the s p r i n g and dashpot r e p l a c e the shear stiffness  and v i s c o u s damping of the element r e s p e c t i v e l y .  movement of the lumped masses a t  the t o p nodes of  The  the s t a t i c l y  10  T  100 Q= y  in  FIG. 5-2  FIG. 5-3  ft.  —  pcf 50-psf  Gmax=3000psf  SINGLE FINITE ELEMENT MODEL  SHEAR DEFORMATION  RELATIONSHIP  71  b i a s e d f i n i t e element the  Newmark  analogy  s h o u l d be v e r y s i m i l a r t o t h e of  a  rigid  T h e r e f o r e a v a l i d comparison The  to  t h e San  b l o c k on a i n c l i n e d p l a n e .  connected  Fernando  to  a  rigid  earthquake motion  base  nodes  for different  element method Newmark  a r e compared  (1965)  displacements coefficient motion.  i n Figure were  with 5-4.  computed  for  shear  values  i n the f i n i t e  displacements The  Newmark  a  range  of the  obtained  by  prediction  of  of  resistance  "N" v a l u e s u s i n g t h e s c a l e d San Fernando earthquake  The  resistance  static  is  scaled to a  maximum a c c e l e r a t i o n of .50g. The maximum d i s p l a c e m e n t top  in  can be made.  s i n g l e f i n i t e element  subjected  sliding  procedure  i s described  i n Chapter  2.2.  The  c o e f f i c i e n t c o r r e s p o n d s t o the y i e l d a c c e l e r a t i o n of  the f i n i t e element.  The y i e l d a c c e l e r a t i o n  i s related  to i t s  s t a t i c shear s t r e s s and was computed as f o l l o w s : N  was  defined  acceleration sliding  by  Newmark  acting  as  being  i n t h e proper  the c o e f f i c i e n t of direction  ( o r i n t h i s case y i e l d i n g of t h e f i n i t e  to  cause  element),  hence:  ( Q  where  Q  y  y  -  Q  s  t  )L= NW  y i e l d s t r e n g t h (50psf)  L  Length of element  N  Newmark's c o e f f i c i e n t  W  weight of the t o p h a l f of t h e element (2.5X10.X100  Q s t  s t a t i c shear  (10ft)  lb/ft ) 3  stress  (5-1 )  10.  •  •  •  1  -  _  LEGEND »SFE • Newmark  D  ^ 5  -  -  D  •  —  • •  _ •  • •  •  • • • B  0.0  1  1  1  .10  .20  .30  N _ A  F1G.5-4  Max Resistance Coeff Max Earthquake Acc  SINGLE FINITE ELEMENT NEWMARK DISPLACEMENT COMPARISON  r.  .40  73  Typical displacements  time  histories  computed  using  of  the  earthquake  induced  t h e two procedures a r e shown i n  F i g u r e 5-5 f o r comparison.  10.0  Time  FIG.  5-5  SINGLE  sees  FINITE E L E M E N T  NEWMARK  MODEL  ANALYSIS  DISPLACEMENT  HISTORY  COMPARISON  A s t a t i c b i a s of 43.75 l b s i s used i n t h e f i n i t e element method, c o r r e s p o n d i n g t o t h e Newmark N v a l u e of .025. As d i s p l a c e m e n t accumulates at  t h e end of  expected  the  w i t h time and t h e maximum v a l u e s o c c u r s  the shaking  period.  a n a l y s i s presented herein p r e d i c t s  I t may be seen t h a t t h e  displacements  that  are i n  74  good agreement w i t h the s i m p l e r Newmark r i g i d p l a s t i c model.  5.3  Dynamic Response of C l a y S t r u c t u r e s The  the  f i n i t e element  clay  slope  and F i g u r e 5-7  s t r u c t u r e and c l a y dam  respectively.  on a r i g i d rock base. that  r e p r e s e n t a t i o n s of the c r o s s s e c t i o n s of  The  are shown on F i g u r e 5-6  s t r u c t u r e s are assumed t o r e s t  The dynamic response  procedure  requires  the p r e - e a r t h q u a k e s t r e s s e s be determined b e f o r e h a n d .  s t a t i c a n a l y s i s i s performed u s i n g the SOILSTRESS  developed  by Byrne  same f i n i t e element g r i d  finite  (1981).  element  The  program  For each s t r u c t u r e , the  was  used  for  nature  of  static  both  the  static  and  dynamic a n a l y s i s . The  non-linear  behavior i s hyperbolic  modelled  assuming  relationship.  the  Duncan  one  iterative  this  problem  is  strains the  solved  used  According until  stress-strain  and  Chang  (1970)  A secant modulus as determined by the  h y p e r b o l i c r e l a t i o n s h i p i s used i n a approach.  shear  to  agreement  is  step  linear  method, obtained  elastic  the  linear  between  the  t o compute secant moduli i n each s o i l element  and  s t r a i n s d e v e l o p e d u s i n g these assumed secant m o d u l i . S t a t i c and dynamic s o i l parameters used f o r the a n a l y s i s of  the  c l a y s t r u c t u r e s i n l a t e r s e c t i o n s of t h i s c h a p t e r are  in  Table  I.  For  c l a y s , Duncan 1980 roughly initial  equal  to  saturated suggested E 600Su.  undrained ;  normally consolidated  (initial  Assuming  a  shown  Young's  modulus)  Poisson r a t i o  s t a t i c shear modulus Gj i s e q u a l t o  200SU  .  .50,  is the  Laboratory  80  160ft  Scale  FIG. 5-6  FINITE  E L E M E N T GRID  O F C L A Y SLOPE  STRUCTURE  I  I  0  80  1  160 ft  Scale  FIG.  5-7  FINITE ELEMENT GRID O F CLAY D A M  -si  CO  7 7  and  f i e l d s t u d i e s have shown t h a t f o r t h e same  same  conditions,  maximum dynamic shear modulus  than  the s t a t i c  equivalent  Accordingly, value  a  G  factor  m  a  under t h e  to  T h i s i s i n agreement  The  3 0 0 0 S u .  ( t h e r e f o r e making  G  structures.  with  m  a  x  shear  strength  t h e Seed  S  u  to 5 .  et. G  m  a  i s assumed  c o n s t a n t ) throughout t h e c l a y  a l . , ranges  x  constant  slope  and  In t h e dynamic a n a l y s i s , t h e i n i t i a l value of  bulk modulus i s c a l c u l a t e d from t h e G Poisson  4  of  f i n d i n g s , where l a b o r a t o r y r e s u l t s have shown  1 0 0 0  i s greater  x  t h e dynamic maximum shear modulus i s a s s i g n e d t h e  l O O O S u .  ( 1 9 6 9 )  dam  G ; by  soil  value  ma%  and an  assumed  r a t i o equal t o . 4 5 .  TABLE I S t a t i c and Dynamic S o i l  Analysi s  STATIC  Shear Strength lb/ft 2  2 9 4 0 .  Unit Wt. lb/ft  3  Properties  P o i sson Ratio  1 2 5 .  . 4 9  Shear modulus parameters  Failure Ratio  n= 0 Kg= 2 7 8 .  0 . 9  ( G i / S ^  2 0 0 . )  .01  DYNAMIC  elas-plas 2 9 4 0 .  1 2 5 .  G  max  =  . 4 5  /  S  u  1 0 0 0 .  1.0  hyperbolic In t h e n o n - l i n e a r method, t h e s t a t i c s t r e s s e s o b t a i n e d the  SOILSTRESS  program  a r e assumed  s t a t e f o r t h e dynamic a n a l y s i s . static  t o be t h e i n i t i a l  from  stress  F i n i t e elements a r e a s s i g n e d  a  b i a s on t h e i r shear s t r e s s - shear s t r a i n c u r v e equal t o  t h e i r s t a t i c shear s t r e s s .  The i n i t i a l  shear  modulus  i s the  78  G  m a x  value. It  has  been observed  i n the f i e l d t h a t s e i s m i c l o a d i n g of  e a r t h s t r u c t u r e s cause l o n g i t u d i n a l The  presence  of  a  tensile  vertical  crack  allows  tensile  f r e e movement i n a  d i r e c t i o n away and p e r p e n d i c u l a r t o the c r a c k p l a n e . horizontal  movement  corresponds  to  free  approach i n which t e n s i l e c r a c k s are modelled isotropic  assumptions  b e h a v i o r may  of  where Young's modulus movement  and  reduced  i s l i m i t e d by  the  T e n s i l e crack  i n an a n i s o t r o p i c manner, .  in  the  direction  of  by s o f t e n i n g of the s t i f f n e s s moduli G and B.  location  throughout would  The  within  allow  horizontal  movement  element.  at  carried  As w e l l , r e d u c i n g G t o s i g n i f i c a n t l y low v a l u e s  actually  model a l i q u i f i e d  analysis  was  c o n d i t i o n , which i s not the  incorporate carried  out  tensile by  crack  reducing  behavior  initial  principal assigned  static  stress a  value  o  stress 3  would  intent. into  the  o n l y the bulk modulus  v a l u e s , i n elements where t e n s i l e s t r e s s o c c u r . the  it  the e n t i r e  element.  to  the  Therefore  not be s t r i c t l y c o r r e c t t o reduce G throughout  be  While a  an element, shear s t r e s s can be  the i n t a c t remainder of the  An attempt  free  For an i s o t r o p i c m a t e r i a l , t h i s e f f e c t can  v e r t i c a l t e n s i l e c r a c k may crack  free  i s kept at the o r i g i n a l v a l u e i n the d i r e c t i o n of  the c r a c k p l a n e . achieved  is  This  volume change.  the p r e s e n t a n a l y s i s .  o n l y be p r o p e r l y modelled  cracks.  Starting  from  s t a t e , i f d u r i n g s h a k i n g the minor  becomes  tensile,  the  bulk  modulus  of 1 p e r c e n t of i t s i n i t i a l v a l u e .  modulus of the f a i l e d element i s l e f t at the reduced the remainder of the earthquake  loading.  is  The  bulk  value  for  79  5.3.1  The  Comparison w i t h t h e Newmark and the M a k d i s i - S e e d A n a l y s i s non-linear  dynamic response of a c l a y s l o p e s t r u c t u r e  and dam s u b j e c t e d t o compared  against  t h e San  Fernando  r e s u l t s obtained  both Newmark, and  earthquake  motion  from p r o c e d u r e s developed by  Makdisi-Seed.  Newmark's r e s i s t a n c e c o e f f i c i e n t "N" i s t h e v a l u e horizontal  seismic  of the  c o e f f i c i e n t determined from a p s e u d o - s t a t i c  conventional slope s t a b i l i t y of  is  a n a l y s i s which w i l l g i v e  a  factor  s a f e t y equal t o u n i t y .  F a i l u r e s u r f a c e , and s o i l p r o p e r t i e s  used i n the s l o p e s t a b i l i t y  a n a l y s e s of t h e c l a y s t r u c t u r e s a r e  shown (or for  i n F i g u r e 5-8 and F i g u r e 5-9.  The r e s i s t a a c e c o e f f i c i e n t  y i e l d a c c e l e r a t i o n ) v a l u e s a r e c o i n c i d e n t l y equal t o the  two  clay  structures.  Computations of  u s i n g Newmark's procedure a r e performed u s i n g the earthquake  scaled  to  O.lOg  displacements San  Fernando  maximum a c c e l e r a t i o n v a l u e s r a n g i n g  0.4 t o 1.Og., f o r t h e determined  yield  acceleration  value  from of  0.1Og. Refinements throughout t h e Makdisi-Seed  to slope  allow and  slide  mass  as  a n a l y s i s a r e made as f o l l o w s .  s t r u c t u r e s when s u b j e c t e d t o evaluated  for variation  using  the  San  the dynamic f i n i t e  in acceleration  prescribed  The response of t h e  Fernando  earthquake  element program QUAD4.  a n a l y s i s employs t h e use of shear modulus r e d u c t i o n and ratios  curves  which  by the  The  The  damping  a r e b u i l t i n t o t h e computer program.  attempt i s made t o modify these c u r v e s .  is  s o i l elements  No are  FIG. 5-8  SLOPE STABILITY ANALYSIS  ON  CtAY  SLOPE  STRUCTURE  4. Center of Failure C i r c l e  FIG. 5 - 9 SLOPE STABILITY ANALYSIS  ON  CLAY  DAM CO  82  are a s s i g n e d the same G  v a l u e of lOOOSu.  m a x  the average a c c e l e r a t i o n of the s l i d e dynamic  response  analysis  is  d i p l a c e m e n t of the s l i d e mass. weighted  The time h i s t o r y of  mass  obtained  from  to  estimate  permanent  used  The average a c c e l e r a t i o n  the  is  the  average of the a c c e l e r a t i o n v a l u e s of the nodes w i t h i n  the s l i d e mass as g i v e n by, a ( t ) = Zntj . X ; ( t ) Im j where  mj and Xj are the mass and a c c e l e r a t i o n at node i .  summation i s taken over a l l the nodes w i t h i n the s l i d e The  displacement  for  both  the by  and  approaches  are  acceleration  i n excess of the y i e l d a c c e l e r a t i o n  for  determined  Newmark  integrating  the d u r a t i o n of the earthquake In  the  finite  and  mass.  Makdisi-Seed the  effective  value  (O.lOg)  motion.  element method each mode w i t h i n the domain  d i s p l a c e s a d i f f e r e n t amount. Newmark  The  Makdisi-Seed  In  order  methods  an  to  compare  average  of  with the  the final  d i s p l a c e m e n t s of the nodes w i t h i n the s l i d e mass i s used.  5.3.1.1  C l a y Slope  Comparison  The r e s u l t s from the Newmark linear  finite  element  type  analysis  computed  using  the  Typical  in  Figure  5-11.  The  the  displacement  non-  time  Newmark a n a l y s i s and n o n - l i n e a r  f i n i t e element method under the same earthquake shown  and  of the c l a y s l o p e are i n good  agreement as shown i n F i g u r e 5-10. histories  methods  non-linear  conditions finite  are  element  10B-  •  LEGEND °  A  5h  • N FE • Newmark M-S  _  A  A •  A  a  A •  .1.0  .20  N _ A FIG 5-10  .30  Max Resistance Coeff Max Earthquake Acc  CLAY SLOPE NON-LINEAR FINITE ELEMENTMAKDISI  SEED-NEWMARK  DISPLACEMENT  COMPARISON  84  j r~*-  N ^NFE  0  J  5  10  Time  FIG 5-11  sees  CLAY SLOPE NON-LINEAR FINITE ELEMENT NEWMARK ANALYSIS DISPLACEMENT HISTORY COMPARISON  15  85  displacement  history  shows  an  some  uphill  Newmark a n a l y s i s i g n o r e s , and my be one displacements  are  larger.  t y p i c a l displacement Figure  5-12,  is  For  of  pseudo  or  slip  Newmark  sliding  plane.  of  The  grid  the  in  entire  circular  arc  The y i e l d a c c e l e r a t i o n o b t a i n e d  s t a t i c a n a l y s i s i s f o r the f a i l u r e s u r f a c e w i t h  the minimum f a c t o r of s a f e t y . for  the  the c l a y s l o p e s t r u c t u r e , the  translation  f a i l u r e does not seem a p p a r e n t . a  reason why  p a t t e r n as shown by the d i s p l a c e d  embankment a l o n g a h o r i z o n t a l  from  motion which the  T h e r e f o r e the y i e l d  acceleration  any o t h e r f a i l u r e s u r f a c e ( i n t h i s case of a s l i d i n g  would be somewhat h i g h e r v a l u e .  A  higher  yield  block)  acceleration  v a l u e would r e s u l t i n lower Newmark p r e d i c t e d d i s p l a c e m e n t s . From  the  formulation  of  the  dynamic  a n a l y s i s , the f o r c e s on a s t r u c t u r e caused  by  finite the  element  rigid  base  acceleration  i s r e p r e s e n t e d by i n e r t i a l o a d s a t the f r e e nodes,  (Chapter  Inertia forces  3).  for  any  given  row  of  elements  should  be g r e a t e s t f o r the row of elements j u s t above the  base.  With  sliding  this  along  reasonable.  the  results  indicate that elements  bottom  row  the  of  in  located  by Seed  (1979).  of  non-linear  the  general, at  the  tensile top  v e r t i c a l l y downwards from t h e r e .  predicted  elements  As w e l l , h o r i z o n t a l s l i d i n g  have been observed The  consideration,  of  behavior  ( F i g u r e 5-12) embankment  finite stresses  is  of is  slopes  element a n a l y s i s first  of the embankment, and This  rigid  in  f i e l d o b s e r v a t i o n s , where v e r t i c a l t e n s i l e c r a c k s  occur  progress  agreement are  in  with  known to  Displacements Magnificied  30 times  FIG. 5-12  DISPLACED GRID OF CLAY SLOPE STRUCTURE  00  CO  87  occur al.,  in  seismicly  1973).  loaded  W i t h these t e n s i l e c r a c k s , t h e r e i s  additional  deformations.  t e n s i l e f a i l e d elements structure  s l o p e s and embankments, (Seed e t .  do  As in  observed  the  upper  in  potiental  Figure  section  5-12,  of  the  D i s p l a c e m e n t , as noted  previously,  i s o t r o p i c assumptions  the p r e s e n t a n a l y s i s may observed  present  may  account  for  the  and p r e d i c t e d b e h a v i o r .  be  mainly  elements.  i n the treatment of t e n s i l e c r a c k s i n  noted t h a t w i t h the development of t e n s i l e bias  slope  is  a t t r i b u t e d t o the shear movement i n the bottom row of  between  the  not r e s u l t i n a d d i t i o n a l r e l a t i v e d i s p l a c e m e n t of  these elements.  The  for  reduced.  Therefore  lack  of  agreement  Though, i t s h o u l d be cracks, in  any  static  an a n i s o t r o p i c  a n a l y s i s , a d d i t i o n a l d i s p l a c e m e n t s where t e n s i l e c r a c k s  develop  may  not be as s i g n i f i c a n t , as t h e r e i s an accompanying r e d u c t i o n  in  statis bias.  O v e r a l l , d e f o r m a t i o n s may  not be  substantially  d i f f e r e n t than t h a t p r e d i c t e d i n the p r e s e n t a n a l y s i s . another i m p o r t a n t c o n s i d e r a t i o n . structure  during  the  Static gravity  earthquake  dynamically softened s o i l m a t e r i a l . earthquake While  loadings  may  explain  5.3.1.2 It  is  C l a y Dam observed  p r e d i c t s markedly  carried  the by  displacements.  the observed f i e l d d i s p l a c e m e n t s , the  s t a t i c e f f e c t has not been f o r m u l a t e d i n t o t h i s r e a s o n , agreement may  is  of  Under these c o n d i t i o n s , the  l o a d i n g w i l l cause a d d i t i o n a l s t a t i c  this  loads  There i s  the  analysis.  For  not be p o s s i b l e .  Comparison i n F i g u r e 5-13  t h a t the n o n - l i n e a r method  lower d i s p l a c e m e n t s f o r the c l a y dam  than does  D  LEGEND  A  • N FE n Newmark A M-S •  A  •  A •  A •  : .10  .20  N _ A  FIG 5-13  i_2  2  i  CLAY  .30  Max Resistance Coeff Max Earthquake Acc  DAM  NON-LINEAR MAKDISI  FINITE ELEMENT -  SEED-NEWMARK  DISPLACEMENT  COMPARISON  89  the  Newmark  displacement  or  Makdisi-Seed  history  i n Figure  element method p r e d i c t s t h a t significant  method.  during  5-14,  the  dam  negative displacements.  can be e x p l a i n e d i n t h i s way. shaking,  outward  As  indicated  the  by  the  non-linear f i n i t e  sliding  mass  undergoes  The n e g a t i v e d i s p l a c e m e n t s  I t i s reasonable t o  assume  that  movement of a s l i d e mass on the o t h e r  embankment s l o p e s h o u l d o c c u r .  The outward movement of a  slide  mass on one s l o p e w i l l a l l o w movement of the o t h e r s l i d e mass i n its  inward  direction,  as  resistance  d i r e c t i o n i s g r e a t l y reduced due t o outwarding s l i d e mass. for  overall  reduced.  outward  to  temporary  movement failure  loading  of the  I f t h i s i s t h e c a s e , t h e r e i s a tendency displacement  of  the  slide  mass  t o be  dam  after  A t y p i c a l d i s p l a c e m e n t p a t t e r n of the c l a y  earthquake  i n that  Figure  5-15, shows net outward movement of  the each dam s l o p e . In a d d i t i o n . , lower d i s p l a c e m e n t e s t i m a t e s may as  the  o v e r a l l s t a t i c b i a s on the c l a y dam i s z e r o .  of the shear s t r e s s e s on the dam s l o p e s and  be  opposite.  are  expected The t o t a l  essentially  equal  Permanent d i s p l a c e m e n t s s h o u l d be lower than i n  the c l a y s l o p e s t r u c t u r e , where t h e r e i s a net s t a t i c b i a s . Newmark type e s t i m a t e s conservative  under  of  permanent  these  displacement  conditions.  are  The  overly  Displacements  for  s y m m e t r i c a l r e s i s t a n c e (Newmark, 1965) a r e an o r d e r of magnitude lower  than  f o r unsymmetrical  resistance.  Comparison  with  s y m m e t r i c a l d i s p l a c e m e n t s may be more a p p l i c a b l e i n t h e c l a y dam where net s t a t i c b i a s i s z e r o .  VI  A1/  1  ft  w  \f\rJ\ \A AA ^ 1  V  •v  1  a  0  v •  5  mQx= -  10 Time sees  FIG. 5-14 CLAY DAM DISPLACEMENT TIME HISTORY  7 5  9  15  Displacements Magnificied  30 times  FIG. 5-15  DISPLACED  GRID O F  CLAY D A M  92  5.3.2  As  Comparison of H y p e r b o l i c and E l a s t i c - P l a s t i c Shear s t r e s s - s t r a i n Laws  outlined  i n Chapter  h y p e r b o l i c and h y s t e r e t i c law,  while  chapter,  used  is a  f o r comparative  The  h y p e r b o l i c curve relationship,  manner.  simple  relationship.  3,  in The  approximation  elastic-plastic  to  R= f  earlier the  shear  in  this  hyperbolic  curve a p p r o x i m a t i o n  i s shown i n F i g u r e 5-16.  (where  structure  elastic-plastic  purposes  U s i n g the  t o the  hyperbolic  1.0) t h e dynamic a n a l y s i s w i l l c l o s e r  model and p r e d i c t the d i s p l a c e m e n t s . slope  shear, s o i l behaves i n a  For t h i s reason, t h e  clay  i s r e - a n a l y s e d , and comparison of the s i m p l e r  and a c t u a l r e l a t i o n s h i p i s made. The  hyperbolic  displacement two  and  elastic-plastic  a r e shown on F i g u r e 5-17.  of  As might be expected the  s t r e s s s t r a i n law r e s u l t s agree q u i t e c l o s e l y f o r the l a r g e  d i s t u r b a n c e s where t h e i r s t i f f n e s s a r e s i m i l a r d i s t u r b a n c e s the d i s p l a c e m e n t s  shear  strains.  shown on F i g u r e 5-18.  , while f o r small  p r e d i c t e d by the h y p e r b o l i c  law a r e g r e a t e r which r e f l e c t s t h i s small  predictions  laws  lesser  T y p i c a l displacement  shear  stiffness  at  time h i s t o r i e s a r e  Shear strain  FIG. 5-16 ELASTIC PLASTIC APPROXIMATION OF HYPERBOLIC CURVE  10  LEGEND O Hyperbolic • Elastic-Plastic  o  O  .10  .20 . N A  FIG. 5-17  Max Resistance Coeff Max Earthquake Acc  HYPERBOLIC  ELASTIC PLASTIC  DISPLACEMENT COMPARISON CLAY  SLOPE  .30  95  6.0 Hyperbolic  ^ \ V  /Elastic Plastic  _ 3.0 c to E  a  0  5  max= -  10 Time sees  FIG. 5-18 HYPERBOLIC  ELASTIC PLASTIC  DISPLACEMENT (CLAY SLOPE)  HISTORY COMPARISON  7 5  9  15  96  CHAPTER 6 SUMMARY AND CONCLUSIONS 6.1  Summary A t w o - d i m e n s i o n a l f i n i t e element  predicting  the  structures to behavior approach  of  stress  seismic the  and  permanent  loading  soil  method  of  displacements  i s presented.  i s modelled  analysis  by  of e a r t h  The  inelastic  an i n c r e m e n t a l l i n e a r  i n which t h e tangent shear modulus i s v a r i e d  l e v e l of t h e shear s t r a i n .  for  with  the  The shear s t r e s s s t r a i n r e l a t i o n was  modelled by h y p e r b o l i c l o a d i n g and u n l o a d i n g c u r v e s l e a d i n g t o a Masing  type  of  energy  d i s s i p a t i o n , and by the more c l a s s i c a l  e l a s t i c - p l a s t i c c o n s t i t u t i v e law. The e q u a t i o n s of motion of t h e s t r u c t u r e a r e the  Newmark  step-by-step  time i n t e g r a t i o n p r o c e d u r e .  time s t e p t h e tangent shear and Hysteretic  solved  bulk  modulus  using  At each  are evaluated.  damping as a r e s u l t of the s t r e s s - s t r a i n l o a d i n g and  u n l o a d i n g c u r v e s i s i n h e r e n t i n t h e model.  V i s c o u s damping  may  a l s o be i n c l u d e d . The number of analysed  non-linear soil using  dynamic response a n a l y s i s was a p p l i e d t o a  structures. Newmark  type  The  same  approaches  structures  were  so t h a t a comparison  c o u l d be made on t h e p r e d i c t i o n of permanent d i s p l a c e m e n t s . order  also  In  t h a t comparison be m e a n i n g f u l , the n o n - l i n e a r a n a l y s i s i s  based on an e l a s t i c - p l a s t i c shear law. In o r d e r t o c l o s e r model and p r e d i c t d i s p l a c e m e n t s a c l a y s l o p e was r e - a n a l y s e d u s i n g the  97  h y p e r b o l i c shear  relationship.  W h i l e the Newmark methods g i v e s a s i n g l e v a l u e e s t i m a t e permanent freedom  displacements, analysis  is  t h e move  desirable  rigorous as  the  of  multi-degree  of  distribution  of  d i s p l a c e m e n t s w i t h i n t h e s t r u c t u r e s can be o b t a i n e d .  6.2  Conclusions The c o n c l u s i o n s reached 1)  method  The  Newmark  from t h i s r e s e a r c h a r e as f o l l o w s :  methods and the n o n - l i n e a r f i n i t e element  for predicting  single  value  estimates  of  permanent  d i s p l a c e m e n t s of t h e c l a y s l o p e s t r u c t u r e a r e i n good agreement. However,  the n o n - l i n e a r  method  h o r i z o n t a l f a i l u r e s u r f a c e which assumed  circular  predicts  movement  i s d i f f e r e n t than  the  a r c f a i l u r e and more i n keeping w i t h  along  a  Newmark observed  movements. 2)  For p r e d i c t i o n  of  displacement  of  dam  s l o p e s , the  Newmark methods a r e o v e r l y c o n s e r v a t i v e . 3)  The  simpler  elastic-plastic  shear law and t h e a c t u a l  h y p e r b o l i c r e l a t i o n s h i p g i v e s d i s p l a c e m e n t e s t i m a t e s t h a t a r e of the same o r d e r .  6.3  Suggestions  f o r Further  Research  I n c o r p o r a t i n g a pore p r e s s u r e model i n t o t h e a n a l y s i s allow  the  evaluation  of  t h e very  important  problem  will of  98  liquefaction potential. used  to  The e f f e c t i v e s t r e s s  p r e d i c t displacements  analysis  can  be  f i e l d s of s a t u r a t e d c o h e s i o n l e s s  soil structures. T e n s i l e c r a c k b e h a v i o r may o n l y be p r o p e r l y modelled anisotropic  manner,  where  Young's  modulus  i s reduced  in  an  i n the  d i r e c t i o n of f r e e movement (away and p e r p e n d i c u l a r t o the  crack  p l a n e ) and i s kept a t the o r i g i n a l v a l u e i n the d i r e c t i o n of the crack plane. Static  gravity  d u r i n g earthquake  l o a d s c a r r i e d by d y n a m i c a l l y s o f t e n e d  l o a d i n g may cause  The e f f e c t s h o u l d be modelled If  possible,  displacements.  and i n c l u d e d i n the a n a l y s i s .  correlation  s h o u l d be c a r r i e d o u t .  additional  soil  study  w i t h observed  field  data  99  BIBLIOGRAPHY 1.  Byrne P.M., and Janzen W., 1981, S o i l s t r e s s : A Computer Program f o r N o n l i n e a r a n a l y s i s of S t r e s s e s and Deformations i n S o i l , S o i l Mechanics S e r i e s NO. 52, Department of C i v i l E n g i n e e r i n g , U n i v e r s i t y of B r i t i s h Columbia.  2.  Clough, R.W., and P e n z i e n , 1975, Dynamics of S t r u c t u r e s ,McGraw-Hi11, Kogakusha, L t d .  3.  Cook R.D., 1974, Concepts and A p p l i c a t i o n s of F i n i t e Element A n a l y s i s , John W i l e y & Sons.  4.  D i b a j , M., and P e n z i e n , J . , 1969, N o n l i n e a r S e i s m i c Response of E a r t h S t r u c t u r e s , Earthquake E n g i n e e r i n g Research Center U n i v e r s i t y of C a l i f o r n i a , B e r k e l e y , C a l i f o r n i a , January 1969.  5.  Duncan, J.M., Byrne, P.M., Wong, R.S., and Mabry, P. 1980, S t r e n t h , S t r e s s - S t r a i n and Bulk Modulus Parameters f o r F i n i t e Element A n a l y s e s , Report No. UCB/GT/80-10 t o N a t i o n a l Science Foundation.  6.  Duncan, J.M., and Chang, C.Y., 1970, Non-Linear A n a l y s i s of S t r e s s and S t r a i n i n S o i l s , J o u r n a l of the S o i l Mechanics and Foundations D i v i s i o n ASCE, V o l . 96, No. SM5, pp. 1629-1653.  7.  F i n n , L.W.D., M a r t i n , G.R., and Lee, M.K.W., 1979, C o p a r i s o n of Dynamic A n a l y s e s f o r S a t u r a t e d Sands , pp. 473-491.  8.  F i n n , L.W.D., and M i l l e r R.I.S., 1971, Dynamic A n a l y s i s of Plane Non-Linear E a r t h S t r u c t u r e s ,  9.  F i n n , L.W.D., Byrne P.M., and M a r t i n G.R, 1978, S e i s m i c Response and L i q u e f a c t i o n of Sands , JSMFD, ASCE, V o l . No. GT8, August 1978 pp. 841-856.  10.  Goodman, R.E., and Seed, H.B., 1966 Earthquake-Induced Displacements of Sand Embankments , J o u r n a l of the S o i l Mechanics and Foundations D i v i s i o n , ASCE, V o l . 92, No. SM2, March 1966.  11.  H a r d i n , B.O., and D r n e v i c h , V.P., 1972, Shear Modulus and Damping i n S o i l s : Design E q u a t i o n s and Curves , JSMFD, ASCE, V o l . 98, No. SM7, pp. 667-692.  12.  I d r i s s , I.M., Dobry, R., and S i n g h , R.D., 1979, N o n l i n e a r Behavior of S o f t C l a y s D u r i n g C y c l i c Loading , J o u r n a l of G e o t e c h n i c a l E n g i n e e r i n g D i v i s i o n , ASCE, V o l . No. GT12, Dec. 1979, pp 14271447.  100  13.  I d r i s s , I,M., and Seed, H.B., 1974, S e i s m i c Response by V a r i a b l e Damping F i n i t e Elements , J o u r n a l of the G e o t e c h n i c a l E n g i n e e r i n g D i v i s i o n , ASCE, V o l . No. GT1, January 1974, pp. 1-13.  14.  Lee, M.K., 1977 Dynamic E f f e c t S t r e s s Response A n a l y s i s , Ph.D. T h e s i s The F a c u l t y of Graduate S t u d i e s , Department of C i v i l E n g i n e e r i n g , The U n i v e r s i t y of B r i t i s h Columbia.  15.  M a k d i s i , F . I . , 1976, Performance and A n a l y s i s of E a r t h Dams d u r i n g Strong Earthquakes , D i s s e r t a t i o n f o r Ph.D. i n E n g i n e e r i n g , Graduate D i v i s i o n , U n i v e r i t y of C a l i f o r n i a , B e r k e l e y , C a l i f o r n i a , 1976.  16.  M a k d i s i , F . I . , and Seed, H.B., 1978, S i m p l i f i e d Procedure f o r E s t i m a t i n g Dam and Embankment Earthquake Induced P i s p l a c e m e n t s , J o u r n a l of t h e G e o t e c h n i c a l E n g i n e e r i n g D i v i s i o n , ASCE, V o l . 104, No. GT7, 1978, pp. 849-867.  17.  M a s s i n g , G., 1926, Eiqenspanningen and V e r f e s t i g u n g beim Messing , P r o c . 2nd I n t e r n a t i o n a l Cong. A p p l i e d Mech., 1926 pp. 332-335.  18.  Newmark, N.M., 1965, E f f e c t s of Earthquakes on Dams and Embankments , Geotechnique, London, E n g l a n d , V o l XV, No. 2, 1965.  19.  Salgado F.M, 1981, S e i s m i c Response of R e t a i n i n g S t r u c t u r e s , M a s t e r s T h e s i s , The F a c u l t y of Graduate S t u d i e s , Department of C i v i l E n g i n e e r i n g , The U n i v e r s i t y of B r i t i s h Columbia, 1981.  20.  Sarma, S.K., 1975, S e i s m i c S t a b i l i t y of E a r t h Dams and Embankments , Geotechnique, Lond, England, V o l 25, No. 4, 1975, pp. 743-761.  21.  Seed, H.B., 1966, A Method f o r Earthquake R e s i s t a n t Design of E a r t h Dams , J o u r n a l of t h e G e o t e c h n i c a l E n g i n e e r i n g D i v i s i o n , ASCE. no. SM1 1966.  22.  Seed, H.B., 1979, C o n s i d e r a t i o n s i n t h e Earthquake R e s i s t a n t Design of E a r t h and Rock F i l l Dams , G e o t e c h n i q u e , London, E n g l a n d , No. 3, 1979, pp. 215-263.  23.  Seed, H.B., and I d r i s s , I.M., 1970, S o i l M o d u l i and Damping F a c t o r f o r Dynamic Response A n a l y s i s , Report No. EERC 70-10, Earthquake E n g i n e e r i n g Research C e n t e r , Berkeley, C a l i f o r n i a .  24.  Seed, H.B., M a k d i s i , F . I . , L e e , K.L., and I d r i s s , I.M., 1973, A n a l y s i s of t h e S l i d e s i n the San Fernando Dams d u r i n g t h e Earthquake of F e b r u a r y 9, 1971 , Earthquake E n g i n e e r i n g Research Center Report No. EERC 73-2  101 U n i v e r i t y of C a l i f o r n i a , B e r k e l e y , June. 25.  Seed, H.B., M a k d i s i , F . I . , and DeAlba, P., 1975, Performance of E a r t h Dams D u r i n g Earthquakes , J o u r n a l of the G e o t e c h n i c a l E n g i n e e r i n g D i v i s i o n , ASCE, V o l . 104, No. GT7, 1975.  26.  W i l s o n , E.L., and Clough, R.W., 1962, Dynamic Response by Step-by-Step M a t r i x A n a l y s i s , Symposium on Use of Computers i n C i v i l E n g i n e e r i n g , L i b s o n , P o r t u g a l , October 1 962.  27.  Z i e n k i e w i c z O.C., 1979, The F i n i t e Elemenet Method , Tata M c G r a w - H i l l P u b l i s h i n g Company L i m i t e d .  102  APPENDIX I F o r m u l a t i o n of the Element S t i f f n e s s M a t r i x f o r Plane I s o p a r a m e t r i c Q u a d r i l a t e r a l Elements  Consider a q u a d r i l a t e r a l having eight namely u that  and v  has  F i g u r e I - 1 a may  sides  but  a  distortion  N,  N  0  2  0  N  N  2  3  0  An element  of  a  'parent'  I-1b.  A mapping f u n c t i o n i s e x p r e s s e d  0  freedom  i s o t h e r w i s e of a r b i t r a r y shape  be c o n s i d e r e d as  r e c t a n g u l a r element, F i g u r e  0  of  at each of the f o u r c o r n e r nodes i .  straight  N ,  degrees  0 N  y i  0  N ,  0  3  as,  x  2  (1-1)  N ,  yo  where, (1  N i =  -  -  7?),  N  2  =  (  1 +  4  3  relates  4  a  unit  square i n the £TJ c o o r d i n a t e s to  p o i n t s i n the q u a d r i l a t e r a l i n xy shape x  are  determined  i,yi,x ,y ,....y<,. 2  2  (1-2)  4  4  mapping  1)  N = ( 1 - £) ( 1 + T?)  N = (l+jMl+n)-,  This  j) ( 1 4  by  the  coordinates eight  whose nodal  size  and  coordinates  u  2  1  x  FIG I - l a  2  X  QUADRILATERAL  ElEMENT  V  0,1  i.o  FIG I-lb UNIT SQUARE  €  1 04  The  axes £7} a r e  orthogonal merely  for a  rectangular  dimensionless  coordinates.  The  i n general  not  orthogonal.  They  are  element, i n which case they a r e  forms  of  rectangular  centroidal  2 by 2 r e c t a n g u l a r element, f o r which we may  w r i t e £= x and 77= y, i s a c o n v e n i e n t  s p e c i a l case t o r e f e r  when  s t u d y i n g t h e f o l l o w i n g development. Displacements  within  the  element a r e d e f i n e d by t h e same  i n t e r p o l a t i o n f u n c t i o n s as used t o d e f i n e the element shape hence t h e name i s o p a r a m e t r i c . {f}=  {u v}= [N]{u, v, u  2  where [N] i s t h e r e c t a n g u l a r indicated  i n the equation  displacements. to  and  Thus, ... v„} = [N]{d} matrix  of  (1-3)  the N  equations  as  (1-1 ), and {d} i s t h e v e c t o r of nodal  The d i s p l a c e m e n t s  u and v a r e d i r e c t e d  parallel  x and y, and not p a r a l l e l t o t h e £ and 17 axes. Note  that  f o r any  point  i n t h e u n i t i n £77 c o o r d i n a t e s ,  equation  1-1 g i v e s a p o i n t ( x , y ) i n t h e  equation  1-3 w i l l g i v e t h e d i s p l a c e m e n t s  nodal  displacements.  Thus  (x,y) i n t h e r e a l c o o r d i n a t e s  the  real  coordinates  and  u and v i n terms of the  displacements  of every  i s d e f i n e d i n terms of  the  point nodal  di splacements. Steps  taken  to  formulate  the  element  stiffness  a c c o r d i n g t o t h e approach o u t l i n e d by Cooke (1975) a r e in  the  remainder  of  t h i s appendix.  will  be  detailed  Because i t i s i m p o s s i b l y  t e d i o u s t o w r i t e t h e shape f u n c t i o n s i n terms of x formulation  matrix  and  y, the  c a r r i e d out i n terms of t h e i s o p a r a m e t r i c  105  coordinates  Equations  £77.  (1-3)  ( 1 - 1 ) ,  are  rewritten  for  convenience i n the form, x= IN, Xj  (1-4)  4  U= I N ; U i  v=  ZN;v  A r e l a t i o n s h i p between d e r i v a t i v e s i n the two c o o r d i n a t e s , from the c h a i n r u l e of d i f f e r e n t i a t i o n ,  <k  (  ( ).,  < ).y  i s as  follows:  )•*  (1-5)  =[J]  -  where commas denote p a r t i a l d i f f e r e n t i a t i o n . (1-4)  [ J ] i s the J a c o b i a n  N  [J]=  N,  M  Ni..  N , 2  matrix, xi  yi  3.£  x  y  N ,n  x  N  2  From e q u a t i o n  3  2  3  x„  y  2  (1-6)  3  y.  where, f o r example, N  = " (1  -  (1-7)  T?)  U s i n g the i n v e r s e r e l a t i o n [J]" , differentials 1  U,y •  V,y <  -  =  0  0  1 J  0  0  2  2  2  be w r i t t e n  [J*]=  as  system a s ,  Jl 1 J12 J  (1-4) and l e t t i n g  i n the xy system may  d i f f e r e n t i a l s i n the  •<  of e q u a t i o n  0  0  Ji 1 J 1 2  0  0  J  2  1 J 22  u. U  A  £  (1-8)  106  The  s t r a i n - d i s p l a c e m e n t r e l a t i o n may 1 0 (e}  =  €y  V  and,  •  0  —  0  0  J  from e q u a t i o n  0  0  0  1  be w r i t t e n  U,y  (1-9)  •  V,*  1 1 0  ^  V,y  (1-4), -  —  0  «  Ui Vi  0  u v u v  0 0  Ni.K  —  Substitution thereafter  equation  substitution  y i e l d s the r e l a t i o n {e}= The  of  stationary  Cooke (1975), and integration  of  elasticity  expressed  in  £TJ  the  3 3  (1-10)  (1-9),  i n t o the r e s u l t ,  matrix  using  and  and  is  obtained  j f ^ B ] [D] [B]dxdy. [B]  coordinates.  is The  as  established  elasticity  (3-20), and the  above,  constitutive  (3-21), (3-22).  patterns  For the  i n terms of  the  present a  shear  for further  s t i f f n e s s m a t r i c e s are  by i n t e g r a t i o n i n the $77 c o o r d i n a t e system u s i n g  the  m a t r i x [D] i s  and b u l k , p a r t i a l s t i f f n e s s m a t r i c e s , see Chap 3.3.1 element ' p a r t i a l '  by  by  Where [D] i s  1  G, B model, e q u a t i o n s  The  the  p o t i e n t a l energy i s w e l l d e s c r i b e d  a n a l y s i s the s t i f f n e s s m a t r i x i s expressed  details.  and  [B]{d}.  the e x p r e s s i o n matrix,  (1-10)  2  equation  equation  Z i e n k i e w i c z (1979)  determined from e q u a t i o n for  of  into  d e r i v a t i o n of the element s t i f f n e s s  principle  the  (1-8)  2  --A - I  i=3  of  as  obtained following  107  transformation [k] =  relationship.  //[B] [D][B]dxdy  The determinant yields  [B] [D][B]det[J]d£d7j  T  area  of  dxdy  [J]  from  (1-11)  T  is a  area  magnification  factor  that  and which i s a f u n c t i o n of  d£d7j,  p o s i t i o n w i t h the element. Having expressed integration  can  a l l m a t r i c e s i n terms  be  carried  in  of  £  £ 7 ? coordinates.  and  rj, the  However, i n  g e n e r a l the i n t e g r a t i o n cannot be performed e x a c t l y due complexity  of  the  polynomials  to  the  £ and r\ t h a t appear i n the  in  denominator of [ J * ] . Hence, n u m e r i c a l  integration  is  commonly  used. The  Gauss  method  of n u m e r i c a l e v a l u a t i o n of t h e d e f i n i t e  i n t e g r a l i s used i n the p r e s e n t a n a l y s i s . quadrature  formula  In two dimensions the  i s o b t a i n e d by i n t e g r a t i n g f i r s t w i t h r e s p e c t  t o one c o o r d i n a t e and then w i t h r e s p e c t t o the o t h e r ,  [k] =  f U , 7 ? ) d £ d T j = /ZW, f ( £ , 7 7 ) d r ? LWi  where  f ( £ , 7 j )  f ( £ , 7 ? )  drj  =  LEW,  Wj  f  (1-12)  ( £ , 7 ? )  r e p r e s e n t s the e x p r e s s i o n  [B] [D] [B]det [ J ] .  The i n t e g r a l i s e v a l u a t e d as t h e summation of v a l u e s of the function  at  selected  a p p r o p r i a t e 'weight' points  W .  points,  each m u l t i p l i e d by an  The Gauss method l o c a t e s t h e  sampling  so t h a t f o r a g i v e n number of them, g r e a t e s t a c c u r a c y i s  obtained. respect  sampling  Sampling to  the  points  center  are  located  of the element.  symmetrically  with  In g e n e r a l , the Gauss  108  quadrature using  'n' p o i n t s  p o l y n o m i a l of degree 2n-1 function degree  f(£,T?)  is  is  exact  or l e s s .  effectively  if  the  is  In u s i n g n p o i n t s , the replaced  by  a  given  a p o l y n o m i a l of  2n-1.  A sequence of s o l u t i o n s t o a problem may using  integral  successively  be  f i n e r meshes of elements.  be e x p e c t e d t o converge t o the c o r r e c t  result  generated  by  The  sequence  may  if  the  element d i s p l a c e m e n t f i e l d s s a t i s f y the c r i t e r i a o f ; continuity  of  assumed  invariance,  d i s p l a c e m e n t s w i t h i n elements, r i g i d body modes,  c o n s t a n t s t r a i n b e h a v i o r , and  interelement  compatibility.  The  v a l i d i t y of the i s o p a r a m e t r i c  element i s examined w i t h regard  to  convergence by Cooke, w i t h the element d i s p l a c e m e n t f i e l d s shown to  satisfy  the  above c r i t e r i a .  i n t e g r a t e d elements, there is,  numerically  addition consideration.  That  convergence t o the c o r r e c t s o l u t i o n i s p o s s i b l e only i f the  numerical i n t e g r a t i o n exact.  This  numerically  adequate  the is  element relaxed  stiffness and  to  is  evaluate  essentially  the  expression.  As  refined  if  numerical  the element a r e a e x a c t l y .  by n o t i n g t h a t , from the  derivation,  be  is  re-stated;  formation  same a  as  mesh  of  a  integration is  refined  f o r an element assumes the  integration This  stiffness  is  statement  stationary  potential matrix  is  of  a  strain  energy  and  a  constant  strain  strain  energy  c o n d i t i o n comes t o p r e v a i l i n each element, expression  can  integral  i n t e g r a t e d elements y i e l d convergence toward c o r r e c t  be e x p l a i n e d  energy  of  condition  r e s u l t s as the mesh  can  i s one  But as w e l l , w i t h  the  form, (1-13)  109  Hence, i n the l i m i t correctly  assumed  the  strain  i f the  energy  of  the  structure  is  volume of each element i s c o r r e c t l y  assessed. E x a m i n a t i o n of the J a c o b i a n determinant y i e l d s of  Gauss  elements. and  points  needed  to  For a plane q u a d r a t i c element d e t [ J ] c o n t a i n terms £  3  can  be accepted.  the  +.57735...  Thus e q u a t i o n 1  1  (1 . )(1 .)f Ui ,i? ) 2  7jt  minimum  1/ 3.=  (1-11) can be w r i t t e n ,  (1.)(l.)f(^ ,T? )  [k] =  3  I t can be shown (Cooke 1975) t h a t W f o r  each Gauss p o i n t i s 1.0 and the Gauss p o i n t v a l u e s a r e +  where  number  o b t a i n the area of the p a r t i c u l a r  T j , hence a 2 by 2 Gauss r u l e ( f o u r p o i n t s ) i s  that  the  = .57735..., and  + (l.)(1.)fU ,7?,) 2  (1.)(1-.)fU ,»?2)  +  2  £, 2  TJ = 2  -.57735...  +  (1-14)  

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-0062733/manifest

Comment

Related Items