Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

A simple anisotropic bounding surface plasticity model for cyclic response of clays Seidalinov, Gaziz 2012

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

Item Metadata

Download

Media
24-ubc_2012_fall_seidalinov_gaziz.pdf [ 3.93MB ]
Metadata
JSON: 24-1.0073338.json
JSON-LD: 24-1.0073338-ld.json
RDF/XML (Pretty): 24-1.0073338-rdf.xml
RDF/JSON: 24-1.0073338-rdf.json
Turtle: 24-1.0073338-turtle.txt
N-Triples: 24-1.0073338-rdf-ntriples.txt
Original Record: 24-1.0073338-source.json
Full Text
24-1.0073338-fulltext.txt
Citation
24-1.0073338.ris

Full Text

A SIMPLE ANISOTROPIC BOUNDING SURFACE PLASTICITY MODEL FOR CYCLIC RESPONSE OF CLAYS by Gaziz Seidalinov  A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in The Faculty of Graduate Studies (Civil Engineering)  THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) November 2012 c Gaziz Seidalinov 2012 ⃝  Abstract Natural clays are anisotropic in their in-situ state and have an undisturbed shear strength in excess of the remoulded strength. In addition, most of the structures founded on clay deposits must be designed to withstand cyclic loads such as seismic ground motions or ocean waves. When subjected to earthquake or wind induced cyclic loadings, clay exhibits complex response. A realistic modeling of clay response under irregular cyclic loading requires an appropriate description of the stress–strain relationship. This thesis extends the formulation of a Simple ANIsotropic CLAY plasticity (SANICLAY) model by incorporation of a bounding surface formulation for successful simulations of clay response under cyclic loading. The most important aspects of the proposed formulation are the position of a projection center and the ability to capture continuous stiffness degradation. The proposed projection center is established in the instant of any stress reversal and it realistically reflects the experimentally observed plastic strains. A damage parameter is also adopted to better simulate the continuous stiffness reduction during the course of applied cyclic loading. The proposed model is developed with the aim of maintaining the simplicity and yet including an adequate level of sophistication for successful modeling of the key features of clay response. The model formulation is presented in detail followed by its qualitative and quantitative comparison with experimentally observed clay response. The presented model validation demonstrates the capabilities of the model in capturing a number of important characteristic features of clay response in both monotonic and cyclic loadings. Limitations and recommendations for future work are discussed as well.  ii  Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  ii  Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  iii  List of Tables  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  iv  List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  v  List of Symbols  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  vi  Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  x  1 Introduction  1  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  2 Experimental and Theoretical Modeling of Clay Response 2.1  2.2  . . . . . . . .  4  . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  4  2.1.1  Discussion on fast and slow undrained cyclic loadings . . . . . . . . .  4  2.1.2  Bacchus (1969) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  5  2.1.3  Sangrey et al. (1969) . . . . . . . . . . . . . . . . . . . . . . . . . . .  6  2.1.4  Takahashi et al. (1980)  7  2.1.5  Zergoun (1991); Zergoun and Vaid (1994)  Experimental modeling  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  8  Theoretical modeling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  10  2.2.1  Classical elastoplastic models . . . . . . . . . . . . . . . . . . . . . .  11  2.2.2  Plastic strains inside the yield surface  . . . . . . . . . . . . . . . . .  13  2.2.3  Development of the bounding surface formulation . . . . . . . . . . .  14  2.2.4  Bounding surface models for monotonic response of clays . . . . . . .  22  2.2.5  Bounding surface models for cyclic response of clays  26  . . . . . . . . .  iii  Table of Contents 2.3  . . . . . . . . . . . . . . . . . . . . . .  32  3 Development of the Proposed Model . . . . . . . . . . . . . . . . . . . . . .  34  3.1  3.2  3.3  Missing features in reviewed models  Reference Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  34  3.1.1  Formulation in the triaxial space . . . . . . . . . . . . . . . . . . . .  34  3.1.2  Formulation in the multiaxial space  . . . . . . . . . . . . . . . . . .  38  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  42  3.2.1  Formulation in the triaxial space . . . . . . . . . . . . . . . . . . . .  42  3.2.2  Formulation in the multiaxial space  . . . . . . . . . . . . . . . . . .  50  Sensitivity analysis on new model constants . . . . . . . . . . . . . . . . . .  55  3.3.1  Parameter h0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  56  3.3.2  Parameter ad . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  58  Proposed Model  4 Qualitative Model Performance  . . . . . . . . . . . . . . . . . . . . . . . . .  4.1  The advantage of the proposed projection center  4.2  Undrained cyclic triaxial loading  4.3  Circular stress path in the deviatoric plane  4.4  61  . . . . . . . . . . . . . . .  61  . . . . . . . . . . . . . . . . . . . . . . . .  65  . . . . . . . . . . . . . . . . . .  67  4.3.1  Response to cyclic circular stress path loading . . . . . . . . . . . . .  68  4.3.2  Effect of parameter C in circular stress path loading  71  . . . . . . . . .  Generation of secant shear modulus and damping ratio variation with shear strain amplitude  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  4.5  Generation of cyclic stress ratio versus number of cycles  . . . . . . . . . . .  4.6  Generation of peak axial strain and residual pore water pressure versus number of cycles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  5 Quantitative Model Performance  . . . . . . . . . . . . . . . . . . . . . . . .  72 77  80 85  5.1  Experiments required for model calibration  . . . . . . . . . . . . . . . . . .  85  5.2  Reconstituted Georgia kaolin clay . . . . . . . . . . . . . . . . . . . . . . . .  87  5.2.1  Model calibration  . . . . . . . . . . . . . . . . . . . . . . . . . . . .  88  5.2.2  Model performance  . . . . . . . . . . . . . . . . . . . . . . . . . . .  92  Structured Cloverdale clay  . . . . . . . . . . . . . . . . . . . . . . . . . . .  98  5.3  iv  Table of Contents  5.4  5.3.1  Model calibration  5.3.2  Model performance  . . . . . . . . . . . . . . . . . . . . . . . . . . . .  98  . . . . . . . . . . . . . . . . . . . . . . . . . . . 102  Reconstituted Ariake clay (fast cyclic loading) . . . . . . . . . . . . . . . . . 111 5.4.1  Model calibration  5.4.2  Model performance  . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 . . . . . . . . . . . . . . . . . . . . . . . . . . . 112  6 Summary, Conclusions, and Recommendations for Future Research  . . 115  6.1  Summary and conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115  6.2  Recommendations for further research  . . . . . . . . . . . . . . . . . . . . . 117  Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118  v  List of Tables 3.1  Model parameters for sensitivity analysis . . . . . . . . . . . . . . . . . . . .  56  3.2  Initial internal variables for sensitivity analysis . . . . . . . . . . . . . . . . .  56  4.1  Summary of model parameters and initial internal variables for simulations of undrained cyclic triaxial loading . . . . . . . . . . . . . . . . . . . . . . . . .  4.2  65  Summary of model parameters for simulations of circular stress path in the deviatoric plane . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  68  5.1  Model constants for the proposed SANICLAY-D-B model . . . . . . . . . . .  87  5.2  Initial internal variables for the proposed SANICLAY-D-B model . . . . . .  87  5.3  Model parameters for Georgia kaolin Clay . . . . . . . . . . . . . . . . . . .  91  5.4  Initial internal variables for Georgia kaolin clay . . . . . . . . . . . . . . . .  92  5.5  Model parameters for Cloverdale clay . . . . . . . . . . . . . . . . . . . . . . 103  5.6  Initial internal variables for Cloverdale clay . . . . . . . . . . . . . . . . . . . 103  5.7  Model parameters for Ariake clay . . . . . . . . . . . . . . . . . . . . . . . . 114  5.8  Initial internal variables for Ariake clay . . . . . . . . . . . . . . . . . . . . . 114  vi  List of Figures 2.1  Comparison of accumulation of pore water pressure in fast (dash lines) and slow (solid lines) undrained cyclic triaxial loading tests with drainage allowed at the end of loading for a range of applied stress amplitudes, τcy /2Suc , from c Zergoun and Vaid (1994) ⃝Canadian Science Publishing or its licensors . . .  10  2.2  Radial mapping rule from Dafalias (1986b) . . . . . . . . . . . . . . . . . . .  18  2.3  Radial mapping rule from Dafalias and Herrmann (1986) . . . . . . . . . . .  19  2.4  Radial mapping rule from Anandarajah and Dafalias (1986) . . . . . . . . .  21  3.1  Reference model schematic diagram in the triaxial plane . . . . . . . . . . .  36  3.2  Proposed model schematic diagram in the triaxial plane . . . . . . . . . . . .  44  3.3  Proposed projection center . . . . . . . . . . . . . . . . . . . . . . . . . . . .  47  3.4  Multiaxial representation of: (a–b) the bounding surface with N -cone, and (c–d) plastic potential with critical state failure envelope . . . . . . . . . . .  52  3.5  Parameter h0 sensitivity in undrained monotonic triaxial loading . . . . . . .  57  3.6  Parameter h0 sensitivity in undrained cyclic triaxial loading . . . . . . . . .  59  3.7  Parameter ad sensitivity in undrained cyclic triaxial loading . . . . . . . . .  60  4.1  Typical stress–strain path in undrained cyclic triaxial loading test after Sheu (1984) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  62  4.2  Projection center on rotational variable α in the triaxial plane . . . . . . . .  63  4.3  Stress–strain path simulations in undrained cyclic triaxial loading with the projection center on rotational variable α: (a) Cp = 0, and (b) Cp = 0.6 . . .  4.4  64  Stress–strain path simulations in undrained cyclic triaxial loading with the projection center proposed in this thesis . . . . . . . . . . . . . . . . . . . .  64  vii  List of Figures 4.5  SANICLAY family models performance in undrained cyclic triaxial loading: (a–b) SANICLAY, and (c–d) SANICLAY-D . . . . . . . . . . . . . . . . . .  66  4.6  Comparison of model performance with bounding surface effect . . . . . . .  67  4.7  Isotropic unloading, constant-p triaxial compression, and circular stress path: (a) deviatoric view, (b) shifted view . . . . . . . . . . . . . . . . . . . . . . .  4.8  69  Observed strains in the strain space (deviatoric view on the left, shifted view on the right) induced by circular stress path (clockwise) in the deviatoric plane 70  4.9  Effect of parameter C on induced strain path during circular stress path loading 71  4.10 Definition of equivalent shear modulus and the quantities used for the determination of damping ratio . . . . . . . . . . . . . . . . . . . . . . . . . . . .  73  4.11 Experimental data after Vucetic and Dobry (1988): (a) degradation of Gsec and (b) variation of damping ratio with γ  . . . . . . . . . . . . . . . . . . .  74  4.12 Simulations of Gsec /Gmax degradation with γ: (a) ad = 0 and (b) ad = 60 . .  75  4.13 Simulations of damping ratio variation with γ: (a) ad = 0 and (b) ad = 60 .  77  4.14 Experimental data after Zergoun and Vaid (1994) on cyclic stress ratio CSR required to reach εa = 3% in a number of cycles . . . . . . . . . . . . . . . .  78  4.15 Generation of CSR required to reach εa = 3 % . . . . . . . . . . . . . . . . .  79  4.16 Development of: (a) peak axial strain εa , and (b) residual pore water pressure u versus number of cycles, experimental data from Zergoun and Vaid (1994)  81  4.17 Generation of: (a) peak axial strain εa , and (b) residual pore water pressure u versus number of cycles with different values of h0 and fixed ad = 40 . . . .  83  4.18 Generation of: (a) peak axial strain εa , and (b) residual pore water pressure u versus number of cycles with fixed h0 = 50 and varied values of ad . . . . . 5.1  Calibration of Cc and Cr in e– log p plane based on isotropic consolidation test after Sheu (1984) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  5.2  84  88  Calibration of Mc and Me in the triaxial plane based on undrained triaxial compression and extension tests starting from normally consolidated state, experimental results after Sheu (1984) . . . . . . . . . . . . . . . . . . . . .  89  viii  List of Figures 5.3  Calibration of C based on undrained triaxial compression test on normally consolidated clay, experimental results after Sheu (1984) . . . . . . . . . . .  5.4  91  Initial hardening parameter calibrated to h0 = 50 based on undrained cyclic triaxial loading, experimental results after Sheu (1984) . . . . . . . . . . . .  92  5.5  Simulation of isotropic consolidation test, experimental results after Sheu (1984) 93  5.6  Simulation of CIUC and CIUE tests on normally consolidated clay, experimental results after Sheu (1984) . . . . . . . . . . . . . . . . . . . . . . . . .  5.7  Undrained cyclic triaxial loading stress paths with the digitized experimental results after Sheu (1984) (left) and model simulations (right) . . . . . . . . .  5.8  96  Undrained cyclic triaxial loading stress–strain paths with the digitized experimental results after Sheu (1984) (left) and model simulations (right) . . . .  5.9  94  97  Calibration of Cc and Cr in e– log p plane based on oedometer test after Zergoun (1991) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .  99  5.10 Calibration of ϕc in the triaxial plane based on undrained triaxial compression test on normally consolidated and overconsolidated specimens, experimental results after Zergoun (1991) . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 5.11 Calibration of: (a) Si , and (b) ki based on oedometer test after Zergoun (1991)101 5.12 Calibration of C based on CIUC test on initially normally consolidated clay, experimental data after Zergoun (1991) . . . . . . . . . . . . . . . . . . . . . 102 5.13 Simulation of oedometer test on Cloverdale clay, experimental results after Zergoun (1991) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104 5.14 Simulations of CIUC and CIUE tests on normally and overconsolidated clay, experimental results after Zergoun (1991) . . . . . . . . . . . . . . . . . . . . 105 5.15 Undrained cyclic triaxial loading stress path with the digitized experimental results after Zergoun (1991); Zergoun and Vaid (1994) (left), and model simulations (right) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106 5.16 Undrained cyclic triaxial loading stress–strain path with the digitized experimental results after Zergoun (1991); Zergoun and Vaid (1994) (left), and model simulations (right) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107  ix  List of Figures 5.17 Simulations of peak εa versus number of cycles in undrained cyclic triaxial loading for the range of available cyclic stress amplitudes . . . . . . . . . . . 110 5.18 Simulations of peak εa versus number of cycles in undrained cyclic triaxial loading for the range of available cyclic stress amplitudes with the original SANICLAY model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111 5.19 Calibration of h0 and ad based on peak εa versus number of cycles in undrained cyclic triaxial loading, experimental results after Yasuhara et al. (1992) . . . 113 5.20 Simulation of cyclic stress ratio CSR versus number of cycles required to reach εDA = 5% . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114  x  List of Symbols List of symbols is applicable for Chapters 3–6. ad  rate of d evolution  A  destructuration distributing model parameter  b  similarity ratio of f = 0 to F = 0  C  rate of α evolution  Cc  compression line slope in e– log p plane  Cp  model parameter controlling position of projection center  Cr  unloading–reloading line slope in e– log p plane  Ce  fourth-order elastic modulus tensor  Cep  fourth-order elastoplastic modulus tensor  d  damage parameter  e  void ratio  ein  initial e  e  deviatoric strain tensor  ep  plastic deviatoric strain tensor  f =0  loading surface  F =0  bounding surface  G=0  plastic potential  G  shear modulus  Gs  specific gravity of soil particles  Gsec  secant shear modulus  Gmax  maximum tangential shear modulus  h  bounding surface hardening function  xi  List of Symbols h0  initial value of h  H  Heaviside step function  I  identity tensor  ki  model parameter, controlling rate of destructuration  K  bulk modulus  Kp  plastic modulus  ¯p K  plastic modulus for stress states on bounding surface  K0  earth pressure coefficient at rest  L  loading index  M  critical state ratio in p–q plane  Mc  M in compression  Me  M in extension  n  unit tensor along r − α  nc  unit tensor along rc − α  N  peak stress ratio on bounding surface  p  mean effective stress  p¯  image p on bounding surface  pa  atmospheric pressure  pc  p at projection center  p0  size of bounding surface along p-axis  p0,d  destructured value of p0  pα  size of plastic potential along p-axis  q  deviatoric stress  q¯  image q on bounding surface  qb  closest qc projection onto bounding surface  qc  q at projection center  qcy  amplitude of uniform cyclic loading  qn  set of n internal variables  qα  qc projection on α  r  distance from projection center to image stress xii  List of Symbols r  back stress-ratio tensor  s  elastic nucleus factor  s  deviatoric stress tensor  sc  s at projection center  ¯s  s on bounding surface  sb  projection of sc on bounding surface along nc  Si  initial structuration factor  sα  projection of sc on α  Su  undrained shear strength  Suc  undrained shear strength in compression  u  pore water pressure  x  model parameter required for α saturation  α  deviatoric back-stress ratio  α  deviatoric back-stress ratio tensor  εa  axial strain  εpd  destructuration plastic strain  εDA  double amplitude axial strain  εpq  plastic deviatoric strain  εr  radial strain  εpv  plastic volumetric strain  ε  strain tensor  εp  plastic strain tensor  δ  distance from current stress to image stress  δ  δ generalized to multiaxial stress space  η  stress ratio  ηK0  η for K0 -loading  Θ  Lode angle  γ  shear strain  κ  slope of unloading–reloading line in e– ln p plane  λ  slope of compression line in e– ln p plane xiii  List of Symbols ν  Poisson’s ratio  σa  axial stress  σe  equivalent σc  σr  radial stress  σc  consolidation stress  σ  stress tensor  ¯ σ  σ for stress states on bounding surface  σ1  major principal stress  σ2  intermediate principal stress  σ3  minor principal stress  τ  shear stress  τcy  shear stress amplitude of uniform cyclic loading  ϕ  friction angle  ϕc  friction angle in compression  ∇F  bounding surface gradient  ∇f  loading surface gradient  ⟨⟩  Macauley brackets  ˙  time derivative or rate  ¯  association with stress state on bounding surface  :  inner product of two tensors  ⊗  outer product of two tensors  tr  trace operator  xiv  Acknowledgements I offer my sincerest gratitude to my supervisor, Dr. Mahdi Taiebat, who has supported me with his guidance throughout my research. His encouragement and high level of competence inspired me to explore the subject. I am truly indebted and thankful to him for the time he devoted to discussions on the subject matter and for his moral support. His sometimes critical comments helped me to become a better researcher. I attribute the level of this thesis to his effort. I would also like to thank Professor Yannis Dafalias for his help in formulating the motion of the projection center and for reviewing this thesis. Thanks to the Department of Civil Engineering faculty and staff members, the University of British Columbia, especially to the library staff for the resources provided. I would also like to thank all the people who have become my friends during my studies at UBC. I would like to thank my wife, Ainur, for her patience, love, and constant support. Thanks will not be enough for my parents, Ilyas and Karlygash, from whom I have been constantly away. Special thanks to my lovely brother, Arystan, who has been taking care of them and tried to fill my missing spot. Financial support provided by Bolashak International Scholarship of the President of the Republic of Kazakhstan and partial financial support from the Natural Sciences and Engineering Research Council of Canada (NSERC) are gratefully acknowledged.  xv  1. Introduction All stresses considered in this thesis are effective stresses. A constitutive model is a mathematical description of the stress–strain relationship for a material. The complexity of a constitutive model depends on the type of material to be simulated. Unlike man-made materials (e.g., concrete, steel), soil behavior is governed by the presence of water in addition to solid phase of soil particles. The complex nature of clays cannot be simulated satisfactory by simple models (e.g., von Mises, Mohr–Coulomb), as their behavior is known to be dependent on confining pressure and overconsolidation ratio (OCR), and they exhibit a significant degree of anisotropy in the in-situ state. The dependence on confining pressure and to some extent on OCR can be simulated by the Modified Cam-Clay model (MCC) (Roscoe and Burland, 1966), which is based on Critical State Soil Mechanics (CSSM) theory developed at Cambridge University (Roscoe and Poorooshasb, 1963; Roscoe et al., 1963, 1958). The Simple ANIsotropic CLAY plasticity (SANICLAY) model developed by Dafalias (1986a), Dafalias et al. (2006) allows successive simulations of induced anisotropy by incorporation of rotational and distortional hardening of a so-called yield surface, which encloses elastic domain. Destructuration theory incorporated into SANICLAY is a further model extension (SANICLAY-D), proposed by Taiebat, Dafalias and Peek (2010), enabling one to simulate structure loss in sensitive clays. Experimental evidence on clays suggests that inelastic strains occur within the yield surface. This issue becomes fundamentally important for the loadings that involve stress paths inside the yield surface. Undrained cyclic loading on clays is accompanied by the reduction of mean stress, resulting in a stress path mainly contained within the yield surface. CSSM theory suggests that the reduction of mean stress is associated with plastic strains. As the classical elastoplastic models assume purely elastic strains for the stress states inside the yield surface, no change in mean effective stress can be simulated. A mechanism allowing  1  Chapter 1. Introduction reduction of mean stress for undrained cyclic loading is required. Iwan (1967) and Mr´oz (1967, 1969) independently, introduced similar models with several surfaces, enclosing each other, of progressively larger size. As the stress point reaches a larger surface, the stiffness reduces, resulting in increased hysteretic work. A bounding surface framework, first presented for metals by Dafalias and Popov (1975) and then incorporated for soils by Dafalias (1986b), Dafalias and Herrmann (1986), and Anandarajah and Dafalias (1986), is another approach to simulate plastic deformations within the yield surface. Unlike in the models by Iwan and Mr´oz, stiffness change within the yield surface is continuous in the bounding surface formulation. In addition, the formulation does not require storing in memory a finite number of surfaces for better simulation of stiffness reduction for the stress states inside the yield surface. The bounding surface framework attracted the attention of researchers owing to its simplicity and has been incorporated into different models to simulate soil response. In the bounding surface formulation, the yield surface is replaced by a bounding surface, which plays a similar role as the yield surface in that it represents consolidation history. The bounding surface allows a certain degree of plastic deformation for the stress states inside it. The model proposed in this thesis (SANICLAY-D-B) incorporates the bounding surface formulation into the SANICLAY-D model, which is chosen to account only for isotropic destructuration sufficient to simulate structure loss. The SANICLAY-D-B model is a continuation of the work presented by Taiebat, Dafalias and Kaynia (2010). The proposed model includes a modified projection center to avoid overdamping and a damage parameter for the reduction of stiffness required to simulate the cyclic softening observed in experiments. The plastic modulus within the bounding surface depends on the relative distance between the current stress point and the image point, mapped from a projection center through the stress point on to the bounding surface. In addition to the distance dependence, the plastic modulus is also controlled by a hardening parameter, either constant or function. In the proposed formulation, the hardening parameter is a function of a newly introduced damage parameter. The projection center is located at the point of stress reversal at the instant of the reversal and further evolves to maintain its relative position within the evolving bounding surface. In general, the projection center can be anywhere within the bounding surface but never outside it. The number of constants that should be calibrated is 10 for reconstituted 2  Chapter 1. Introduction and 11 for structured clays, including two new model parameters related to the bounding surface. The thesis is presented in a manner that allows the readers, who are not familiar with the subject, to understand the proposed model formulation, its advantages and predictive capabilities. All stresses considered are effective stresses which govern the response of soil. In Chapter 2, a review of experimentally observed clay response in undrained cyclic and theoretical modeling of inelastic strains within the yield surface is covered. Theoretical modeling part also includes a brief introduction to classical elastoplastic models. The proposed SANICLAY-D-B formulation followed by a sensitivity analysis of new model parameters is presented to the reader in Chapter 3. A qualitative model performance is demonstrated in Chapter 4, consisting of: • a comparison of different projection center approaches reported in the literature and proposed for the current formulation, • a comparison of the proposed model with the previous SANICLAY family models in undrained cyclic triaxial loading and circular stress path in the deviatoric plane, • a model simulation of secant shear modulus and damping ratio variation with shear strain amplitude, • a model simulation of the cyclic stress ratio required to reach a certain strain level in a number of cycles, • a model simulation of attained axial strains and residual pore water pressure in a number of cycles. In Chapter 5, model calibration and quantitative performance are illustrated for three clays: Georgia kaolin clay (reconstituted), Cloverdale clay (structured) in undrained slow cyclic triaxial loading, and Ariake clay (reconstituted) in fast undrained cyclic triaxial loading. Finally, a summary, conclusions, and recommendations for future research are presented in Chapter 6.  3  2. Experimental and Theoretical Modeling of Clay Response In this chapter, a relatively detailed literature review is conducted on the experimentally observed response of clays in cyclic loading and a number of the proposed constitutive models for simulating that response. A bounding surface framework incorporated into the reviewed models appears to be an attractive approach for simulations of undrained cyclic loading. Finally, the features that need to be reflected in cyclic loading simulations are discussed.  2.1 2.1.1  Experimental modeling Discussion on fast and slow undrained cyclic loadings  Imposed undrained cyclic loading on clays can be divided into two categories: slow and fast undrained cyclic loadings. Slow undrained cyclic loading tests are those in which the frequency of applied cyclic loading is slow enough to achieve equalization of pore water pressure in the whole specimen; otherwise, the tests are referred to as fast undrained cyclic loading ones in this study. It is important to understand advantages and disadvantages of both. The review of the experimental data is conducted on the tests performed in triaxial cell. During triaxial testing pore water pressure is usually measured on top and bottom platens of tested specimens. A measured pore water pressure is then deducted from applied total stresses to come up with effective stresses, governing the response of a tested material. In this way, positive pore water pressure is associated with reduction of effective stresses, while negative pore water pressure with increase of effective stresses. An element test requires the measured data on specimen boundaries to be representative of the entire specimen. The main advantage of slow undrained cyclic loading is that equalization of pore water pressure 4  2.1. Experimental modeling throughout the specimen is allowed, which ensures measurements of approximately average pressure in the entire specimen. Thus, interpreted effective stresses from slow undrained cyclic loading are more reliable. This is especially important for clays, where dissipation of pore water pressure requires a considerable amount of time. Fast cyclic loading tests miss this important feature and become less reliable in terms of representative effective stress measurements for the soil element. The rate effect of applied loading also plays an important role. The rate effect becomes more pronounced with increased plasticity of soil. Many investigators rely on the supposition that measured response is a result of the rate of applied loading only; however, in fast undrained cyclic loading, measured response is also a result of measuring un-equalized pore water pressure on the boundaries of the specimen, which is not representative of the material only. Of course, strain rate effect is present, but it is very difficult or even impossible to differentiate soil response due to the rate effect from the measured response due to the pore water pressure inequalities within the specimen. Fast undrained cyclic loading tests invoke much higher measured strength than that in slow undrained cyclic loading. It is true that the increase of strength is due to the rate effect, but unreliable pore water pressure measurements are likely to erroneously increase measured strength as well. Thus, the strength in fast undrained cyclic loading is overestimated, rendering results that cannot be relied on. Although slow undrained cyclic loading does not account for the rate effect and therefore underestimates the strength of the soil, it gives conservative results. The emphasis is placed on slow undrained cyclic triaxial loading tests throughout the discussion of available experimental data.  2.1.2  Bacchus (1969)  Undrained cyclic strain-controlled triaxial tests with constant strain amplitude were performed by Bacchus (1969) on commercially available Webb Mark 4 clay with plasticity index, PI = 26 %. Tested specimens consist of 28% of clay. Bacchus (1969) was well aware of a need of pore water pressure equalization within cyclically loaded specimen; however, he adopted 0.2 Hz frequency as the appropriate one. From reported results on undrained cyclic triaxial loading on normally consolidated spec5  2.1. Experimental modeling imens the main findings could be summarized as follows: • in all applied cyclic strain levels, during the first initial cycles there was a decrease in mean stress p; the rate of reduction of p decreased with cyclic loading; • larger strain amplitude tests involved larger shear stresses to reach desired strains during first few cycles; however, with an increasing number of cycles, lower shear stresses were attained for larger strain amplitudes in a given number of cycles; this was an indication of more rapid stiffness loss for high strain amplitude tests; • small strain amplitudes stress–strain loops exhibited a subtle change in shape, and stress–strain loops became “S-shaped” with an increasing number of cycles; this behavior might be attributed to the presence of a significant amount of granular material. Undrained cyclic triaxial tests were also performed on samples isotropically consolidated to a range of preconsolidation stresses and with different overconsolidation ratios (OCRs). Specimens with larger OCRs exhibited an initial increase in p and reduction with further cycles. The extent of p increase, followed by its decrease, was larger for higher OCRs. In general, stress–strain loops for overconsolidated specimens were similar to those obtained for normally consolidated specimens with the difference that less reduction in stresses and energy was observed during initial cycles.  2.1.3  Sangrey et al. (1969)  Sangrey et al. (1969) paid a close attention to reliable pore water pressure measurements in undrained cyclic triaxial tests. They performed a series of one-way cyclic loading tests with constant stress amplitude on undisturbed Newfield clay. To ensure pore water pressure equalization within the specimens all tests were conducted with an axial strain rate of about 0.0002 %/min. Their main finding was that for each OCR there existed a “critical cyclic stress level” (CCSL). The tests with cyclic stress amplitudes below CCSL attained equilibrium, and resultant strains were finite, whereas applied loading with stress amplitude above CCSL caused ultimate collapse with infinite axial strains. While normally consolidated specimens 6  2.1. Experimental modeling subjected to undrained cyclic triaxial loading exhibited stress path migration towards the origin until the equilibrium was reached, overconsolidated samples showed an increase of p within the initial cycle. The increase of p was then followed by its reduction that was smaller for higher OCRs.  2.1.4  Takahashi et al. (1980)  Takahashi et al. (1980) discussed observations from undrained cyclic triaxial tests with different frequencies on a low-plasticity sandy clay, Lower Cromer Till. The loading was applied in constant stress amplitude tests with sinusoidal waveform. In addition to conventional pore water pressure transducers, piezometer probes at mid-height of tested specimens were used. The period of applied cyclic loading was restricted to a minimum of 10 min (of about 0.0017 Hz maximum frequency), ensuring reliable measurements of pore water pressure at the mid-height of tested specimens. The following observations were made in examination of undrained cyclic triaxial tests with 30-min period and applied magnitude of (σa − σr ) /2Suc = 0.75, where Suc is an undrained shear strength in compression: • a migration of the stress path towards the origin induced by the rise of pore water pressure within a constant volume specimen; • the stress path changed its shape as it approached the critical stress ratio; loops developed in extension owing to a more dialative response than in compression; • as the stress path reached the critical stress ratio, axial strains development accelerated and larger magnitude of axial strains was observed in extension; • the stress path traveled beyond the critical stress ratio in extension as the rate of applied loading was increased; Takahashi et al. (1980) further aimed to investigate what influenced the rate of stress path migration towards the origin. From their findings they concluded that the main factors were: cyclic stress ratio (CSR), consolidation history, and rate of applied loading. A larger CSR initially induced more rapid reduction of p. However, once the critical stress ratio was 7  2.1. Experimental modeling reached, p reduction occurred at a decreasing rate. A lower CSR, apparently, approached the critical stress ratio at lower p. As a result, overall reduction of p for lower CSRs was larger once the stress path was stabilized. Specimens sheared from high OCRs exhibited an increase in p, i.e., a dialative response during initial cycles, which was then followed by reduction of p with the magnitude of reduction inversely proportional to OCR. Takahashi et al. (1980) indicated that direct evidence on the effect of cyclic frequency was not available as the results were obtained by measuring pore water pressure at the boundaries; however, it is worth to elaborate on reported results of observed stress paths with different cyclic loading frequencies. In general, reduction of p during cyclic loading occurred to a lesser degree as loading frequency was increased. The first series of tests were performed with cycle periods of 100 and 900 s, even though Takahashi et al. (1980) mentioned that the minimum cyclic period should be 10 min to ensure reliable pore water pressure measurements. Thus, the reported low reduction of p for smaller periods was partly due to not stabilized pore water pressures within the tested specimens. In the second series, soil response during cyclic tests with periods of 10, 30, 240, and 480 min were studied. The results clearly showed that as the loading period was increased the reduction of p was larger. Takahashi et al. (1980) ensured that pore water pressures at the peripheries of the specimens were measured accurately and radial non-uniformities were unlikely. They admitted that there might had been local irregularities of pore water pressure within the specimens.  2.1.5  Zergoun (1991); Zergoun and Vaid (1994)  Zergoun (1991), Zergoun and Vaid (1994) justified the need of slow undrained cyclic triaxial loading for clays. They mentioned that non-uniformities within a specimen were dependent on the degree of end restraints and the rate of loading. As the center of specimen was believed to be the most representative of the whole specimen, the only way to obtain reliable pore water pressure measurements was to perform slow undrained cyclic triaxial loading tests with 0.05%/hr for strain-controlled and 60 kPa/hr for stress-controlled tests. Zergoun (1991), Zergoun and Vaid (1994) compared slow and fast undrained cyclic triaxial loading tests with constant stress amplitude on structured Cloverdale clay with PI = 24%. When comparing results for a particular stress amplitude in slow and fast cyclic loading the 8  2.1. Experimental modeling following observations were made: • the effect of stress amplitude was similar; • a much smaller number of cycles was required to induce a particular peak axial strain in slow cyclic loading; this was accompanied by more extensive p reduction of the stress path; • axial strains accumulation accelerated once the stress path nearly reached equilibrium for slow cyclic loading; in fast cyclic loading such acceleration occurred well before the equilibrium was attained; • larger values of pore water pressure were measured during slow cyclic loading; in Fig. 2.1 accumulation of pore water pressure reached a plateau in slow cyclic loading; referring to Fig. 2.1, fast cyclic loadings were accompanied by a pore water pressure increase during the rest periods at which applied stresses were held constant, while in slow cyclic loadings pore water pressure was unchanged; this indicated that equilibrium of pore water pressure was not achieved during fast cyclic loading and it did so during the rest periods involving further deformations; even though pore water pressure increased in rest periods after fast cyclic loading, residual pore water pressure was below that measured in slow cyclic loading for a given stress amplitude. From their findings Zergoun (1991), Zergoun and Vaid (1994) confirmed that during fast undrained cyclic loading tests the effective stresses could not be adequately measured. In the tests with different stress amplitudes, Zergoun (1991), Zergoun and Vaid (1994) found a threshold stress amplitude. Axial strains decelerated if the applied stress amplitude was below the threshold and accelerate if it was above. Larger axial strains developed in extension for all stress amplitudes with a decreasing difference between strains in compression and extension as the amplitude of applied loading was decreased. For stress amplitudes above the threshold, the following three stages of deformations were traced: • primary stage: axial strains develop at a decreasing rate; • secondary stage: strains develop at an approximately constant rate; 9  2.1. Experimental modeling  Figure 2.1: Comparison of accumulation of pore water pressure in fast (dash lines) and slow (solid lines) undrained cyclic triaxial loading tests with drainage allowed at the end of loading for a range of applied stress amplitudes, τcy /2Suc , from Zergoun and Vaid (1994) c ⃝Canadian Science Publishing or its licensors • tertiary stage: strains develop at an increasing rate. For very high stress amplitudes the secondary stage was negligible, and the decreasing rate of strains during the primary stage was followed by an increasing rate in the tertiary stage. Two more important observations should be considered. The initial direction of applied cyclic loading played a significant role. Two undrained tests with the same stress amplitude, one starting from compression and the other from extension, were compared. The cyclic test starting from extension exhibited more severe deformations than the test starting from compression. This implied the importance of an anisotropy for tested clay. The effect of a step increase in amplitude of cyclic loading was studied as well. In 10 cycles, the amplitude of loading was increased from τcy /Suc = 0.63 to 0.67, where τcy is an amplitude of applied shear stress. The induced strains, residual strains, and pore water pressure in a number of cycles simply corresponded to transformation of the response to the one that would had been obtained for τcy /Suc = 0.67 with no change of amplitude of applied loading.  10  2.2. Theoretical modeling  2.2  Theoretical modeling  Before going into detail of theoretical modeling of undrained cyclic loading, a review of general framework used in classical elastoplastic models is proposed to the reader.  2.2.1  Classical elastoplastic models  Classical elastoplastic models assume the existence of a yield surface in a six dimensional stress space. Its general equation can be expressed as f (σ, qn ) = 0  (2.1)  where σ is a stress tensor; and qn are n internal variables controlling the geometry of the yield surface. This yield surface encloses the region inside which any stress increment that does not cross the yield surface induces only elastic strains. Once a stress increment crosses the yield surface a plastic strain increment is induced in addition to the elastic one. A basic kinematic assumption is a decomposition of a total strain tensor ε into elastic εe and plastic εp parts. The rate of ε can be decomposed accordingly: ε˙ = ε˙ e + ε˙ p  (2.2)  where a superposed dot stands for the material time derivative, i.e., rate. An input into a constitutive model could be either a strain increment or stress increment. General constitutive relations between stress increment and elastic strain increment are then given for the corresponding inputs, respectively, σ˙ = Ee : ε˙e ;  (2.3)  ε˙ e = Ce : σ˙  (2.4)  where operator : is an inner product of two adjacent tensors; Ee and Ce are fourth order tensors of elastic moduli and compliances, respectively. Rates of plastic strain and internal variables are assumed to be given by Eqs. 2.5 and 11  2.2. Theoretical modeling 2.6. ε˙ p = ⟨L⟩Np  (2.5)  q˙ n = ⟨L⟩¯ qn  (2.6)  where L is a plastic scalar-valued multiplier, or loading index to be defined in the sequel of this section; ⟨⟩ are Macauley brackets such that ⟨L⟩ = L if L > 0, and ⟨L⟩ = 0 if L ≤ 0; Np ¯ n are tensorial and scalar/tensorial, respectively, functions of state variables σ and qn . and q It is then common to define Np = ∂g/∂σ, where ∂g/∂σ is a gradient of a plastic potential given by: g (σ, qn ) = 0  (2.7)  An associative flow rule considers the yield surface and plastic potential coinciding with each other, while a non-associative flow rule assumes non-coincident yield surface and plastic potential. Then Eqs. 2.3 and 2.4 are re-written on the basis of Eqs. 2.2 and 2.5: ( ) ∂g σ˙ = E : ε˙ − ⟨L⟩ ∂σ e  ε˙ = Ce : σ˙ + ⟨L⟩  ∂g ∂σ  (2.8)  (2.9)  The consistency condition requires that the stress point cannot lie outside the domain enclosed by the yield surface. Thus, once a stress increment crosses the yield surface, internal variables are updated to maintain an updated stress point on the yield surface. The consistency condition is expressed in mathematical form as ∂f ∂f f˙ = : σ˙ + : q˙ n ∂σ ∂qn  (2.10)  Substituting Eq. 2.6 into Eq. 2.10 gives the loading index 1 L= Kp  (  ∂f : σ˙ ∂σ  ) (2.11)  12  2.2. Theoretical modeling ¯ n is a plastic modulus. where Kp = − (∂f /∂qn ) : q The constitutive relationships in Eqs. 2.8 and 2.9 can be re-written on the basis of Eq. 2.11 so that the the constitutive matrices are: [  ] (Ee : ∂g/∂σ) ⊗ (∂f /∂σ) : Ee σ˙ = E ε˙ = E − H (L) : ε˙ Kp + ∂f /∂σ : Ee : ∂g/∂σ ep  e  ] (∂g/∂σ) ⊗ (∂g/∂σ) ε˙ = C : σ˙ = C + : σ˙ Kp  (2.12)  [  ep  e  (2.13)  where H (L) is the Heaviside step function, H (L) = 1 for L > 0 and H (L) = 0 for L ≤ 0; Eep and Cep represent an elastoplastic modulus and compliances matrices. Eqs. 2.12 and 2.13 are complete sets of constitutive relations for inputs of strain and stress increments, respectively, in a classical elastoplastic model. It is sufficient to provide the yield surface f = 0, plastic potential g = 0 or any other appropriate function for Np , and evolution laws for plastic strain rate and rate of internal variables shown generally in Eqs. 2.5 and 2.6  2.2.2  Plastic strains inside the yield surface  The observed experimental undrained cyclic response of clays in triaxial tests illustrated that the behavior of clay is not elastic within the yield surface. From what the question follows if it is necessary to account for inelastic strains within the yield surface in constitutive modeling. The answer depends on what type of loading is required to be simulated. For monotonic loading, concerned with peak and ultimate strength of clay or the stress states mainly on the yield surface, classical elastoplastic models with incorporated certain complexities are able to predict soil response with fair accuracy. If monotonic loading is concerned with simulated soil response for the stress states within the yield surface, model allowing inelastic strains within the yield surface should be considered. For earthquake loadings, the imposed stress path always involves loading and unloading processes. Even if the loading stress path is on the yield surface, unloading is directed towards the interior of the yield surface which is then followed by reloading. The necessity to simulate plastic deformations within the  13  2.2. Theoretical modeling yield surface becomes pronounced as the number of such loading-unloading-reloading cycles increases. Earthquake-induced deformations within the soil often involve a considerable number of these cycles, which does require simulation of inelastic soil response inside the yield surface for successive model predictions. As has been discussed in section 2.1.1, experimental testing requires that slow undrained cyclic loading should be performed for accurate pore water pressure measurements. In slow undrained cyclic loading tests rate effect is still present, but to a much lower degree than in fast undrained cyclic loading. Therefore, literature review is conducted below on rateindependent models with incorporated bounding surface formulation. The bounding surface framework is an attractive tool to model undrained cyclic loading imposed on clays used by many researchers. The bounding surface formulation can be incorporated into any classical elastoplastic model by adjusting a plastic modulus Kp and leaving all other relationships unchanged. First, development of the bounding surface formulations is discussed. Then, rate-independent elastoplastic models with incorporated bounding surface formulation for monotonic and cyclic loading simulations are reviewed. Not all of the reviewed models are referred to by the authors proposing them as bounding surface plasticity models; however, the concept adopted for plastic strains simulation within the yield surface is similar to that used in the bounding surface formulation. To keep a common definition of the yield surface, inside which only elastic strains are allowed, bounding surface term is used for the yield surface with inelastic strains allowed inside it. As in various publications, the bounding surface is referred to in different ways, a conventional yield surface term is to be used in such cases.  2.2.3  Development of the bounding surface formulation  Attempts to simulate inelastic response within the conventional yield surface were made in the early 1950s; slip theory proposed by Batdford and Budiansky (1949) was adopted by Koiter (1953) for this purpose. Further development of models describing soil behavior in multiaxial stress space was conducted separately by Iwan (1967) and Mr´oz (1967, 1969). They adopted an idea of family of surfaces of progressively larger size, one enclosing the other, with constant work-hardening moduli for each surface. The inner-most surface enclosed a 14  2.2. Theoretical modeling truly elastic region, while the outer-most served as the conventional yield surface. The change of the work-hardening modulus was linear for the stress states on a single surface until a surface of larger size was reached; at this point the work-hardening modulus changed and corresponded to that of a larger surface at the current stress point until next surface of larger size was attained. In this way, piecewise constant work-hardening moduli allowed simulation of inelastic strains inside the outer-most surface. The main disadvantage of this approach is that a number of surfaces have to be stored in the model memory. The need of a simpler approach inspired development of the bounding surface formulation proposed by Dafalias and Popov (1975). Rules of plastic modulus variation were proposed based on uniaxial cyclic loading for metals and generalized to multiaxial stress space. Since its introduction, the bounding surface formulation has been used by various investigators. In a series of publications by Dafalias (1986a), Dafalias and Herrmann (1986), and Anandarajah and Dafalias (1986), the bounding surface formulation is reviewed for simulations of soil response. Because the current work is dedicated to cyclic response of clays, the publications mentioned will be discussed in detail. The bounding surface formulation introduces inelastic strains for the stress states inside the conventional yield surface through plastic modulus. All other relations assumed for the bounding surface are similar to those used in classical elastoplastic models. The conventional yield surface is replaced by the bounding surface. Within the bounding surface F = 0, plastic strains depend on the distance from the actual stress point σij to the corresponding image point σ ¯ij on the bounding surface, where the latter is determined by a mapping rule. The degree of plasticity increases as this distance becomes smaller and reaches its limit once the stress point is on the bounding surface. A loading surface f = 0 passes through σij and is enclosed by the bounding surface so that both surfaces coincide once a current stress point is on the bounding surface. The bounding surface concept requires conformance of the following: 1. the mapping rule has to be such that a unique σ ¯ij is defined for any stress state within the bounding surface, and ∂f /∂σij = ∂F/∂ σ ¯ij , where ∂f /∂σij and ∂F/∂ σ ¯ij are gradients of the loading surface and the bounding surface,respectively, for σij = σ ¯ij to guarantee that the loading surface never intersects the bounding surface; the mapping 15  2.2. Theoretical modeling rule is not invertible in a sense that an infinite number of the stress states within the bounding surface could correspond to a given σ ¯ij . 2. plastic strains are allowed for the stress states within the bounding surface; thus, as σ ¯ij must remain on the bounding surface, the plastic modulus for the stress states on the bounding surface is obtained from the consistency condition, F˙ = 0: ¯ p = − ∂F : q¯˙n K ∂ q¯n  (2.14)  ¯ p through Euclidian distance: 3. The actual plastic modulus Kp is related to K δ = [(¯ σij − σij ) (¯ σij − σij )]1/2  (2.15)  ¯ p for δ > 0, and Kp → K ¯ p for δ → 0. such that Kp > K Dafalias (1986b) reviewed two different mapping rules; each of the two mapping rules deserves to be discussed in more detail. Original bounding surface mapping rule This mapping rule was originally proposed by Dafalias (1975). The loading and bounding surfaces, having kinematic hardening rules that can be independent, were expressed respectively as follows: f (σij − αij , qn ) = 0  (2.16)  F (¯ σij − βij , qn ) = 0  (2.17)  where αij and βij are the corresponding kinematic hardening internal variables. Adopted kinematic and other hardening rules for the two surfaces should be such that surfaces never intersect. Then the mapping rule, with the assumption that bounding and loading surfaces had similar shapes, was given by: σ ¯ij = b (σij − αij ) + βij  (2.18) 16  2.2. Theoretical modeling where b is a similarity ratio between the two surfaces. The proposed plastic modulus for any stress state was: ¯ p + h (|δ in |) Kp = K  δ: n ⟨(δ in − δ) : n⟩  (2.19)  where δ = (¯ σij − σij ); δ in is equal to δ at the moment of initiation of new loading process; and n is a loading direction. The δ in was updated when (δ in : n − δ : n) ≤ 0. Notice that the first term was always larger than the second for a loading process and vice versa for unloading as the terms changed signs. The use of tensorial quantities was advocated by the case of neutral loading simulations (see Dafalias, 1986b, for details). With the original bounding surface mapping rule, Dafalias (1986b) was aware that the problem of “overshooting” exists. This was the case for short unloading with new δ in after the loading process, followed by reloading with updated δ in . This would result in a very small denominator in Eq. 2.19, i.e., a large plastic modulus, as opposed to the case where no such small perturbation had occurred and stress–strain curve in reloading would overshoot the one in loading. However, the overshooting was not infinite, as it was bounded by the bounding surface. Radial mapping rule Dafalias (1986b) referred to the work of Hashiguchi and Ueno (1977) as the pioneering work for the radial mapping rule, where the origin had been used as a projection center. The proposed radial mapping rule, shown in Fig. 2.2, required back-stress ratio α, which was the projection center, to be always located within or on the bounding surface. The image point on the bounding surface was a projection from α through the current stress point. The projection center also served as the center of homology between the loading and bounding surfaces, such that ∂f /∂σij points in the same direction as ∂F/∂ σ ¯ij . This ensures that the two surfaces never intersect. The radial mapping rule was given by σ ¯ij = b (σij − αij ) + αij  (2.20)  where the value of b was determined by plugging Eq. 2.20 into F = 0, 1 ≤ b ≤ ∞. Because of the two surfaces being homologous to each other no consistency condition was required 17  Ljj = dF/d&ij, define indirectly a loading surface / = 0 shown by a dashed line in Fig. 1, which is homologous to F = 0 with a,, being the center of 2.2. Theoretical modeling  ELASTIC NUCLEUS LOADING SURFACE J  ki  -BOUNDING SURFACE F ( £ f M , , „ ) = !>  Figure 2.2: Radial mapping rule from Dafalias (1986b)  FIG. 1.—Schematic Diagram of Bounding Surface and Related Concepts for the stress states within the bounding surface.  977  The plastic modulus equation for any stress state was proposed as follows: ¯p + H ˆ Kp = K  δ ¯ p + H⟨ ˆ b − s⟩−1 =K ⟨r − sδ⟩ b−1  (2.21)  ˆ is a positive shape where r is the distance from the projection center to the image point, H hardening function, s is an elastic factor, and r/δ = b/ (b − 1). When δ ≥ r/s, i.e., for the stress states within the elastic nucleus, Kp = ∞, rendering purely elastic strains. For s = 1, the elastic domain degenerated into a single point, whereas for s → ∞, the elastic domain coincided with the bounding surface, rendering a conventional elastoplastic model. A particular radial mapping rule for the bounding surface formulations was presented by Dafalias and Herrmann (1986), as a continuation of previous work by Dafalias (1986b), for isotropic cohesive soils within the Critical State Soil Mechanics (CSSM) theory. In addition 18  2.2. Theoretical modeling to the radial mapping, the authors included Lode angle dependency in the proposed model, which enhanced predictive capabilities for more general stress states. For the formulation, I was a first invariant of the stress tensor, and J was a second invariant of the deviatoric stress tensor. The proposed projection center in Fig. 2.3 was positioned on the hydrostatic line, represented by I-axis, and was related to the size of the bounding surface along this axis by: Ic = Cp I0  (2.22)  where Ic is the position of the projection center on the I-axis, Cp is a new model constant in the range of 0 ≤ Cp < 1, and I0 is the size of the bounding surface. Similar to Eq. 2.20, image stress for a given projection center was: I¯ = b (I − Cp I0 ) + Cp I0  (2.23)  The composite bounding surface, consisting of two ellipses and one hyperbola with continuous tangents, as shown in Fig. 2.3, was used for better predictions at high OCRs. In addition, the use of R > 1 parameter, which determined the position of point I1 in Fig. 2.3 and related two major axes of the ellipse 1 as I1 = I0 /R, was considered. Lode angle was denoted by α. The plastic modulus for any stress state was the same as the one in Eq. 2.21. Dafalias ˆ and considered the and Herrmann (1986) elaborated on the shape hardening function H following equation: ˆ = 1 + ein g ∗2 pa [z m h(α) + (1 − z m ) h0 ] H λ−κ  (2.24)  where ein is an initial void ratio; λ and κ are slopes of the normal compression and the ( )1.2 rebound lines, respectively, in the e-ln p plane; g ∗ = (∂F/∂ p¯)2 + (∂F/∂ q¯)2 was used to account for a jump of a bounding surface gradient in transition points between the three surfaces comprising the bounding surface; pa is an atmospheric pressure; z = J/J1 ; m = 0.02, rendering z m ≈ 1 to maintain dependency on h(α) rather than h0 ; h(α) is a shape hardening function of h, dependent on α, and it required determination of hc and he for triaxial compression and extension; h0 = (hc + he ) /2.  19  2.2. Theoretical modeling  ELLIPSE 1 HYPERBOLA  ELLIPSE 2 ELASTIC  Figure 2.3: Radial mapping rule from Dafalias and Herrmann (1986)  FIG. 1.—Schematic Illustration of Bounding Surface and Radial Mapping Rule in Stress Invariants Space  Effects of hc and C parameters in triaxial compression were shown in detail. While the  it will be the explained sequel. Eq. 30 has specific former controlled degree in of the plastic deformation for the the following stress states withinforms. the bounding For ellipse 1: surface, the latter was successfully used to control contracting or dialative responses for  F=(I-I0stress )[l+^I (31) 0)+(R-lfU) overconsolidated states. Calibration of the =0. parameters and model performance in undrained triaxial compression suggested that the projection center along the Formonotonic the hyperbola: 2  hydrostatic axis was an attractive feature of the model for monotonic loading simulations LJjL = 0. I- h (32) N/ RJ \N R N R\ for high OCR states.  For ellipse Anandarajah and2:Dafalias (1986) further incorporated inherent and induced anisotropy 2 F = (7 - TIsurface 2QI0] + p/ = 0 The orientation of the bounding(33a) 0)(I - (X into the bounding soil+ plasticity model. surface was T(Z + TF) T2 controlled by evolving anisotropy, represented by the second-order tensor δija in Fig. 2.4, (33b) £ = - - = - — = — ; p Z(Z + 2TF') Z + 2TF' during the loading process. The projection center was determined as Cp I0a , where Cp was N i N RA 2 and Herrmann (1986) and I0a was the same model parameter proposed by= Dafalias (1 + y VT + 7 ) (33c)a size of :; Z " N VT+f' R the bounding surface; the radialp mapping rule was also7 identical, given the projection center. The dependence on e occurs through loie ), and the dependence on a through the parameters N(a), R(<x), and A(a) according to = g(a,c)Q (34«) Due Q(a) to the bounding surface orientation being controlled by evolved anisotropy, the c Qe and the shape-hardening function were adjusted accordingly, with the latter plastic modulus c=~ (34b) C g(a,c) = — > , 1 + c - (1 - c) sin 3a  (34c) 20  1270  2.2. Theoretical modeling  FIG. 1.—Schematic Bounding Surface and Geometrical InterpretaFigure 2.4: RadialIllustration mappingofrule from Anandarajah and Dafalias (1986) tion of Stress Variables  given by:  was accounted for in a rather simple manner. Additionally, the concept of elastic nucleus was introduced in order to model the cyclic stabili[ ( ) ( )] 1 + ein bymclays. zation m a 2 a 2 ˆ =exhibited ¯ ¯ H p [z h(α) + (1 − z ) h ] 9 ∂F/∂ I + 3 ∂F/∂ J (2.25) a 0 Some important modifications are incorporated in applying the conλ−κ cept of bounding surface for anisotropic clays. Referring to Fig. 1, the bounding surface is assumed to be oriented at any stage along a line (/"The model was tested in comparison with experimental in undrained triaxial comaxis) whose direction changes during the loading, data representing the evolution of anisotropy. A symmetric second-order tensor 8°;, such that pression and teststheonorientation isotropically and bounding K0 -consolidated 8|8| extension = 3, defines of the surface specimens. axis (/"-axis),The use of in the stress with space. the concept of bounding the bounding surface theAccording projectiontocenter on evolving I a -axis surface, improvedevery model perforpoint in the stress space, such as point B in Fig. 1, is associated with a corresponding point on0 -consolidated the boundingsamples surfacefor viathe a suitable map-within the mance during undrained image shearing on K stress states ping rule. The current theory, as in the isotropic formulation, employs bounding asurface. radial mapping rule (Fig. 1) where point C is obtained by radially projecting point B from point A, which is assumed to be on the /"-axis. The Herrmann and forAnandarajah In the projection radial mapping adopted by Dafaliasofand centerrule A is the counterpart point Ic of the(1986) isotropic mulation, in figure of Part II (11). The scalar between and Dafalias (1986) shown for isotropic and1anisotropic cohesive soils, distance respectively, the projection the origin and point.A is taken as CI"0 where C is a model parameter, referred to as the on projection parameter. The gradient of occurred the bounding center imposed a restriction plastic strains. Purely elastic strains upon unloading surface at the image point C, defines the loading/unloading direction, until loading index = ((∂F/∂ σ ¯ijof ) σ˙ plastic became positive the definition of the plastic as well as Lthe direction rate for theasactual stress point ij ) /Kp strain B. The plastic modulus Kp at the actual stress point a,y, is assumed to strains increment is ε˙ij = ¯ijmodulus , where Macauley brackets rendered nullimage for negative L. be a function of ⟨L⟩ the ∂F/∂ plasticσ Kp (bounding modulus) at the point might <r,y, and scalar distance 8 between the simulations, image and the This restriction be the appropriate for monotonic loading whileactual unsatisfactory stress points. 1294  21  2.2. Theoretical modeling results are expected for undrained cyclic loading simulations.  2.2.4  Bounding surface models for monotonic response of clays  Although the study in this thesis is dedicated to simulations of cyclic clay behavior, it is also of interest to investigate incorporated bounding surface formulation proposed for monotonic loading simulations, as some important points and observations can be inferred. Whittle (1993) The bounding surface formulation was incorporated into the MIT-E3 model by Whittle (1993). The model was able to predict initial and induced anisotropy through a rotational hardening rule. A projection center was always positioned at the origin of the stress space for adopted radial mapping rule. An elastoplastic modulus HI was proposed based on K0 -loading on normally consolidated clays to simulate parallel virgin consolidation line in compression plane for any radial consolidation path, strain-hardening, peak strength, strain-softening, and critical stress states. For the stress states within the bounding surface the actual elastoplastic modulus H was related to HI and value of H at first yield (first loading for stress states within the bounding surface) which is denoted by H0 : H = HI + H0 g 2  (2.26)  where g2 is a mapping function describing the relative position of the current stress state. MIT-E3 model required calibration of 15 parameters, including initial void ratio. Model simulations were then compared with experimental results in undrained compression (CK0 UC) and undrained extension (CK0 UE) tests on K0 -consolidated clay. For OCR > 4, the model overestimated shear strength in compression and underestimated in extension. Jung and Yune (2011) Jung and Yune (2011) incorporated the bounding surface formulation into the ABSM model previously proposed by Banerjee and Yousif (1986). The main goal was to simulate smallstrain behavior of a particular clay. This was achieved by the use of a nonlinear anisotropic 22  2.2. Theoretical modeling elastic model obtained from singular value decomposition of a compliance matrix. Further, they thoroughly investigated different projection center locations in adopted radial mapping rule. From stress probes in different directions, they concluded that the projection center at the last stress reversal rendered the most appropriate plastic strains. The plastic modulus for any stress state was modified by inclusion of an extra distance parameter to further enhance simulations at small strains. (  [ ∗  H = A + A0 (1 − h)  δ/δ0 1 − δ/δ0  )µ  ( +h  δ/δ1 1 − δ/δ1  )µ ] (2.27)  where A∗ is a plastic modulus at the image stress, A0 is a reference modulus at the apex of the bounding surface, h is a weighting factor, definition of δ is the same as in previously described models, δ0 is equal to δ at initiation of a new loading process, δ1 is a distance from the current stress point to the image stress point where the soil became stiffer prior to shearing. Computed stress–strain curves with 12 calibrated model parameters gave a good match for small strains observed in triaxial tests, while underestimated shear stiffness at larger strains. Suebsuk et al. (2011) The bounding surface formulation was incorporated into Modified Cam-Clay model with destructuration (MSCC) by Suebsuk et al. (2011). The origin of the stress space was adopted as the projection center in a radial mapping rule. Suebsuk et al. (2011) utilized ratio dependence of the plastic modulus via introduction of the image stress ratio, which was a function of current and image stress states: α = [(p/pj ) + (q/qj )]  (2.28)  where (p, q) and (pj , qj ) are current and image stress states, respectively. The plastic modulus was then given by the following equation: (1 − α)2 H = Hj + (hHj,i ) α  (2.29)  23  2.2. Theoretical modeling where h is a model parameter, Hj is a plastic modulus at image point, Hj,i is equal to Hj at initiation of the loading process. The proposed model required calibration of 12 parameters: 6 parameters were the same as in the original Modified Cam Clay model, 5 parameters were required for destructuration simulations, and parameter h was needed for the bounding surface formulation. The model predictions were compared with undrained monotonic triaxial tests. Volumetric strains in the e-ln p plane induced by destructuration in consolidation process were closely matched by the model. Jiang et al. (2012); Ling et al. (2002) The model proposed by Ling et al. (2002) is an extension of the anisotropic model by Dafalias (1987), where rotational hardening was derived based on the critical state concept and deviatoric back-stress ratio α as a dissipating energy measurement. The associative flow rule and inherent anisotropy were used in the model; the bounding surface was a single ellipsoid in the multiaxial stress space. Isotropic, rotational, and distortional hardening rules were adopted for the bounding surface evolution. Distortional hardening was required to control the evolution of R. The projection center was always located on a back-stress ratio, representing stressinduced anisotropy, in the same way as proposed by Anandarajah and Dafalias (1986), while the plastic modulus was modified to reduce the number of model parameters as follows: ¯ p + 1 + ein pa Rij Rij Kp = K λ−κ  (  δ ⟨r − sδ⟩  )W (2.30)  where Rij = ∂F/∂ σ ¯ij and W is a model parameter replacing hc , he , z, and m, and definition of δ is the same as in previously described models. In total, 12 model constants were required for simulations. Model calibration and performance were illustrated as compared with experimental data on isotropically and anisotropically consolidated kaolin clay, San Fransisco Bay Mud, and Boston Blue clay to a range of OCRs. Overall model performance was enhanced for overconsolidated states; however, because of the associative flow rule, the model was not able to predict the softening observed in CK0 UC tests on normally consolidated specimens. 24  2.2. Theoretical modeling Jiang et al. (2012) recently proposed an extension of the model discussed above. The nonassociative flow rule proposed by Jiang and Ling (2010) was adopted in a new formulation. This flow rule considered the bounding surface and the plastic potential possessing the same back-stress ratio, controlling rotation of the two surfaces, and having different peak deviatoric stresses. For the plastic potential, such peak deviatoric stress represents the critical stress ratio at which infinite shear strains occur without change in stresses, while peak deviatoric stress is used for simulations of the shear strength of soil. The projection center was not changed from the preceding model, while the plastic modulus was slightly modified in that the W power was omitted and the hardening parameter, h, was included in the second term of Eq. 2.30. With the newly proposed non-associative flow rule, softening during CK0 UC stress path can be successfully simulated. Huang et al. (2011) Huang et al. (2011) proposed an anisotropic model with incorporated destructuration and the bounding surface formulation. The analytical description of the bounding surface was adopted from Ling et al. (2002), as was evolution of anisotropy from Wheeler et al. (2003). An exponential decay for a scalar valued variable, accounting for destructuration, was proposed with the rate of destructuration parameter and its initial value. The radial mapping rule was utilized in the bounding surface formulation with the projection center fixed at the origin. The proposed plastic modulus considered degradation via introduced cumulative plastic deviatoric strain: [( ˆ p + ζpa Hp = H  ∂F ∂p  )2  ( +  ∂F ∂q  )2 ] [(  ψ = ψ0 exp (−ξεps ) where εps =  ∫  δ0 δ0 − δ  ]  )ψ −1  (2.31)  (2.32)  ε˙ps is a cumulative plastic deviatoric strain; ζ, ψo , and ξ are model parame-  ters; the ψ function was used to simulate degradation of stiffness with accumulated deviatoric strain; and definitions of δ and δ0 are the same as in previously described models. A total of 13 model parameters were required for model simulations. Simulations were per25  2.2. Theoretical modeling formed for undrained triaxial compression on normally consolidated and overconsolidated clay. While shearing of normally consolidated specimens was simulated relatively well, lower shear strength was obtained for overconsolidated specimens than observed in experiments. The destructuration effect was captured to some extent in the compression plane.  2.2.5  Bounding surface models for cyclic response of clays  The review of the bounding surface models that have been developed and evaluated for clay response in undrained cyclic loading is conducted below. Only some models are selected and presented in this discussion, as they illustrate some important features for cyclic loading simulations. Other models available in the literature are either somehow related to the ones discussed below or found to be irrelevant for elaboration in the current work (e.g., viscoplasticity models as the time effect is not being considered). Zienkiewicz et al. (1985) A simple model was proposed by Zienkiewicz et al. (1985) which incorporated the bounding surface formulation. The projection center fixed at the origin was considered in the radial mapping rule. Zienkiewicz et al. (1985) distinguished the plastic modulus in unloading and reloading for the stress states within the bounding surface. Unloading was defined for the stress increments pointing towards the interior of the loading surface, passing through the current stress point and the origin, and reloading for the stress increments pointing outwards. Very simple plastic moduli for unloading and reloading were considered. The former was given by:  [ HU =  HLcs  HL =  HLcs  and the latter by  δ0 δ  [  δ0 δ  ]γ1 (2.33) ]γ (2.34)  where HLcs is a plastic modulus at the image stress, γ and γ1 are model constants, and definitions of δ and δ0 are the same as in previously described models. In this way the plastic modulus was larger than HLcs for the stress states inside the bounding surface and was equal to HLcs for the stress states on the bounding surface. 26  2.2. Theoretical modeling The proposed model simulations were then compared with actual experimental data from undrained cyclic triaxial tests. It was clear that the model predicts considerable overdamping as opposed to experimental data and a way to model plastic strains in unloading was required. Based on observations, Zienkiewicz et al. (1985) suggested the use of cumulative plastic deviatoric strain, as the proposed model was not able to simulate accumulation of plastic strains for clays once the stress path reached equilibrium in undrained cyclic triaxial loading. This was also done in the hope of achieving further reduction of stresses to simulate liquefaction for sands; however, it was still not accomplished with inclusion of cumulative plastic deviatoric strains. Zienkiewicz et al. (1985) expressed the need to simulate plastic strains in unloading to avoid overdamping. Al-Tabaa and Wood (1989); Stallebrass and Taylor (1997) Al-Tabaa and Wood (1989) extended the Modified Cam-Clay model by adding a kinematically hardening yield locus within the Cam-Clay locus. The adopted translation of the yield locus guaranteed that the outward normal of the inner yield locus was in the same direction as the outward normal of the Cam-Clay locus. This avoided the possibility that the two loci would intersect. The response inside the yield locus was assumed to be elastic; once the stress point was on the yield locus inelastic strains took place. While the inner yield locus was allowed translating within the Cam-Clay locus, the latter was allowed only expanding or shrinking, as in the Modified Cam-Clay model. Expansion/shrinkage of both loci occurred any time plastic volumetric strains took place. The degree of plasticity for the stress states on the yield locus and inside the Cam-Clay locus was dependent on the distance from the current stress point to its conjugate point on the Cam-Clay locus. This enabled smooth transition of soil response from the initial stress point on the inner yield locus to the stress point on the Cam-Clay locus. This approach of smooth transition is somewhat similar to that proposed for soils by Dafalias (1986b); the difference is that, instead of keeping in memory an additional memory parameter or having a projection center, the memory of the kinematic yield surface was required. The model was attractive in that it required only 2 more parameters in addition to those required for the Modified Cam Clay model. With the 6 model parameters close match of 27  2.2. Theoretical modeling volumetric and shear strains accumulation was achieved during one-way undrained cyclic triaxial loading consisting of 6 cycles. Stallebrass and Taylor (1997) further incorporated an additional kinematic hardening surface, history surface, which was always inside the outer surface and enclosed the inner yield locus. The history surface was required to account for an effect of recent stress history and yield at small strains. The translation laws, ensuring that the three surfaces never intersect, required that the centers of the surfaces moved in the direction of a vector from the current stress state to its conjugate point on the next surface. The surfaces translations occurred anytime plastic deformations took place. The proposed 3-SKH model required 8 model parameters. The model performance was illustrated for a single element simulations in monotonic and two cycles of undrained cyclic triaxial loading with overpredicted stiffness reduction as opposed to observed experiments. Stallebrass and Taylor (1997) implemented the proposed model into finite element code and simulated the available centrifuge test of a stresses distribution underneath a foundation. The analysis showed relatively good model predictive capabilities under a low level of loads, while for increased loads the model underpredicted displacements underneath the foundation. Wathugala and Desai (1993) The Hierarchical Single Surface (HiSS) model approach by Desai et al. (1986) for isotropic material with associative flow rule was used by Wathugala and Desai (1993) as a basis model; the bounding surface formulation was incorporated for plastic strains generation inside the conventional yield surface. Wathugala and Desai (1993) suggested that other complexity features (e.g., non-associative flow rule and anisotropy) could be included in the proposed model if desired. The model considered virgin loading, unloading, and reloading separately. Virgin loading corresponded to the stress states on the conventional yield surface and increments of a stress directed outwards of the conventional yield surface. For virgin loading, the usual approach of finding the plastic modulus from the consistency condition was used. A reference surface, passing through the current stress point, was used to distinguish between unloading and reloading in the same way as loading surface is used in other models. 28  2.2. Theoretical modeling This surface was similar in shape to the conventional yield surface and also passed through the origin; both surfaces were homologous to each other with the center of homology lying at the origin. Based on laboratory tests, Wathugala and Desai (1993) proposed unloading to be essentially elastic, while for reloading the plastic modulus was a function of the plastic modulus for virgin loading such that the former approached the latter: H  RL  =  HIVL 1  +  r1 HIVL 2  ( )r αps 2 1− αr  (2.35)  where HIVL and HIVL are plastic moduli at image stress and intersection of the conventional 1 2 yield surface with the hydrostatic axis, respectively, r1 and r2 are model constants, αps and αr are distances from the origin to the current stress state and to the image stress point, respectively. The αps /αr was used to impose distance dependency of the plastic modulus for the reloading process. The model required calibration of 12 parameters based on isotropic consolidation and undrained triaxial compression/extension tests. Simulations of undrained cyclic triaxial test, consisting of 4 cycles of two-way loading, suggested that the model was able to capture generation of plastic deviatoric strains within a considered number of cycles, while it overestimated damping as a unloading was simulated as elastic. The model was then implemented into finite elements code for pile installation, consolidation after installation, application of load tests, and final cyclic loading simulation; the whole procedure of model implementation, calibration, and simulation was described by Desai et al. (1993). The results obtained gave a relatively good match to field observations of pile segment displacement, shear transfer, and horizontal stress. Yu et al. (2007) A Clay And Sand Model (CASM) proposed by Yu (1998) is a unified model for simulations of soil behavior; CASM is based on the Modified Cam-Clay model with 2 extra parameters required to control the shape of the conventional yield surface. One of the parameters controled the position of the point where the conventional yield surface intersected with the CSL, while the other controled a “drop-like” shape.  29  2.2. Theoretical modeling Incorporation of the bounding surface formulation into CASM model required 5 extra parameters, rendering in total 12 parameters. The radial mapping rule assumed the projection center at the origin of the stress space. The proposed model treated the plastic modulus for virgin loading, unloading, and reloading separately. The description of the virgin loading plastic modulus is omitted and interested readers are referred to the original paper. The plastic modulus for unloading process was given by: H = HU  1 1−γ  (2.36)  where HU is an unloading hardening parameter, representing the limit for H in unloading, γ = p/pj = q/qj is a non-dimensional ratio between current and image stresses. Thus, for an unloading process with initial stress state on the bounding surface, HU ≤ H ≤ ∞. In reloading the plastic modulus was given by: ( H = Hj + HR  1−γ γ  )  (  1 + εpq  )k  (2.37)  where Hj is equal to H at the image point, HR is a reloading hardening parameter, εpq is a cumulative plastic deviatoric strain, and k is a model parameter controlling the rate of a shakedown. In reloading, H = ∞ for the stress states far away from the projected image point, and H → Hj as γ → 1. The εpq allowed for stiffness degradation as plastic deviatoric strains accumulated throughout loading. The model was further compared with bubble-type models proposed by Al-Tabaa and Wood (1989); Stallebrass and Taylor (1997) in simulations of one-way and two-way undrained cyclic triaxial loading with available experimental results. In terms of deviatoric strains, the proposed model illustrated better predictive capabilities and was satisfactory for volumetric strains. For two-way cyclic loading simulations, the model predictions of deviatoric strains were also satisfactory. Li and Meissner (2002) Li and Meissner (2002) proposed a two-surface plasticity model with an associative flow  30  2.2. Theoretical modeling rule and incorporated bounding surface formulation for simulation of clay cyclic undrained behavior. The bounding surface represented the consolidation history of clay, while the loading surface passing through current stress point was related to each loading process. Li and Meissner (2002) utilized the kinematic hardening rule of the bounding surface that allowed its translation in the stress space to simulate induced anisotropy. Expansion and contraction of the bounding surface were also permitted. The elastic domain was reduced to a single point at the last stress reversal. The plastic modulus within the bounding surface was dependent on the relative configurations of the loading and bounding surfaces and varied from its local to prescribed value on the bounding surface as the stress point approached the latter. The authors extended Masing’s rule of stress reversal to a multiaxial case. They defined a memory center as the point where loading changed its direction so that it pointed towards the interior of the loading surface, i.e., the memory center was the point at the last stress reversal. The loading and bounding surfaces tangentially touched each other at the memory center and never intersected. The loading surface was homologous to the bounding surface and passed through the memory center and current stress point. Once the loading surface coincided with the bounding surface, the plastic modulus reached its prescribed value on the bounding surface, and the two surfaces evolved identically. If the stress path was reversed and directed towards the interior of the loading surface, the memory center was updated. A new loading surface evolved, starting from the point at the memory center, and the bounding surface was translated to the same point with its size unchanged so that the two surfaces were tangential. For the stress reversals occurring inside the bounding surface the latter indeed translated from its previous position to the memory center, whereas for the stress reversals on the bounding surface the latter stayed in the same position. To simulate cyclic degradation of clay stiffness, Li and Meissner (2002) used a varying power parameter γ enabling reduction of the plastic modulus: ( ) ¯ m + HM − H ¯m hm = H  (  δ δ0  )γ (2.38)  ¯ m is a plastic modulus for the stress points on the bounding surface, HM is a local where H  31  2.3. Missing features in reviewed models value of the plastic modulus, and definitions of δ and δ0 are the same as in previously described models. The γ was constant for monotonic loading, while for cyclic loading it increased hyperbolically with cumulative shear strain. It should be noted that the proposed HM required specification of the applied dynamic shear stress amplitude qd and monotonic shear stress at failure qf , which must be known prior to cyclic loading. This makes the model not applicable for simulation of irregular cyclic loading, where amplitude of shear stress cannot be specified in advance. The model required 12 parameters, including qd and qf , for cyclic loading simulations. Li and Meissner (2002) did not indicate if the clay tested in monotonic and undrained cyclic triaxial loadings was the same and prepared under similar conditions; however, a different set of parameters was used for monotonic and each cyclic loadings. Monotonic loading simulations closely matched experimentally observed clay responses. One-way and two-way undrained cyclic triaxial tests were well simulated in terms of reduction of p and accumulation of shear strains.  2.3  Missing features in reviewed models  The presented literature review on experimental undrained cyclic response of clay suggests plastic strains should be reflected in model simulations for the stress states inside the conventional yield surface. The stress path migration towards the origin was accompanied by plastic deviatoric/shear strains development. When the stress path almost stopped migrating towards the origin, plastic deviatoric/shear strains kept developing with an increasing rate (Takahashi et al., 1980; Zergoun, 1991; Zergoun and Vaid, 1994). The bounding surface formulation is used by different researchers to simulate the cyclic response observed in the experiments. However, Zienkiewicz et al. (1985) noticed that model simulations with the plastic modulus possessing only the distance dependence lacked the feature of simulating increased rate of deviatoric/shear strain accumulation once the stress path stopped migrating towards the origin. They proposed cumulative plastic deviatoric strain for continuous reduction of plastic modulus throughout applied cyclic loading. A similar approach for reducing the plastic modulus during simulations of undrained cyclic triaxial loading was used by Yu 32  2.3. Missing features in reviewed models (1998) for the isotropic model, and by Li and Meissner (2002) for the kinematic hardening two-surface model. Huang et al. (2011) used the same idea for monotonic loading simulations in an anisotropic model with destructuration. In addition to the observed continuous development of plastic deviatoric/shear strains during undrained cyclic triaxial tests, a certain cyclic stress amplitude existed below which the rate of deviatoric/shear strains development decreased, This trend was observed in all experiments (Bacchus, 1969; Sangrey et al., 1969; Takahashi et al., 1980; Zergoun and Vaid, 1994). Another important conclusion made by Zienkiewicz et al. (1985) was that with the projection center fixed at the origin simulated unloading exhibited purely elastic strains, while plasticity was in effect suggested by experiments. This resulted in overdamping produced in simulations of undrained cyclic triaxial loading. The projection center at the last stress reversal appears to solve a problem of overdamping in simulations of clay response. Jung and Yune (2011) concluded, based on stress probe experiments, that the projection center at the point of the last stress reversal predicted plastic strain increments closer to those observed in experiments; however, their focus was on application of the bounding surface formulation for simulations at small strains in monotonic loading rather than cyclic loading. Li and Meissner (2002), based on Masing’s rule generalization, adopted the projection center fixed at the last stress reversal and were able to obtain relatively good predictions of one-way and two-way undrained cyclic triaxial loading. The main disadvantage of their model was that it required the knowledge of the applied stress amplitude prior to cyclic loading. Some of the models reviewed treated the plastic modulus in loading, unloading, and reloading separately (e.g., Wathugala and Desai, 1993; Yu, 1998; Zienkiewicz et al., 1985). This typically required a larger number of parameters (e.g., Yu, 1998; Zienkiewicz et al., 1985) and treatment of unloading involving only elastic strains (e.g., Wathugala and Desai, 1993). The goal of this thesis is to incorporate the bounding surface formulation into an existing model for better predictions of clay response in undrained cyclic loading imposed by earthquakes. In the development of the model, the position of the projection center and simulation of cyclic softening will be addressed based on the observed experimental evidence and the current trend in available constitutive models. The complexities incorporated 33  2.3. Missing features in reviewed models into the developed model should be adequate to maintain its simplicity for further model implementation.  34  3. Development of the Proposed Model A proposed model is developed by incorporation of the bounding surface formulation into an existing reference model. First, the reference model itself is described with some simplifications in the triaxial plane and then generalized to the multiaxial space. The proposed model is then presented in the same manner. The bounding surface constants in the proposed model are studied in a sensitivity analysis.  3.1  Reference Model  A Simple ANIsotropic CLAY plasticity model with destructuration (SANICLAY-D) proposed by Taiebat, Dafalias and Peek (2010) is chosen as a reference model. The SANICLAYD model incorporates induced anisotropy and destructuration mechanism. Some simplifications for SANICLAY-D, which are believed to be sufficient for successive simulations in monotonic loading, are applied.  3.1.1  Formulation in the triaxial space  The formulation of the reference model in the triaxial space is presented in terms of the following stress and strain quantities: p=  1 (σa + 2σr ) , 3  εv = (εa + 2εr ) ,  q = (σa − σr )  (3.1)  2 (εa − εr ) 3  (3.2)  εq =  35  3.1. Reference Model where subscripts a and r denote the axial and radial directions, respectively, and subscripts v and q denote volumetric and deviatoric components, respectively. Elastic relations With the total strain rate additive decomposition, ε˙ = ε˙ e + ε˙ p , into elastic (subscript e) and plastic (subscript p) parts, the isotropic hypoelastic relations ε˙ev =  p˙ , K  ε˙eq =  q˙ 3G  (3.3)  are adopted. K and G are bulk and shear moduli expressed as: K=  p (1 + e) , κ  G=  3K (1 + e) 2 (1 + ν)  (3.4)  where e is a void ratio and κ is a slope of the rebound line in the e– ln p plane. Flow rule and yield surface The volumetric and deviatoric plastic strain rates are given by the flow rule ε˙pv = ⟨L⟩  ∂g , ∂p  ε˙pq = ⟨L⟩  ∂g ∂q  (3.5)  where g is the plastic potential and L is a plastic multiplier or a loading index to be defined in the sequel of the formulation. ⟨L⟩ are Macauley brackets rendering ⟨L⟩ = L when L > 0 and ⟨L⟩ = 0 when L ≤ 0. The flow rule requires the plastic potential g = 0, which is proposed by Dafalias (1986a): ) ( g = (q − pα)2 − M 2 − α2 p (pα − p)  (3.6)  where pα is the value of p at q = pα, adjusted so that g = 0 passes through a given (p, q), α is a non-dimensional anisotropic variable or rotational hardening variable, M is a critical stress ratio, Mc and Me are critical stress ratios in compression and extension, respectively. M = Mc when the stress ratio η = q/p > α, and M = Me when the stress ratio η = q/p < α.  36  3.1. Reference Model One must have |α| < M for real-valued p and q in Eq. 3.40. The plastic potential is of the shape of a rotated and distorted ellipse, as illustrated in Fig. 3.1.  q  Mc  N  (p, q)  ∇g  α  p  Yield surface (f=0)  Me N  Plastic potential (g=0)  Figure 3.1: Reference model schematic diagram in the triaxial plane On the basis of Eqs. 3.5 and 3.40, the flow rule is given by ε˙pv = ⟨L⟩  ( ) ∂g = ⟨L⟩p M 2 − η 2 , ∂p  ε˙pq = ⟨L⟩  ∂g = ⟨L⟩2p (η − α) ∂q  (3.7)  The yield surface was assumed similar to the plastic potential by Dafalias (1986a), Dafalias et al. (2006), and Taiebat, Dafalias and Peek (2010) with M replaced by N and α by β, where N represented peak q-stress, and β was an anisotropy variable for the yield surface. Jiang and Ling (2010) proposed the yield surface possessing the same anisotropy variable as the plastic potential. This approach is attractive because it eliminates storing the extra tensor. Following Jiang and Ling (2010), the expression for the yield surface is  37  3.1. Reference Model then given by  ( ) f = (q − pα)2 − N 2 − α2 p (p0 − p)  (3.8)  where p0 is the size of the yield surface and represents the isotropic hardening variable, and α is the mentioned anisotropic variable. Clearly, one must have |α| < N for real–valued p, q in Eq. 3.8. The yield surface is also shown in Fig. 3.1. Rate evolution equations for the internal variables Evolution laws for internal variables p0 and α are necessary in the sequel of the formulation. The p0 evolution is given by p˙0 = ⟨L⟩¯ p0  (3.9)  p¯0 is defined to account for the destructuration mechanism as p¯0 = Si p¯0d + S¯i p0d  (3.10)  where Si ≥ 1 is an isotropic structuration factor, and p0d is a destructured value of p0 (or the value of p0 when Si = 1). p¯0d and S¯i are given, respectively, by ( p¯0d =  1+e λ−κ  ( S¯i = −ki  )  1+e λ−κ  ( p0d  ∂g ∂p  ) (3.11)  ) (Si − 1) ε¯pd  (3.12)  where λ is a slope of the normal compression line in the e– ln p plane, ki is a rate of destructuration. Degradation of Si is considered through destructuration plastic strain rate ε˙pd defined as ε˙pd  =  ⟨L⟩¯ εpd  √ p2 = (1 − A)ε˙p2 v + Aε˙q  (3.13)  A is a parameter distributing destructuration between plastic volumetric strain and plastic deviatoric strain increments; by default A = 0.5.  38  3.1. Reference Model The evolution of rotational hardening α is expressed as ( α˙ = ⟨L⟩¯ α = ⟨L⟩  1+e λ−κ  )  ( C  p¯ p0  )2  ( ) ∂g |η − xα| αb − α ∂p  (3.14)  where the definition of α ¯ is self evident, C is a model constant controlling the rate of rotational hardening, and αb = min (N, Me ) is a bounding value of α to guarantee real values of p, q in Eqs. 3.40 and 3.8. Loading index and plastic modulus For the completion of the model, loading index L and plastic modulus Kp are obtained from the consistency condition f˙ = 0 applied to Eq. 3.8, which in combination with Eqs. 3.9 and 3.14 yields 1 L= Kp  (  ( Kp = −  3.1.2  ∂f ∂f p˙ + q˙ ∂p ∂q  )  ∂f ∂f p¯0 + α ¯ ∂p0 ∂α  (3.15) ) (3.16)  Formulation in the multiaxial space  In multiaxial stress and strain spaces second-order tensors are denoted by boldface symbols. Triaxial stress quantities are generalized to multiaxial space, given a stress tensor σ that can be decomposed into its hydrostatic pI and deviatoric s components, as follows: p=  trσ , 3  s = σ − pI  (3.17)  where tr is a trace operator and I is an identity tensor; Triaxial strain quantities are similarly generalized to multiaxial space, given that a strain tensor ε can be decomposed into its volumetric εv I/3 and deviatoric e parts: εv = trε,  e=ε−  εv I 3  (3.18)  A systematic generalization of the triaxial constitutive relations is based on the observa39  3.1. Reference Model tion that under triaxial conditions the following holds true: x = [(3/2) x : x]1/2  (3.19)  where x is any deviatoric stress tensor and x is its triaxial counterpart; the symbol : denotes the trace of the product of adjacent tensors. On the basis of Eq. 3.19: q = [(3/2) s : s]1/2  (3.20)  εq = [(2/3) e : e]1/2  (3.21)  α = [(3/2) α : α]1/2  (3.22)  Elastic relations The strain tensor and its rate are decomposed into elastic and plastic counterparts, and the same is applied for εv and e. The multiaxial generalization of Eq. 3.3 is given by ε˙ e =  s˙ p˙ + I 2G 3K  (3.23)  Flow rule and yield surface The generalization of a flow rule is as follows: ε˙ p = ⟨L⟩  ∂g , ∂σ  (3.24)  where g is the plastic potential and L is the loading index. The Eq. 3.40 for the plastic potential generalizes to ( ) 3 3 2 g = (s − pα) : (s − pα) − M − α : α p (pα − p) = 0 2 2  (3.25)  40  3.1. Reference Model M is chosen to be Lode angle dependent: M=  2m Me ;m= (1 + m) − (1 − m) cos 3Θ Mc  cos 3Θ =  √ 6 tr n3 ; n =  (3.26a)  r−α  (3.26b)  [(r − α) : (r − α)]1/2  where Θ is a Lode angle in the range from 0 to π/3 corresponding to the stress-ratio definition of compression and extension, respectively, where effective stress-ratio (r − α) rather than r is used to define Θ. Given stress tensor σ, Eq. 3.25 yields the following solution for pα : pα =  (3/2) (s − pα) : (s − pα) +p (M 2 − (3/2) α : α) p  (3.27)  The gradient ∂g/∂σ in Eq. 3.24 is calculated based on Eqs. 3.25 and 3.26, as follows: ( ) ∂g 1 3 ∂g ∂Θ 2 = 3 (s − pα) + p M − r : r I + ∂σ ¯ 3 2 ∂Θ ∂ σ ¯  ∂g = 6M 2 p (p0 − p) ∂Θ  (  (1 − m) sin 3Θ (1 + m) − (1 − m) cos 3Θ  )  [ ] −3 n3 − (trn3 )n − 13 I (1 + tr(n2 α) − trn3 tr(nα)) ∂Θ = ∂σ p [(3/2) (r − α) : (r − α)]1/2 sin 3Θ ( tr  ∂Θ ∂σ  )  (3.28a)  (3.28b)  (3.28c)  ( ) ( ) 3 ∂g 3 [tr (n2 α) − trn3 tr (nα)] 2 = p M − r: r + 2 ∂Θ [(3/2) (s − pα) : (s − pα)]1/2 (1 − 6tr2 n3 )1/2 (3.28d)  Similar to the plastic potential, the yield surface is generalized to ( ) 3 3 2 f = (s − pα) : (s − pα) − N − α : α p (p0 − p) = 0 2 2  (3.29)  41  3.1. Reference Model N is assumed to be Lode angle independent, which significantly simplifies the generalization of the yield surface gradient ∂f /∂σ: ( ) ∂f 1 3 2 = 3 (s − pα) + p N − r : r I ∂σ ¯ 3 2  (3.30)  Rate evolution equations for the internal variables The p0 evolution in multiaxial space is given by p˙0 = ⟨L⟩¯ p0  (3.31)  p¯0 is defined to account for destructuration mechanism as in the triaxial plane p¯0 = Si p¯0d + S¯i p0d  (3.32)  p¯0d and S¯i are generalized, respectively, to ( p¯0d =  1+e λ−κ (  S¯i = −ki  )  ( p0d tr  1+e λ−κ  ∂g ∂σ  ) (3.33)  ) (Si − 1) ε¯pd  (3.34)  where Si is considered through ε˙pd and defined in multiaxial space as ε˙pd  =  ⟨L⟩¯ εpd  √ ˙ p : e˙ p = (1 − A)ε˙p2 v + A (2/3) e  (3.35)  e˙ p = ε˙p − (1/3) ε˙pv I is a plastic deviatoric strain tensor. The evolution law for the plastic potential and the yield surface rotational hardening α is generalized to α = ⟨L⟩α ¯ ( α ¯=  1+e λ−κ  )  ( C  p p0  (  )2 tr  ∂g ∂σ ¯  ) [  ]1/2 ( b ) 3 (r − xα) : (r − xα) α −α 2  (3.36a)  (3.36b)  42  3.2. Proposed Model √ αb =  3 r − xα min (N, Me ) nx ; nx = 1 2 [(r − xα) : (r − xα)] 2  (3.36c)  Loading index and plastic modulus For the completion of the model, loading index L and plastic modulus Kp are obtained from the consistency condition f˙ = 0 applied to Eq. 3.29, which in conjunction with Eqs. 3.31 and 3.36 yields the following: 1 L= Kp ( Kp = −  3.2  (  ∂f : σ˙ ∂σ  )  ) ∂f ∂f ¯ p¯0 + :α ∂p0 ∂α  (3.37)  (3.38)  Proposed Model  The bounding surface formulation is incorporated in this section into the reference model described in section 3.1. The resultant model (SANICLAY-D-B) is first presented in the triaxial plane and then generalized to multiaxial space. It is pertinent to say that the bounding surface approach determines only the plastic modulus different from the one in the reference model, while other constitutive relations are unchanged. The plastic modulus then requires consideration of extra complexities associated with the bounding surface formulation. It might be necessary to incorporate additional complexities into the proposed for better model predictions.  3.2.1  Formulation in the triaxial space  Similar to section 3.1, triaxial stress and strain quantities are used in formulation of the proposed model.  43  3.2. Proposed Model Bounding surface, plastic potential and loading surface In the SANICLAY-D-B model a bounding surface replaces the yield surface that served as a consolidation history in SANICLAY-D. The bounding surface is then given by ( ) F = (¯ q − p¯α)2 − N 2 − α2 p¯ (p0 − p¯) = 0  (3.39)  where p0 is a size of the bounding surface, (¯ p, q¯) are image stresses on the bounding surface (bar over symbols implies relevance to the bounding surface) to be explained in the sequel of the formulation, and α and N are the same as in the reference model. The plastic potential passes through the image stresses (¯ p, q¯) and is expressed as ( ) G = (¯ q − p¯α)2 − M 2 − α2 p¯ (pα − p¯)  (3.40)  In all the relations involved in the reference model, f = 0 and g = 0 are replaced by F = 0 and G = 0, respectively, to render the proposed SANICLAY-D-B model. A loading surface f = 0 is a surface passing through current stresses (p, q) and is homologous to the bounding surface with a projection center (pc , qc ) being the center of homology (to be defined in the sequel of the formulation), as shown in Fig. 3.2. The analytical expression of the loading surface is not required for SANICLAY-D-B, but it is helpful to visualize the loading surface in order to understand the loading, unloading and reloading processes explained later. Image stress The bounding surface formulation postulates that for any stress (p, q) a unique image stress (¯ p, q¯) exists on the bounding surface (Dafalias, 1986b). The image stress (¯ p, q¯) is obtained by a radial mapping of a stress point (p, q) from a projection center (pc , qc ) onto the bounding surface as shown in Fig. 3.2. Given a projection center (pc , qc ), image stresses are expressed as p¯ = pc + b (p − pc )  (3.41)  44  3.2. Proposed Model  q  Mc  N ∇F  (¯ p, q¯)  ∇f  α  (p, q)  (pc , qc )  p  Loading surface (f=0) Bounding surface (F=0)  Me  N  Plastic potential (G=0)  Figure 3.2: Proposed model schematic diagram in the triaxial plane  q¯ = qc + b (q − qc )  (3.42)  where b is a similarity ratio between the loading surface f = 0 and the bounding surface F = 0. The value b can be obtained by plugging Eqs. 3.41 and 3.42 into Eq. 3.39. This is solved term by term as follows: (¯ q − p¯α)2 = b2 A1 + bB1 + C1  (3.43)  45  3.2. Proposed Model where A1 = [(q − pα) − (qc − pc α)]2  (3.44)  B1 = 2q [(q − pα) − (qc − pc α)] (qc − pc α)  (3.45)  C1 = (qc − pc α)2  (3.46)  Similarly, p¯(p0 − p¯) = −b2 A2 − bB2 − C2  (3.47)  where A2 = (p − pc )2  (3.48)  B2 = (2pc − p0 ) (p − pc )  (3.49)  C2 = pc (pc − p0 )  (3.50)  Then a quadratic equation is obtained from Eqs. 3.43 and 3.47 Ab2 + Bb + C = 0  (3.51)  where ( ) A = A1 / N 2 − α2 + A2 ( ) B = B1 / N 2 − α 2 + B2 ( ) C = C1 / N 2 − α 2 + C2  (3.52) (3.53) (3.54)  The positive root corresponds to the image stress (¯ p, q¯) projected through the stress point (p, q) from the projection center (pc , qc ). Projection center Determination of the image stress on the bounding surface requires the projection center. The projection center at the last point of stress reversal can be used to simulate plastic 46  3.2. Proposed Model strains in unloading, required for better undrained cyclic loading simulations. The issue of simulations plastic strains in unloading during cyclic loading was discussed by Zienkiewicz et al. (1985). Jung and Yune (2011) showed from stress probing tests that simulations with the projection at the last stress reversal are favorable for monotonic loadings. The model proposed in this thesis, thus, considers projection center (pc , qc ) established at the last point of stress reversal. Geometrically, stress reversal is considered in the proposed model any time stress increments (p, ˙ q) ˙ are pointing into the interior of the loading surface (see Fig. 3.2). The projection center (pc , qc ) is updated at the moment of the stress reversal. However, if the projection center (pc , qc ) is simply chosen to be a fixed point until the next stress reversal occurs, then there could be a case such that the projection center stays outside the bounding surface. This is true for the reference model which adopts isotropic and rotational hardening of the yield surface (the bounding surface in the proposed model) relatively to the origin. The projection center outside of the bounding surface violates the concept of the bounding surface suggesting the uniqueness of the image stress for a given projection center and stress state. Thus, the idea of the projection center at the last stress reversal should be modified for the proposed model to maintain the desired effect of plastic strains upon unloading and, at the same time, to keep the well-established bounding surface concept. The proposed idea is to maintain the relative position of the projection center while the bounding surface evolves. Consider (pc , qc ) as an updated projection center at the moment of stress reversal. If one requires that the relative position of the projection center is kept unaltered owing to only a change of the size of the bounding surface, then the projection center (pc , qc ) should evolve accordingly: p˙c =  p˙0 pc p0  (3.55)  q˙c =  p˙0 qc p0  (3.56)  The relative position of the projection center (pc , qc ) should also be maintained upon rotation and distortion of the bounding surface. For this to be in effect, the following ratio 47  3.2. Proposed Model should be kept constant:  qc − qα =X qb − qα  (3.57)  where, referring to Fig. 3.3, qα = pc α; qb is the closest to the projection center point on the bounding surface for a given pc and X is a constant ratio. q  N  qb  Projection center  qc α qα  pc  p  Bounding surface (F=0) N  Figure 3.3: Proposed projection center The denominator of Eq. 3.57 can be obtained from F = 0 in Eq. 3.39 as qb − qα =  [( 2 ) ]1/2 N − α2 pc (p0 − pc )  (3.58)  Then the rate of qc due to rotation and distortion of the bounding surface is obtained based on Eq. 3.57, which yields [ q˙c = pc − X  pc (p0 − pc ) α [(N 2 − α2 ) pc (p0 − pc )]1/2  ] α˙  (3.59)  Notice that the first term of Eq. 3.59 considers the change of qc due to rotation of the bounding surface without distortion, while the second term adds the effect of distortion 48  3.2. Proposed Model induced by the bounding surface rotation. The two terms in conjunction with each other guarantee that the ratio X is maintained constant and the projection center is never outside the bounding surface. The rate of the projection center evolution due to isotropic, rotational and distortional hardening is expressed by combining Eqs. 3.55, 3.56, and 3.59: p˙c =  p˙0 pc p0  (3.60)  [ ] p˙0 pc (p0 − pc ) α q˙c = qc + pc − X α˙ p0 [(N 2 − α2 ) pc (p0 − pc )]1/2  (3.61)  Loading index determining the stress reversal The loading index L is the same as that given in the reference model by Eq. 3.15 and is used in the proposed SANICLAY-D-B model to define the point of the stress reversal. The stress reversal is assumed to take place whenever L = 0. The projection center is then updated to the point of the stress reversal. The use of L = 0 rather than L < 0 is required to update the projection center to the point of initiation of shearing preceded by unloading in isotropic consolidation. The neutral loading in the bounding surface formulation corresponds to the stress path sliding along the loading surface, which renders L = 0 and as a result zero plastic strains. With the proposed update of the projection center the possibility of the neutral loading is eliminated. However, the stress path with L → 0+ would correspond to almost neutral loading which does not involve update of the projection center. Plastic modulus For the stress states on the bounding surface, the plastic modulus is obtained from the consistency condition F˙ = 0 as ( ¯p = − K  ∂F ∂F p¯0 + α ¯ ∂p0 ∂α  ) (3.62)  49  3.2. Proposed Model Then, following Dafalias (1986b), the plastic modulus Kp for any stress state (p, q) is assumed ¯ p at the image point (¯ to be a function of plastic modulus K p, q¯) and the Euclidian distances between (p, q) and (¯ p, q¯) denoted by δ and between (pc , qc ) and (¯ p, q¯) denoted by r: 3 h p30 ¯ p + h p0 δ = K ¯p + Kp = K ⟨r − sδ⟩ ⟨b/(b − 1) − s⟩  (3.63)  where h is a positive hardening function, r/δ = b/(b − 1) with 1 ≤ b ≤ ∞ defined previously, and s is an indirect measure of the radius of the elastic nucleus, which is homologous to the bounding surface with the projection center as the center of homology. The presence ¯ p , i.e., stress to the of p0 is required to maintain similar units for the second term as for K third power. For s = ∞, the elastic nucleus coincides with the bounding surface, implying that the stress state inside the bounding surface considers purely elastic deformations. The elastic nucleus degenerates to a point with s = 1. Splitting the total plastic modulus into its image part and an additional term is a key feature of the bounding surface formulation. The second term in Eq. 3.63 controls Kp in the following way. Within the bounding surface where b >> 1, the term takes large values ¯ p , also resulting in high Kp ; as b → 1, the term assumes smaller values and relative to K ¯ p . Thus, the plastic modulus for the stress states within the bounding surface Kp → K decreases from infinite for b >> 1 to finite values, approaching the plastic modulus at image stress as b → 1. This allows a smooth transition from purely elastic strains to elastoplastic strains on the bounding surface. The default value of the elastic nucleus is chosen to be s = 1 for the current formulation. One could argue that it is necessary to have an elastic region at the initiation of a new loading process following reversal of stresses. However, this could be readily achieved by the model to a certain extent, since the projection center at the last stress reversal renders b → ∞ at initiation of a new loading process, resulting in Kp → ∞, which is equivalent to an initially purely elastic response. Hardening function Zienkiewicz et al. (1985) proposed a cumulative plastic deviatoric strain to induce reduction of plastic modulus Kp , as the use of only the distance dependent Kp was not sufficient for  50  3.2. Proposed Model undrained cyclic triaxial loading simulations. Similarly, in the proposed SANICLAY-D-B model hardening function h in Eq. 3.63 is assumed to be a decaying function of its initial value h0 and new internal variable d, simulating damage, h=  h0 1+d  (3.64)  As d starts evolving, h decreases from initial value h0 , leading to reduction of Kp . The proposed d evolution has a linear relationship with an increment of the plastic deviatoric strain ε˙pq , d˙ = ad ε˙pq  (3.65)  where ad is a new model parameter that controls the rate of d evolution. Internal variable d evolves any time plastic deviatoric strains ε˙pq occur.  3.2.2  Formulation in the multiaxial space  Generalization of the proposed model is performed following the observation given by Eq. 3.19: q¯ = [(3/2) ¯ s: ¯ s]1/2  (3.66)  where, as for triaxial quantities, the bar over bold symbols implies the relevance to F = 0. Bounding surface and plastic potential The analytical expression for the bounding surface, given by Eq. 3.39, is generalized to ( ) 3 3 2 F = (¯ s − p¯α) : (¯ s − p¯α) − N − α : α p¯ (p0 − p¯) = 0 2 2  (3.67)  where N is assumed to be independent of Θ, as in the reference model. The plastic potential form in multiaxial space is given by ( ) 3 3 2 G = (¯ s − p¯α) : (¯ s − p¯α) − M − α : α p¯ (pα − p¯) = 0 2 2  (3.68)  where M is Lode angle dependent and is expressed by Eqs 3.26a and 3.26b in the reference 51  3.2. Proposed Model model with ¯r substituted for r. The bounding surface and plastic potential, corresponding to the triaxial stress state shown in Fig. 3.2, are plotted in the stress space in Figs. 3.4(a–b) and 3.4(c–d), respectively. The bounding surface is a rotated and distorted ellipsoid with the anisotropic tensor α being a central axis of the ellipsoid. Observe that N -cone is formed in the multiaxial space, as the parameter N is chosen to be independent of Lode angle in the proposed formulation. The N -cone is fixed in the stress space with the space diagonal being a central axis of the cone. In the triaxial plane the gradients of the bounding surface at the intersection of the bounding surface with the stress ratio N in compression and extension are orthogonal to the p-axis. Similarly, the gradients of the bounding surface at the intersection of the bounding surface with N -cone are orthogonal to the stress space diagonal (σ11 = σ22 = σ33 ) in the stress space. The plastic potential passes through the image point on the bounding surface (see Fig. 3.2). Similar to the bounding surface, rotation of the plastic potential is controlled by the anisotropic tensor α, which is also a central axis of the plastic potential. In the stress space, the shape of the plastic potential is not similar to the shape of the bounding surface as the critical state failure envelope, chosen to be Lode angle dependent, controls the shape of the plastic potential. Critical state failure envelope is fixed in the stress space with the stress space diagonal being its central axis. The gradients of the plastic potential at the intersection of the plastic potential with the critical state failure envelope are orthogonal to the stress space diagonal. Image stress Given projection center tensor σ c can be decomposed into its hydrostatic pc I and deviatoric sc counterparts, multiaxial generalizations of Eqs. 3.41 and 3.42 are given by p¯ = pc + b (p − pc )  (3.69)  ¯ s = sc + b (s − sc )  (3.70)  The value of b is determined by plugging Eqs. 3.69 and 3.70 into Eq. 3.67, which in 52  3.2. Proposed Model  (a)  (b)  (c)  (d)  Figure 3.4: Multiaxial representation of: (a–b) the bounding surface with N -cone, and (c–d) plastic potential with critical state failure envelope  53  3.2. Proposed Model multiaxial space is solved term by term as follows: 3 (¯ s − p¯α) : (¯ s − p¯α) = b2 A1 + bB1 + C1 2  (3.71)  3 3 (s − sc ) : (s − sc ) − 3 (p − pc ) (s − sc ) : α + (p − pc )2 α : α 2 2  (3.72)  where A1 =  B1 = 3 (s − sc ) : sc − 3 (pc s + (p − 2pc ) sc ) : α + 3pc (p − pc ) α : α  (3.73)  3 3 C1 = sc : sc − 3pc sc : α + p2c α : α 2 2  (3.74)  Eq. 3.47 is unchanged in multiaxial generalization. Then a generalized form of Eq. 3.51 is given by: Ab2 + Bb + C = 0  (3.75)  where ( ) A = A1 / N 2 − (3/2) α : α + A2 ( ) B = B1 / N 2 − (3/2) α : α + B2 ( ) C = C1 / N 2 − (3/2) α : α + C2  (3.76) (3.77) (3.78)  Similar to the model formulation in triaxial plane, positive root is of interest for determination of the image stress given by p¯ and ¯s. Projection center Projection center evolution due to only the change in size of the bounding surface is given in multiaxial space by p˙0 pc p0 p˙0 s˙ c = sc p0  p˙c =  (3.79) (3.80)  54  3.2. Proposed Model and pc evolution is unchanged in generalization. The multiaxial stress space equation equivalent to Eq. 3.57 is (sc − sα ) : (sc − sα ) =X (sb − sα ) : (sb − sα )  (3.81)  where sb and sα = pc α are deviatoric tensors generalizing qb and qα , respectively. sb is determined along unit tensor nc , which is expressed as nc = (rc − α) / [(rc − α) : (rc − α)]1/2  (3.82)  where rc = sc /pc . Thus, it is necessary to find a scalar-valued multiplier ab to solve for sb sb = sα + nc ab  (3.83)  ab can be determined by plugging Eq. 3.84 into Eq. 3.67 yielding [  (N 2 − (3/2) α : α) pc (p0 − pc ) ab = (3/2) nc : nc  ]1/2 (3.84)  Given ab , Eq. 3.81 can be used to give the rate of ¯sc due to rotational and distortional hardening of the bounding surface: [ s˙ c = pc − X  pc (p0 − pc ) [(3/2) α : α]1/2 [(N 2 − (3/2) α : α) pc (p0 − pc )]1/2  ] ˙ α  (3.85)  The complete evolution of the projection center in multiaxial stress space is then given by a combination of Eqs 3.80 and 3.85 for sc : p˙c =  p˙0 pc p0  [ ] p˙0 pc (p0 − pc ) [(3/2) α : α]1/2 ˙ α s˙ c = sc + pc − X p0 [(N 2 − (3/2) α : α) pc (p0 − pc )]1/2  (3.86) (3.87)  55  3.3. Sensitivity analysis on new model constants Loading index and plastic modulus The loading index L in multiaxial space is the same as that given in the reference model by Eq. 3.37. A plastic modulus corresponding to the stress states on the bounding surface is generalized to  ( ¯p = − K  ) ∂F ∂F ¯ p¯0 + α ∂p0 ∂α  (3.88)  ¯ p , the plastic modulus Kp for any stress state Given the multiaxial generalization of b and K is of the same form given in Eq. 3.63: 3 h p30 ¯ p + h p0 δ = K ¯p + Kp = K ⟨r − δ⟩ ⟨b/(b − 1) − s⟩  (3.89)  Hardening function For the generalization of hardening function h, only d evolution must be modified for multiaxial space: h=  h0 1+d  d˙ = ad [(2/3)˙ep : e˙ p ]1/2  (3.90)  (3.91)  where e˙ p is a plastic deviatoric strain increment tensor. SANICLAY-D-B model generalization is complete at this stage.  3.3  Sensitivity analysis on new model constants  The proposed model has two constants corresponding to the bounding surface. Constant h0 is an initial value of h used to determine the plastic modulus within the bounding surface. Constant ad regulates the rate of evolution of the newly introduced damage parameter d. Sensitivity analyses on h0 and ad are performed in this section with the set of model constants and initial internal variables shown in Tables 3.1 and 3.2, respectively. Destructuration is not considered for illustration purposes.  56  3.3. Sensitivity analysis on new model constants  Table 3.1: Model parameters for sensitivity analysis Model constant Designation Elasticity κ ν Critical state Mc Me N λ Rotational hardening C x Destructuration (no effect) ki Bounding surface h0 ad  Value 0.03 0.2 1 1 1 0.15 5 1.7 0 * *  * These constants are to be analyzed for sensitivity.  Table 3.2: Initial internal variables for sensitivity analysis Model internal variable Designation Size of the bounding surface p0 Initial damage parameter d Structuration factor Si Initial void ratio ein Orientation of plastic potential and bounding surface α Initial stress state (p, q)  3.3.1  Value 200 kPa 0 1 0.7 0 (200, 0) kPa  Parameter h0  The mapping rule employed in the model formulation involves a plastic modulus dependent on a relative distance from the current stress state to the image stress. Eq. 3.89 shows that the plastic modulus is also controlled by hardening function h, which is determined by parameters h0 and d. In this section stress and stress–strain path variations with altered values of h0 are illustrated; ad = 0 is adopted so that h = h0 throughout imposed loading. Simulations of undrained triaxial compression starting from an isotropically consolidated state (CIUC) test on overconsolidated clay are performed. Undrained monotonic triaxial loading is simulated for overconsolidation ratio OCR = 1.5 and 3 with the results shown in Figs 3.5(a–b) and 3.5(c–d), respectively. The projection center is updated to the points of stress reversals. The corresponding image point during the simulations starting from OCR = 1.5 renders a gradient dictating contractive response. For lower h0 values, the 57  3.3. Sensitivity analysis on new model constants second term of Eq. 3.89 renders lower Kp , and more plasticity is involved in the process, i.e., a higher degree of contraction and reduction of p. With increased h0 , Kp is increased and less contraction is simulated. The projection center for simulations starting from OCR = 3 dictates the gradient of the plastic potential imposing negative plastic volumetric strains, i.e., dialative response. The dialation results in increased p, as observed in simulations shown in Fig. 3.5(c–d). A lower h0 provides a stress path that is more favorable to increasing p than that with a higher h0 . 0.8  0.4  q/p0,in  q/p0,in  240 0.6 120  240 120  0.6  60 30 h0=15  0.2 0 0  0.8  OCR=1.5  0.2  0.4  0.6 0.8 p/p0,in  0.4 0.2  1  0 0  1.2  60 30 h0=15 2  4  10  6  8  10  (b) OCR=3  240 120  0.6 240 120 q/p0,in  q/po,in  8  0.8  0.8  0.4 60 30 h0=15  0.2 0 0  6  a  (a)  0.6  ε (%)  0.2  0.4  0.6 0.8 p/po,in (c)  0.4 60 30 h0=15  0.2  1  1.2  0 0  2  4  εa (%)  (d)  Figure 3.5: Parameter h0 sensitivity in undrained monotonic triaxial loading (ad = 0, p0,in : initial size of the bounding surface): (a–b) OCR = 1.5, and (c–d) OCR = 3 For illustration of h0 in cyclic loading, 6 cycles of undrained cyclic triaxial loading are simulated with four different h0 values starting from the normally consolidated state and are shown in Fig. 3.6. As in monotonic loading, a lower plastic modulus is involved for lower h0 values, resulting in a more contractive undrained stress path and larger axial strains. 58  3.3. Sensitivity analysis on new model constants The stress path migration towards the origin is not observed once OCR = 2 is reached since the projection center invokes image stresses at the critical state ratio where no plastic volumetric strains are involved. On the basis of the volumetric counterpart of Eq. 3.23, zero plastic volumetric strain correspond to p˙ = 0. The consequence of the imposed stress path and constant h0 is that the stress–strain curve traces an unchanged path. To have a systematic calibration approach, h0 can be determined with ad = 0 based on initial cycles. These cycles should not involve considerable accumulation of plastic deviatoric strains, since this would require consideration of damage parameter d as well. Model users should be aware of small changes in the simulated response during initial cycles once ad ̸= 0 is utilized and consider this while calibrating h0 .  3.3.2  Parameter ad  To demonstrate the effect of the rate ad , loading involving considerable accumulation of plastic deviatoric strains is required. Therefore, only cyclic loading simulations are considered for sensitivity analysis on ad . 6 cycles of undrained cyclic triaxial loading starting from the normally consolidated state are imposed in simulations with different values of ad . h0 = 400 is considered for the simulations to achieve a wider range of h reduction during the course of applied loading. Within the initial cycles little difference is observed in the results shown in Fig. 3.7, since accumulation of plastic deviatoric strain εpq is negligible. As the loading proceeds, damage parameter d evolves, with the rate determined by parameter ad . For ad = 0, the stress path reaches the state where it traces the same path in subsequent cycles, as also shown in h0 sensitivity simulations, which results in the stress–strain tracing an unchanged path as well. As ad is increased from zero, damage parameter d evolves during shearing, leading to continuous reduction of h. Observe that for ad > 0 the stress–strain curves keep evolving, with increasing axial strains in both sides even when the corresponding stress paths has reached equilibrium. The rate of axial strains development is controlled by the parameter ad . This illustrates the capability of the current model to simulate cyclic softening. Parameter h0 should be determined prior to calibration of ad . If only monotonic loading is of interest, then the damage parameter effect might not be required and one can set ad = 0. 59  0.4  0.4  0.2  0.2 q/p0,in  q/p0,in  3.3. Sensitivity analysis on new model constants  0  0 −0.2  −0.2 h =15  −0.4 0  0  0.2  0.4  0.6 0.8 p/p0,in  1  −0.4 −4 −3 −2 −1  1.2  0.4  0.4  0.2  0.2  0 −0.2 −0.4 0  4  2  3  4  2  3  4  −0.2 h0=30  0.2  0.4  0.6 0.8 p/p0,in  1  −0.4 −4 −3 −2 −1  1.2  0 1 εa (%)  (d)  0.4  0.4  0.2  0.2 q/p0,in  0,in  q/p  3  0  (c)  0 −0.2 −0.4 0  2  (b)  q/p0,in  q/p0,in  (a)  0 1 εa (%)  0 −0.2  h0=60  0.2  0.4  0.6 0.8 p/p0,in  (e)  1  1.2  −0.4 −4 −3 −2 −1  0 1 εa (%)  (f)  Figure 3.6: Parameter h0 sensitivity in undrained cyclic triaxial loading (ad = 0, p0,in : initial size of the bounding surface): (a–b) h0 = 15, (c–d) h0 = 30, (e–f) h0 = 60  60  0.4  0.4  0.2  0.2 q/p0,in  q/p0,in  3.3. Sensitivity analysis on new model constants  0 −0.2 −0.4 0  0 −0.2  ad=0  0.2  0.4  0.6 0.8 p/p0,in  1  −0.4 −4 −3 −2 −1  1.2  2  3  4  2  3  4  2  3  4  (b)  0.4  0.4  0.2  0.2 q/po,in  q/po,in  (a)  0 1 εa (%)  0  0 −0.2  −0.2 a =15  −0.4 0  d  0.2  0.4  0.6 0.8 p/po,in  1  −0.4 −4 −3 −2 −1  1.2  (d)  0.4  0.4  0.2  0.2 q/p0,in  q/p  0,in  (c)  0 −0.2 −0.4 0  0 1 εa (%)  0 −0.2  ad=30  0.2  0.4  0.6 0.8 p/p0,in  (e)  1  1.2  −0.4 −4 −3 −2 −1  0 1 εa (%)  (f)  Figure 3.7: Parameter ad sensitivity in undrained cyclic triaxial loading (h0 = 400, p0,in : initial size of the bounding surface): (a–b) ad = 0, (c–d) ad = 15, (e–f) ad = 30  61  4. Qualitative Model Performance SANICLAY-D-B model formulated in the preceding Chapter is qualitatively compared with typical experimental response of clay. The advantage of the proposed projection center is demonstrated. Then the model is compared with the previous SANICLAY family models in undrained cyclic triaxial loading and circular stress path in the deviatoric plane. The proposed model performance is also illustrated in generation of the plots often used in practice for determination of dynamic characteristics of clay.  4.1  The advantage of the proposed projection center  In this section the projection center on anisotropic variable α is qualitatively compared with the projection center proposed in this thesis. Undrained cyclic triaxial loading simulations with the stress amplitude qcy /p0,in = 0.39, where qcy is an amplitude of applied cyclic loading and p0,in is an initial size of the bounding surface, are performed for this purpose. For a reference to the reader, typical experimental results from undrained cyclic triaxial test performed by Sheu (1984) are shown in Fig. 4.1. The projection center on rotational variable α was initially proposed by Anandarajah and Dafalias (1986) and used by many researchers for monotonic and cyclic loading simulations. It requires the use of the following statement for the projection center in terms of triaxial quantities: pc = Cp p0  (4.1)  qc = p c α  (4.2)  where Cp is a projection center parameter relating the position of the projection center on  62  4.1. The advantage of the proposed projection center  0.6  q /p  experiment (15 cycles)  cy  =0.39  0,in  0.4  q/p  0,in  0.2 0 −0.2 −0.4 −0.6 −10 −8 −6 −4 −2  0 2 εa (%)  4  6  8  10  Figure 4.1: Typical stress–strain path in undrained cyclic triaxial loading test after Sheu (1984) α to the size of the bounding surface p0 . The projection center is illustrated in the triaxial plane in Fig. 4.2. Generalization of Eqs. 4.1 and 4.2 to multiaxial space is given by σ c = Cp p0 (α + I)  (4.3)  The results of 15 cycles of simulations of undrained cyclic triaxial loading with the projection fixed on anisotropic variable α are shown in Fig. 4.3. The simulations start from an isotropically normally consolidated state with h0 = 50 and different values of Cp . The values of ad are adjusted to obtain approximately similar size of the stress–strain loops and are shown on each plot. Other model constants and internal variables are the same as in section 3.3. Parameter Cp = 0, used for simulations shown in Fig.4.3(a), renders the projection center fixed at the origin. In this case the loading surface passes through the origin and the stress point. Observe that the simulated stress–strain curve involves elastic strains in unloading. Recall that unloading is defined for the stress increments pointing towards the interior of the loading surface. No plastic strains are expected in unloading according to Eqs. 3.24 and 3.37 since the projection center fixed on anisotropic variable α invokes the stress increments pointing towards the interior of the loading surface. The plastic strains resume once the stress increments start pointing outwards the loading surface. In comparison with Fig. 4.1, simulated stress-strain path renders highly overdamped response. With 63  4.1. The advantage of the proposed projection center  q Mc  N ∇F (¯ p, q¯)  ∇f α  (p, q)  (pc , qc ) Cp p0 p0  p  Loading surface (f=0)  N Me  Bounding surface (F=0) Plastic potential (G=0)  Figure 4.2: Projection center on rotational variable α in the triaxial plane Cp = 0.6 simulations shown in Fig.4.3(b) the projection center is shifted to p = 0.6 p0 on anisotropic variable α. In this case, for a given stress state , the loading surface is of smaller size than that in simulations with Cp = 0 owing to the center of homology shifted from the origin. This results into reduced elastic part of the stress-strain path. However, overdamping is still present while comparing the simulated response with typical experimental results shown in Fig. 4.1. Similar simulation of undrained cyclic triaxial loading is performed with the projection center proposed in this thesis and the results are shown in Fig. 4.4. The value of ad = 25 is chosen to obtain stress–strain loops of comparable size with the loops in simulations results obtained with the projection center on anisotropic variable α. All other parameters are the same as before. Recall from the model formulation that the projection center is updated to the point of the stress reversal at the instant of stress increment pointing towards the interior of the loading surface. Then new loading surface starts evolving with the updated projection center from a point at the last stress reversal and stress increments point outwards the loading surface. Thus, the projection center at the last stress reversal guarantees that 64  4.1. The advantage of the proposed projection center  0.6  q /p  Projection center 0.4 on α  cy  0.6  =0.39  cy  =0.39  0,in  0.2 q/po,in  0.2 q/po,in  q /p  Projection center 0.4 on α  0,in  0 −0.2  −0.6 −10 −8 −6 −4 −2  −0.2  Cp=0 h0=50 ad=75  −0.4 0 2 εa (%)  4  6  8  0 Cp=0.6 h0=50 ad=28  −0.4 −0.6 −10 −8 −6 −4 −2  10  (a)  0 2 εa (%)  4  6  8  10  (b)  Figure 4.3: Stress–strain path simulations in undrained cyclic triaxial loading with the projection center on rotational variable α: (a) Cp = 0, and (b) Cp = 0.6 plastic strains will occur almost immediately after the reversal, although with a relatively smaller magnitude that increases as the loading proceeds. The results obtained suggest that the problem of overdamping is eliminated with the projection center proposed in this thesis.  0.6 Proposed projection 0.4 center  qcy/p0,in=0.39  q/p  o,in  0.2 0 −0.2 h =50 0 a =25  −0.4 −0.6 −10 −8 −6 −4 −2  d  0 2 εa (%)  4  6  8  10  Figure 4.4: Stress–strain path simulations in undrained cyclic triaxial loading with the projection center proposed in this thesis  65  4.2. Undrained cyclic triaxial loading  4.2  Undrained cyclic triaxial loading  SANICLAY developed by Dafalias (1986a); Dafalias et al. (2006) is a model of Cam-Clay class with incorporated stress-induced anisotropy. Its extension, SANICLAY-D, proposed by Taiebat, Dafalias and Peek (2010) can be successfully used to simulate clay destructuration. So far, SANICLAY family models have been used only for monotonic loading simulations. The proposed SANICLAY-D-B model is a further extension, enabling one to simulate undrained cyclic loading. Its qualitative performance is demonstrated in this section in undrained cyclic triaxial loading with the constant stress amplitude. Simulations are first performed with the previous SANICLAY family models. A summary of model simulations and corresponding model abbreviations are given in Table 4.1. Using SANICLAY-D-B, one can adjust model parameters shown in the table to reproduce the results that would correspond to previous SANICLAY family models. The model parameters and initial internal variables not shown in the table are given previously in Tables 3.1 and 3.2, respectively. Table 4.1: Summary of model parameters undrained cyclic triaxial loading Abbreviation SANICLAY SANICLAY-D SANICLAY-D-B SANICLAY-D-B (with damage)  and initial internal variables for simulations of ki 0 2 0 0  Si 1 3 1 1  h0 10100 10100 60 60  ad 0 0 0 25  The first two simulation results of initial 6 cycles of undrained cyclic triaxial loading are obtained with the original SANICLAY and SANICLAY-D models, shown in Figs. 4.5(a–b) and 4.5(c–d), respectively. During the initial quarter cycle the stress state is always on the bounding surface so that the plastic strains simulated by both models. Observe that SANICLAY-D model invokes larger plastic strains owing to the destructuration effect. Upon unloading, the stress path lies within the bounding surface rendering only elastic strains for both models until the stress path reaches the bounding surface again in extension. Note that if no rotational hardening had been considered, as in the Modified Cam-Clay model, the stress path would not have reached the bounding surface in extension. Once the stress path reaches the bounding surface, plastic strains resume in the simulations. As in the first quarter 66  4.2. Undrained cyclic triaxial loading cycle, SANICLAY-D model exhibits more plastic strains induced by the destructuration mechanism. For the subsequent cycles the stress path lies entirely within the bounding surface, leading to purely elastic response for both SANICLAY and SANICLAY-D models.  0.4  0.4  0.2  0.2  0  SANICLAY ki = 0 Si = 1 h0 = 10100 ad = 0  −0.2 −0.4 0  0.2  0.4  0.6 0.8 p/p0,in  1  q/po,in  q/po,in  Neither anisotropy nor destructuration induce further reduction of p.  0 −0.2 −0.4 −2  1.2  −1  0.4  0.4  0.2  0.2  0  SANICLAY-D ki = 2 Si = 3 h0 = 10100 ad = 0  −0.2 −0.4 0  0.2  0.4  0.6 0.8 p/p0,in  (c)  1  2  1  2  (b)  1  1.2  q/po,in  q/po,in  (a)  0 εa (%)  0 −0.2 −0.4 −2  −1  0 εa (%)  (d)  Figure 4.5: SANICLAY family models performance in undrained cyclic triaxial loading: (a–b) SANICLAY, and (c–d) SANICLAY-D Simulation results obtained with the proposed SANICLAY-D-B model are shown in Fig. 4.6. Destructuration is not considered in order to illustrate only the effect of the bounding surface. The results shown in Fig. 4.6(a–b) are obtained without damage effect. Notice the reduction of p within the bounding surface owing to allowed plastic strains. However, as the stress path reaches OCR = 2, the image stress projected from the projection center through the stress point falls onto the critical state ratio. As a result, neither compression nor dialation are involved, and the stress path traces similar path in each cycle. The h = h0 renders the plastic modulus purely dependent on the relative distance from the stress 67  4.3. Circular stress path in the deviatoric plane point to the image stress. As it was explained in section 3.3, this is accompanied by tracing unchanging stress–strain path. The damage effect in SANICLAY-D-B improves the model performance through reduction of hardening function h as can be seen in Fig. 4.6(c-d). Even though the stress path traces the same route at OCR = 2, the corresponding stress–strain path exhibits continuously reducing stiffness. In this way, SANICLAY-D-B with damage effect can be successfully used to simulate continuously evolving plastic strains in undrained  0.4  0.4  0.2  0.2  0 SANICLAY-D-B ki = 0 Si = 1 h0 = 60 ad = 0  −0.2 −0.4 0  0.2  0.4  0.6 0.8 p/p0,in  1  q/po,in  q/po,in  cyclic loading.  0 −0.2 −0.4 −10 −8 −6 −4 −2  1.2  0.4  0.4  0.2  0.2  0  SANICLAY-D-B ki = 0 Si = 1 h0 = 60 ad = 25  −0.2 −0.4 0  0.2  0.4  0.6 0.8 p/p0,in  (c)  4  6  8  10  4  6  8  10  (b)  1  1.2  q/po,in  q/po,in  (a)  0 2 εa (%)  0 −0.2 −0.4 −10 −8 −6 −4 −2  0 2 εa (%)  (d)  Figure 4.6: Comparison of model performance with bounding surface effect: SANICLAY-B and (c–d) SANICLAY-B with damage  4.3  (a–b)  Circular stress path in the deviatoric plane  Undrained cyclic triaxial loading tests are commonly used for model calibration to subsequently simulate soil response induced by earthquake loadings. However, in reality, earth68  4.3. Circular stress path in the deviatoric plane quakes involve a more complicated stress path than imposed in triaxial tests. A qualitative model simulation of strains induced by a circular stress path in the deviatoric plane is illustrated in this section.  4.3.1  Response to cyclic circular stress path loading  Model parameters and internal variables from section 3.3 are considered for the model simulations. Different combinations of h0 and ad , shown in Table 4.2, are used to reproduce the original SANICLAY and the proposed SANICLAY-D-B models with and without damage effect. Each simulation consists of three stages starting from normally consolidated state. In the first stage, isotropic unloading from p = 200 kPa to 150 kPa is induced, corresponding to the stress path A-B in Fig. 4.7. During the second stage, stress path B-C in Fig. 4.7, q is increased from 0 to 50 kPa in triaxial compression with constant p. The applied q is such that the stress point at the end of the second stage and subsequently applied cyclic stress path lie inside the bounding surface. For the last and main stage, the circular stress path in the deviatoric plane is applied clockwise, during which p and q are kept constant. To induce considerable plastic strains within the bounding surface, three cycles of the circular stress path are simulated. Table 4.2: Summary of model parameters for simulations of circular stress path in the deviatoric plane Abbreviation h0 ad SANICLAY 10100 0 SANICLAY-B 100 0 SANICLAY-B (with damage) 100 60 The results of the three simulations in the strain space are shown in Fig. 4.8 with the deviatoric view on the left and the shifted view on the right for better illustration of volumetric strains. Observe that the strain path followed in simulations with the original SANICLAY model in Fig. 4.8(a–b) traces unchanging circular strain path within the three cycles as the simulated response is completely elastic. The results of simulations with the effect of the bounding surface with and without damage effect are shown in Figs. 4.8(c–d) and 4.8(e–f), respectively. The two sets of results 69  4.3. Circular stress path in the deviatoric plane σ /σ 11  c  Space diagonal  Imposed stress path  C  σ /σ 11  C  A  c  B  σ /σ 33  A,B  σ /σ  σ22/σc  33  c  σ /σ 22  (a)  c  c  (b)  Figure 4.7: Isotropic unloading, constant-p triaxial compression, and circular stress path: (a) deviatoric view, (b) shifted view differ in that more plastic strains are involved with the considered damage effect. Notice that this time plastic strains are induced throughout the whole process of the applied loading. The imposed image stress dictates positive plastic volumetric strains as can be seen from the strain path in the shifted view. In the deviatoric view the strain path is slightly elongated along ε33 -axis during the first cycle. Such response is due to the distance dependence of the plastic modulus Kp . At the initiation of the circular stress path loading projection center is updated to the point C shown in Fig. 4.7. Then the relative distance from the current stress point to the image stress reduces during the first-half cycle, rendering increasing plastic strains. Shortly after the first half-cycle is complete, the projection center is updated and simulated plastic strains become initially negligible. Again the relative distance reduces as the loading proceeds and plastic strains increase until the projection center is updated shortly after the completion of the first cycle. In this way larger degree of plastic strains approximately after the completion of first half-cycle and first cycle results into the strain path somewhat elongated along ε33 -axis. Notice that update of the projection center in subsequent cycles occurs at later stages and the elongation of the strain paths is rotated clockwise. One could argue that it might be unrealistic for clays to exhibit larger plastic strains at particular points during circular stress path. However, the circular stress path imposes the  70  4.3. Circular stress path in the deviatoric plane ε  (%)  11  SANICLAY: h0=1e10, ad=0  ε  11  (%)  C  Space diagonal C  ε  33  ε22 (%)  ε  33  (%)  ε  (%)  22  (a)  (%)  (b) ε  11  (%)  SANICLAY−B: h0=100, ad=0  ε  11  (%) Space diagonal  C  ε33 (%)  C  ε22 (%)  ε  33  ε  (%)  22  (c)  (%)  (d)  ε11 (%) SANICLAY−B: h =100, a =60 0  d  ε  11  (%) Space diagonal  C  ε33 (%)  C  ε22 (%)  ε  33  (e)  ε  (%)  22  (%)  (f)  Figure 4.8: Observed strains in the strain space (deviatoric view on the left, shifted view on the right) induced by the circular stress path (clockwise) in the deviatoric plane: (a– b) original SANICLAY model, (c–d) proposed SANICLAY-D-B without damage, and (e–f) proposed SANICLAY-D-B with damage effect 71  4.3. Circular stress path in the deviatoric plane direction of the stress increments close to being orthogonal to the bounding surface gradient, which renders small values of L on the basis of Eq. 3.37. According to Eq. 3.24, the resulting plastic strains are also small. Therefore, even though plastic strains are abruptly reduced to zero upon the update of the projection center, the overall plastic strain magnitudes are not large enough to cause major deviations from expected clay response. It is more important that, unlike the original SANICLAY model, the bounding surface formulation enables one to simulate plastic strains for the stress states inside the bounding surface.  4.3.2  Effect of parameter C in circular stress path loading  Apparently the rate of the rotational hardening of the bounding surface can control the induced strain path to some extent during the applied circular stress path loading. Simulation results of one cycle of circular stress path loading similar to that discussed above but with different values of C, and fixed h0 = 50, and ad = 0 (i.e., no damage effect) are shown in Fig. 4.9. It should be noted that similar effect of C would apply for the stress states on the bounding surface. ε11 (%) SANICLAY−B: h0=50, ad=0  15 9 C=3  ε22 (%)  ε33 (%)  Figure 4.9: Effect of parameter C on induced strain path during circular stress path loading  72  4.4. Generation of secant shear modulus and damping ratio variation with shear strain amplitude Observe that, for lower values of C, the bounding surface does not rotate significantly and it has to expand more to accommodate the imposed stress path. The less rotated and more expanded bounding surface gives more freedom for the stress path to be directed towards the interior of the loading surface and the projection center is updated at earlier stages. Then the updated projection center invokes the gradients of plastic potential mainly on the left side of the bounding surface in the deviatoric view. As a result, the induced increments of strain tensor, given by a Eq. 3.24, and strain path are pointing mainly along ε33 -axis. The resulting strains at the end of a cycle are above initial strains. As value of C increases the update of the projection center occurs at later stages and the resulting strain path becomes rounded with the strain path below initial strains at the end of a cycle.  4.4  Generation of secant shear modulus and damping ratio variation with shear strain amplitude  Equivalent linear analysis is widely used in practice for site response analysis, where equivalent shear modulus and damping ratio for a given strain amplitude are required for representation of soil properties to simulate non-linear soil response. Variation of the representative soil properties with applied shear strain amplitude in cyclic loading is obtained from laboratory tests. In this section, the effect of the model constants h0 and ad on simulated variation of representative soil properties with a shear strain amplitude is studied. Simulations of undrained cyclic loading in Direct Simple Shear (DSS) test are performed with the constraints corresponding to an element at the center of a DSS device. Model constants and initial internal variables are given in Tables 3.1 and 3.2, respectively. In order to come up with representative soil properties, a hysteresis loop at 5th cycle is considered following Seed et al. (1986). A schematic diagram of a typical hysteresis loop is shown in Fig. 4.10. The hysteresis loop is represented by a closed loop as explained by Hardin and Drnevich (1972). The shear modulus is measured as a secant shear modulus Gsec determined by the extreme points on the hysteresis loop, while the damping is proportional to the area inside the loop (Seed et al., 1986). It is typical to normalize the secant shear  73  4.4. Generation of secant shear modulus and damping ratio variation with shear strain amplitude modulus Gsec by a maximum tangential shear modulus Gmax for illustration of results. The damping ratio can be obtained from commonly used equation Damping ratio (%) =  AL 100 4πAT  (4.4)  where AL is an area of a hysteresis loop; and AT is a stored equivalent elastic energy.  Figure 4.10: Definition of equivalent shear modulus and the quantities used for the determination of damping ratio Typical experimental data from Vucetic and Dobry (1991) on variation of secant shear modulus Gsec and damping ratio with shear strain γ is shown in Fig. 4.11 to the reader for qualitative comparison with the simulated results. The curves are obtained from the experiments on soils with a range of plasticity index (PI). The simulated variation of secant shear modulus Gsec with shear strain amplitude γ is shown in Fig. 4.12. Values of h0 = 50, 200, and 800 are adopted with ad = 0 and ad = 60 for the results shown in Figs. 4.12(a) and 4.12(b), respectively. Gsec starts reducing at γ = 0.1% for the simulations obtained with the original SANICLAY model (dash line). With the proposed model adopting ad = 0, Gsec starts reducing already at γ = 0.01%. Further reduction of Gsec is controlled by the parameter h0 , giving a lower Gsec for a given γ as h0 is decreased. Observe that as γ → 10% obtained values of Gsec tend to converge. This is true for any values of h0 adopted in the simulations as the stress path has to reach the bounding surface, where the plastic modulus and the degree of plastic strains reach their 74  4.4. Generation of secant shear modulus and damping ratio variation with shear strain amplitude 1  0.6  G  sec  /G  max  0.8  0.4 0.2 0 −4 10  PI=200 PI=100 PI=50 PI=30 PI=15 −3  10  −2  −1  10  10  0  10  1  10  (a)  60  Damping ratio (%)  50 40  PI=200 PI=100 PI=50 PI=30 PI=15  30 20 10 0 −4 10  −3  10  −2  −1  10 10 Shear strain (%)  0  10  1  10  (b)  Figure 4.11: Experimental data after Vucetic and Dobry (1988): (a) degradation of Gsec and (b) variation of damping ratio with γ limiting values. For the results obtained with ad = 60, Gsec at lower values of γ does not differ significantly from those ones obtained with ad = 0 as the damage parameter d is still not evolved sufficiently. The effect of d becomes more pronounced as γ is increased and the resulting Gsec variation curves shift down. The curves tend to converge for large values of γ as for simulations with ad = 0; however, larger γ needs to be achieved for the curves to converge. In comparison with the experimentally obtained Gsec degradation shown in Fig. 4.11(a), the results produced by SANICLAY-D-B model with damage effect are more plausible than those obtained with the original SANICLAY model. However, the proposed  75  4.4. Generation of secant shear modulus and damping ratio variation with shear strain amplitude 1 ad=0 SANICLAY  0.6  G  sec  /G  max  0.8  0.4 h =800 0  0.2  h =200 0  h0=50 0 −4 10  −3  10  −2  −1  10  10  0  1  10  10  (a)  1 a =60 d SANICLAY  0.6  G  sec  /G  max  0.8  0.4 h0=800 0.2  h0=200 h =50 0  0 −4 10  −3  10  −2  −1  10 10 Shear strain, γ, (%)  0  10  1  10  (b)  Figure 4.12: Simulations of Gsec /Gmax degradation with γ: (a) ad = 0 and (b) ad = 60 SANICLAY-D-B model is not able to replicate Gsec degradation at strains γ < 0.01% as mainly elastic strains are predicted by the model. It is interesting to note that the variation of Gsec could potentially serve as a basis for calibration of the bounding surface parameters h0 and ad . The former could be calibrated based on strain levels up to γ = 0.03%, and the latter based on larger strains. The variation of damping ratio with shear strain γ is illustrated in Fig. 4.13 with the similar values of h0 as before. The values of ad = 0 and ad = 60 correspond to the results shown in Figs. 4.13(a) and 4.13(b), respectively. SANICLAY model (dash line) predicts only elastic strains for γ < 0.4% and, as a result, no damping variation since the stress path 76  4.4. Generation of secant shear modulus and damping ratio variation with shear strain amplitude apparently does not reach the bounding surface at the 5th cycle. The damping ratio starts increasing abruptly at approximately γ = 0.4% and reaches about 60% of damping ratio at γ = 10%. The proposed model simulates the damping ratio that starts increasing at γ = (0.01–0.1) %, depending on the adopted value of h0 . The damping ratio initially evolves faster for lower values of h0 in the simulated results obtained with both ad = 0 and ad = 60. The rate of damping ratio evolution then drops down slightly for lower values of h0 , and the rate increases for higher values of h0 . The reason is a reduction of damping, represented by AL in Fig. 4.10, with decreasing value of h0 . This renders lower damping ratio on the basis of Eq. 4.4. The results obtained with ad = 60 show more significant drop down of the rate of the increase of damping ratio at larger γ as damage parameter d is considerably evolved at 5th cycles. Clearly, SANICLAY-D-B model with the damage parameter effect can better simulate experimentally obtained variation of damping ratio with shear strain γ, shown in Fig. 4.11(b), than the original SANICLAY model and SANICLAY-D-B model without the damage effect. Similar to Gsec degradation, the model cannot predict variation of damping ratio for γ < 0.01% owing to simulated elastic strains only. Results obtained above are based on the model simulations with the fixed set of model parameters, and their change may result in a different outcome. The parameters are typically calibrated based on experimental data, and each set corresponds to a particular soil. Some of the parameters can be related to particular soil properties. For example, parameter λ can be correlated to PI (Wood, 1990): λ ≈ 0.6PI  (4.5)  This renders approximate PI = 25% for λ = 0.15 used in the simulations. Referring to the above, the simulations can be approximately compared with the curves corresponding to PI = 30% in Fig. 4.11. It then follows that the model is able to successfully simulate the variation of the representative soil properties for γ > 1%. Satisfactory model performance could be considered for γ = (0.1–1) %, and much less satisfactory for γ < 0.01%. For practical purpose, strains around well-designed geotechnical structures, such as retaining walls, foundations, and tunnels, are generally small (Clayton, 2011), and might be considered of about γ = (0.01–1)%. However, for deformation analyses of dams and natural slopes, it  77  4.5. Generation of cyclic stress ratio versus number of cycles 60  Damping ratio (%)  h0=800 50  h =200  40  h0=50  SANICLAY  0  30 20 10 ad=0 0 −4 10  −3  10  −2  −1  10  10  0  10  1  10  (a)  60  Damping ratio (%)  50 40  h0=800  SANICLAY  h0=200 h0=50  30 20 10 ad=60 0 −4 10  −3  10  −2  −1  10 10 Shear strain, γ, (%)  0  10  1  10  (b)  Figure 4.13: Simulations of damping ratio variation with γ: (a) ad = 0 and (b) ad = 60 is important to consider larger strains (γ > 1%). Thus, the proposed model is adequate for most geotechnical problems.  4.5  Generation of cyclic stress ratio versus number of cycles  The effect of model parameters h0 and ad is considered in this section for generation of commonly used in practice plots of cyclic stress ratio (CSR) required to reach a particular  78  4.5. Generation of cyclic stress ratio versus number of cycles strain in a number of cycles. Simulations of undrained cyclic triaxial loading are performed for this purpose. CSR = τcy /σc = qcy /2σc is typically used for sands, and CSR = τcy /Su = qcy /2Su for clays, where σc is a consolidation stress and Su is an undrained shear strength used for normalisation. The latter CSR is adopted for simulations presented in this section with Su = 64 kPa determined from undrained triaxial compression test simulation starting from isotropically normally consolidated state. Plots are generated for a given strain, representing adopted failure criterion, and results are obtained for a particular CSR and required number of cycles to reach this strain. In this section plots are generated for CSR required to reach peak εa = 3%. Experimental data after Zergoun and Vaid (1994) is shown in Fig. 4.14 to the reader for qualitative comparison with the simulated results. 1.2  cy  Cyclic stress ratio, q /2 Su  experiment 1 0.8 0.6  Fitted curve  0.4 0.2 0 0  20  40 60 Number of cycles  80  100  Figure 4.14: Experimental data after Zergoun and Vaid (1994) on cyclic stress ratio CSR required to reach εa = 3% in a number of cycles Simulations are performed with the model parameters used in previous sections, except parameters h0 and ad . Each plot in Fig. 4.15 corresponds to simulations with values of h0 = 50, 200, and 800, and fixed ad . The range of CSR = (0.4–1) is considered in simulations. Only results within 100 cycles are used for illustration of qualitative model performance. The dash line corresponds to the original SANICLAY model simulations in Fig. 4.15. Observe that εa = 3% can be reached in one cycle only with CSR = 1, corresponding to the stress level equivalent to undrained shear strength Su. The εa = 3% cannot be attained with CSR < 1, rendering number of cycles → ∞ as suggested by the horizontal dash line. 79  4.5. Generation of cyclic stress ratio versus number of cycles  SANICLAY 1  cy  Cyclic stress ratio, q /2 Su  1.2  0.8 0.6 0.4  h =800 0  h =200  0.2  0  h0=50  a =0 0 0  d  20  40  60  80  100  (a)  SANICLAY 1  cy  Cyclic stress ratio, q /2 Su  1.2  0.8 0.6 0.4  h0=800  0.2  h0=200 h =50  a =60 d 0 0  0  20  40 60 Number of cycles  80  100  (b)  Figure 4.15: Generation of CSR required to reach εa = 3 % SANICLAY-D-B simulation results shown in Fig. 4.15(a) are obtained with ad = 0, rendering h = h0 throughout cyclic loading. Notice that simulations with different values of h0 fall very close to each other for CSR = 1 and 0.9. Only simulations with h0 = 50 attain εa = 3% within 100 cycles for CSR < 0.9, and none of the simulations fall within 100 cycles for CSR ≤ 0.6. In Fig. 4.15(b), simulations are performed with ad = 60, rendering h decreasing from its initial values of h0 as damage parameter d evolves. This time fewer cycles are required to reach εa = 3% for all simulations with the same h0 values as in Fig. 4.15(a). It should be noted that the simulated damage effect becomes more pronounced with larger number of cycles 80  4.6. Generation of peak axial strain and residual pore water pressure versus number of cycles The results simulated with the original SANICLAY model are clearly far from the experimental data shown in Fig. 4.14, while simulations with the proposed SANICLAY-D-B model significantly improve the results. One could argue that the use of the damage parameter d is not required as certain values of h0 and ad = 0 provide the results (Fig. 4.15(a)) that are close to those observed in experiments (Fig. 4.14). However, CSR versus number of cycles plot is generated only for a particular strain level, which neglects the strain history from before and after the strain level. Thus, different combination of h0 and ad can provide one with approximately similar CSR versus number of cycles plot, while at any other strain levels the simulated response would be completely different. In other words, a close qualitative match to experimental data obtained with, say, h0 = 800 and ad = 0 can be achieved with much higher h0 and ad ̸= 0, which might be required for better simulations at strain levels other than εa = 3%. From the above simulations, the effect of h0 and ad on the generation of CSR versus number of cycles can be summarized: • as h0 is decreased, simulated soil response becomes softer and a lower plastic modulus Kp is involved for the stress states within the bounding surface; this results into smaller number of cycles required to bring the simulated response to a certain strain level; • with increased ad , the damage parameter d evolves more rapidly as plastic deviatoric strains accumulate throughout cyclic loading; increased damage parameter d reduces h from h0 , and the plastic modulus Kp decreases, leading to larger plastic strains within the bounding surface.  4.6  Generation of peak axial strain and residual pore water pressure versus number of cycles  In this section effect of model parameters h0 and ad are demonstrated for generation of peak axial strain εa and residual pore water pressure u versus number of cycles. Simulations of undrained cyclic triaxial loading with qcy /2Su = 0.5 of applied cyclic stress amplitude are performed with the model parameters used in previous sections, except parameters h0 and 81  4.6. Generation of peak axial strain and residual pore water pressure versus number of cycles ad . Typical experimental results from undrained cyclic triaxial loading tests with τcy /Su = qcy /2Su = 0.69 reported by Zergoun and Vaid (1994) are shown in Fig. 4.16 as a reference for the reader. 10 experiment  8 qcy/2Su=0.69  Peak εa (%)  6 4 2 0 −2 −4 −6 −8 −10 0 10  1  2  10  10  (a)  1 experiment  q /2Su=0.69 cy  Residual u/σ  e  0.8 0.6 0.4 0.2 0 0 10  1  10 Number of cycles  2  10  (b)  Figure 4.16: Development of: (a) peak axial strain εa , and (b) residual pore water pressure u versus number of cycles, experimental data from Zergoun and Vaid (1994) The results obtained with the original SANICLAY model are shown with dash line in Fig.4.17. Observe that plastic strains are induced only during the initial loading as the stress state lies on the bounding surface and during the following reloading once the stress path hits the bounding surface in extension. In subsequent cycles, no plastic strains are simulated as suggested by the horizontal dash line. This is accompanied by residual pore water pressure 82  4.6. Generation of peak axial strain and residual pore water pressure versus number of cycles u accumulation only during simulated plastic strains, i.e., no further build-up of residual u. For the simulation results obtained with different values of h0 and fixed ad = 40 shown in Fig. 4.17, parameter ad ̸= 0 is chosen to induce accumulation of plastic strains for the stress path reaching equilibrium at OCR = 2. Observe that for lower values of h0 fewer number of cycles are required to induce a certain strain level, shown in Fig. 4.17(a). Accumulation of residual u in Fig. 4.17(b) can be interpreted to indicate when the stress path reaches equilibrium, where no further development of residual u is simulated. Clearly, the lower h0 values, the faster such equilibrium is reached since the plastic strains induce larger εpv resulting into faster p reduction. Note that no further peak εa would have been possible once OCR = 2 is reached, should ad = 0 had been utilized. In comparison with experimental data shown in Fig. 4.16, the simulations with the original SANICLAY model are not able to realistically replicate the observed response of clay. The proposed SANICLAY-D-B model with the evolving damage parameter d reproduces the observed response qualitatively well, although residual u is about twice smaller than observed in experiments owing to the equilibrium reached in simulations by the stress path at OCR = 2. In Fig. 4.18, the results are obtained with fixed h0 = 50 and different values of ad . Simulations results obtained with ad = 0, render h = h0 throughout the applied cyclic loading. The resultant peak εa develops until the stress path reaches the equilibrium at OCR = 2. For ad ̸= 0, the rate of damage parameter d evolution is larger with larger ad , which invokes faster reduction of plastic modulus Kp and fewer cycles are required to reach a certain strain level in Fig. 4.18(a). Observe in Fig. 4.18(b) that unlike simulations with different h0 values, variation of ad does not involve considerable difference in accumulation of residual u, because effect of ad is more pronounced with increasing number of cycles. The proposed SANICLAY-D-B model simulation with no effect from damage parameter d cannot replicate experimental results shown in Fig. 4.16, while the simulations with evolving damage parameter d can reproduce accumulation of peak εa . However, the proposed model predicts about twice smaller residual u than that observed in the experiments, because the stress path migration is not allowed once the stress path reaches OCR = 2. The results of qualitative model performance illustrated in this section demonstrate the 83  4.6. Generation of peak axial strain and residual pore water pressure versus number of cycles ability of the model to simulate cyclic response of clays by means of new model parameters h0 and ad . Parameter h0 is necessary to simulate undrained cyclic response within initial few cycles and thus accumulation of residual u, while parameter ad enables one to simulate reduction of stiffness with increasing number of cycles observed in experiments.  10 q /2Su=0.5  8  cy  Peak εa (%)  6  h =25 0  4  50  200  100  400 SANICLAY  2 0 −2 −4 −6  a =40 d  −8  −10 0 10  1  2  10  10  (a)  1 ad=40  q /2Su=0.5 cy  Residual u/σ  c  0.8 0.6  h0=25 50  0.4 0.2 0 0 10  100 200  400  SANICLAY  1  10 Number of cycles  2  10  (b)  Figure 4.17: Generation of: (a) peak axial strain εa , and (b) residual pore water pressure u versus number of cycles with different values of h0 and fixed ad = 40  84  4.6. Generation of peak axial strain and residual pore water pressure versus number of cycles  10 qcy/2Su=0.5  8  Peak εa (%)  6  a =80 d  40  10  20  4  5  2 0  2 0 −2 −4 −6  h0=50  −8  −10 0 10  1  2  10  10  (a)  1  h0=50  qcy/2Su=0.5  Residual u/σ  e  0.8 0.6 0.4  ad=80 0  0.2 0 0 10  1  10 Number of cycles  2  10  (b)  Figure 4.18: Generation of: (a) peak axial strain εa , and (b) residual pore water pressure u versus number of cycles with fixed h0 = 50 and varied values of ad  85  5. Quantitative Model Performance This Chapter illustrates the proposed SANICLAY-D-B model performance in comparison with experimental data in monotonic and cyclic loading and is intended to serve as a guidance for model parameters calibration. Experiments required for calibration of model parameters are discussed. Then the model is tested quantitatively against experimental data on three different clays.  5.1  Experiments required for model calibration  Abbreviations of the tests and their description used in this section are given below: CIUC undrained triaxial compression test on isotropically consolidated clay CIUE undrained triaxial extension test on isotropically consolidated clay CK0 UC undrained triaxial compression test on K0 -consolidated clay CK0 UE undrained triaxial extension test on K0 -consolidated clay The proposed SANICLAY-D-B model requires 11 model constants, whose roles are summarized in Table 5.3. Constants h0 and ad are novel in regards to SANICLAY-D model proposed by Taiebat, Dafalias and Peek (2010). Initial internal variables required for the proposed model with their definition are shown in Table 5.2. The reader is referred to the works by Dafalias et al. (2006) and Taiebat, Dafalias and Peek (2010) for comprehensive description of anisotropy and destructuration related parameters calibration, respectively. A brief description of the calibration procedure for the required constants is given as follows: • Calibration of λ and κ is based on, preferably isotropic, consolidation tests. The use of data from the K0 -consolidation test is also acceptable. 86  5.1. Experiments required for model calibration • Calibration of ν is based on initial K0 -unloading stress path preceded by K0 -consolidation test. Typical range of ν = 0.2–0.3 for clays can be assumed, if the required stress path is not available. • Calibration of Mc and Me is based on triaxial compression and extension tests, respectively. • Calibration of N is based on CK0 UC test on normally consolidated clay. If not available, N can be calibrated along with C from CIUC or CIUE test on normally consolidated clay. • Calibration of C is based on the stress path from CK0 UE test. If not available, C can be calibrated along with N from CIUC or CIUE test on normally consolidated clay. C = 3–20 is the typical range for most clays (Dafalias et al., 2006). • Calibration of x is based on the closed-form relation given by Eq. 27 in Dafalias et al. (2006). • Calibration of ki for structured clays is based on consolidation test to stresses, ensuring complete loss of structure. ki = 0 should be used for unstructured clays. • Calibration of h0 is based on stress and stress–strain paths from initial cycles of undrained cyclic loading, preferably on normally consolidated clay. If monotonic loading is of interest, the calibration can be performed in comparison with undrained monotonic loading test. At this stage ad = 0 should be used. Isotropic consolidation prior to cyclic loading guarantees that damage parameter d does not evolve. However, for K0 -consolidated clay, calibrated h0 should be adjusted after ad is known to render h that is equal to the first calibrated h0 prior to cyclic loading. • Calibration of ad is based on similar test data used for h0 calibration, but with much larger number of cycles. ad = 0 could be used for monotonic loading, if cyclic loading is not of interest. Below, model calibration and performance are shown for reconstituted Georgia kaolin clay (experimental data after Sheu (1984)), structured Cloverdale clay (experimental data 87  5.2. Reconstituted Georgia kaolin clay  Table 5.1: Model constants Model constant Designation Elasticity κ ν Critical state Mc Me N λ Rotational hardening C x Destructuration ki Bounding surface h0 ad  for the proposed SANICLAY-D-B model Description of its role Compressibility of overconsolidated clay Poisson’s ratio Critical stress ratio in triaxial compression Critical stress ratio in triaxial extension Shape of the bounding surface Compressibility of normally consolidated clay Rate of evolution of anisotropy Saturation limit of anisotropy, constant η = q/p path Rate of destructuration Initial hardening parameter Rate of damage evolution  Table 5.2: Initial internal variables for the proposed SANICLAY-D-B model Designation Description of its role p0 Size of the bounding surface Si Structuration factor ein Initial void ratio α Orientation of plastic potential and bounding surface (p, q) Initial stress state in the triaxial plane  after Zergoun (1991); Zergoun and Vaid (1994)) in slow cyclic loading, and reconstituted Ariake clay (experimental data after Yasuhara et al. (1992)) in fast cyclic loading.  5.2  Reconstituted Georgia kaolin clay  Experimental data after Sheu (1984) on Georgia kaolin clay is used for model calibration. The soil tested was composed of 62% of clay and 38% of silt. The plasticity index for this clay was 20%. Kaolin clay was mixed with salty water for preparation of Georgia kaolin clay in order to reduce the effect of anisotropy through formation of a flocculated structure. The mixture was poured into the consolidometer with two-way drainage and incremental vertical pressures were applied with the final pressure of 205 kPa. Then block sample was trimmed into triaxial specimens.  88  5.2. Reconstituted Georgia kaolin clay  5.2.1  Model calibration  Model calibration is based on available experimental results. Parameters for destructuration are set to the values corresponding to already destructured clay as the clay considered was reconstituted. Details of the calibration of model parameters for this clay are presented below. Parameters λ and κ. In Fig. 5.1 the response of Georgia kaolin clay in isotropic consolidation with three unloading-reloading sequences is presented in the e– log p plane. The slopes of the normal compression and the rebound lines are Cc and Cr , respectively, in the e– log p plane, whereas the required λ and κ are measured in the e– ln p plane. The relationships between the slopes are as follows: λ = Cc / ln 10 and κ = Cr / ln 10. Then the values Cc = 0.28 and Cr = 0.085 measured in Fig. 5.1 render λ = 0.121 and κ = 0.037. 1  experiment  0.9 Cr=0.085  e  0.8  Cc=0.28  0.7  0.6  0.5 1 10  2  10  3  10  p (kPa)  Figure 5.1: Calibration of Cc and Cr in e– log p plane based on isotropic consolidation test after Sheu (1984) Parameters Mc and Me . CIUC and CIUE stress paths on initially normally consolidated clay to a range of preconsolidation stresses σc = 276, 345,414, and 552 kPa are shown in Fig. 5.2. From direct measurements at the ultimate stress states, critical state ratios 89  5.2. Reconstituted Georgia kaolin clay are estimated as Mc = 0.87 and Me = 0.86. For later calibration of x, friction angle in compression ϕc might be required to come up with the K0 value; in this case, Mc is related to ϕc by Mc =  6 sin ϕc 3 − sin ϕc  (5.1)  400 Compression σc (kPa):  Mc=0.87  414 345 276  q (kPa)  200  0 Extension σc (kPa): 552 414 345 276  −200  Me=0.86 −400  0  200  400 p (kPa)  600  800  Figure 5.2: Calibration of Mc and Me in the triaxial plane based on undrained triaxial compression and extension tests starting from normally consolidated state, experimental results after Sheu (1984) Parameter N . Due to unavailability of CK0 UC test, the value for N must be assumed. Looking ahead, N = Mc or N = Me would require parameter C to be very small relatively to its typical values for clay in the range of C = 3–20 mentioned previously. In addition, cyclic loading simulations would suggest that the value should be considerably smaller for better simulations, especially for high cyclic stress amplitudes. For these reasons, the decision is made to adopt N = 0.8. 90  5.2. Reconstituted Georgia kaolin clay Parameter ν. Due to unavailability of K0 -unloading preceded by K0 -consolidation test for Georgia kaolin clay, parameter ν = 0.2 is assumed. This value has been calibrated for other clays by Taiebat, Dafalias and Peek (2010), and Dafalias et al. (2006). Parameter x. The parameter x for successful K0 -loading simulation can be obtained using the following equation proposed by Dafalias et al. (2006): x=  3 BϵηK 0  +  3 ηK 0  2ηK0 ϵ(1 − κ/λ) 2 (1 + ν) κ ; B=− 2 2 9 (1 − 2ν) λ + [2 (1 − κ/λ) − BMc ] ϵηK0 − Mc  (5.2)  where ϵ = 3/2 for K0 -loading path; K0 is a measured value of the earth pressure coefficient at rest; and ηK0 = 3 (1 − K0 ) / (1 + 2 K0 ). Should K0 is not available, Jacky’s equation, K0 = 1 − sin ϕ, can be used. For Georgia kaolin clay, ϕ = ϕc = 22.3◦ from Eq. 5.1 yields K0 = 0.62, ηK0 = 0.51 and x = 1.69. Parameter C. Calibration of C requires all other constants calibrated in advance (Dafalias et al., 2006). In the proposed model novel bounding surface parameters h0 and ad are not required to be known for calibration of C. This is true since parameters h0 and ad are in effect for overconsolidated states, while C calibration is based on the stress path corresponding to normally consolidated state. From the available data, CIUC test is the most appropriate for C calibration, ensuring rotation of the bounding surface. The C = 3 is chosen as an optimum value from trial-and-error simulations shown in Fig. 5.3. Bounding surface parameters h0 and ad . Calibration of h0 and ad requires experimental data involving the stress path corresponding to overconsolidated states. As cyclic loading is of concern, the parameters should be calibrated based on cyclic loading test. To be consistent, calibration of the two parameters should be performed from the same test results. Parameter h0 is responsible for the initial part of the cyclic loading, when accumulation of damage parameter d is still not significant to induce considerable reduction of h. Therefore, initial cycles from available undrained cyclic triaxial test are used for h0 calibration. h0 = 50 is calibrated from trial-and-error simulations with ad = 0 shown in Fig. 5.4. Once h0 is known, ad calibration is performed based on the same undrained cyclic triaxial test used for h0 calibration, but considering all available 15 cycles. For Georgia kaolin clay, ad = 7 is calibrated in a trial-and-error procedure.  91  5.2. Reconstituted Georgia kaolin clay  300 experiment C=9 200 q (kPa)  C=6  C=3  100  0 0  100  200 p (kPa)  300  Figure 5.3: Calibration of C based on undrained triaxial compression test on normally consolidated clay, experimental results after Sheu (1984) Calibrated constants and initial internal variables for Georgia kaolin clay are shown in Tables 5.3 and 5.4, respectively. In Table 5.4, ein = 1.5 and p = 1 kPa are used for simulations of undrained cyclic triaxial loading. For simulations of consolidation, CIUC and CIUE tests, different values of ein and p are used, which are to be specified for the corresponding tests simulation. Table 5.3: Model parameters for Georgia kaolin Clay Model constant Designation Elasticity κ ν Critical state Mc Me N λ Rotational hardening C x Destructuration (no effect) ki Bounding surface h0 ad  Value 0.037 0.2 0.87 0.86 0.8 0.121 3 1.69 0 50 7  92  5.2. Reconstituted Georgia kaolin clay  200  200  qcy=140.7 kPa  experiment simulation  experiment simulation  100 q (kPa)  q (kPa)  100  0  0  −100  −100  −200 0  −200 −3  100  200  300 p (kPa)  400  (a)  500  h0=50  −2  −1  0 εa (%)  1  2  3  (b)  Figure 5.4: Initial hardening parameter calibrated to h0 = 50 based on undrained cyclic triaxial loading, experimental results after Sheu (1984) Table 5.4: Initial internal variables for Georgia kaolin clay Model internal variable Designation Size of the bounding surface p0 Structuration factor Si Initial void ratio ein Orientation of plastic potential and bounding surface α Initial stress state (p, q)  5.2.2  Value 1 kPa 1 1.5 0 (1, 0) kPa  Model performance  Model constants calibrated previously (Table 5.3) are used in this section to simulate isotropic consolidation, CIUC, CIUE, and undrained cyclic triaxial tests on normally consolidated Georgia kaolin clay. Isotropic consolidation test The results of isotropic consolidation test simulation in comparison with experimental data are shown in Fig. 5.5. p0 = 75 kPa and ein = 0.88 are adopted for the simulation. Notice that plastic volumetric strains are involved upon reloading, even before the past maximum preconsolidation stress is reached. The magnitude of the plastic volumetric strains in reloading is in agreement with experimental results, suggesting that h0 = 50 calibrated based on  93  5.2. Reconstituted Georgia kaolin clay cyclic loading can be successfully used for monotonic loading simulations.  1  experiment simulation  0.9  e  0.8  0.7  0.6  0.5 1 10  2  10  3  10  p (kPa)  Figure 5.5: Simulation of isotropic consolidation test, experimental results after Sheu (1984)  Undrained triaxial compression and extension tests Simulated stress and stress–strain paths of CIUC and CIUE loading on clay normally consolidated to a range of preconsolidation stresses σc = 276, 345, 414 and 552 kPa are shown in Fig. 5.6 against available experimental data. For the stress path in compression, simulations closely match experimental results, while in extension the “hook-type” stress path cannot be reproduced by the model because of the evolution of anisotropic variable α being a function of tr (∂F/∂σ) according to Eq. 3.36b, i.e., the stress path freezes as the critical stress ratio is reached (Dafalias and Taiebat, 2012; Taiebat and Dafalias, 2013). In terms of the stress–strain response, simulated path ultimately matches experimental data, although experimental results suggest larger plastic axial strains at earlier stages. Note that there is no effect from parameters h0 and ad for CIUC and CIUE on normally consolidated clay.  94  5.2. Reconstituted Georgia kaolin clay  400 Compression σ (kPa): c  414 345 276  q (kPa)  200  0  Extension σ (kPa): c  −200  −400  552 414 345 276 0  200  400 p (kPa)  600  800  (a) 400 Compression σc (kPa):  q (kPa)  200  414 345 276  0  Extension σc (kPa):  −200  −400 −20  552 414 345 276 −10  0 εa (%)  10  20  (b)  Figure 5.6: Simulation of CIUC and CIUE tests on normally consolidated clay, experimental results after Sheu (1984) 95  5.2. Reconstituted Georgia kaolin clay Undrained cyclic triaxial loading The main focus of the current study is on simulations of cyclic loading. For convenience, stress and stress–strain paths are shown separately in Figs. 5.7 and 5.8, respectively, with the digitized experimental data on the left and simulations on the right. The results of 22 cycles for qcy = 121.4 kPa, 15 cycles for qcy = 136 and 140.7 kPa, and 3 cycles for qcy = 165.5 kPa are available from the experiments. Simulations are not exhibiting the same amount of p reduction as suggested in experiments owing to the balance reached by the simulated stress path at OCR = 2 (see section 3.3). In the experiments the stress path falls well beyond the critical stress ratio. A similar response was observed by Takahashi et al. (1980). They found that if the drainage was allowed after the undrained loading, then the stress ratio finally reduced to the critical stress ratio. This suggests that unreliable pore water pressure measurements were the reason of the stress ratio exceeding the critical stress ratio. For the stress–strain path, the model is able to reproduce cyclic softening observed in the experiments for the whole range of cyclic stress amplitudes. It is no surprise that simulations for qcy = 140.7 kPa closely match experiments as h0 and ad calibration is based on this particular cyclic stress amplitude. For a slightly lower stress amplitude, qcy = 136 kPa, the simulated stress–strain path is also matching well the experiment. For stress amplitudes much different from 140.7 kPa, simulations slightly deviate from experiments, although the results can be considered satisfactory for simulations with the single set of model constants. Observe that the simulated strain levels are less than experimentally observed ones for 165.5 kPa cyclic stress amplitude and larger for qcy = 121.4 kPa.  96  5.2. Reconstituted Georgia kaolin clay 200  200  qcy=121.4 kPa  simulation (22 cycles)  100 q (kPa)  q (kPa)  100  0  −100  −200 0  qcy=121.4 kPa  experiment (1−6,22 cycles)  0  −100  100  200  300 p (kPa)  400  −200 0  500  100  200  (a) 200  0  −100  simulation (15 cycles)  0  −100  100  200  300 p (kPa)  400  −200 0  500  100  200  (c) 200  200  qcy=140.7 kPa  −100  0  −100  100  200  300 p (kPa)  400  −200 0  500  100  200  200  qcy=165.5 kPa  300 p (kPa)  400  500  (f) experiment (3 cycles)  100  qcy=165.5 kPa  experiment (3 cycles)  100 q (kPa)  q (kPa)  500  simulation (15 cycles)  (e)  0  −100  −200 0  400  100 q (kPa)  q (kPa)  qcy=140.7 kPa  experiment (1−6,15 cycles)  0  200  300 p (kPa)  (d)  100  −200 0  500  100 q (kPa)  q (kPa)  qcy=136 kPa  experiment (1−6,15 cycles)  100  −200 0  400  (b) 200  qcy=136 kPa  300 p (kPa)  0  −100  100  200  300 p (kPa)  (g)  400  500  −200 0  100  200  300 p (kPa)  400  500  (h)  Figure 5.7: Undrained cyclic triaxial loading stress paths with the digitized experimental results after Sheu (1984) (left) and model simulations (right) for cyclic stress amplitudes: (a–b) qcy = 121.4 kPa, (c–d) qcy = 136 kPa, (e–f) qcy = 140.7 kPa, and (g–h) qcy = 165.5 kPa 97  5.2. Reconstituted Georgia kaolin clay 200  200  experiment (1−6,19−22 cycles) qcy=121.4 kPa  100 q (kPa)  q (kPa)  100  0  −100  −200 −15  0  −100  −10  −5  ε (%)  0  5  −200 −15  10  −10  a  q (kPa)  q (kPa)  10  100  0  −100  0  −100  −10  −5  ε (%)  0  5  −200 −15  10  −10  (c)a  −5  ε (%)  0  5  10  (d)a  200  200 experiment (1−7,13−15 cycles)  qcy=140.7 kPa  simulation (15 cycles)  100 q (kPa)  100 q (kPa)  5  qcy=136 kPa  simulation (15 cycles)  100  0  −100  0  −100  −10  −5  ε (%)  0  5  −200 −15  10  −10  (e)a  −5  ε (%)  0  5  10  (f)a 200  200  qcy=165.5 kPa  simulation (3 cycles)  experiment (3 cycles)  100 q (kPa)  100 q (kPa)  0  200 experiment (15 cycles)  0  0  −100  −100  −200 −15  ε (%)  (b)  200  −200 −15  −5  a  (a)  −200 −15  qcy=121.4 kPa  simulation (22 cycles)  −10  −5  εa (%)  (g)  0  5  10  −200 −15  −10  −5  εa (%)  0  5  10  (h)  Figure 5.8: Undrained cyclic triaxial loading stress–strain paths with the digitized experimental results after Sheu (1984) (left) and model simulations (right) for cyclic stress amplitudes: (a–b) qcy = 121.4 kPa, (c–d) qcy = 136 kPa, (e–f) qcy = 140.7 kPa, and (g–h) qcy = 165.5 kPa 98  5.3. Structured Cloverdale clay  5.3  Structured Cloverdale clay  Experimental data after Zergoun (1991), Zergoun and Vaid (1994) on Cloverdale clay are used for model calibration. Soil tested was sensitive undisturbed marine clay consisting of 49% of clay and 45% of silt with PI = 24%. The sensitivity was a result of surface infiltration following uplift above the sea level. The clay was block sampled from depth of (2.0–2.6) m below ground level. Then the triaxial specimens were trimmed from the same horizon and hydrostatically loaded under an effective stress of 200 kPa for 48 h. The consolidation stress beyond apparent preconsolidation stress destructured clay to some extend inducing reduced sensitivity. Prior to experiments, the specimens were left undrained under constant total stress for 24 h to allow pore water pressure generation due to arrest of secondary consolidation.  5.3.1  Model calibration  Calibration of the model parameters given below is based on available experimental data. The clay considered requires calibration of the parameters ki and Si . Parameters λ and κ. The response of Cloverdale clay in an oedometer test during onedimensional compression and swelling is shown in Fig. 5.9. Lacking lateral stress measurements in this test, the data are plotted in the e– log σa plane, where σa is an axial (vertical) stress. Cc and Cr measured in Fig. 5.9 are related to λ and κ as before, giving: λ = Cc / ln 10 = 0.21 and κ = Cr / ln 10 = 0.03. Note that while this indirect estimation of λ is accurate, determination of κ is only approximated, given the relative increase of the p/σa ratio in swelling due to the increase of K0 for the overconsolidated state (Wood, 1990), i.e., rate of p and σa change is no longer constant as in the loading process. Therefore, κ calibration based on isotropic consolidation test, if available, is preferred (Papadimitriou et al., 2005). Parameters Mc and Me . The response of Cloverdale clay in CIUC tests on normally consolidated and overconsolidated specimens is shown in Fig. 5.10. The data is presented in the MIT stress plane, normalized by the value of the equivalent consolidation stress σe , which is an equivalent stress on the normal consolidation line in the e– log σc plane corresponding to 99  5.3. Structured Cloverdale clay  1.6  experiment  1.4 Cc=0.48 1.2 e  Cr=0.07 1  0.8  0.6 0 10  1  10  2  10 σa (kPa)  3  10  Figure 5.9: Calibration of Cc and Cr in e– log p plane based on oedometer test after Zergoun (1991) the same void ratio as the overconsolidated clay specimen (Zergoun, 1991). Given these tests, the friction angle in compression is estimated as ϕc = 32◦ , rendering Mc = 1.29 according to Eq. 5.1. Me = 1.27 is chosen for better undrained cyclic triaxial loading simulations. Parameter N . Due to unavailability of a CK0 UC test the value of N = 1 is assumed. Similar to the N preference for Georgia kaolin clay in section 5.2.1, this value is chosen owing to its appropriateness for calibration of the constant C and better undrained cyclic triaxial loading simulations. Parameter ν. Due to unavailability of K0 -unloading stress path following K0 -consolidation for Cloverdale clay the parameter ν = 0.2 is assumed. Parameter x. Utilizing Eq. 5.2 along with Jacky’s equation, ϕ = 32◦ yields K0 = 0.47, ηK0 = 0.82 and x = 1.73. Parameters Si and ki . The best experiment for calibration of the initial internal variable Si and the value of destructuration parameter ki is an isotropic compression test on an intact sample until complete loss of structure occurs (Taiebat, Dafalias and Peek, 2010). The only  100  5.3. Structured Cloverdale clay  (σa − σr) /2 σe  0.6 Compression, OCR: 1.0 1.3 2.0 3.8  0.4  φc=32°  0.2  0  0  0.2  0.4  0.6 (σa + σr) /2 σe  0.8  1  1.2  Figure 5.10: Calibration of ϕc in the triaxial plane based on undrained triaxial compression test on normally consolidated and overconsolidated specimens, experimental results after Zergoun (1991) available consolidation test is oedometric consolidation, which is used for the calibration of the parameters. Parameter Si can be estimated from the compression plane as the ratio of the axial stress (or mean stress, in the case of isotropic loading) at the end of the initial stiff regime and the corresponding stress on the Intrinsic Compression Curve (ICC) at the same void ratio (Taiebat, Dafalias and Peek, 2010). The location of the ICC can be approximated by the dotted line in Fig. 5.11(a) and the ratio yields Si = 90/24 = 3.7 for Cloverdale clay. The parameter ki is calibrated to ki = 2 from trial-and-error simulations shown in Fig. 5.11(b). The initial value of p0 = 94 kPa has been used to successfully simulate the initial stiff (elastic) part of the response. New model constants that still require calibration are set to h0 = 1e10 and ad = 0 to simulate stiff response for overconsolidated state, since the calibration of ki involves normally consolidated state. Parameter C. CIUC is the most appropriate test from the available data for calibration of C, involving rotation of the bounding surface. From trial-and-error calibration procedure C = 3 is selected as an optimum value. Simulation results with different C values for the CIUC test on normally consolidated clay to p = 200 kPa are shown in Fig. 5.12. Bounding surface parameters h0 and ad . Model parameters h0 = 550 and ad = 68 are  101  5.3. Structured Cloverdale clay  1.6  experiment 24  90  1.4 Si=90/24=3.7 1.2 e  ICC (approx.) 1  0.8  0.6 0 10  1  10  2  3  10 σa (kPa)  10  (a)  1.6  experiment  1.4  k =1.6 i  ki=2.0 1.2  k =2.4  e  i  1  0.8  0.6 0 10  1  10  2  10 σa (kPa)  3  10  (b)  Figure 5.11: Calibration of: (a) Si , and (b) ki based on oedometer test after Zergoun (1991)  102  5.3. Structured Cloverdale clay 0.6  (σa − σr) /2 σe  experiment  0.4  C=9 C=6  0.2 C=3  0  0  0.2  0.4  0.6 (σa + σr) /2 σe  0.8  1  Figure 5.12: Calibration of C based on CIUC test on initially normally consolidated clay, experimental data after Zergoun (1991) calibrated based on the available experimental cyclic loading with τcy /Suc = 0.75 cyclic stress amplitude. The 11 calibrated model constants are shown in Table 5.5, and initial internal variables in Table 5.6. These parameters are to be used in simulations of the available experiments.  5.3.2  Model performance  Model performance for Cloverdale clay is illustrated in this section. Simulations are conducted for one-dimensional consolidation, CIUC on normally consolidated and overconsolidated clay, and undrained cyclic triaxial tests on normally consolidated specimens. Oedometer test In Fig. 5.13 experimental results and model simulations are shown for an oedometer test on Cloverdale clay. Destructuration observed in laboratory testing can be well simulated by the model. In the reloading following initial unloading, plastic volumetric strains for the stress states before the preconsolidation stress is attained are predicted by the model in agreement  103  5.3. Structured Cloverdale clay  Table 5.5: Model parameters for Cloverdale clay Model constant Designation Elasticity κ ν Critical state Mc Me N λ Rotational hardening C x Destructuration ki Bounding surface h0 ad  Table 5.6: Initial internal variables for Cloverdale clay Model internal variable Designation Size of the bounding surface p0 Structuration factor Si Initial void ratio ein Orientation of plastic potential and bounding surface α Initial stress state (p, q)  Value 0.03 0.2 1.29 1.27 1 0.21 3 1.73 2 550 68  Value 94 kPa 3.7 1.48 0 (1, 0) kPa  with the experiment. For unloading, simulated response closely matches experiment for σa ≥ 20 kPa; however, the dialative nature of clay in K0 -unloading cannot be predicted for σa < 20 kPa, resulting in slightly underpredicted void ratio at the end of unloading. One might anticipate a bounding surface providing the desired dialation as the image point upon K0 -unloading should be located so that the negative plastic volumetric strains are predicted. It is indeed so and the proposed model better simulates dialative response than the previous SANICLAY family models, but the degree of plastic strains with calibrated h0 and ad does not allow reaching dialation observed in the experiment. During parameters calibration, it was noticed that the parameter C is the factor that can control the degree of dialation through rotational hardening. For lower values of C, less rotation of the bounding surface renders more dialative response simulated for K0 -unloading. Thus, adjusted rotational hardening in conjunction with the bounding surface formulation can potentially enhance model predictions in K0 -loading.  104  5.3. Structured Cloverdale clay  1.6  experiment simulation  1.4  e  1.2  1  0.8  0.6 0 10  1  10  2  10 σa (kPa)  3  10  Figure 5.13: Simulation of oedometer test on Cloverdale clay, experimental results after Zergoun (1991) Undrained triaxial compression and extension tests Isotropic compression to p = 200 kPa is imposed prior to simulations of undrained triaxial compression and extension tests. Simulations of CIUC and CIUE tests on normally consolidated and overconsolidated specimens are performed with the results shown in Fig. 5.14. Simulated response approximately follows the experimental results in terms of the stress path for OCR ≤ 2. As OCR increases the stress path deviates from experiment in that the dialative response observed in actual tests is not captured. Similarly, the stress–strain path is closely matching for OCR ≤ 2 throughout the whole process of undrained loading. For OCR = 3.8, initial stress–strain simulation exhibits larger plastic strains than desired and then falls close to the experiment. Undrained cyclic triaxial loading Experimental data, where slow undrained cyclic triaxial loading was performed by Zergoun (1991) to maintain the uniformity of pore water pressure, is used in this section. Prior to 105  5.3. Structured Cloverdale clay  0.6 Compression, OCR: 1.0 0.4  2.0 3.8  a  r  (σ − σ ) /2 σ  e  1.3  0.2  0  0  0.2  0.4  0.6 (σ + σ ) /2 σ a  r  0.8  1  1.2  e  (a) 0.6  2.0 3.8  0.4  a  r  (σ − σ ) /2 σ  e  Compression, OCR: 1.0 1.3  0.2  0  0  2  4  6  8  10  εa (%)  (b)  Figure 5.14: Simulations of CIUC and CIUE tests on normally and overconsolidated clay, experimental results after Zergoun (1991)  106  5.3. Structured Cloverdale clay simulations of undrained cyclic triaxial loading, isotropic consolidation to p = 200 kPa is simulated. The stress and stress–strain path in undrained cyclic triaxial loading are shown separately in Figs. 5.15 and 5.16, respectively, with the experimental results on the left and model simulations on the right. Unlike undrained cyclic triaxial loading tests on Georgia kaolin clay, the slow cyclic loading imposed on Cloverdale clay renders the stress path that does not travel beyond the critical stress ratio, as shown in Fig.5.15. The simulated stress path that freezes at OCR = 2 (see section 3.3) still deviates from the observed experiments, but to a smaller degree than the simulated stress path for Georgia kaolin clay. In terms of the stress–strain simulations, τcy /Suc = 0.75 cyclic stress amplitude that is used for h0 and ad calibration is apparently well simulated. For τcy /Suc = 0.70, the stress–strain response closely matches the experiment in terms of attained peak strains. Simulation of τcy /Suc = 0.79 cyclic loading provides strains that are much smaller than the ones observed during the experiments. Finally, simulations in terms of peak εa versus number of cycles for the range of available cyclic stress amplitudes in undrained cyclic triaxial loading are shown in Fig. 5.17(b). Experimental data is digitized and is shown in Fig. 5.17(a) for comparison . Considering that earthquake loading is of interest, a M7.5 earthquake corresponds to approximately 30 equivalent number of cycles imposed on fine-grained material according to Idriss and Boulanger (2008). Therefore, the results of simulations and experimental data are plotted within 30 cycles. Evolution of peak εa with number of cycles is closely simulated for τcy /Suc = 0.75 and 0.69. Simulated τcy /Suc = 0.79 cyclic loading requires about 5 extra cycles to reach the experimentally observed peak εa = 10%. For lower cyclic stress amplitudes, the picture becomes less satisfactory, simulating a softer clay response than desired. Cyclic stress amplitude of τcy /Suc = 0.62 reaches approximately εa = 10% in 30 cycles as simulated, while in experiments only εa = 2% is attained in 30 cycles. With a τcy /Suc = 0.57, peak εa = 2% is simulated in 30 cycles, while less than 1% is reached in the experiments. A decreased rate of peak εa with a number of cycles is observed in the tests for τcy /Suc = 0.54 approaching a plateau, whereas strains keep accumulating in simulation. Such plateau, however, is attained for lower stress amplitudes, although at a strain levels above that observed in the experiments. It should be noted that the data is plotted in logarithmic scale for the vertical axis 107  5.3. Structured Cloverdale clay  0.4  τ /Su =0.70  0.4  experiment (5 cycles)  c  0.2 0 −0.2 −0.4 0  0.2  0.4 0.6 0.8 (σa + σr) /2 σe  1  τ /Su =0.70 cy  (σa − σr) /2 σe  (σa − σr) /2 σe  cy  0.2 0 −0.2 −0.4 0  1.2  0.2  0.4 0.6 0.8 (σa + σr) /2 σe  (a) 0.4  τ /Su =0.75  0.4  experiment (13 cycles)  c  0.2 0 −0.2 −0.4 0  0.2  0.4 0.6 0.8 (σa + σr) /2 σe  1  τ /Su =0.75 cy  τ /Su =0.79  0.2 0 −0.2 −0.4 0  1.2  0.2  0.4 0.6 0.8 (σa + σr) /2 σe  0.2 0 −0.2 −0.4 0  0.4  experiment (4 cycles)  c  0.2  0.4 0.6 0.8 (σa + σr) /2 σe  (e)  1  1.2  (d)  1  1.2  τ /Su =0.79 cy  (σa − σr) /2 σe  (σa − σr) /2 σe  cy  1.2  simulation (13 cycles)  c  (c) 0.4  1  (b)  (σa − σr) /2 σe  (σa − σr) /2 σe  cy  simulation (5 cycles)  c  simulation (4 cycles)  c  0.2 0 −0.2 −0.4 0  0.2  0.4 0.6 0.8 (σa + σr) /2 σe  1  1.2  (f)  Figure 5.15: Undrained cyclic triaxial loading stress path with the digitized experimental results after Zergoun (1991); Zergoun and Vaid (1994) (left), and model simulations (right) for cyclic stress amplitudes: (a-b) τcy /Suc = 0.70, (c-d) τcy /Suc = 0.75, and (e-f) τcy /Suc = 0.79  108  5.3. Structured Cloverdale clay  0.4  experiment (5 cycles)  (σa − σr) /2 σe  (σa + σr) /2 σe  0.4 0.2 0 −0.2 −0.4 −3  −2  −1  0 εa (%)  1  2  −0.2  (σa − σr) /2 σe  (σa + σr) /2 σe  0.4  0.2 0 −0.2 0 2 εa (%)  4  6  (σa − σr) /2 σe  (σa + σr) /2 σe  2  3  τ /Su =0.75  simulation (13 cycles)  cy  c  0 −0.2  0.4  0 −0.2  (e)  1  0 2 εa (%)  4  6  8 10  (d)  0.2  0 2 εa (%)  0 εa (%)  −0.4 −10 −8 −6 −4 −2  8 10  experiment (4 cycles)  −0.4 −10 −8 −6 −4 −2  −1  0.2  (c) 0.4  −2  (b)  experiment (13 cycles)  −0.4 −10 −8 −6 −4 −2  c  0  (a) 0.4  cy  0.2  −0.4 −3  3  τ /Su =0.70  simulation (5 cycles)  4  6  8 10  τ /Su =0.79  simulation (4 cycles)  cy  c  0.2 0 −0.2 −0.4 −10 −8 −6 −4 −2  0 2 εa (%)  4  6  8 10  (f)  Figure 5.16: Undrained cyclic triaxial loading stress path with the digitized experimental results after Zergoun (1991); Zergoun and Vaid (1994) (left), and model simulations (right) for cyclic stress amplitudes: (a-b) τcy /Suc = 0.70, (c-d) τcy /Suc = 0.75, and (e-f) τcy /Suc = 0.79  109  5.3. Structured Cloverdale clay and therefore deviation from experimental data is somewhat exaggerated for lower cyclic stress amplitudes. Observe that at lower cyclic stress amplitudes the simulated curves (Fig.5.17(b)) lie above the ones experimentally obtained (Fig.5.17(a)) from initiation of cyclic loading. This suggests that either too large plastic strains are simulated at low strain levels or elasticity model constnats/relation used is/are not appropriate at low level of strains. Simulations results, shown in Fig. 5.18, obtained with the original SANICLAY model reveal that approximately the same peak εa are attained during initial cycles for the lowest cyclic stress amplitude. This eliminates the possibility of too large plastic strains induced by the model, implying that either adopted elasticity constants κ and ν or elasticity relation in the proposed model might be not appropriate. Changing parameter κ would render inappropriate simulation of the consolidation process and would affect all other simulations, thus requiring re-calibration of some parameters. Parameter ν can potentially enhance the simulated elastic response, but this would require its reduction to values far below the typical range for clays. Simulations with reduced ν = 0.05 are shown in Fig. 5.17(c). Note that the overall simulated response does not change noticeably, except for the lowest τcy /Suc , where elastic strains are governing strain during initial cycles. Thus, only simulations of very low τcy /Suc can be controlled to some extent by elastic parameters. Undrained cyclic triaxial loading tests are also simulated with the original SANICLAY model for the whole range of available cyclic stress amplitudes with the results shown in Fig. 5.18. Model parameters and internal variables are the same as those used for the proposed SANICLAY-D-B model simulations. Parameters h0 and ad are not required for the original SANICLAY model. Clearly, no accumulation of peak εa can be simulated as the stress path falls inside the conventional yield surface, rendering purely elastic strains after the first half of cycle for all τcy /Suc . From the above simulated response of structured clay, the proposed SANICLAY-D-B model is considered to give adequate results within 30 cycles of applied cyclic loading, which is equivalent to a M7.5 earthquake (Idriss and Boulanger, 2008). Larger number of cycles might be considered with the reduced degree of satisfaction of the proposed model performance for low τcy /Suc . 110  5.3. Structured Cloverdale clay 1  10  Experiment  0.79 0.69  0.57  0.75  0  Peak εa (%)  0.62  10  −1  10  −2  10  0  compression extension 5  0.54 10  0.47  0.36  15 20 Number of cycles  τcy/Suc=0.27 25  30  (a)  1  10  0.79 0.75  Peak εa (%)  0  10  −1  0.69 0.62 0.57 0.54 0.47 0.36 τcy/Suc=0.27  10  −2  10  0  compession extension 5  SANICLAY−D−B 10  15 20 Number of cycles  25  30  (b)  1  10  0.79 0.75  Peak ε a (%)  0  10  −1  10  0.69 0.62  −2  10  0  0.57 0.54 0.47 0.36 τcy/Suc=0.27  compression extension 5  SANICLAY−D−B (ν=0.05) 10  15 20 Number of cycles  25  30  (c)  Figure 5.17: Simulations of peak εa versus number of cycles in undrained cyclic triaxial loading for the range of available cyclic stress amplitudes: (a) experimental results from Zergoun and Vaid (1994), and model simulations with (b) calibrated model parameters, (c) calibrated model parameters and ν = 0.05 111  5.4. Reconstituted Ariake clay (fast cyclic loading) 1  10  compression extension  Peak εa (%)  0  0.79  10  SANICLAY  0.75  −1  10  0.69 0.62 0.57  0.54 0.47 0.36 τ /Su =0.27 cy c  −2  10  0  5  10  15 20 Number of cycles  25  30  Figure 5.18: Simulations of peak εa versus number of cycles in undrained cyclic triaxial loading for the range of available cyclic stress amplitudes with the original SANICLAY model  5.4  Reconstituted Ariake clay (fast cyclic loading)  In order to illustrate if the model is able to reproduce the clay response in fast cyclic loading, although with less reliable experimental data (see section 2.1.1 for details), model parameters are calibrated approximately based on available data by Yasuhara et al. (1992). Two reconstituted marine clays, namely two Ariake clays, were used with the following index properties: PI = 69, Gs = 2.65; and PI = 72, Gs = 2.58. A clay slurry was consolidated under 55 kPa of vertical stress in a large consolidation vessel. Triaxial specimens were cut from the block samples formed in the consolidation vessel. The specimens were isotropically consolidated to σc = 200 kPa with the drainage allowed through filter paper surrounding the specimen in the triaxial cell. Then isotropically consolidated specimens were subjected to undrained cyclic triaxial loading with constant stress amplitude.  5.4.1  Model calibration  Because few experimental results are available, most of the model parameters are either approximated from known empirical relations or assumed as typical values for clays. A brief description of an approximate calibration procedure is given below.  112  5.4. Reconstituted Ariake clay (fast cyclic loading) Parameters λ and κ. Given average PI ∼ 71% and Gs ∼ 2.62, Eq. 4.5 yields λ = 0.41. The κ = 0.05 is chosen simply as 1/8 λ. Parameters Mc and Me . Yasuhara et al. (1992) obtained Mc = 1.68, corresponding to ϕ = ϕc = 41◦ from Eq. 5.1, based on monotonic undrained triaxial compression test with a strain rate of 0.05%/min. Such values for critical stress ratio M and friction angle ϕ are questionable and perhaps are erroneous owing to the high rate of applied monotonic loading resulting in non-equalized pore water pressure distribution within the tested specimen (see section 2.1.1). However, all undrained cyclic triaxial tests were performed with high frequency loading of 1 Hz. Therefore, it was decided to adopt the reported Mc = 1.68 and assumed value of Me = 1.65 for simulations of high frequency cyclic loading. Parameter N . For Ariake clay, N = 1.68 is assumed. Parameter x. Based on Eq. 5.2 and the adopted value of Mc , x = 1.76 is determined. Parameter ν. For Ariake clay, ν = 0.2 is chosen as this value is representative for most clays. Parameter C. C = 15 is calibrated from available CIUC test. Bounding surface parameters h0 and ad . For calibration of parameters h0 and ad , plots of peak εa versus number of cycles generated for qcy /2σc = 0.3 and qcy /2σc = 0.25 in undrained cyclic triaxial loading tests are used. This is an example of the parameters calibration based on the history of peak εa accumulation with number of cycles rather than already shown for Georgia kaolin and Cloverdale clay calibration of the parameters based on stress and stress– strain paths. In trial-and-error simulations the most appropriate parameters are h0 = 1600 and ad = 80 with the final simulation shown in Fig. 5.19. Calibrated model parameters are shown in Tables 5.7 and 5.8.  5.4.2  Model performance  With the calibration of parameters completed, it is of interest to see how well the model can replicate the experimental data in cyclic stress ratio CSR versus number of cycles required to reach double amplitude axial strain εDA = 5%. The comparison of the experimental data with the simulated results is shown in Fig. 5.20. CSR = qcy /2σc is adopted in accordance with the experiments. For CSR = 0.3, model simulates εDA = 5% reached in 10 cycles, while 113  a  Peak ε (%)  5.4. Reconstituted Ariake clay (fast cyclic loading)  5 4 3 2 1 0 −1 −2 −3 −4 −5 0  experiment simulation  0.3 qcy/2σc=0.225  5  10  15 20 25 Number of cycles  30  35  40  Figure 5.19: Calibration of h0 and ad based on peak εa versus number of cycles in undrained cyclic triaxial loading, experimental results after Yasuhara et al. (1992) in experiments the same strain level is reached in 6 cycles. As CSR decreases, simulated results fall closer to experiments. The two fitted curves for simulation and experimental data points intersect at CSR = 0.24. For further reduced CSR, approximately the same number of cycles is required to reach εDA = 5%. In general, reasonable results are obtained with the calibrated model parameters. It is pertinent to say that the use of CSR versus number of cycles for calibration purposes is not suggested as it considers only one particular strain level that a given CSR reaches in a number of cycles and neglects the history of strains from before and after that particular strain level. As has been pointed out previously, an equivalent number of uniform cyclic loading representing a M7.5 earthquake is about 30 (Idriss and Boulanger, 2008). Considering this and the model simulations discussed previously, calibrated parameters can be used with satisfactory prediction of clay cyclic behavior for an earthquake magnitude of about M7.5 and possibly larger with reduced degree of satisfaction. Note that the model calibration is performed for 1 Hz frequency of cyclic loading, and the same model parameters cannot be used for loadings with different frequencies as the rate-independent model is not able to differentiate between the applied frequencies.  114  5.4. Reconstituted Ariake clay (fast cyclic loading)  Table 5.7: Model parameters for Ariake clay Designation κ ν Critical state Mc Me N λ Rotational hardening C x Destructuration (no effect) ki Bounding surface h0 ad Model constant Elasticity  Table 5.8: Initial internal variables for Ariake clay Model internal variable Designation Size of the bounding surface p0 Structuration factor Si Initial void ratio ein Orientation of plastic potential and bounding surface α Initial stress state (p, q)  Value 0.05 0.2 1.68 1.65 1.68 0.41 15 1.76 0 1600 80  Value 37 kPa 1 2 0 (1, 0) kPa  0.4  experiment simulation  cy  Cyclic stress ratio, q /2 σ  c  0.5  0.3 0.2 Fitted curves 0.1 0 0  10  20 30 Number of cycles  40  50  Figure 5.20: [Simulation of cyclic stress ratio CSR versus number of cycles required to reach εDA = 5%, experimental results after Yasuhara et al. (1992) and fitted power function y = a CSR−b 115  6. Summary, Conclusions, and Recommendations for Future Research 6.1  Summary and conclusions  The bounding surface formulation is incorporated into the SANICLAY-D model (Taiebat, Dafalias and Peek, 2010) to simulate plastic strains within the conventional yield surface during undrained cyclic loading. In addition to the effect of the bounding surface, damage parameter d is proposed to further reduce hardening function h required for plastic modulus Kp reduction. An evolving projection center is adopted in a way that allows it to keep its relative position within the bounding surface established at the moment of stress reversal. The effect of newly introduced initial hardening parameter h0 and rate of damage evolution ad is illustrated in the sensitivity analysis. The advantage of the proposed projection center is demonstrated in the qualitative model performance. The model is compared with previous SANICLAY family models in undrained cyclic triaxial loading, demonstrating the ability of the model to simulate plastic strains for the stress states inside the bounding surface and the advantage of damage parameter d. The bounding surface effect in generalized loading clearly illustrates the ability of the model to simulate plastic strains for circular stress path loading in the deviatoric plane with the stress path entirely lying inside the bounding surface. For better understanding of a role of newly introduced model parameters h0 and ad from practical point of view, the effect of the parameters on the generation of the following plots is illustrated: • secant shear modulus Gsec and damping ratio variation with shear strain amplitude; 116  6.1. Summary and conclusions • cyclic stress ratio CSR versus number of cycles required to reach a certain strain level; • peak axial strain εa and residual pore water pressure u versus number of cycles Calibration and model performance in comparison with experimental data is conducted for reconstituted and structured clays. The calibration process of the model parameters is shown and the model performance in monotonic and cyclic loadings is demonstrated. Destructuration of structured clay in monotonic loading can be well simulated by the destructuration theory proposed by Taiebat, Dafalias and Peek (2010). It was reported by Idriss and Boulanger (2008) that 30 cycles of uniform cyclic loading are representative of a M7.5 earthquake for cyclic loading tests on clays. Thus, the number of cycles considered in simulations of cyclic loading is approximately 30. For reconstituted Georgia kaolin clay, simulations of undrained cyclic triaxial loading closely match experimentally observed strains for higher cyclic stress amplitudes, whereas for lower cyclic stress amplitudes the model slightly overpredicts the observed strains. Reasonable simulation results are obtained for structured Cloverdale clay in undrained cyclic triaxial loading with wide range of available cyclic stress amplitudes, although for lower cyclic stress amplitudes the strains are also overpredicted. Finally, model simulation of Ariake clay in fast undrained cyclic triaxial loading illustrates that the new model parameters calibrated based on sufficient experimental data can be successfully used to reproduce experimentally obtained cyclic stress ratio CSR versus number of cycles required to reach a certain strain level. The proposed model seems to have an adequate level of complexity for capturing many important characteristic features of clay response in both monotonic and cyclic loading. In the literature reviewed here (as well as those not discussed in the current work), none of the rate-independent bounding surface models demonstrated the ability to predict well a wide range of cyclic stress amplitudes. Some models were missing features of damage effect proposed in this work, others were not able to adequately simulate the damping, and yet others used model parameters that were not feasible to be calibrated for complex loading such as earthquake (e.g., cyclic stress amplitude parameter).  117  6.2. Recommendations for further research  6.2  Recommendations for further research  The proposed SANICLAY-D-B model provides adequately detailed framework for predictions of stress-strain response of clays under irregular cyclic loading that may be induced by earthquake. The model has been validated in a number of element tests. The validation process can be continued for various experiments on different clays. Creation and collection of a more comprehensive experimental database on laboratory element tests on clays would provide the basis for further research. The model should also be properly implemented into a finite element or finite difference framework. This would allow a detailed study on various problems of seismic wave propagation through clay deposits. In particular, the tool may be used for analysis of various foundations, dams, tunnels, and slopes subjected to earthquake loading. The model has been designed with particular attention on maintaining the simplicity of the formulation and the model calibration. However, if necessary, one may be able to address some of the shortcomings of the proposed model as discussed below. This would of course require a higher level of sophistication in the formulation and subsequently in the calibration process of the model. Although the model has significantly improved from its previous versions, it still may not sufficiently reproduce variations of shear modulus and damping ratio observed in experiments at very small strains. This shortcoming may be due to the very simple elastic relation adopted in the model and may be solved by adopting a more sophisticated nonlinear elasticity formulation that provides hysteretic stress–strain relationships in small strain levels (e.g., Jung and Yune, 2011). Also the present model may overpredict the plastic strains at small cyclic stress amplitudes with increasing number of cycles. The hardening function h may be the key issue for addressing some of the shortcoming of the model in wide range of cyclic stress amplitude and a more extensive study on this function may provide a further enhanced formulation, if needed. Finally the time dependency and creep have not been addressed in the presented formulation and may be considered in the future steps of the work.  118  Bibliography Al-Tabaa, A. and Wood, D. (1989), Experimentally based bubble model for clay., in ‘Proceedings of the Third International Symposium NUMOG III’, Elsevier Applied Science Publishers Limited, Niagara Falls, Canada, pp. 91–99. Anandarajah, A. and Dafalias, Y. (1986), ‘Bounding surface plasticity. III: Application to anisotropic cohesive soils.’, Journal of Engineering Mechanics 112(12), 1292–1318. Bacchus, D. (1969), Cyclic deformation of a clay, Ph.D. dissertation, The University of Auckland, Auckland, New Zealand. Banerjee, P. and Yousif, N. (1986), ‘A plasticity model for the mechanical behaviour of anisotropically consolidated clay’, International Journal for Numerical and Analytical Methods in Geomechanics 10(5), 521–541. Batdford, S. and Budiansky, B. (1949), A mathematical theory of plasticity based on the concept of slip, Technical Report TN 1871, NACA. Clayton, C. (2011), ‘Stiffness at small strain: Research and practice’, Geotechnique 61(1), 5– 37. Dafalias, Y. (1975), On Cyclic and Anisotropic Plasticity: I) A General Model including Material Behavior under Stress Reversals. II) Anisotropic Hardening for Initially Orthotropic Materials, Ph.D. dissertation, University of California, Berkeley, United States. Dafalias, Y. (1986a), ‘An anisotropic critical state soil plasticity model’, Mechanics Research Communications 13(6), 341–347. Dafalias, Y. (1986b), ‘Bounding surface plasticity. I: Mathematical foundation and hypoplasticity.’, Journal of Engineering Mechanics 112(9), 966–987. 119  Bibliography Dafalias, Y. (1987), An anisotropic critical state clay palsticity model, in C. Desai, E. Krempl, P. Kiousis and T. Kundu, eds, ‘Constitutive laws for engineering materials: theory and application’, pp. 513–521. Dafalias, Y. and Herrmann, L. (1986), ‘Bounding surface plasticity. II: Application to isotropic cohesive soils.’, Journal of Engineering Mechanics 112(12), 1263–1291. Dafalias, Y., Manzari, M. and Papadimitriou, A. (2006), ‘SANICLAY: simple anisotropic clay plasticity model’, International Journal for Numerical and Analytical Methods in Geomechanics 30(12), 1231–1257. Dafalias, Y. and Popov, E. (1975), ‘A model of nonlinearly hardening materials for complex loading.’, Acta Mechanics 21, 173–192. Dafalias, Y. and Taiebat, M. (2012), Advantages and disadvantages of some rotational hardening rules in clay plasticity, in ‘Third International Workshop on Modern Trends in Geomechanics IW-MTG3, 4–5 September 2012’. Desai, C., Somasundaram, S. and Frantziskonis, G. (1986), ‘A hierarchical approach for constitutive modelling of geologic materials’, International Journal for Numerical and Analytical Methods in Geomechanics 10(3), 225–257. Desai, C., Wathugala, G. and Matlock, H. (1993), ‘Constitutive model for cyclic behavior of clays. II: applications’, Journal of Geotechnical Engineering 119(4), 730–748. Hardin, B. and Drnevich, V. (1972), ‘Shear modulus and damping in soils: Measurement and parameter effects.’, American Society of Civil Engineers, Journal of the Soil Mechanics and Foundations Division 98(SM6), 603–624. Hashiguchi, K. and Ueno, M. (1977), ‘Elasto-plastic constitutive laws of granular materials’, Proceedings of the International Conference on Soil Mechanics and Foundation Engineering = Comptes Rendus du Congres International de Mecanique des Sols et des Travaux de Fondations (9), 73–82.  120  Bibliography Huang, M., Liu, Y. and Sheng, D. (2011), ‘Simulation of yielding and stress–stain behavior of shanghai soft clay’, Computers and Geotechnics 38(3), 341–353. Idriss, I. and Boulanger, R. (2008), Soil Liquefaction during Earthquakes, Earthquake Engineering Research Institute (EERI), Oakland, California, USA. Iwan, W. (1967), ‘On a class of models for the yielding behavior of continuous and composite systems’, Journal of Applied Mechanics 34(3), 612–617. Jiang, J. and Ling, H. (2010), ‘A framework of an anisotropic elastoplastic model for clays’, Mechanics Research Communications 37(4), 394–398. Jiang, J., Ling, H. and Kaliakin, V. (2012), ‘An associative and non-associative anisotropic bounding surface nodel for clay’, Journal of Applied Mechanics, Transactions ASME 79(3). Jung, Y. and Yune, C. (2011), ‘Reappraisal of the anisotropic bounding surface model of small-strain behavior for clays’, KSCE Journal of Civil Engineering 15(3), 463–472. Koiter, W. (1953), ‘Stress-strain relations, uniqueness and variational theorems for elasticplastic materials with a singular yield surface’, Quarterly of Applied Mathematics 11, 350. Li, T. and Meissner, H. (2002), ‘Two-surface plasticity model for cyclic undrained behavior of clays’, Journal of Geotechnical and Geoenvironmental Engineering 128(7), 613–626. Ling, H., Yue, D., Kaliakin, V. and Themelis, N. (2002), ‘Anisotropic elastoplastic bounding surface model for cohesive soils’, Journal of Engineering Mechanics 128(7), 748–758. Mr´oz, Z. (1967), ‘On the description of anisotropic workhardening’, Journal of the Mechanics and Physics of Solids 15(3), 163–175. Mr´oz, Z. (1969), ‘An attempt to describe the behavior of metals under cyclic loads using a more general workhardening model’, Acta Mechanica 7, 199–212. Papadimitriou, A., Manzari, M. and Dafalias, Y. (2005), Calibration of a simple anisotropic plasticity model for soft clays, number 128, pp. 415–424.  121  Bibliography Roscoe, K. and Burland, J. (1966), On the generalised stress-strain behavior of “wet” clay, in Heymann and Leckie, eds, ‘Engineering Plasticity’, Cambridge University Press, pp. 53– 609. Roscoe, K. and Poorooshasb, H. (1963), ‘A theoretical and experimental study of strains in triaxial compression tests on normally consolidated clays’, Geotechnique 13(1), 12–38. Roscoe, K., Schofield, A. and Thurairajah, A. (1963), ‘Yielding of clays in states wetter than critical’, Geotechnique 13(3), 211–240. Roscoe, K., Schofield, A. and Wroth, C. (1958), ‘On the yielding of soils’, Geotechnique 8(1), 22–53. Sangrey, D., Henkel, D. and Esrig, M. (1969), ‘The effective stress response of a saturated clay soil to repeated loading’, Canadian Geotechnical Journal 6(3), 241–252. Seed, H., Wong, R., Idriss, I. and Tokimatsu, K. (1986), ‘Moduli and damping factors for dynamic analyses of cohesionless soils’, Journal of Geotechnical Engineering 112(11), 1016– 1032. Sheu, W. (1984), Modeling of stress-strain-strength behavior of a clay under cyclic loading (Soil Dynamics), Ph.D. dissertation, Univeristy of Colorado, Boulder, Colorado, USA. Stallebrass, S. and Taylor, R. (1997), ‘The development and evaluation of a constitutive model for the prediction of ground movements in overconsolidated clay’, Geotechnique 47(2), 235–253. Suebsuk, J., Horpibulsuk, S. and Liu, M. (2011), ‘A critical state model for overconsolidated structured clays’, Computers and Geotechnics 38(5), 648–658. Taiebat, M. and Dafalias, Y. (2013), Rotational hardening and uniqueness of critical state line clay plasticity, in ‘Second International Symposium on Constitutive Modeling of Geomaterials: Advances and New Applications, 15–16 October 2012’. Taiebat, M., Dafalias, Y. and Kaynia, A. (2010), Bounding surface model for natural anisotropic clays, in P. Papanastasiou, E. Papamichos, A. Zervos and M. Stavropoulou, 122  Bibliography eds, ‘Ninth HSTAM International Congress on Mechanics, Proceedings of Vardoulakis mini-symposia’, Limassol, Cyprus, pp. 43–47. Taiebat, M., Dafalias, Y. and Peek, R. (2010), ‘A destructuration theory and its application to SANICLAY model’, International Journal for Numerical and Analytical Methods in Geomechanics 34(10), 1009–1040. Takahashi, M., Hight, D. and Vaughan, P. (1980), Effective stress changes observed during undrained cyclic triaxial tests on clay, in G. Pande and O. Zienkiewicz, eds, ‘Soils under Cyclic and Transient Loading’, Vol. 1, 2, pp. 201–209. Vucetic, M. and Dobry, R. (1988), ‘Degradation of marine clays under cyclic loading.’, Journal of Geotechnical Engineering 114(2), 133–149. Vucetic, M. and Dobry, R. (1991), ‘Effect of soil plasticity on cyclic response’, Journal of Geotechnical Engineering 117(1), 89–107. Wathugala, G. and Desai, C. (1993), ‘Constitutive model for cyclic behavior of clays. I: theory’, Journal of Geotechnical Engineering 119(4), 714–729. Wheeler, S., N¨a¨at¨anen, A., Karstunen, M. and Lojander, M. (2003), ‘An anisotropic elastoplastic model for soft clays’, Canadian Geotechnical Journal 40(2), 403–418. Whittle, A. (1993), ‘Evaluation of a constitutive model for overconsolidated clays’, Geotechnique 43(2), 289–313. Wood, D. (1990), Soil behaviour and critical state soil mechanics, Cambridge University Press. Yasuhara, K., Hirao, K. and FL Hyde, A. (1992), ‘Effects of cyclic loading on undrained strength and compressibility of clay’, Soils and Foundations 32(1), 100–116. Yu, H. (1998), ‘CASM: A unified state parameter model for clay and sand’, International Journal for Numerical and Analytical Methods in Geomechanics 22(8), 621–653. Yu, H., Khong, C. and Wang, J. (2007), ‘A unified plasticity model for cyclic behaviour of clay and sand’, Mechanics Research Communications 34(2), 97–114. 123  Bibliography Zergoun, M. (1991), Effective stress response of clay to undrained cyclic loading, Ph.D. dissertation, The University of British Columbia, Vancouver, British Columbia, Canada. Zergoun, M. and Vaid, Y. (1994), ‘Effective stress response of clay to undrained cyclic loading’, Canadian Geotechnical Journal 31(5), 714–727. Zienkiewicz, O., Leung, K. and Pastor, M. (1985), ‘Simple model for transient soil loading in earthquake analysis. I: Basic model and its application.’, International Journal for Numerical and Analytical Methods in Geomechanics 9(5), 453–476.  124  

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.24.1-0073338/manifest

Comment

Related Items