Numge 2010

cover_NUMGE2010_D.
qxd 29-04-2010 12:44 Pagina 1
Numerical Methods in Geotechnical Engineering

contains 153 scientific papers presented at the 7th European Conference on BENZ & NORDAL
Numerical Methods in Geotechnical Engineering, NUMGE 2010, held at editors
Norwegian University of Science and Technology (NTNU) in Trondheim,
Norway, 2–4 June 2010.
IN GEOTECHNICAL ENGINEERING
NUMERICAL METHODS
The contributions cover topics from emerging research to engineering
practice, and are organized into the following sections:
- Constitutive modeling;
- Computer codes and algorithms;
- Discontinuum and particulate modeling;
- Large deformation – large strain analysis;
- Flow and consolidation;
- Unsaturated soil mechanics;
- Artificial intelligence;
- Reliability and probability analysis;
- Dynamic problems and Geohazards;
- Slopes and cuts;
- Embankments, shallow foundations, and settlements;
- Piles;
- Deep excavations and retaining walls;
- Tunnels and caverns;
- Ground improvement modeling;
- Offshore geotechnical engineering;
NUMERICAL METHODS
- Numerical methods and Eurocode.
IN GEOTECHNICAL ENGINEERING
The book is of special interest to academics and engineers in geotechnical
engineering.
2010
an informa business
2010
EDITED BY:
THOMAS BENZ
STEINAR NORDAL
NUMERICAL METHODS IN GEOTECHNICAL ENGINEERING
PROCEEDINGS OF THE SEVENTH EUROPEAN CONFERENCE ON NUMERICAL METHODS IN
GEOTECHNICAL ENGINEERING, TRONDHEIM, NORWAY, 2–4 JUNE 2010
Numerical Methods in Geotechnical

Engineering
Edited by
Thomas Benz & Steinar Nordal
Department of Civil and Transport Engineering, Norwegian University of
Science and Technology (NTNU), Trondheim, Norway
CRC Press
Taylor & Francis Group
6000 Broken Sound Parkway NW, Suite 300
Boca Raton, FL 33487-2742
© 2000 by Taylor & Francis Group, LLC
CRC Press is an imprint of Taylor & Francis Group, an Informa business
No claim to original U.S. Government works

Version Date: 20140602
International Standard Book Number-13: 978-0-203-84236-2 (eBook - PDF)
This book contains information obtained from authentic and highly regarded sources. Reasonable efforts have been
made to publish reliable data and information, but the author and publisher cannot assume responsibility for the
validity of all materials or the consequences of their use. The authors and publishers have attempted to trace the
copyright holders of all material reproduced in this publication and apologize to copyright holders if permission to
publish in this form has not been obtained. If any copyright material has not been acknowledged please write and let
us know so we may rectify in any future reprint.
Except as permitted under U.S. Copyright Law, no part of this book may be reprinted, reproduced, transmitted, or
utilized in any form by any electronic, mechanical, or other means, now known or hereafter invented, including pho-
tocopying, microfilming, and recording, or in any information storage or retrieval system, without written permis-
sion from the publishers.
For permission to photocopy or use material electronically from this work, please access www.copyright.com (http://
www.copyright.com/) or contact the Copyright Clearance Center, Inc. (CCC), 222 Rosewood Drive, Danvers, MA
01923, 978-750-8400. CCC is a not-for-profit organization that provides licenses and registration for a variety of
users. For organizations that have been granted a photocopy license by the CCC, a separate system of payment has
been arranged.
Trademark Notice: Product or corporate names may be trademarks or registered trademarks, and are used only for
identification and explanation without intent to infringe.
Visit the Taylor & Francis Web site at
http://www.taylorandfrancis.com
and the CRC Press Web site at
http://www.crcpress.com
Numerical Methods in Geotechnical Engineering – Benz & Nordal (eds)
© 2010 Taylor & Francis Group, London, ISBN 978-0-415-59239-0
Table of Contents
Preface XIII
Scientific Committee (ERTC 7) XV
Constitutive modelling
A non-associated creep model for structured anisotropic clay (n-SAC) 3

G. Grimstad & S.A. Degago
A state dependent constitutive model for sand-structure interfaces 9
A. Lashkari
Adaptive integration of hypoplasticity 15
W. Fellin, M. Mittendorfer & A. Ostermann
An anisotropic bubble model for soft clays 21
N. Sivasithamparam, D. Kamrat-Pietraszewska & M. Karstunen
An anisotropic model for structured soils 27
G. Belokas & M. Kavvadas
An examination of strain space versus stress space for the formulation of
elastoplastic constitutive models 33
K.C. Ellison, K. Soga & B. Simpson
Anisotropic small strain stiffness within the multilaminate framework 39
B. Schädlich & H.F. Schweiger
Application of discontinuity layout optimization to problems involving non-associative friction 45
A.F. Babiker, C.C. Smith & M. Gilbert
Associated plasticity for nonassociated frictional materials 51
K. Krabbenhoft, A.V. Lyamin & S.W. Sloan
Comparison of methods for calculation of settlements of soft clay 57
H.P. Jostad & S.A. Degago
Effect of yield surface shape on the simulated elasto-plastic response of cohesive soils 63
A.G. Papadimitriou, A.D. Vranna, Y.F. Dafalias & M.T. Manzari
Impact of input data on soil model calibration using Genetic Algorithms 69
D. Taborda, A. Pedro, P.A.L.F. Coelho & D. Antunes
Influence of destructuration of soft clay on time-dependant settlements 75
D.F.T. Nash
Modeling liquefaction behavior of sands by means of hypoplastic model 81
A.B. Tsegaye, F. Molenkamp, R.B.J. Brinkgreve, P.G. Bonnier, R. de Jager & V. Galavi
Modeling of creep mechanism and damage of rock salt 89
B. Leuger, K. Staudtmeister, S. Yıldırım & D. Zapf
Modeling static liquefaction within multilaminate framework 95
A.B. Tsegaye, V. Galavi, R.B.J. Brinkgreve, R. de Jager, F. Molenkamp & P.G. Bonnier
On the differences between the Drucker-Prager criterion and exact implementation of
the Mohr-Coulomb criterion in FEM calculations 101
J. Clausen, L. Andersen & L. Damkilde
V
Simulation of mechanical behaviour of Toyoura sand using Severn Trent constitutive model 107
S. Miliziano, G.M. Rotisciani & F.M. Soccodato
Soil parameter identification for cyclic loading 113
A. Papon, Z.-Y. Yin, K. Moreau, Y. Riou & P.-Y. Hicher
Study of tensorial damage in a porous geomaterial 119
M. Mozayan Kharazi, C. Arson & B. Gatmiri
Time- and stress-compressibility of clays during primary consolidation 125
S.A. Degago, H.P. Jostad, M. Olsson, G. Grimstad & S. Nordal
Uncertainty and sensitivity analysis of laboratory test simulations using an elastoplastic model 131
F. Lopez-Caballero & A. Modaressi-Farahmand-Razavi
Validation of empirical formulas to derive model parameters for sands 137
R.B.J. Brinkgreve, E. Engin & H.K. Engin
Computer codes and algorithms

3D parallel computing FEA in offshore foundation design 145
L. Andresen, H. Sturm, M. Vöge & K. Skau
70-line 3D finite deformation elastoplastic finite-element code 151
W.M. Coombs, R.S. Crouch & C.E. Augarde
A simple time stepping algorithm for material point method 157
W.T. Sołowski & D. Sheng
Analysis of the stability of sheet pile walls using Discontinuity Layout Optimization 163
S.D. Clarke, C.C. Smith & M. Gilbert
Application of Discontinuity Layout Optimization to geotechnical limit analysis problems 169
M. Gilbert, C.C. Smith, I.W. Haslam & T.J. Pritchard
Enhancing solution procedures of a new numerical scheme for dynamic analysis
of soil-structure interaction problems 175
M.H. Bazyar & Ch. Song
Numerical bearing capacity computation and load-displacement behavior of shallow
foundations by stress level based ZEL method 181
M. Jahanandish, M. Veiskarami & A. Ghahramani
Simple quality indicators for FE analysis based on stress maps for Gauss points 187
C. Vulpe, N. Droniuc, E. Bourgeois & Ph. Mestat
The upper bound limit analysis of bearing capacity problems using the finite element method 193
A.I.M. AL-Janabi, A.A.R. Orabi & A.Y.A. Baqir
Discontinuum and particulate modelling

A numerical simulation on centrifuge liquefaction model using microscopic fluid
coupling scheme with Discrete Element Method 201
Y. Shimizu & Y. Inagawa
Discrete element modeling of low strength rock 207
N.B. Yenigül & M. Alvarez Grima
Effect of drying on a granular slope physical model analysed by Discrete Element Method (DEM) 213
F. Gabrieli, S. Cola, P. Simonini & F. Calvetti
Isotropic compression of cohesive-frictional particles with rolling resistance 219
S. Luding
Size effects on a virtual calibration chamber 225
J. Butlanska, M. Arroyo & A. Gens
VI
Large deformation – large strain analysis
A Coupled Eulerian-Lagrangian approach to solve geotechnical problems involving

large deformations 233
S. Henke, G. Qiu & J. Grabe
Advances in meshless methods with application to geotechnics 239
C.E. Heaney, C.E. Augarde, A.J. Deeks, W.M. Coombs & R.S. Crouch
An ALE Finite Element Method for cohesionless soil at large strains:
Computational aspects and applications 245
D. Aubram, F. Rackwitz & S.A. Savidis
Analysis of dynamic penetration of objects into soil layers 251
J.P. Carter & M. Nazem
Large deformation analysis of the installation of Dynamic Anchor 255
H. Sturm & L. Andresen
Modelling of installation effects of driven piles using hypoplasticity 261
H.D. Pham, H.K. Engin, R.B.J. Brinkgreve, & A.F. van Tol
Flow and consolidation

A multiscale approach for the consideration of spatial groundwater flow in the stability
analysis of a large excavation pit 269
H. Montenegro & R. Kauther
A numerical model for the electrokinetic treatment of natural soils with calcite 275
F. Cattaneo, C. Jommi & G. Musso
Analysis of artificial ground freezing in the Pari-Duomo platform tunnel of the Naples metro 281
S. Papakonstantinou, E. Pimentel & G. Anagnostou
Large scale hydraulic conductivity of the soil deposits of the Venezia Lagoon from
numerical back-analysis 285
E. Giacomini, F. Colleselli, F. Cattaneo, C. Jommi & G. Mayerle
Numerical analyses of granulometric stability of moraine dam cores 291
F. Federico & A. Montanaro
Numerical prediction of time-dependent rock swelling based on an example of
a major tunnel project in Ontario/Canada 297
A. Kirsch & T. Marcher
Some features of the coupled consolidation models used for the evaluation of the dissipation test 303
E. Imre & P. Rózsa
Steady state seepage flow through zoned earth structures affected by permeability defects 311
F. Federico, F. Calzoletti & A. Montanaro
Unsaturated soil mechanics
A comparison between numerical integration algorithms for unsaturated soils constitutive models 319
F. Cattaneo, G. Della Vecchia, C. Jommi & G. Maffioli
Comparison of stress update algorithms for partially saturated soil models 325
M. Hofmann, G. Hofstetter & A. Ostermann
Modelling of the hysteretic soil–water retention curve for unsaturated soils 331
A. Tsiampousi, L. Zdravkovic & D.M. Potts
Numerical integration and analysis of equilibrium in unsaturated multiphase media 337
R. Tamagnini, M. Mavroulidou & M.J. Gunn
VII
Artificial intelligence
A genetic algorithm for solving slope stability problems: From Bishop to a free slip plane 345
R. van der Meij & J.B. Sellmeijer
Simulation of the mechanical behavior of railway ballast by intelligent computing 351
M.A. Shahin
Three dimensional site characterization model of Suurpelto (Finland) using support vector machine 355
A. Pijush Samui & T. Länsivaara
Reliability and probability analysis

Evaluation of soil variability and its consequences 363
M. Huber, P.A. Vermeer & A. Bárdossy
Inverse modelling including spatial variability applied to the construction of a road embankment 369
A. Hommels, F. Molenkamp, M. Huber & P.A. Vermeer
Reliability analysis of piping in embankment dam 375
A. Noorzad & M. Rohaninejad
Spatial variability of soil parameters in an analysis of a strip footing using hypoplastic model 383
R. Suchomel & D. Mašín
Validating models against experience in foundation engineering, using the ROC curve 389
A.M.J. Mens & A.F. van Tol
Dynamic problems and Geohazards
A 2.5D finite element model for simulation of unbounded domains under dynamic loading 397
P. Alves Costa, R. Calçada, J. Couto Marques & A. Silva Cardoso
A comparison of different approaches for the modelling of shallow foundations in
seismic soil-structure interaction problems 405
S. Grange, D. Salciarini, P. Kotronis & C. Tamagnini
A finite element approach for dynamic seepage flows 411
R. Stucchi, A. Cividini & G. Gioda
A method to solve Biot’s u-U formulation for soil dynamics applications using
the ABAQUS/explicit platform 417
F.J. Ye, S.H. Goh & F.H. Lee
Alternative formulations for cyclic nonlinear elastic models: Parametric study and
comparative analyses 423
D. Taborda, L. Zdravkovic, S. Kontoe & D.M. Potts
Analysis of the effect of pile length in a pile group on the transfer and impedance functions
in soil-pile interaction models 429
A. Mahboubi & K. Panaghi
Dynamic fragmentation in rock avalanches: A numerical model of micromechanical behaviour 435
K.L. Rait & E.T. Bowman
Evaluation of the efficiency of a model of rockfall protection structures based
on real-scale experiments 441
F. Bourrier, Ph. Gotteland, A. Heymann & S. Lambert
Evaluation of viscous damping due to solid-fluid interaction in a poroelastic layer subjected
to shear dynamic actions 447
J. Grazina, P.L. Pinto & D. Taborda
Non linear numerical modeling of slopes stability under seismic loading – reinforcement effect 453
F. Hage Chehade, M. Sadek & I. Shahrour
VIII
Numerical analysis of blast impact on sealings of neighbouring structures 459
W. Krajewski, O. Reul & L. te Kamp
Numerical analysis of the seismic behavior of vertical shaft 465
S. Jeong, Y. Kim, S. Lee, J. Jang & Y. Lee
Numerical and experimental study of the detection of underground heterogeneities 471
P. Alfonsi, E. Bourgeois, F. Rocher-Lacoste, L. Lenti, & M. Froumentin
Numerical modelling of impacts on granular materials with a combined
discrete – continuum approach 477
A. Breugnot, Ph. Gotteland & P. Villard
Numerical simulations of the dynamic impact force of fluidized debris flows onto structures 483
F. Federico & A. Amoruso
Three dimensional analysis of seismic performance of an earthfill dam in Ethiopia 489
B.G. Tensay & W. Wu
Slopes and cuts
Effect of updated geometry in analyses of progressive failure 497

A.S. Gylland & H.P. Jostad
Evaluation of the effective width method for strip footings on slopes under undrained loading 503
K. Georgiadis & E. Skoufaki
Failure induced pore pressure by simple procedure in LEM 509
T. Länsivaara
Investigation of soil property sensitivity in progressive failure 515
A.S. Gylland, M.S. Sayd, H.P. Jostad & S. Bernander
Short-term slope stability calculation according to Eurocode 7 521
V. Thakur, S. Nordal & S. Hove
Embankments, shallow foundations, and settlements

3D settlement analysis using GIS and FEM: A case study in Sliedrecht area, the Netherlands 529
N.B. Yenigül & A.S. Elkadi
A comparison of 1D, 2D, and 3D settlement analyses of the Tower of Pisa 535
A.J. Klettke & L. Edgers
Analysis of a full scale failure test on old railway embankment 541
J. Mansikkamäki & T. Länsivaara
Analysis of ground movements induced by diaphragm wall installation 547
B. Garitte, M. Arroyo & A. Gens
Bearing capacity of a surface footing founded on a layered clay subsoil 553
Z. Bournta & L. Zdravkovic
Finite element analysis of the main embankment at Empingham dam 557
A. Grammatikopoulou, N. Kovacevic, L. Zdravkovic & D.M. Potts
Forecasting of the stability of the tailing dam in permafrost region on the basis
of numerical methods 563
A.B. Lolaev, A.P. Akopov, A.Kh. Oganesian, M.N. Sumin & V.V. Butygin
Numerical modeling of the mechanical response of recycled materials in embankments 569
M.M. Villani, X. Liu, A. Scarpas & A. D’Andrea
Rail track structural analysis using three-dimensional numerical models 575
A. Paixão & E. Fortunato
IX
Three dimensional analyses of ring foundations 581
M. Laman, A. Yildiz, M. Ornek & A. Demir
Piles
A back analysis of vertical load tests on bored piles in granular soil 589
L. Tosini, A. Cividini & G.Gioda
A numerical study on the effects of time on the axial load capacity of piles in soft clays 595
K.P. Giannopoulos, L. Zdravkovic & D.M. Potts
Analysis of foundation solution of new building in built-up area 601
Ž. Arbanas, V. Jagodnik & S. Dugonjić
Collapse of thin-walled model piles during hard driving 607
J. Bergan, S. Øren Holo & S. Nordal
Dynamic analysis of large diameter piles Statnamic load test 613
K.J. Bakker, F.J.M. Hoefsloot & E. de Jong
Finite difference analysis of pile on sloping ground under passive loading 619
K. Muthukkumaran & M. Gokul Khrishnan
Ground displacements due to pile driving in Gothenburg clay 625
T. Edstam & A. Kullingsjö
Lateral loading of pile foundations due to embankment construction 631
A. Feddema, J. Breedeveld & A.F. van Tol
Modelling of piled rafts with different pile models 637
S.W. Lee, W.W.L. Cheang, W.M. Swolfs & R.B.J. Brinkgreve
Modelling performance of jack-in piles 643
S. Jie & S.-A. Tan
Numerical analyses of axial load capacity of rock socketed piles in Turkey 649
M. Kirkit, H. Kılıç & C. Akgüner
Numerical simulation of low-strain integrity tests on model piles 655
J. Fischer, C. Missal, M. Breustedt & J. Stahlmann
Response of pile groups in clays under lateral loading based on 3-D numerical experiments 661
E.M. Comodromos, M.C. Papadopoulou & I.K. Rentzeperis
Selection of the proper hammer in pile driving and estimation of the total driving time 667
A. Afshani, A. Fakher & M. Palassi
Settlement analysis of a large piled raft foundation 673
M. Wehnert, T. Benz, P. Gollub & T. Cubaleski
Study of a complex deep foundation system using 3D Finite Element analysis 679
F. Tschuchnigg & H.F. Schweiger
The influence of pile displacement on soil plug capacity of open-ended pipe pile in sand 685
L. Sa, L. Grande, H. Jianchuan & L. Guohui
Deep excavations and retaining walls

3D modelling of a deep excavation in a sloping site for the assessment of induced
ground movements 693
O.J. Gastebled & S. Baghery
Analysis of an excavation in asymmetrical soil conditions: The Marquês station 699
A. Pedro, J. Almeida e Sousa, D. Taborda & P. França
Comparison of finite element predictions with results from a centrifuge test representing
a double anchor wall in sand 705
P.J. Bourne-Webb, D.M. Potts & D. König
X
Crane monopile foundation analysis 711
A. Mar
Influence of excavation and wall geometry on the base stability of excavations in soft clays 717
T. Akhlaghi, H. Norouzi & P. Hamidi
Numerical modelling of a steel sheet-pile quay wall for the harbour of Ravenna, Italy 723
D. Segato, V.M.E. Fruzzetti, P. Ruggeri, E. Sakellariadi & G. Scarpelli
Numerical modelling of spatial passive earth pressure in sand 729
M. Achmus, S. Ghassoun & K. Abdel-Rahman
Practical numerical modelling for very high reinforced earth walls 735
A. Mar, D.M. Tonks & D.A. Gorman
Short term three dimensional back-analysis of the One New Change basement in London 741
R. Fuentes, A. Pillai & M. Devriendt
Tunnels and caverns
3D analysis of a micropile umbrella for stabilizing the tunnel face of a NATM tunnel 749
F. Schmidt, C. Sagaseta & H. Konietzky
Analysis and design of a two span arch cut & cover structure 755
S. Kumar, T. Suckling, L. Macdonald & H.C. Yeow
Analysis of a bolt-reinforced tunnel face using a homogenized model 761
E. Bourgeois & E. Seyedi Hosseininia
Class A prediction of the effects induced by the Metro C construction on a preexisting
building, in Rome 767
F. Buselli, A. Logarzo, S. Miliziano & A. Zechini
Estimated settlements during the Brescia Metrobus tunnel excavation 773
A. Sanzeni, L. Zinelli & F. Colleselli
Numerical investigation of the face stability of shallow tunnels in sand 779
A. Kirsch
Numerical modeling of a bolt-reinforced tunnel in a fractured ground 785
E. Seyedi Hosseininia, E. Bourgeois & A. Pouya
On the effects of modelling gap closure and assumed soil behavior on the FE predictions of
ground movements induced by tunneling in soft clay 789
C. Miriano & C. Tamagnini
Role of numerical modelling in the current practice of tunnel and cavern design
for hydroelectric projects 795
C. Vibert, G. Colombet & O.J. Gastebled
Some modeling techniques for deep tunnels in rock with FE-continuum models 801
T. Marcher
Stress-strain behaviour of a soft-rock pillar acted upon vertical loads 807
F. Federico, S. Screpanti & G. Rastiello
Tunnel face stability with groundwater flow 813
P.M. Ströhle & P.A. Vermeer
Viscoplastic models for the analysis of tunnel reinforcement in squeezing rock conditions 819
G. Barla, D. Debernardi & D. Sterpi
Ground improvement modelling

3D FEM analysis of soil improving resin injections underneath a mediaeval tower in Italy 827
M. Gabassi, A. Pasquetto, G. Vinco & F. Mansueto
XI
A numerical study of factors governing the performance of stone columns supporting rigid
footings on soft clay 833
M.M. Killeen & B.A. McCabe
Calibration and verification of numerical model of ground improved by dynamic replacement 839
S. Kwiecien
Identification and quantification of the mechanical response of soil-wall structures
in soft ground improvement 845
X. Liu, Y. Zhao, A. Scarpas & A. de Bondt
Modelling embankments on floating stone columns 851
D. Kamrat-Pietraszewska & M. Karstunen
Numerical investigation of the mechanical behaviour of Vibro Replacement stone columns
in soft soils 857
T. Meier, E. Nacke, I. Herle & W. Wehr
Numerical modelling of consolidation around stone columns 863
J. Castro & C. Sagaseta
Numerical modeling of inertial soil-inclusion interaction 869
X. Zhang, Ph. Gotteland, P. Foray, S. Lambert & A. Hatem
Performance of geogrid-encased stone columns as a reinforcement of soft ground 875
M. Elsawy, K. Lesny & W. Richwien
Offshore geotechnical engineering

A new elasto-plastic spring element for cyclic loading of piles using the p-y-curve concept 883
O. Hededal & R. Klinkvort
Behaviour of cyclic laterally loaded large diameter monopiles in saturated sand 889
H. Ercan Taşan, F. Rackwitz & S.A. Savidis
Caisson movement caused by wave slamming—a comparison of ABAQUS and FLAC analyses 895
L. Andersen, H.F. Burcharth, T. Lykke Andersen & A.H. Augustesen
Comparison of calculation approaches for monopiles for offshore wind turbines 901
A.H. Augustesen, S.P.H. Sørensen, L.B. Ibsen, L. Andersen, M. Møller & K.T. Brødbæk
Effects of diameter on initial stiffness of p-y curves for large-diameter piles in sand 907
S.P.H. Sørensen, L.B. Ibsen & A.H. Augustesen
Numerical investigations for the pile foundation of an offshore wind turbine under
transient lateral load 913
P. Cuéllar, M. Pastor, P. Mira, J.A. Fernández-Merodo, M. Baeßler & W. Rücker
Numerical study of piping limits for installation of large diameter buckets in layered sand 921
L.B. Ibsen & C.L. Thilsted
Shallow circular foundations under undrained general combined loading in three-dimensional space 927
B. Bienen
Undrained ultimate capacity of suction anchors using an advanced constitutive model 933
S. Panayides & M. Rouainia
Numerical methods and Eurocode

Embedded cantilever retaining wall ULS design by FEA in accordance with EN 1997-1 941
A.S. Lees & S. Perdikou
Ultimate Limit State Design to Eurocode 7 using numerical methods 947
C.C. Smith & M. Gilbert
Author index 953
XII
Preface
These proceedings present 154 scientific papers written for the 7th European Conference on Numerical Methods
in Geotechnical Engineering, NUMGE 2010, held at Norwegian University of Science and Technology (NTNU)
in Trondheim, Norway from 2nd to 4th June 2010.
NUMGE 2010 is the seventh conference in a series of conferences organized by the ERTC7 (Numerical
Methods in Geotechnical Engineering) under the auspices of the International Society for Soil Mechanics
and Geotechnical Engineering (ISSMGE). The first conference in this series was held in Germany in 1986 in
Stuttgart and was followed by conferences every fourth year, every time in a new country in Europe: Spain in
1990 (Santander), United Kingdom in 1994 (Manchester), Italy in 1998 (Udine); France in 2002 (Paris) and
Austria in 2006 (Graz).
Following the traditions of the preceding conferences, NUMGE 2010 provides a forum for exchange of
ideas and discussion on topics related to geotechnical numerical modeling. Both senior and young researchers,
scientists and engineers from Europe and overseas countries have met at NUMGE 2010 to share and exchange
their knowledge.
The papers for NUMGE 2010 cover topics from emerging research to engineering practice. For the proceed-
ings the contributions are organized into the following sections:
Computer codes and algorithms
Discontinuum and particulate modelling
Large deformation – large strain analysis
Flow and consolidation
Unsaturated soil mechanics
Artificial intelligence
Reliability and probability analysis
Dynamic problems and Geohazards
Slopes and cuts
Embankments, shallow foundations, and settlements
Piles
Deep excavations and retaining walls
Tunnels and caverns
Ground improvement modelling
Offshore geotechnical engineering
Numerical methods and Eurocode
The editors would like to thank all authors for their contributions, for their cooperation during the review
process and for participating in the conference. Each paper has been reviewed by a minimum of two reviewers and
the editors are grateful for help from the reviewers in achieving quality. The national representatives in ERTC7
are thanked for promoting the conference in their respective home countries. Special thanks go to Professor
Cesar Sagaseta for keeping up the work within ERTC7.
This conference is jointly organized by NTNU, NGI/ICG, and SINTEF. These institutions and all conference
sponsors are gratefully acknowledged for their generous support. Sincere thanks go to the staff at the Geotech-
nical Division at NTNU and at the Conference Secretariat, NTNU Videre for all help in organizing NUMGE 2010.
Thomas Benz and Steinar Nordal
XIII
Scientific Committee (ERTC 7)
Chairman
C. Sagaseta, Spain
Core Members
I. Vanicek – ISSMGE Vice President Europe
P. Mestat, France
S. Nordal, Norway
M. Pastor, Spain
J. Pestana, U.S.A.
D. Potts, U.K.
H. Schweiger, Austria
S. Sloan, Australia
National Representatives
S. Aleynikov, Russia
K. Bagi, Hungary
R. Brinkgreve, The Netherlands
I. Bojtár, Hungary
A. Bolle, Belgium
H. Burd, U.K.
A. Cividini, Italy
G. Dounias, Greece
T. Edstam, Sweden
P. Fritz, Switzerland
M. Gryczmanski, Poland
O. Hededal, Denmark
I. Herle, Czech Republic
F. Kopf, Austria
T. Länsivaara, Finland
J.C. Marques, Portugal
T. Schanz, Germany
H. Walter, Austria
Local Organizing Committee at Geotechnical Division NTNU, Norway

S. Nordal, Conference chairman
T. Benz, Editor in chief conference proceedings
A. Bihs
A. Emdal
L. Grande
M. Skjåk Bræk
P. Paniagua Lopez
Reviewers not being member of ERTC 7 or the local organizing committee

L. Andresen
C. Athanasiu
M.G. Bæverfjord
B. Bostrøm
S. Degago
G. Eiksund
XV
G. Grimstad
A. Gylland
F. Hage Chehade
H.P. Jostad
S. Kirkebo
M. Leoni
R. van der Meij
R. Schwab
C. Tamagnini
V. Thakur
D. Unteregger
B.V. Vangelsten
M. de Vries
M. Wehnert
Conference secretariat NTNU Videre

A. Bye, Coordinator
XVI
A non-associated creep model for structured anisotropic clay (n-SAC)
G. Grimstad
Norwegian Geotechnical Institute, NGI, Oslo, Norway
S.A. Degago
Norwegian University of Science and Technology, NTNU, Trondheim, Norway
ABSTRACT: Elastoplastic models, based on experiments on reconstituted clays, tend to adapt an associated
flow rule. This assumption is then included in models for natural clay. However, laboratory experiments indicate
that the idea of associated flow in natural clay is insufficient. Therefore a new model, abbreviated as n-SAC,
is proposed in this paper. The model incorporates creep, using the time resistance concept, with a single creep
parameter determined from oedometer tests. Two different cap surfaces are defined in the model, i.e. the refer-
ence surface (or alternatively surface of equivalent stress measure, peq ) and the potential surface, Q. Different
(kinematic) rotational hardening rules for the two surfaces are defined along with two hardening rules for the
size of the reference surface. The two size hardening rules consists of one for the decrease in compressibility
for equivalent reconstituted material and one for loss of unstable structure. A fully implicit backward Euler
implementation scheme for the n-SAC model is used for the simulations shown in this paper.
1 BACKGROUND big consequence of wrongly predicting the horizontal

to vertical compression ratio is “bad” predictions of
Grimstad et. al. (2008) proposed a model based on settlements close to fillings.
the S-CLAY1S model (Karstunen et. al. 2005). The
model proposed in this paper has some of the elements
found in this extended S-CLAY1S model. However, 2 MODEL FORMULATION
the model differs considerably in the formulation for
flow direction. In S-CLAY1 the yield surface rota- 2.1 The reference and plastic potential surface
tion is determined from the earth pressure coefficient
Grimstad et. al. (2008) presented rate of the plastic
under virgin loading, K0NC , and the friction angle,
multiplier as a function of an equivalent stress ratio,
ϕ, (Wheeler et. al. 2003). However, Grimstad (2009)
as given by equation (1). This expression is derived
showed that this imposes unnatural limitations on the
from the time resistance concept (Janbu, 1969). The
input parameters, i.e. a “narrow” band of “good” com-
full derivation can be found in Grimstad (2009).
binations of ϕ and K0NC . The n-SAC model proposed
in this paper allows a wider range of input parameters.
Prediction of the strain behavior under various stress
paths, based on experimental evidence from e.g. Feng
(1991), indicates that the non-associated flow rule
to be reasonable. Feng conducted a series of experi- where ζi = intrinsic viscoplastic compressibility coef-
ments, on natural anisotropic clays, where both volume ficient; rsi = intrinsic time resistance number; peq
and axial deformation were measured in isotropic equivalent effective stress; pmi = intrinsic reference
compression. A decrease in both volume and height stress; τ = reference time; x = amount of unstable
were measured when considerably permanent defor- structure; and
mations were experienced. This leads to a positive
ratio between the permanent vertical and horizontal
strain increments in isotropic compression test. How-
ever, for friction angles less than about 40◦ and K0NC
from Jaky’s formula (suggested used in Wheeler et. al. where αK0NC = rotation of the potential surface in K0NC
2003) the potential surface used in S-CLAY1 gives loading (virgin oedometer loading); ηK0NC = mobi-
a negative ratio. The conclusion is that the rotational lization in K0NC loading; and MfC = critical state line
hardening rule for the potential surface of S-CLAY1 in compression loading.
and/or the shape of the potential surface needs to be The equivalent stress, peq , is calculated from the
adjusted in order to describe general soil behavior. One shape of the reference surface, which in Grimstad et. al.
3
(2008) is assumed to be associated to the potential (2002), Dafalias et. al. (2006) etc.. Wheeler et. al.
surface, Q. The shapes of these surfaces are identical (2003) argues that the rotation is dependent on the
to that of the elastoplastic Anisotropic Modified Cam deviatoric part of the plastic strain and not only the vol-
Clay Model (Dafalias, 1986), later used in S-CLAY1S umetric part as suggested by Dafalias (1986). Dafalias
(Karstunen et. al. 2005) and SANICLAY (Dafalias et. al. (2006) states three requirements for the rota-
et. al. 2006). However, unlike Dafalias (1986) and tional rules. The rotational rule proposed by Wheeler
Karstunen et. al. (2005), the n-SAC model takes a et. al. (2003), unlike Dafalias (1986), fulfills all three
similar approach as Dafalias et. al. (2006) where requirements for certain limits of input parameters
non-associated yield and potential surfaces are used. (Grimstad 2009). In the n-SAC model, the rotation
This non-association allows simulation of “softening” (rotational vector αd ) of the potential surface is depen-
response in undrained shearing without including dent solely on the volumetric strain, while the rotation
deviatoric strain dependent destructuration or spe- of the reference surface (βd ) depends generally both
cial features in the rotational hardening rule of the on volumetric and deviatoric strain. The two rotational
yield surface. The equivalent stress is calculated from hardening rules are given in equation (7) and (8).
equation (3) while the plastic potential is given in
equation (4).
where p = mean stress; σ d = deviatoric stress vec-

tor; βd = deviatoric rotational vector; M = Lode angle
dependent peak of the reference surface of in p -q
space
ζi , MC , αK0NC and ηK0NC are internal model param-

eters determined via the user input parameters given
where Mf = Lode angle dependent critical state line later in this paper.
in p -q space; αd = deviatoric rotational vector. Note Note that peq is used rather than Q for the rotation
that the precise definitions of the deviatoric vectors of the reference surface. This is introduced along with
are given in the appendix of this paper together with the Macaulay brackets to fulfill the three requirements
the Mohr-Coulomb criterion. posted by Dafalias et. al. (2006). The dependency on
In addition to introducing the non-associated flow the deviatoric strain in the rotational rule of the refer-
rule, a “creep” limit is introduced by setting the param- ence surface is included to ensure a unique rotation at
eter tmax or alternatively OCRmax below which creep “critical state”.This quality has previously been argued
should not occur. Equation (1) is then modified to that for by e.g. Wheeler et. al. (2003) and it is also used by
of equation (5). e.g. Pestana and Whittle (1999).
Modeling loss of “unstable structure” (destructura-
tion) was in a general form proposed by Gens and
Nova (1993). The same formulation was later used
in a simpler form, in e.g. Karstunen et. al. (2005)
and Grimstad et. al. (2008), in which loss of attrac-
where < > is the Macaulay brackets tion was not included. In the multilaminate framework,
Alternatively tmax could be expressed in term of Cudny and Vermeer (2004) limited destructuration to
a maximum creep induced over-consolidation ratio, be dependent on solely “normal” strain. The equivalent
OCRmax , as: approach, in standard continuum models, is visco-
plastic volume strain dependent destructuration. In
such formulations, volumetric softening can be pro-
hibited by selecting a destructuration parameter under
where OCR is defined according to the reference a certain limit. In n-SAC such requirement is auto-
time τ. matically obtained when ω = 0 in equation (9) and
rsmin > 0. However, among others Grimstad et. al.
(2008) argues that the destructuration should also be
2.2 The various hardening rules
dependent on the visco-plastic deviatoric strain, such
The precise form of the rotational hardening rules has that the “true” softening response in undrained shear-
widely been discussed by e.g. Karstunen and Wheeler ing of natural clays can be modeled. Hence equation
4
(9) is proposed as a destructuration rule for the n-SAC
model.
Isotropic hardening of the intrinsic reference stress,

pmi , is given by equation (10).
Figure 1. concept of destructuration – effect on time resis-

where ζi = a hardening parameter determined via tance number.
ref
{Eoed }i , Eref , ν and pref as given in equation (11).
The mobilization under virgin loading in oedome-

ter, ηK0NC is given by:
where λi and κ are alternative input parameters to

the model, which is in accordance with the notation
commonly used for the modified cam clay model
(MCCM), Roscoe & Burland (1968)
The value for the internal parameters βK0NC and Table 1 gives a summary of the user input parameters
αK0NC , determined from steady state rotation under to the model. In addition to these parameters the initial
K0NC loading, are given by equation (12) and (13). state variables must be given. This includes the stress
vector, σ , the initial value of the intrinsic reference
stress, pmi0 and the initial rotation vectors α0 and β0 .
Note that pmi0 can be found from the OCR and α0
and β0 from αK0NC and βK0NC . The initial structure is
defined in equation (17).
where rsi = the intrinsic time resistance number and

rsmin = the minimum measured time resistance num-
ber, see sketch in Figure 1.
In practice 11 external parameters are then left to
be determined. This can be accomplished from only
two laboratory tests. First an incremental oedometer
tests to determine rsmin and rsi has to be executed (rsi
can also be found from an incremental oedometer test
on a remolded sample). The normalized oedometer
stiffness in the NC region for a remolded sample,
ref
Equation (14)–(16) gives more details on the calcula- {E oed }i /pref , and the normalized stress dependent
tion of some of the other internal parameters. isotropic Young’s modulus, Eref /pref , should also be
established from this test.
The second test is one undrained triaxial compres-
sions test at a OCR in the interval of 1.0-OCRmax to
determine ω (deviatoric destructuration parameter), ϕp
(friction angle at peak of the reference curve) and ϕcs
where ϕcs is the critical state friction angle
(the critical state friction angle). Note that the normal-
ized stiffness, Eref /pref , could also be found from this
undrained test.
The Poisson ratio, ν, is set by default to 0.15. How-
ever, if one for instance has radial stress measurement
where ϕp is the friction angle at the peak of the in the oedometer test it can be determined from the
reference curve. stress path in the OC region. The reference time, τ, is
5
usually 1 day, since it is common to determine OCR Table 1. Model input parameters.
for the 24 h load duration in the incremental oedometer
ref
test. ν K0NC Eref /pref {E oed }i /pref rsmin rsi
0.15 0.55 160 16 100 800

3 IMPLEMENTATION SCHEME
ω ϕp ϕcs τ tmax
The n-SAC model is implemented in an implicit back-
ward Euler integration scheme. Procedures for implicit
0.3 25◦ 33◦ 1 day 1e3yr
backward Euler implementation of general viscoplac-
tic models may be found in for instance de Borst and
Heeres (2002). In this particular case the following
residuals are defined (equation (18) to (23)):
where Dn+1 is the mean stress dependent isotropic

elasticity matrix under the assumption of a constant
Poisson ratio at the end of the calculation step
Figure 2. Result of test 1 in p − q space.
Typically the scheme will converge in a few iterations.

However, to improve performance sub-stepping will
be initiated if the number of iteration becomes higher
than 50 or if the estimated reciprocal condition number
of ∂rin+1 /∂vin+1 becomes less than 1E-12.
The 21 state variables, including the effective stresses,

are stored in a vector v:
4 NUMERICAL PERFORMANCE
In order to speed up the implementation process the In order to test the performance of the model and
number of state variables could be reduced by making implementation three different tests where ran with
use of the dependencies, i.e. αd,xx + αd,yy + αd,zz = 0, the input parameters found in Table 1. The initial state
βd,xx + βd,yy + βd,zz = 0. variables were generated from an initial vertical stress
To find the new state, a standard iterative Newton- of 72.3 kPa and an OCR of 1.383.
Raphson scheme is used as given in equation (25) Test 1 contained 181 radial strain paths of 30 steps
to (27). The iteration is ran until r∗n+1 T · r∗n+1 < TOL2 . with time increments of 0.1 day for each step given
Where r∗n+1 is a normalized version of the residual vec- in εv − εq space under the condition that dε2v + dε2q =
tors, rn+1 . The normalization is done in such a way that 1e-6. The result of test 1 is given in p − q space in
the tolerance check is irrespective of the magnitude Figure 2.
and dimension of the state variables. TOL = 1E-6 is Test 2 and test 3 are both undrained tests, consist-
used in the particular simulations shown in this paper. ing of 91 paths of 100 steps with time increments
6
Figure 5. Comparing 100 and 10 steps for a simulation of
an undrained triaxial compression test.
peak or at 6% shear strain) are 0.51 and 0.70 respec-

tively. This is in the range reported by Bjerrum (1973)
Figure 3. Result of test 2 in σd1 − σd2 − σd3 space. and by Whittle (1993).
6 CONCLUSIONS AND RECOMMENDATIONS
This paper describes a constitutive model which

includes a variety of effects observed in clay behav-
ior. It has relatively few additional input parameters,
to that of for instance the Anisotropic Creep model,
ACM (Leoni et. al. 2008). The few “extra” input can
easily be calibrated from standard laboratory tests.
The paper presents the implementation scheme and
a variety of model simulations, i.e. constant rate of
Figure 4. Result of test 3 in plane strain deviatoric space. strain tests under various strain paths. The implemen-
tation is under these test conditions is illustrated to be
of 0.001 day for each step under the condition that sufficiently robust.
dε1 − dε3 = 6e-4. Test 2 involves only εz , εy and εz , The model and the implementation will go through
while test 3 were ran under plane strain condition (i.e. further testing, specifically trying to reproduce actual
involving εz , εy and γxy ). Results are shown in Figure 3 measured behavior in laboratory and in the field. Such
and Figure 4 in the deviatoric stress space correspond- cases will involve both settlement predictions of actual
ing to the respective strain space in which increments or hypothetical cases, compared to that of measure-
were applied. ments and along with the response predicted by other
In order to compare accuracy in the integration similar models.
method, the three tests were then repeated with 1/10 The model will hopefully be a candidate for fur-
of the number of steps used for previous simulations. ther studies on the modeling of clay behavior. Further
Result of one such comparison, for undrained triax- development to the n-SAC model could be to improve
ial compression) is shown in 5. This particular case is the behavior at small strains by incorporating non-
showing satisfactory results. linear response at small strain and elastic anisotropy.
5 DISCUSSION ACKNOWLEDGEMENTS
The particular input shown in Table 1 gives an ini- The work presented was partly carried out as a part
tial value of αNC
K0 of 0.1122 (equation (12)). This value of project 5, Geomechanical modeling, at the Inter-
for α gives a maximum ratio of horizontal to verti- national Centre of Geohazards, ICG, a Centre of
vp vp
cal visco-plastic strain (εh /εv ) of 6.7 in an isotropic Excellence (CoE) with funding from The Research
consolidation test. This is close to that measured in Council of Norway. Most of the work was finished
for instance Batiscan clay by Feng (1991). Associated when both the authors were PhD students at NTNU
anisotropic critical state models for clays do not nec- under the supervision of Professor Steinar Nordal.
essary guarantee a positive ratio for this case, or they Nordal is acknowledged for his contributions in dis-
fail to reproduce the measured K0NC − ϕ combination. cussing the content of this paper. The Marie Curie
The suP /suA and the suDSS /suA ratios (found in Figure 4 at Research Training Network “Advanced Modeling of
7
Ground Improvement on Soft Soils (AMGISS)” (Con- APPENDIX
tract No MRTN-CT-2004-512120) supported by the
European Community through the program “Human Definition of the deviatoric stress and rotational vec-
Resource and Mobility” is also acknowledged. tors:
REFERENCES
Bjerrum, L., 1973. Problems of Soil Mechanics and Con-
struction on Soft Clays, State of the Art Report to Session
IV, 8th International Conference on Soil Mechanics and
Foundation Engineering, Moscow, also in NGI report 100
(1974).
Cudny, M. and Vermeer, P. A. 2004. On the modelling of
anisotropy and destructuration of soft clays within the
multi-laminate framework, Computers and Geotechnics
31: 1–22.
Dafalias, Y. F. 1986. An anisotropic critical state soil plas-
ticity model, Mechanics research communications 13(6):
341–347.
Dafalias,Y. F., Manzari, M. T. and Papadimitriou, A. G. 2006.
SANICLAY: simple anisotropic clay plasticity model, Int.
J. Numer. Anal. Meth. Geomech. 30: 1231–1257.
de Borst, R. and Heeres, O. M. 2002,A unified approach to the
implicit integration of standard, non-standard and viscous
plasticity models, Int. J. Numer. Anal. Meth. Geomech. 26:
1059–1070.
Feng, T.W. 1991. Compressibility and permeability of natu-
ral soft clays and surcharging to reduce settlements. PhD
diss., University of Illinois at Urbana-Champaign, Urbana
Illinois.
Gens, A. and Nova, R. 1993. Conceptual bases for a constitu-
tive model for bonded soils and weak rocks, Geotechnical
Engineering of Hard Soils - Soft Rocks, Anagnostopoulos
et. al. (eds) Balkema, Rotterdam.
Grimstad, G. 2009. Development of effective stress based
anisotropic models for soft clays, PhD diss., Norwegian
University of Science andTechnology, NTNU,Trondheim.
Grimstad, G., Degago, S., Nordal, S. and Karstunen, M. 2008. The two Mohr Coulomb criteria in p − q space are
Modelling creep and rate effects using the time resis-
tance concept in a model for anisotropy and destructura-
given for the critical state (potential surface):
tion, Nordic Geotechnical Meeting, Sandefjord, Norway,
195–202.
Janbu, N. 1969. The resistance concept applied to deforma-
tions of soils. Proc. 7th Int. Conf. on Soil Mech. & Found.
Eng, Mexico city 1: 191–196.
and for the peak of the reference surface:
Karstunen, M. and Wheeler, S. 2002. Discussion of “Finite
Strain, Anisotropic Modified Cam Clay Model with Plas-
tic Spin. I: Theory” by George Z. Voyiadjis and Chung R.
Song. Journal of Engineering Mechanics, ASCE 128:
497–498.
Karstunen, M., Krenn, H., Wheeler, S. J., Koskinen, M. and
Zentar, R . 2005. The effect of anisotropy and destructura- where θ α = Modified lode angle (function of σ d and
tion on the behaviour of Murro test embankment. ASCE αd ) and θ β = Modified lode angle (function of σ d and
International Journal of Geomechanics 5(2): 87–97. βd )
Leoni, M., Karstunen, M. & Vermeer, P. A. 2008. Anisotropic The two modified Lode angles are calculated from
creep model for soft soils, Géotechnique 58(3): 215–226. the middle eigenvalues S α and S β of the tensors, sα2 and
Pestana, J. M. and Whittle, A. J. 1999. Formulation of a uni- β
fied constitutive model for clays and sands, Int. J. Numer. s2 , such that: sin (θ α ) = 3/2 · S2α /qα and sin (θ β ) = 3/2 ·
β
Anal. Meth. Geomech. 23: 1215–1243. S2 /qβ , where:
Roscoe, K. H. and Burland, J. B. 1968. On the generalized
stress-strain behavior of wet clay, Engineering plasticity,
535–609, Cambridge university press.
Wheeler, S. J., Näätänen, A., Karstunen, M. and Lojander, M.
2003. An anisotropic elastoplastic model for natural soft
clays, Canadian Geotechnical Journal 40: 403–418.
Whittle, A. J. 1993. Evaluation of a constitutive model for
overconsolidated clays, Géotechnique 43(2): 289–313.
8
A state dependent constitutive model for sand-structure interfaces
A. Lashkari
School of Civil and Environmental Engineering, Shiraz University of Technology, Shiraz, Iran
ABSTRACT: Numerical simulations of geotechnical engineering problems depend strongly on predictive

capacity of constitutive models used for modeling of soil-structure interface behavior. It has been shown that the
mechanical behavior of sand-structure interfaces is significantly influenced by the combined effect of density
and normal stress which is called the effect of interface state. In this paper, a two-surface interface model within
the context of bounding surface plasticity is presented. Several ingredients of the model are defined as direct
functions of interface state. It is shown that the model is capable of distinguishing interfaces in dense state from
loose ones and providing reasonable predictions.
1 INTRODUCTION later introduced by Gajo & Wood (1999). Hereto-

fore, the constitutive model of Manzari & Dafalias
The mechanical behavior of soil-structure interfaces (1997) has been widely developed and applied to var-
is an influential factor in the load-deformation and the ious problems in sands constitutive modeling includ-
bearing capacity of almost all of geotechnical engi- ing liquefaction under cyclic loadings (Dafalias and
neering problems such as retaining structures, rein- Manzari, 2004), inherent anisotropy (Dafalias et al.,
forced soils, piles, under ground and offshore gravity 2004), Non-coaxiality and liquefaction under rota-
structures. A number of experimental techniques have tional shear (Li & Dafalias, 2004; Lashkari &
been suggested in the literature to study the mechani- Latifi, 2007; Lashkari, 2009a), and unsaturated sands
cal behavior of soil-structure interfaces under different (Chiu & Ng, 2003). In the following lines, an interface
stress paths and stiffness boundary conditions. Among model based on the Manzari & Dafalias (1997) plat-
them, Ghionna & Mortara (2002) and Hu & Pu form is suggested. Comparisons between the model
(2004) used direct shear tests, Zeghal & Edil (2002) predictions and experimental results are presented to
employed ring shear apparatus, and Evgin & Fakhar- demonstrate the predictive capability of the model
ian (1996) used simple shear apparatus. On the other for simulation of the state dependent behavior of
hand, considerable developments have been achieved interfaces.
in the field of constitutive modeling of interfaces. In
this regard, Clough & Duncan (1971) suggested an
interface model based on hyperbolic elasticity theory. 2 GENERAL FORMULATION OF THE MODEL
Ghaboussi et al. (1973) suggested an elasto-plastic
cap model for soil-structure interfaces. De Gennaro & 2.1 Definitions of spaces
Frank (2002) considered phase transformation and
In this work, stress vector is defined as
ultimate state in their elasto-plastic interface model.
Ghionna & Mortara (2002) introduced a Cam-Clay
type model for sand-structure interfaces. Recently,
Liu et al. (2006) suggested that the Critical State where σn and τ are, respectively, the compressive nor-
Soil Mechanics concepts can be extended to include mal and tangential components of the stress vector.
the mechanical behavior of interfaces. They also pro- Corresponding to Equation (1), relative displacement
posed a state dependent generalized plasticity inter- vector becomes
face model using the state parameter of Been and
Jefferies (1985). More recently, Lashakari (2010) sug-
gested a generalized plasticity interface model for where v and u are normal and tangential displace-
sand-structure interfaces subjected to rotational shear. ments, respectively.
Following the impressive suggestion of Wood et al. Experimental studies have revealed that the thick-
(1994), Manzari & Dafalias (1997) introduced a ness of the interface zone is 5–10 times of mean grain
two-surface critical state compatible bounding sur- diameter of soil. In this study it is assumed that
face model for state dependent behavior of sands.
A similar constitutive model, Severn-Trent, has been
9
where t is the interface thickness and d50 is mean grain 3 THE MODEL SPECIAL ELEMENTS
diameter of grains in contact with structure. Assuming
that strains are uniformly distributed in the interface 3.1 Yield and plastic strain rate direction vectors
zone, strain vector can be defined as
Similar to the original platform of Manzari & Dafalias
(1997), a narrow wedge-shape yield function is
adopted here
where εn and εt are normal and tangential strains with

respect to the interface plane.
where η = τ/σn is stress ratio. In τ − σn plane, α is
back stress ratio which is the slope of the bisector of
the yield criterion with respect to the positive direction
2.2 Constitutive equation between stress of the σn -axis (Fig. 1). Finally, m indicates the yield
and relative displacement rate vectors surface size. In practice, m = 0.01 M is a reasonable
selection, where M is the slope of critical state line in
In the elasto-plasticity theory, the strain rate vector can τ-σn plane. Based on the yield surface introduced in
be decomposed into elastic and plastic parts Equation (10), {n} becomes
where superscripts “e” and “p” stand for the elastic where s = +1 when η > α , and s = −1 when η < α.
and plastic parts of the strain rate vector, respectively. Analogous to the work of Manzari & Dafalias (1997),
Hereafter, dot sign, . , on each parameter defines the a non-associated flow is introduced through the fol-
rate of corresponding parameter with time. In addition, lowing definition for {R}
one has the following constitutive equation between
the rates of stress and relative displacement vectors in
the elasto-plasticity theory
where Rn plays the role of dilatancy, d, in the proposed

ep
model flow rule.
where [D] is the elasto-plastic stiffness matrix which
is calculated by
3.2 State dependent peak and phase transformation
stress ratios
Liu et al. (2006) have suggested that the Critical State
Soil Mechanics concepts can also be applied to rough
interfaces. To this aim, the application of some proper
In Equation (7), [D]e is the elastic stiffness matrix
where Kn and Kt are respectively the interface normal

and tangential elastic moduli. Similar to sands, it is
assumed that these moduli are pressure dependent
Kn0 and Kt0 are model parameters and pref = 101 kPa
is the atmospheric pressure which plays the role of a
reference stress.
Kp is plastic hardening modulus. {n} and {R} are
two vectors which defines yield and plastic strain rate
directions. Particular definitions for these terms are
presented in sequel. Figure 1. Illustration of the model constitutive surfaces.
10
state parameters into the model formulation is essen- where
tial. State parameters are commonly used in order to
describe the current state of a soil or soil-structure
interface uniquely. Up-to-date reviews on a number of
state parameters can be found in Dafalias and Man-
zari (2004) and Lashkari (2009b). According to the
In above equation, A0 and A1 are model parameters,
latter work, the following constitutive equations are
where it is worthy to note that A0 is usually larger than
suggested here for state dependent peak and phase
A1 . It is observed that when a new tangential loading
transformation stress ratios
starts, the mechanical behavior of both loose and dense
interfaces is significantly contractive. Subsequently,
the mentioned contraction decreases and may turns
into dilation at moderate and large tangential displace-
ments. Considering Equation (19), α = αin < αd at the
Table 1. Model parameters used in model predictions shown

in Figs. 2–7.
where nb and nd are model parameters, and Parameter M-M-G * E-F ** S-R ***
Kt 0 (MPa) 1.57 5.0 3.0

Kn 0 (MPa) 1.85 5.85 3.60
A0 2.0 8.0 5.0
A1 0.7 1.4 0.4
In the above equation, emax , emin , and e are respec- h0 0.38 0.25 0.6
tively the maximum possible, minimum, and current M 0.63 0.638 0.64
amounts of the interface void ratio. In Equations (13) e0 0.787 1.01 0.85
and (14), σnc is normal stress corresponding to the cur- λ 0.0557 0.09 0.074
rent amount of void ratio. The location of critical state nb 0.110 0.12 0.35
line in e-σnc plane is defined by nd 0.055 0.0674 0.10
*: experiments reported by Mortara et al. (2007).

**: experiments reported by Evgin & Fakharian (1996).
***: experiments reported by Shahrour & Rezaie (1997).
In equation (16), e0 and λ are model parameters.
3.3 Dilatancy function and plastic hardening

modulus
The model constitutive surfaces are demonstrated in
Fig. 1. Similar to the yield surface, one can define back
stress ratios corresponding to bounding and dilatancy
surfaces
where αb and αd are respectively the bounding and

dilatancy back stress ratios.
Now, the model plastic hardening modulus is
where h0 is a model parameter and αin is the initial

amount of α when the most recent tangential loading
has started.
Dilatancy function is defined in the following form
Figure 2. The model predictions compared with three con-
stant normal stress tests [experimental data taken from
Mortara et al. 2007].
11
stant normal stress tests [experimental data taken from
Evgin & Fakharian 1996].
Mortara et al. (2007) published the results of a

series of Gioia Tauro sand-steel interface tests car-
ried out with direct shear apparatus. Samples were
prepared by tamping method with initial relative den-
sity ID = 60%. The physical properties of the GT40
fraction of Gioia Tauro sand considered here are:
Figure 3. The model predictions compared with three con- Gs = 2.69, emax = 0.96, and emin = 0.60.
stant normal stiffness tests [experimental data taken from The model of this study is calibrated versus exper-
Mortara et al. 2007]. imental data of Mortara et al. (2007). Amounts of the
model parameters are given in Table 1. The model
predictions on tangential Stress vs. horizontal dis-
start of√a tangential loading and therefore, one has placement and normal displacement vs. horizontal
A = A0 pref /σn . As the loading proceeds, the ratio displacement for three constant normal stress tests
(α − αin )/(sαb − αin ) approaches toward 1 and leads to with σn = 25, 150, and 300 kPa are illustrated against
A = A0 at large tangential loading. The particular def- experimental data in Fig. 2. In addition, for three
inition of A given in Equation (19) enables the model constant normal stiffness tests with the same initial
for proper prediction of volume change behavior. stress conditions and K = 1.0 GPa/m, the model pre-
dictions for tangential stress vs. horizontal displace-
ment, tangential stress vs. normal stress, and normal
4 THE MODEL EVALUATION displacement vs. horizontal displacement are shown in
Fig. 3.
In experimental studies, stiffness boundary condition Evgin & Fakharian (1996) reported the results of a
in the direction normal to the interface plane is usually series of Ottawa sand-rough steel interface tests car-
defined by ried out using modified simple shear apparatus. The
physical properties of Ottawa sand are: Gs = 2.65,
emax = 1.024, and emin = 0.651. Grains are angular and
mainly made of quartz. The initial density of samples
In experiments, two types of stiffness conditions are is 84%. Using parameters presented in Table 1, simi-
usual lar comparisons are shown for evaluation of the model
versus this set of experiments. The model predictions
1. K = 0 (constant normal stress condition) in which are compared with experimental data in Fig. 4 for
σ̇n = 0 and v̇ = 0 three constant normal stress tests σn = 100, 300, and
2. K = constant (constant normal stiffness condition) 500 kPa. For three other tests carried out under con-
in which σ̇n = 0 and v̇ = 0 stant normal stiffness condition K = 800 kPa/m and
12
Figure 6. The model predictions versus experimental
results for two dense (ID = 90%) and loose (ID = 15%)
Hostun sand-steel interfaces subjected to constant normal
stress σn = 100 kPa condition [experimental data taken from
Shahrour & Rezaie 1997].

stant normal stiffness tests [experimental data taken from
Evgin & Fakharian 1996].
σn = 100, 200, and 300 kPa, the model simulations

are depicted against corresponding experiments in
Fig. 5.
Using direct shear apparatus, Shahrour & Rezaie
(1997) studied the mechanical behavior of Hos-
tun sand-steel interfaces with different initial den-
sities (ID0 = 15, and 90%). The physical properties
of the Hostun sand are: Gs = 2.65, emax = 1.00, and
emin = 0.653. In all simulations, the interface thickness
is assumed 7 mm. For two dense and loose interfaces
subjected to constant normal stress σn = 100 kPa con-
dition, the model predictions are compared with exper-
imental data in Fig. 6. Similar comparisons are made
for samples with σn = 300k̇Pa in Fig. 7. For predictions
presented, the model parameters are given in Table 1. Figure 7. The model predictions versus experimental
The model has predicted a peak in tangential strength results for two dense (ID = 90%) and loose (ID = 15%)
for interfaces in dense state. In addition, the model Hostun sand-steel interfaces subjected to constant normal
is correctly predicted that the maximum tangential stress σn = 300 kPa condition [experimental data taken from
strength of loose samples is attained at large horizontal Shahrour & Rezaie 1997].
13
displacements. Moreover, the volume change behav- Evgin, E. & Fakharian, K. 1996. Effect of stress path on the
ior of loose interfaces is always contractive. On the behavior of sand-steel interface. Canadian Geotechnical
other hand, dense interfaces demonstrate contraction Journal 33: 853–865.
initially which turns into dilation in moderate and large Gajo, A. & Wood, D.M. 1999. Severn-Trent sand: a
kinematic-hardening constitutive model: the q-p formu-
horizontal displacements. From both figures, it can be lation. Géotechnique 49(5):595–614.
observed that the model can capture the fundamental Ghaboussi, J. Wilson, E. L. & Isenberg, J. 1973. Finite element
aspects of interfaces behaviors. for rock joints and interfaces. J. Soil Mech. & Found. Div.
ASCE 99 (SM10): 833–848.
Ghionna, V.N. & Mortara, G. 2002. An elastoplastic model
5 CONCLUSIONS for sand-structure interface behavior. Géotechnique 52(1):
41–50.
Within the frameworks of bounding surface plasticity Hu, L. & Pu, J. 2004. Testing and modeling of soil-structure
and Critical State Soil Mechanics, a state dependent interface. ASCE Journal of Geotechnical and Geoenviron-
mental Engineering 130(8): 851–860.
sand-structure interface model was presented. The Lashkari, A. 2009a. A constitutive model for sand liquefac-
constitutive model of Manzari and Dafalias (1997) was tion under rotational shear. Iranian Journal of Science &
selected as platform. Technology, Transaction B, Engineering 33 (B1): 31–48.
New elements for state dependency and dilatancy Lashkari, A. 2009b. On the modeling of the state depen-
were suggested. Employing the data reported by three dency of granular soils. Computers and Geotechnics 36:
independent research teams, the model predictions 1237–1245.
were compared with experimental results under var- Lashkari, A. 2010. Modeling sand-structure interfaces under
ious stiffness boundary conditions. Using a unique rotational shear. Mechanics Research Communications
set of parameters for each type of soil-structure inter- 37: 32–37.
Lashkari, A. & Latifi, M. 2007. A constitutive model for
face, it has been shown that the model is capable of non-coaxial flow of sand. Mechanics Research Commu-
providing reasonable predictions for samples of dif- nications 34: 191–200.
ferent initial states subjected to loading under various Li, X.S. & Dafalias, Y.F. 2004. A constitutive framework
stiffness boundary conditions. for anisotropic sand including non-proportional loading.
Géotechnique 54 (1): 41–55.
Liu, H., Song, E. & Ling, H. I. 2006. Constitutive modeling of
REFERENCES soil-structure interface through the concept of critical state
soil mechanics. Mechanics Research Communications 33:
Been, K. & Jefferies, M. G. 1985. A state parameter for sands. 515–531.
Géotechnique 35(2): 99–112. Manzari, M.T. & Dafalias, Y.F. 1997. A critical state two
Dafalias, Y.F. & Manzari, M.T. 2004. Simple plasticity sand surface plasticity model for sands. Géotechnique 47(2):
model accounting for fabric change effects. ASCE Journal 255–272.
of Engineering Mechanics 130(6): 622–634. Mortara, G. Mangiola,A. & Ghionna,V. N. 2007. Cyclic shear
Dafalias, Y.F. Papadimitriou, A.G. & Li, X.S. 2004. Sand stress degradation and post-cyclic behaviour from sand-
plasticity model accounting for inherent fabric anisotropy. steel interface direct shear tests. Canadian Geotechnical
ASCE Journal of Engineering Mechanics 130(11): Journal 44: 739–752.
1319–1333. Shahrour, I. & Rezaie, F. 1997. An elastoplastic constitu-
De Gennaro, V. & Frank, R. 2002. Elasto-plastic analysis tive relation for the soil-structure interface under cyclic
of the interface behavior between granular media and loading. Computers and Geotechnics 21(1): 21–39.
structure. Computers and Geotechnics 29: 547–572. Wood, D.M. Belkheir, K. & Liu D. F. 1994. Strain soften-
Chiu, C.F. & Ng, C.W.W. 2003. A state dependent elastoplas- ing and state parameter for sand modeling. Géotechnique
tic model for saturated and unsaturated soils. Géotech- 44(2):335–339.
nique 53 (9): 809–829. Zeghal, M. & Edil,T. 2002. Soil structure interaction analysis:
Clough, G.W. & Duncan, J.M. 1971. Finite element analysis modeling the interface. Canadian Geotechnical Journal
of retaining wall behavior. J. Soil Mech. & Found. Div. 39: 620–628.
ASCE 97 (SM12): 1657–1672.
14
Adaptive integration of hypoplasticity
W. Fellin
Division of Geotechnical and Tunnel Engineering, Department of Infrastructure,
University of Innsbruck, Austria
M. Mittendorfer & A. Ostermann

Department of Mathematics, University of Innsbruck, Austria
ABSTRACT: Evaluating the stress response of a constitutive relation of the rate type for a given strain increment
can be seen as a time integration. The question whether explicit or implicit methods should be used for this
integration is controversially discussed in the literature. In our previous paper (Fellin et al. 2009), we have
analysed two adaptive second order methods, constructed by extrapolation of the explicit and a semi-implicit
Euler method, respectively. Here, we compare their numerical behaviour on two geotechnical finite element tests.
As constitutive relation, we use hypoplasticity with the intergranular strain concept.
1 INTRODUCTION 2 ADAPTIVE INTEGRATION
The reliable computation of the stress response for a The constitutive rate equations of hypoplasticity
given strain increment is an important issue in compu-
tational geotechnics. Hypoplasticity (Kolymbas 1985)
is a framework for constitutive models of the rate
type specialised for soil behaviour. For the compar- form a system of ordinary differential equations,
ison of the numerical behaviour of an adaptive second see (Fellin et al. 2009) and references therein. Here, D
order semi-implicit method (Fellin et al. 2009) and denotes the Eulerian stretching, T the effective Cauchy
an adaptive second order explicit method (Fellin and stress and Q the additional state or internal variables.
Ostermann 2002), we choose the hypoplastic model Collecting the components of T, Q, and their deriva-
with the intergranular strain concept (Niemunis and tives with respect to the stretching in a vector y, we
Herle 1997). For single element tests the solutions obtain a nonlinear initial value problem
show a stiff behaviour all over the computational
domain (Mittendorfer 2010). Stiffness in the mathe-
matical sense means that certain implicit integrators
perform much more efficiently than explicit ones. For
adaptive explicit methods applied to stiff problems the Its efficient and reliable numerical solution is an
product of the time step size with the dominant eigen- essential step in solving the equilibrium equations.
value of the linearised system lies near the border of In our recent article (Fellin et al. 2009), we have
the stability domain (Hairer and Wanner 1996). This discussed in detail how (2) can be solved efficiently.
can enforce very small step sizes, already observed in Here, we resume briefly two attractive second order
the single element tests in (Fellin et al. 2009). methods that are both endowed with an error estimate
Motivated by these element tests, we compare here and an adaptive step size strategy.
the performance of the proposed integration schemes
for typical geotechnical problems. For this purpose, 2.1 A second order explicit method
we implement a user subroutine for the finite ele-
ment package Abaqus and we compare the numerical Starting from a numerical approximation yn ≈ y(tn ) at
behaviour on two typical finite element problems from time tn , the explicit Euler method
geotechnics: a biaxial test and a sheet pile wall exam-
ple. For these examples not only the behaviour of the
solution is important but also the structure of the spec-
imen. If the stresses in an integration point remain yields a numerical approximation yn+1 at time
roughly constant, only few time steps are required and tn+1 = tn + τn . Due to its simplicity, the method is still
the implicit method cannot exploit its advantages. much in use for integrating (2). Its main drawbacks,
15
however, are its low accuracy and the lacking error 2.3 Error estimation and step size control
control.
Next, we will treat the problem of step size selection.
A simple combination of two consecutive Euler
Our approach is that of (Fellin and Ostermann 2002;
steps, combined with a local extrapolation procedure,
Fellin et al. 2009; Hairer et al. 1993). The difference
avoids both of these drawbacks without destroying the
of the auxiliary values (4) and (7), respectively,
simplicity of the method. In the following, we briefly
describe this method.
Starting from yn , we first perform an Euler step of
size τn
is an asymptotically correct estimate for the local error
of w. For a user-supplied tolerance TOL, we obtain
an optimal step size τopt . We use this for controlling
as well as two Euler steps of size τn /2 the step size. If the estimated error EST is below the
tolerance TOL, the step is accepted and a new larger
step size is chosen for the next step. If the estimated
error EST is larger than TOL, however, we reject the
step and redo it with a smaller step size.
In order to obtain a reliable error estimate, it is
common to use the maximum norm in (9)
Taylor expansions show that the combination
is a second order approximation to the solution. with the scaling factors

The resulting method (5) is called Richardson
extrapolation of the explicit Euler method. In this
paper, we will use the name ERK2 henceforth.
The parameters ai and ri are used to fine-tune the error
2.2 A second order semi-implicit method estimate concerning the absolute and relative error
For stiff problems, explicit methods like ERK2 become tolerances for each entry of y.
inefficient. As a remedy, we proposed in (Fellin et al.
2009) to replace the Euler steps in the construction of
(5) by semi-implicit Euler steps of the form 3 FINITE ELEMENT EXAMPLES
The performance of the proposed time integration

schemes is shown with finite element calculations car-
ried out with Abaqus. The integration schemes were
In contrast to the fully implicit Euler method, this
implemented together with the constitutive equations
method only requires the solution of linear systems
and the calculation of the consistent tangent stiffness in
of equations.
a user defined material subroutine, the so called Umat.
To obtain a second order integration scheme, we
We use the default convergence criteria and load incre-
define the auxiliary variables v and w as before. Start-
mentation of Abaqus. The material parameters for the
ing from the approximation y ≈ y(tn ), we compute a
calculations are listed in Table 1.
step of size τ
Table 1. Parameters for the extended hypoplastic model.
Parameter Values
as well as two steps of size τn /2
ϕc [◦ ] 33
hs [kPa] 1 × 106
n 0.25
ed0 0.55
ec0 0.95
ei0 1.05
α 0.25
The extrapolated value β 1.50
R 1 × 10−4
mR 5.0
mT 2.0
βr 0.5
is the searched second order approximation which will χ 6.0
be called SIRK2 in the remainder of this paper.
16
The following values for the tolerances are used in large strains and large rotations, so the lateral pressure
all numerical experiments: TOL = 10−3 and remains perpendicular to the edge of the specimen.
The calculations were performed without gravity. In
this case, the problem is symmetric with respect to the
horizontal direction.
Figure 2 is a quilt plot for the void ratio e at the end
with AERRi being the lowest resolution of the compo- of the test. The void ratio increased from the initial
nent yi . We set this value for the state variables to: 0.1 value 0.569 in the dashed area up to maximum 0.708
for the stress, 0.01 for the void ratio and 10−6 for the in the shear bands, which are formed during the test.
intergranular strain. The derivatives of the state vari- Both numerical integrations, ERK2 and SIRK2, give
ables are needed to calculate the consistent tangent the same plot.
stiffness (Fellin and Ostermann 2002). These deriva- The time steps used at the end of the test for
tives are included in the error estimation and step size ERK2 and SIRK2 are shown in Figures 3 and 4.
control, with AERR equal to 0.1. The weighting fac- For ERK2 small time step sizes are required in the
tors ri are set to one for state variables and to 100 for shear bands, whereas in other regions the time steps
the derivatives of them. are much larger. This indicates the adaptivity of our
Abaqus uses the time t as parameter throughout the method, which accounts for the stiff behaviour of
calculation. To distinguish between this time and the the constitutive equations. Due to the comparatively
time steps τ used in the constitutive time integration, large deformations in the shear bands combined with
we shall call changes of loads or boundary conditions the stiff behaviour of the hypoplastic equations the
load steps and increments of them load increments t. explicit integrator requires small time steps to meet
the accuracy requirements. It is worth to note that the
constitutive equations show stiff behaviour outside the
3.1 Biaxial test shear band as well. As the strains are rather small in
these regions, the stress remains roughly constant, and
We start with a biaxial test as standard geotechnical therefore larger steps can be accepted.
benchmark example (Hügel 1995). A soil specimen of The semi-implicit method, however, shows a rather
0.04 m width and 0.14 m height is laterally confined balanced allocation of the step sizes over the whole
with a constant stress of size 400 kN/m2 with plain structure, see Figure 4. There are barely elements in
strain condition in the other horizontal direction. The the shear bands where the step sizes are considerably
specimen is compressed vertically by a prescribed dis- smaller than in other regions. However, they are much
placement u = 0.01 m. The material in the dashed area larger than that used by the ERK2 method in the shear
of Figure 1 is given an initially higher void ratio of band.
0.569, whereas the void ratio is 0.506 elsewhere. In An “exact” solution of the load displacement curve,
this way, an initial imperfection in the dashed area is which was obtained with 200 load increments, and
simulated. The dashed area is of the size 0.02 m by the behaviour of the automatic load incrementation
0.02 m. strategy of Abaqus are shown in Figure 5. The total
The biaxial test is modelled with 8 by 28 linear
plane strain elements. The calculation accounts for
Figure 2. Biaxial compression test without gravity: void

Figure 1. Biaxial test. ratio at the end of the test.
17
Figure 5. Biaxial compression test without gravity: load
displacement curve; the continuous line is an “exact” solu-
tion obtained with 200 load increments, the circles denote
the increments of the automatic load incrementation strategy,
which are the same for both methods ERK2 and SIRK2.
Figure 3. Biaxial compression test without gravity: time

step size τ of ERK2 at the end of the test.
Figure 6. Biaxial compression test without gravity: total

number of time steps.
post peak behaviour is not studied here. The solu-

tion obtained by the automatic load incrementation
matches quite well with the exact solution over the
whole time window.
Figure 4. Biaxial compression test without gravity: time The mean value of the conducted time steps in the
step size τ of SIRK2 at the end of the test.
four integration points of the plain strain per element
per increment are summed over all Abaqus load incre-
loading time is t = 1. The initial and maximum load ments, see Figure 6. The minimum number of required
increments are t = 0.1. steps is nearly equal for both methods. The time inte-
In this test Abaqus requires 20 load increments gration with ERK2 integration needs in the shear band
with the same sequence for both time integration meth- much more time steps in total than the integration with
ods. There is some decreasing of the load increments SIRK2. Switching from ERK2 to SIRK2 in those ele-
around the peak due to the inherent mechanical diffi- ments would make sense. However, the number of
culties of material softening. However, the increments elements where the implicit method needs consider-
are increased by Abaqus in the softening branch due to ably fewer time steps than the explicit method is rather
rapid convergence of the equilibrium iterations. This small.
indicates the good performance of the implemented Significant deformations take part only in small
consistent tangent stiffness operator (Fellin and Oster- regions of the computational domain of the biaxial
mann 2002). Note that the mesh dependence of the test. In large regions of the specimen, the imposed load
18
Table 2. Comparison of computational costs in the biaxial
test without gravity. Abaqus: automatic load incrementation
with tstart = tmax = 0.1.
Int.method No. of Inc. No. of It. CPU [s]
ERK2 20 79 200.3
SIRK2 20 83 621.8
Figure 8. Sheet pile wall example.
Table 3. Comparison of computational costs in the sheet

pile wall test. Abaqus: automatic load incrementation with
tstart = tmax = 0.1.
Int.method No. of Inc. No. of It. CPU [s]
ERK2 21 56 28.3
SIRK2 21 52 98.3
Figure 7. Biaxial compression test without gravity: relative

number of time steps per integration call.
increments are small. For such elements, the time step

size of ERK2 and SIRK2 are similar and the advan-
tages of SIRK2 for stiff equations cannot be exploited,
which is shown by a comparison of the computational
costs in Table 2. Using ERK2 in the whole domain is
about three times faster than using SIRK2.
As just mentioned, the fact that significant deforma-
tions take part in small regions of the computational
domain only has a considerable effect on the efficiency Figure 9. Sheet pile wall example: horizontal displacement
at the end of the calculation; deformations are scaled by a
of the integrators. Figure 7 displays the relative number factor 5.
of conducted times steps, taken by Umat. For instance,
in about 10% of the Umat calls SIRK2 takes between
31 and 40 time steps. To exploit the advantages of section area A = 1.2 × 10−2 m2 , moment of iner-
SIRK2, ERK2 should require considerably more time tia I = 3.84 × 10−5 m4 , Young’s modulus E = 2.1 ×
steps. This, however, is not the case for the consid- 108 kN/m2 .
ered example. Nevertheless, we have shown a good After an initial static step to impose the geostatic
performance of 2SIRK2 for our constitutive equations stress state, the second step is the excavation with fixed
in single element tests in (Fellin et al. 2009). bottom of the pit. The third load step is releasing the
bottom of the pit. For the second and third step the
initial and maximum load increments are t = 0.1, i.e.
3.2 Sheet pile wall example
a minimum of 10 load increments per step is required.
As a further typical geotechnical application we The automatic load incrementation strategy ofAbaqus
choose a 3.0 meter deep excavation with a 4.5 meter uses 10 load increments per load step for both methods
cantilever sheet pile wall. The geometry of the model only, i.e. one increment plus two times 10 increments in
is shown in Figure 8. The soil continuum is modelled total, see Table 3. Although the calculation with ERK2
with plane strain linear triangular finite elements. The requires a little bit more iterations, its computational
left and right boundary are fixed in horizontal direc- effort in terms of CPU time is considerably smaller
tion, the bottom boundary in vertical direction. The than the calculation with SIRK2S.
initial void ratio is set to e = 0.735 and the unit weight Figure 9 is a filled contour plot for the horizon-
to γ = 14.99 kN/m3 . The initial stress state is geostatic tal displacement u1 at the end of the calculation.
with the coefficient of earth pressure K0 = 0.5. The resulting horizontal displacement at the top of
The sheet pile wall is modelled with 2-node the wall is u1 = 3.819 × 10−2 m. Significant displace-
plane linear beam elements with the properties: cross ments occur in the area around the sheet pile wall,
19
large discretisation areas around the explored problem,
i.e. the regions where considerable deformations take
part are comparatively small. Due to this fact, in most
of the elements only few time steps with an arbitrary
adaptive integration method have to be conducted. A
time step of an implicit or semi-implicit method con-
sumes more computing time than a time step of an
explicit integrator. As a consequence, implicit or semi-
implicit methods can exploit their advantages only
in regions where an explicit method needs far more
time steps. Such regions are rather small in typical
geotechnical problems as exemplified here, and adap-
tive explicit methods turn out to be the superior choice
for integrating hypoplasticity with intergranular strain
in geotechnical applications. Switching from ERK2
to SIRK2 in regions with many explicit time steps is
worth to think about. However, as these regions are typ-
ically small and any switch algorithm will take some
extra time, the effect on the overall performance is
assumed to be small.
REFERENCES
Fellin, W., M. Mittendorfer, and A. Ostermann (2009). Adap-
tive integration of constitutive rate equations. Computers
Figure 10. Sheet pile wall example: total number of time and Geotechnics 36, 698–708.
steps in constitutive time integration. Fellin, W. and A. Ostermann (2002). Consistent tangent oper-
ators for constitutive rate equations. International Journal
whereas the displacements in the rest of the structure for Numerical and Analytical Methods in Geomechan-
are rather small. ics 26, 1213–1233.
Hairer, E., S. Nørsett, and G. Wanner (1993). Solving Ordi-
Figure 10 shows the total number time steps in each nary Differential Equations I. Nonstiff Problems (2nd ed.).
element, compare Section 3.1. The minimum number Berlin: Springer.
of required steps is nearly equal for both methods. In Hairer, E. and G. Wanner (1996). Solving Ordinary Dif-
the elements around the sheet pile wall more steps ferential Equations II: Stiff and Differential-Algebraic
have to be conducted. The explicit method needs at Problems. Berlin: Springer.
most twice as much steps than the semi-implicit one. Hügel, H. (1995). Prognose von Bodenverformungen, Vol-
However, there are only few elements where ERK2 ume 136 of Veröffentlichung des Institutes für Boden-
requires significantly more steps than SIRK2. Thus, mechanik und Felsmechanik. Universität Fridericiana in
the explicit method is more efficient for this test. Karlsruhe.
Kolymbas, D. (1985). A generalized hypoelastic constitutive
It is worth to note that hypoplasticity cannot handle law. In Proc. XI Int. Conf. Soil Mechanics and Foundation
stress states with tr T > 0. Such stress states could be Engineering, San Francisco, Volume 5, Rotterdam, pp.
predicted if the integrator chooses too large time steps 2626. Balkema.
in regions where the stresses are near to zero, which Mittendorfer, M. (2010). A modular finite element setting for
is the case directly below the ground surfaces. Our nonlinear constitutive models: design and implementa-
adaptive integrator rejects such steps into the unde- tion. Ph. D. thesis, University of Innsbruck.
fined area automatically, without extra checks for the Niemunis, A. and I. Herle (1997). Hypoplastic model for
admissibility of the solution. cohesionsless soils with elastic strain range. Mechanics
of Cohesive-frictional Materials 2, 279–299.
4 CONCLUSION
For sophisticated geotechnical problems one chooses

usually a finite element approach. Generally, one has
20
An anisotropic bubble model for soft clays
Nallathamby Sivasithamparam
Department of Civil Engineering, University of Strathclyde, UK
Plaxis B.V, Delft, The Netherlands
Daniela Kamrat-Pietraszewska & Minna Karstunen

Department of Civil Engineering, University of Strathclyde, UK
ABSTRACT: This paper describes the principles behind a new anisotropic bubble model for natural soils. The
model is a hierarchical extension of the anisotropic S-CLAY1 model. The kinematic yield surface of S-CLAY1
model is treated as a bounding surface and a bubble surface is introduced within the bounding surface. The
bubble surface is similar in shape to the S-CLAY1 yield surface, and assumes an isotropic elastic behaviour
and an associated flow rule. A translation rule of the bubble is used to control the movement of the bubble. The
implementation of the model is first verified by simulating slow cyclic loading with constant deviator stress
on Kaolin clay, and secondly, simulations of undrained triaxial shear tests (in compression and extension) are
made to highlight the effect of evolution of anisotropy, and finally, simulations of high number of loading cycles
performed to examine ratcheting feature of the model.
1 INTRODUCTION for Bubble S-CLAY1) and its isotropic version, which

is very similar to the Al-Tabbaa model, is referred
During past few decades several modifications have as BMCC. With the introduction of the bubble, the
been proposed to enhance elasto-plastic models devel- model allows the simulation of important features of
oped within the framework of kinematic hardening soil behaviour not realized by the S-CLAY1 model,
plasticity. One of the most successful approaches is to such as non-linearity and plasticity from early stages
introduce one or two kinematic surfaces within a con- of loading and hysteretic behaviour during cyclic load-
ventionally defined yield surface (Mroz et al. 1978, ing. BSCLAY model would be ideal for simulating the
1979). Models of this type are often termed kinematic behaviour of overconsolidated soils and/or the cyclic
hardening “bubble” models (Al-Tabbaa et al. 1987, response of soils.
1989). This paper presents a new constitutive model
that is capable of representing anisotropic and cyclic
behaviour of clay. The model implementation is first 2 DESCRIPTION OF THE BSCLAY MODEL
partially verified by comparisons with the Al-Tabbaa
(1987) model simulations for Kaolin clay considering The bounding surface in the general space of model is
isotropic material. Secondly, simulations of undrained formulated based on S-CLAY1 yield surface, which in
triaxial shear tests in compression and extension high- the simplified case of triaxial space considering cross-
light the effect of evolution of anisotropy on the anisotropic sample simplifies as follows:
predicted soil response.
The proposed constitutive model is developed
within the framework of the critical state theory and
bounding surface plasticity. The model is an extension where M is the value of the stress ratio η = q/p (devi-
of the S-CLAY1 model (Wheeler et al. 1999, 2003). ator stress over mean effective stress) at critical state,
The kinematic yield surface of S-CLAY1 is treated pm defines the size of the bounding surface (see Fig-
as bounding surface, and a bubble surface (kinematic ure 1) and α defines the orientation of the bounding
yield surface) is introduced within bounding surface surface analogously to the S-CLAY1 formulation. For
to enclose a truly elastic region. The anisotropic S- the finite element implementation of the model as a
CLAY1 bounding surface can describe the effect of user-defined model in PLAXIS, the equations have
initial anisotropy caused by one-dimensional deposi- been reformulated in terms of deviatoric stress vec-
tion and K0 -consolidation process, and the subsequent tor, a deviatoric fabric tensor (in vector form) and
evolution of anisotropy due to large strains is described mean effective stress, as explained in detail by Wheeler
by a kinematic hardening law of the S-CLAY1 model. et al. (2003). The bubble (kinematic yield) surface,
The new model is called the BSCLAY model (in short based on the ideas by Al-Tabbaa (1987), encloses the
21
the size of the bounding surface is controlled by the
change of plastic strain as follows:
where λ and κ are the slopes of the normal compression

line and swelling line in the e-lnp space and e is the
void ratio. When the bubble touches the bounding sur-
face the BSCLAY model becomes S-CLAY1 model,
and the bubble is dragged with the bounding surface.
Hardening rules and other details for the mathematical
formulation of the S-CLAY1 model can be found in
Wheeler et al. (1999, 2003).
2.2.2 Translation rules of the bubble surface

The translation rules are formulated based on transla-
tion rules proposed byAl-Tabbaa (1987). Two different
translations rules are needed, one for the case when the
bubble surface moves within the bounding surface, and
Figure 1. Schematic illustration of BSCLAY model. one for the case when the two surfaces are in contact.
The translation rule of the bubble surface is formu-
lated in a manner that guarantees that the two surfaces
truly elastic region, and it has the same shape as the (bubble & the bounding surface) can come in contact
bounding surface, but is smaller in size (see Figure 1). at a common normal but never intersect, similarly to
Although in the simplified form shown in this paper, Al-Tabbaa (1987). The centre of the kinematic yield
M is assumed the same in compression and extension, surface moves always along a vector, β, which joins the
in the FE implementation of the model, Lode angle current stress state, C, to its conjugate point, D, on the
dependency has been taken into account. bounding surface, see Figure 1. When bubble moves
within the bounding surface no rotation is allowed.
The translation rule is divided into two components.
2.1 Flow rule One is associated with the change in size of the bub-
ble surface due to isotropic expansion or contraction
Experimental evidence suggests that the assumption of the bounding surface, the other is associated with
of an associated flow rule is a reasonable approxima- the movement of the bubble surface along the vector
tion of natural clays when combined with an inclined β. When the two surfaces are in contact at the cur-
yield curve (Wheeler et al. 1999, 2003, Karstunen et al. rent stress state, the vector β is equal to zero and the
2005, 2008). Therefore, the plastic strain increment translation rule reduces to:
vector is assumed to be normal to the kinematic yield
surface at the current stress state. Consequently, the
flow rule of the model is associated and the plastic
potential is given by Equation 1.
2.2.3 Hardening modulus

2.2 Hardening rules The hardening modulus is defined in such a way that
The evolution of the bounding surface is described by when the two surfaces are in touch, and the yielding
anisotropic hardening; this means the surface rotates is continuous, the model predicts the same behaviour
as a function of large plastic strains (corresponding to as the S-CLAY1 model. It is initially formulated for
stress states at the bounding surface). The evolution of special case when two surfaces are in contact, and then
the bubble surface is described by a combination of modified for the general case when two surfaces are
kinematic and isotropic hardening, in which the sur- not in contact and the stress state is within the bounding
face translates in the stress space following the current surface.
stress point and changes size simultaneously. When the For general case when the two surfaces are not
bubble is not in touch with the bounding surface, it is in contact the hardening function h0 is given by the
not allowed to rotate. This combination of kinematic following equation:
and isotropic hardening forms the translation rule of
the kinematic yield surface.
2.2.1 Isotropic hardening of the bounding surface

The model adopts the volumetric hardening rule of In order to calculate the plastic strains whenever
the Modified Cam Clay (MCC) model. The change in they occur, whether or not the bubble surface and
22
the bounding surface are in contact, Al-Tabbaa (1987) Table 1. Model parameters for Kaolin clay.
replaced h0 with more general expression:
λ κ ν e0 M
0.187 0.017 0.3 1.2 0.9

where h0 is given by Equation 4 and H is a scalar α0 β µ R ψ
quantity which is function of the stress state.Al-Tabbaa 0.0 0.37 60 0.2 1.5
(1987) assumed, after Hashiguchi (1985), that H is a
function of a measure of the proximity of the bubble
surface to the bounding surface. The measure of the
proximity used in this model is the scalar product of to zero, the model simplifies to the isotropic MCC
the vector β and the vector n normal to the bubble model.
surface at the current stress state (see Fig. 1), divided
by the measure of the size of the bubble surface.
4 COMPARISON OF RESULTS
3 REQUIRED MODEL PARAMETERS The implementation of proposed model was first ver-
ified in a case of isotropic soil by comparing against
The proposed formulation of the model in general the results of Al-Tabbaa (1987) model predictions for
stress space requires values for 8 soil constants and slow cyclic triaxial test at constant deviator stress q.
3 state variables. These are: This simulation was initially started from normally
Soil constants: consolidated state corresponding to one-dimensional
loading, and the deviator stress q was kept constant
κ Initial slope of swelling/recompression line in e-
when cyclic loading cycles (unloading/reloading) were
lnp -space (see Al-Tabbaa 1987)
applied by changing p’. Secondly, the proposed model
ν Poisson’s ratio
was used to simulate soil behaviour under undrained
λ Slope of post yield compression line in e−lnp -
triaxial shearing following isotropic and anisotropic
space
triaxial consolidation. Table 1 summarizes the model
M Stress ratio at critical state (in triaxial compres-
parameters which obtained from Al-Tabbaa (1987),
sion)
and the additional soil constant and state variables were
µ Absolute effectiveness of rotational hardening
determined based on the suggestions by Wheeler et al.
β Relative effectiveness of rotational hardening (cal-
(1999, 2003) for β and Karstunen et al. 2005, 2008 and
culated based on M , see Wheeler et al. 2003)
Zentar et al. (2002) for µ. Given all tests by Al-Tabbaa
R Ratio of the size of the bubble surface to that of
were done for reconstituted Kaolin, the initial value
the bounding surface
for anisotropy (α0 ) has been assumed zero.
Exponent in the hardening function H (see Al-
Figure 2 compares the simulations of the isotropic
Tabbaa 1987)
version of the proposed model (BMCC) with the sim-
State variables: ulations of Al-Tabbaa (1987). The initial values of
p and q are 300 kPa and 80 kPa respectively, and
e0 Initial void ratio
q is kept constant while cyclic changes of p are
pm Initial size of the bounding surface (calculated
applied. As mentioned above, initial anisotropy has
based on vertical preconsolidation stress)
been switched off (α0 = 0) and additionally, the evolu-
α0 Initial inclination of the yield surface (calculated
tion of anisotropy was switched off by setting µ equal
based on M, see Wheeler et al. 2003)
to zero. In reality, anisotropy would have been cre-
The soil constants of the BSCLAY model include ated though the initial K0 consolidation, resulting in a
four parameters from the MCC model (κ, λ, M and theoretical value of α0 = 0.35. The mach between the
Poisson’s ratio ν ) that can be determined from conven- two model predictions is overall very good. Although
tional laboratory tests. Two additional parameters (R, BMCC is very similar to the Al-Tabbaa (1987) model,
) are required for introduction of the bubble surface she used the modified compression and swelling
into the S-CLAY1 model. Al-Tabbaa (1987) explains indices instead of λ and κ, and hence small differ-
how these six model parameters can be obtained from ences would be expected. These results suggest that
simple standard tests or multi-stage test using the triax- the proposed model has been implemented correctly.
ial apparatus. Two additional soil constants (µ and β) A corresponding simulation with the anisotropic ver-
and additional state variable (α0 ) govern the evolution sion of the model (BSCLAY), which has not been
of anisotropy and the initial anisotropy, respectively. included in the paper, suggests that for this type of
Wheeler et al. (1999, 2003) discussed the determina- cycling loading and amplitude, anisotropy does not
tion of these three parameters in detail and generally have significant influence in the volumetric response,
no non-standard tests are needed to get reasonable esti- but nevertheless it has major impact on the predicted
mates for these values. The model is hierarchical, so deviatoric straining. Just like the isotropic BMCC ver-
it is possible to reduce the model to the S-CLAY1 sion of the model, the BSCLAY model seems to be
model, by setting R equal to one. Furthermore, if initial able to reproduce well the soil response under slow
anisotropy is switched off, by setting α0 and µ equal cyclic loading.
23
Figure 2. Slow cyclic isotropic constant q triaxial simulation: a) After Al-Tabbaa (1987) model predictions b) BSCLAY
Model simulation.
Figure 3. Simulation of undrained stress path a) after an isotropic stress history b) after a one dimensional stress history.
In Figure 3, thick solid lines represent the pre- initially isotropic and during the initial isotropic con-
dictions of the anisotropic BSCLAY model and the solidation and unloading it stays isotropic according to
dashed lines represent the equivalent results by the both models, as even BCSLAY model predicts soil to
isotropic BMCC model. In both cases the soil is stay isotropic under isotropic loading. The isotropic
24
Figure 4. High number (100 cycles) of cyclic simulation of BSCLAY model a) q/p’ versus εs and b) q/p’ versus εv .
loading is followed by isotropic unloading, corre- A high number of cyclic loading constant q triaxial
sponding to overconsolidation ratios (OCR) of 1, 1.3, simulation was performed with BSCLAY model after
2, 4 and 8. Due to the initial isotropic compression, as a one dimensional stress history. The simulation was
seen in Figure 3(a), both models predict similar stress initially one-dimensionally compressed to σv = 20 kPa
paths for undrained shearing after the isotropic com- then cyclically (100 cycles) loaded between stress
pression at early stages of the simulations, but once ratios of η = 0.45 and η = 0.23, see Figure 4. Shear
the bounding surface is reached, the prediction devi- strain continues to accumulate with increasing num-
ate. The anisotropic version of the model (BSCLAY) ber of cycles. The ratcheting feature of the model
predicts lower excess pore pressures and higher values may over-predict the shear strain after large number
of deviator stresses at failure than BMCC. of cycles. To avoid ratcheting feature of the model, the
Differences between the two model predictions are size of the bubble, R, could be made to function of
very striking in Figure 3b relating to the simulations number of cycles so that the soil will behave elasti-
of anisotropically consolidated undrained shearing in cally after large number of cycles. However, this will
compression and extension. Again, the soil is assumed require further investigations.
initially isotropic, but during the initial K0 consolida-
tion anisotropy evolves in the case of BSCLAY model,
5 CONCLUSIONS
resulting in an α-value of 0.35 at the start of undrained
shearing. Due to the associated flow rule, K0 - load-
A new constitutive model, BSCLAY, which is a hier-
ing results in different predicted stress paths, both for
archical extension of the S-CLAY1 model, has been
loading and unloading. BMCC gives just like the MCC
developed to simulate cyclic loading of anisotropic
model a very poor K0 prediction, and consequently
clays. The model is based on the principles of bound-
in the cases of high OCR the shearing starts close
ing surface plasticity. A bubble surface is introduced
to failure. Overall, during compression the BSCLAY
within S-CLAY1 model to enhance the performance
model predicts higher undrained strength than BMCC,
of the model to describe soil behaviour in over-
and the predicted undrained strength in extension is
consolidated region and under cyclic loading. The
notably lower than in compression. Once Lode angle
comparisons of the model predictions with the Al-
dependency is included, the difference is even more
Tabbaa (1987) model simulations of Kaolin clay under
significant than in the case of Drucker-Prager assumed
different stress paths, considering slow cycling load-
in this paper. In contrast the isotropic BMCC model
ing and shearing under compression and extension,
predicts almost the same value of undrained shear
revealed the predictive capability of the proposed
strength in compression and extension.
model. Ratcheting feature of the model is also verified.
In order to have a unique critical state, the rotational
hardening law of the BSCLAY model is formulated
in such as way that at reaching critical state the ACKNOWLEDGEMENTS
bounding surface keeps rotating until a unique orien-
tation is reached (see Wheeler et al. 2003 for details). The research was carried out as part of a “GEO-
Because of this feature, the results for triaxial exten- INSTALL” (Modelling Installation Effects in Geotech-
sion have strange looking curvature when approaching nical Engineering), supported by the European
critical state. This may require some modification Community through the programme “Marie Curie
when considering finite element applications, such as Industry-Academia Partnerships and Pathways” (Con-
excavations. tract No PIAP-GA-2009-230638).
25
REFERENCES Mróz, Z., Norris, V.A. & Zienkiewicz, O. C. 1979. Applica-
tion of an anisotropic hardening model in the analysis of
Al-Tabbaa, A. 1987. Permeability and stress-strain response elasto-plastic deformation of soils. Géotechnique 29, No.
of Speswhite kaolin. PhD dissertation. University of 1, 1–37.
Cambridge. Wheeler, S.J., Karstunen, M. & Näätänen, A. 1999.
Al Tabbaa & Wood, D.M. 1989 An experimentally based bub- Anisotropic hardening model for normally consolidated
ble model for clay. In: Proc. 3rd Int. Conf. on Numerical soft clay. In G.N. Pande, S. Pietruszczak & H.F. Schweiger
Models in Geomechanics. Niagara Falls, pp. 91–99. (ed.), Proc. 7th Int. Symp. on Numerical Models in
Hashiguchi, K. 1985. Two- and three-surface models of Geomechanics (NUMOG VII), Graz : 33–40. A.A.
plasticity. Proceedings of 5th International Confer- Balkema.
ence on Numerical Methods in Geomechanics, Nagoya, Wheeler, S.J., Näätänen, A., Karstunen, M. & Lojander, M.
pp. 285–292. 2003. An anisotropic elastoplastic model for soft clays.
Karstunen, M & Koskinen, M. 2008. Plastic anisotropy of Canadian Geotechnical Journal 40: 403–418.
soft reconstituted clays. Canadian Geotechnical Journal Zentar, R., Karstunen, M. & Wheeler, S.J. 2002. Influence
45: 314–328. of anisotropy and destructuration on undrained shearing
Karstunen, M.; Krenn, H.; Wheeler, S.J.; Koskinen, M., Zen- of natural clays. In P. Mestat (ed.), Proc. 5th European
tar, R. 2005. The effect of anisotropy and destructuration Conf. on Numerical Methods in Geotechnical Engineering
on the behaviour of Murro test embankment. International (NUMGE 2002), Paris: 21–26. Presses de l’ENPC.
Journal of Geomechanics (ASCE); 5(2): p. 87–97. Zentar, R., Karstunen, M., Wiltafsky, C., Schweiger, H.F. &
Mróz, Z., Norris, V.A. & Zienkiewicz, O. C. 1978. An Koskinen, M. 2002. Comparison of two approaches for
anisotropic hardening model soils and its application modelling anisotropy of soft clays. In G.N. Pande &
to cyclic loading. International Journal for Numer- S. Pietruszczak (ed.), Proc. 8th Int. Symp. on Numer-
ical and Analytical methods in Geomechanics. 2, ical Models in Geomechanics (NUMOG VIII), Rome:
203–221. 115–121. A.A. Balkema.
26
An anisotropic model for structured soils
G. Belokas & M. Kavvadas

National Technical University of Athens, Greece
ABSTRACT: The paper presents the mathematical formulation of the recently developed constitutive Model
for Structured Soils – 2, which is a two surface anisotropic bounding surface plasticity model. This model is able
to reproduce the engineering effects of the structure inducing mechanisms, including the effect of anisotropy
by incorporating: a) distorted and rotated ellipsoids for the Structure Strength Envelope (bounding surface) and
the Plastic Yield Envelope (elastic region – inner surface) to describe bond and stress induced anisotropy, b)
the Intrinsic Strength Envelope as a reference locus that delimits all possible unbonded states, representing a
lower bound of the bounding surface, c) the Intrinsic Compressibility Framework that describes all structureless
states, d) a damage-type mechanism to model bond degradation and e) a non-associated flow rule depending
on structure. The proposed model is modular, its features can be activated simultaneously or selectively, and the
3-D tensorial formulation facilitates direct implementation in finite elements codes.
1 INTRODUCTION the tensorial mathematical formulation of the model

allows its direct implementation in finite elements
Most natural soils are structured and the main structure codes and 3D analyses.
inducing mechanisms are stress history (e.g. preload- The tensorial formulation includes the stress ten-
ing, densification) and bonding (e.g. cementation, sor σ = s + σI (s: deviatoric stress, σ: mean effective
aging, thixotropy). These mechanisms can have a dom- stress) and the strain tensor ε = e + εI (s: deviatoric
inant role on natural soil behaviour affecting stiffness, strain, σ: volumetric strain), with I the unit tensor. In
dilatancy and strength, as well as their anisotropic the following, all stress is effective. Dots over a symbol
characteristics. Bonding results in components of stiff- denote incremental quantities.
ness, dilatancy, strength and anisotropy that cannot
be accounted solely from stress history. Therefore, it
is essential that modern constitutive modelling takes 2 STRUCTURELESS SOIL BEHAVIOUR
into account these mechanisms and their anisotropic
characteristics. This improves greatly the predictive Following Roscoe et al (1963), Lewin & Burland
capability of a constitutive model. (1970) and Leroueil & Vaughan (1990), structure-
The proposed Model for Structured Soils – 2 less materials can occur after thorough remoulding
(MSS-2) simulates the engineering effects of these (to eliminate all memory of original structure) and
structure inducing mechanisms. More specifically it subsequent radial consolidation, i.e., along a path of
advances present elastoplastic constitutive modelling constant stress ratio ησ = (s:s)0.5 / σ. This is supported
practice of structured soils by incorporating: a) dis- by various experimental evidence (see Belokas & Kav-
torted ellipsoids for the Structure Strength Envelope vadas, under review). It has also been assumed that
(bounding surface) and the Plastic Yield Envelope the behaviour of structureless materials is controlled
(elastic region) to describe bond and stress induced by current effective stress state (σ, s) and specific vol-
anisotropy, b) the Intrinsic Strength Envelope as a ume (v = 1 + e) only (e.g. Leroueil & Vaughan 1990),
reference locus that delimits all possible unbonded since remoulding has erased any pre-existing bonding
states, representing a lower bound of the bounding sur- and stress history effects.
face, c) the Intrinsic Compressibility Framework that During a radial consolidation, structureless materi-
describes all structureless states, d) a new damage-type als move along Intrinsic Compression Curves (ICC η ),
mechanism to model bond degradation and e) a plastic which are herein assumed to be parallel straight lines
dilatancy dependence on bonding and anisotropy. in the lnv–lnσ plane (Fig. 1, Equation 1) with slope
MSS-2 has been developed starting from the described by a Modified Cam-Clay compressibility
Kavvadas & Belokas (2001) model and the original parameter (ρc ). Equation 2 correlates the stress path
MSS model (Kavvadas & Amorosi, 2000). It is formu- orientation with the ICC η and is a variance of the
lated based on a hierarchical approach, which results equation proposed by Belokas & Kavvadas (in press).
in a modular and versatile model with features than Equations 1 and 2 are used for the mathematical
can be used simultaneously or selectively. Moreover, description of the structureless soil states. The lower
27
Figure 1. Structureless states and Intrinsic Compressibility
Framework (from Belokas & Kavvadas, 2010).
bound of all possible ICC η is the Critical State Curve

(CSC), which in Critical State Soil Mechanics sep-
arates contractant from dilatant behaviour, and the
upper bound is the Isotropic Intrinsic Comrpession
Curve.
Figure 2. Influence of structure in plastic dilatancy (based
on Belokas 2009 and Belokas & Kavvadas 2010b).
where σ PTC is given by Equation 4 if ICC η is described

by Equation 1.
3 STRUCTURE AND DILATANCY
The ICC η curves delimit all possible states under radial

compression. In Critical State Soil Mechanics, CSC
serves the purpose of a Phase Transition Curve (PTC). where σ o ≥ σ ∗o (see Fig. 2) is a measure of the available
To the left of PTC the behaviour is dilatant, while to structure and defines the Structure Strength Curve in
the right of PTC the behaviour is contractant (Fig. 2a). v – σ plane.
For structured soils, PTC moves to the right, therefore As loading proceeds, it generally results in a loss of
enlarging the dilatant domain (Fig. 2b). The larger the structure and bonding. Moreover, PTC shifts towards
structure. the larger the dilatant domain and the PTC the CSC and σ o tends to σ ∗o .
moves further to the right.
The distinction between dilatant and contractant
behaviour is controlled by the phase parameter, ψσ , 4 FORMULATION OF MSS-2 MODEL
which is defined by Equation 3 and is a variance of the
state parameter proposed by Been & Jefferies (1985). The formulation of the MSS-2 model is based on rate-
When the current state lies: a) to the left of PTC it independent incremental bounding-surface elasto-
is ψσ < 0 (i.e. dilatant behaviour), b) upon PTC it is plasticity. It incorporates the Critical State concepts
ψσ = 0 (e.g. critical state) and c) to the right of PTC it and the behavioural framework for structured soils
is ψσ > 0 (i.e. contractant behaviour). (e.g. Belokas, 2008). MSS-2 encompasses the impor-
tant aspects of soil behaviour described in paragraphs
2 and 3. Moreover, it is able to model structure and
28
(Kavvadas 1983) and is the bounding surface:
Tensor σ K = σ K I + sK is the centre K of SSE and

the half-axes of the ellipsoid are equal to (α) along the
isotropic axis and (cα) along the deviatoric axes.
The above mathematical formulation can use dif-
ferent values of the c-parameter along the various
deviatoric stress components to model shear strength
anisotropy. The size (α) of the SSE is controlled by
the magnitude of structure, while the orientation vec-
tor (bK ≡ sK /σ K ) of the centre of the SSE is a direct
Figure 3. The characteristic surfaces scheme used in the measure of structure-induced anisotropy.
MSS-2 model. The PYE is, for mathematical simplicity, similar in
shape and has axes parallel to the SSE scaled down by
a factor ξ (<<1):
anisotropy evolution and degradation, as well as their

effect on soil strength, dilatancy and stiffness. It has
been developed in a modular way, i.e. its various
individual features can be activated or de-activated
depending on the amount of information or the type
of the problem to be solved.
The model presented herein has been based on the
Kavvadas & Belokas (2001) model and its complete Tensor σ L = σ L I + sL is the centre L of PYE. Since
formulation is presented in Belokas (2008). most soils behave elastically in a very limited strain
domain (order of strain ε ≈ 0.001% ÷ 0.2%), the
model typically uses a very small PYE (ξ ≈ 0.001)
which can also provide realistic modelling of cyclic
4.1 Characteristic Surfaces
loading.
MSS-2 employs three characteristic surfaces (Fig. 3): The ISE corresponds to an equivalent structureless
the Structure Strength Envelope (SSE), the Plastic state, which has the same specific volume (v) and con-
Yield Envelope (PYE) and the Intrinsic Strength Enve- solidation stress ratio (ησ ) as the natural material. This
lope (ISE), each one serving a different purpose. SSE envelope is also a rotated distorted ellipsoid centred at
delimits all possible states of a structured soil (it is the point K ∗ having coordinates: σ ∗K = σK∗ I + s∗K :
bounding surface), PYE bounds the elastic states (it
is the yield surface), while ISE is a reference surface
that corresponds to an equivalent structureless state.
These surfaces are formulated in a tensorial space
consisting of the isotropic axis (σ) and the deviatoric
hyperaxis (s).
The behaviour within the PYE has been assumed to The long axis of the ISE, oriented along the vec-
be hyperelastic. States on PYE represent the onset of tor b∗K ≡ s∗K /σK∗ , describes inherent anisotropy (con-
yielding, which correspond to limited structure degra- trolled by the consolidation stress ratio ησ ) while the
dation. For loadings directing outwards of PYE, the corresponding orientation of the SSE, bK ≡ sK /σ K ,
PYE moves together with current state (kinematic describes structure-induced anisotropy (which can be
hardening) until it comes in touch with SSE (the different). The orientation of the ISE can vary accord-
structure yield of Kavvadas 1998). This is the onset ing to a kinematic hardening rule, modelling the
of appreciable structure degradation and evolution of variation of inherent anisotropy with consolidation.
structure – induced anisotropy. The ISE represents a The size (α∗ ) of the ISE is controlled by classical Cam-
lower bound of the SSE, when all effects of structure Clay type isotropic hardening depending on the current
are eliminated (typically by intense straining). Thus, specific volume (v), the intrinsic compressibility (ρc )
in structureless soils, the SSE and ISE coincide and and the consolidation stress ratio (ησ ) expressed via
the model reduces to a Cam-Clay type twin-surface the parameter Nη (Equation 2 and Fig. 1):
model with a rotated bounding surface (SSE ≡ ISE)
and an internal bubble (PYE).
The geometrical representation of the SSE in the
stress space σ ≡ (σ, s) is a rotated distorted ellipsoid
29
The difference in size between the SSE and ISE i.e., K moves along a radial path passing through the
(α − α∗ ) is a direct measure of the magnitude of origin. As the ratio bK ≡ sK /σK remains constant, pri-
structure. mary anisotropy does not change. The SSE reduces
to the Modified Cam-Clay yield surface if K lies on
4.2 Hardening rules the isotropic axis (σ K = αI), e.g. during an isotropic
consolidation path.
The isotropic and kinematic hardening rules control
For material states on the SSE:
the evolution of the characteristic surfaces during plas-
tic straining. Upon plastic straining, current stress state
(σ) is always on the PYE.
4.2.1 Isotropic hardening rule

MSS-2 model employs a structure degradation The second term in the above expression causes K to
isotropic hardening: deviate from the radial path during straining, altering
primary anisotropy (bK ), while the constant ψ controls
the rate of its evolution.
Motion of the centre L of the PYE
For material states on the SSE, surfaces PYE and
SSE remain in contact (at the current stress state σ)
and the position of L is dictated by the position (σ K )
of K:
For material states inside the SSE, the motion of

point L is such that the current state on the PYE (point
M on PYE in Fig. 3) moves towards a conjugate point
where Equations 11a and 11b give respectively the M on the SSE (normal vectors at these points are paral-
p
plastic volumetric and deviatoric strain increment, (εv , lel). The geometric similarity of PYE and SSE defines
p
εq ) are the accumulated plastic volumetric and devia- the direction vector β ≡ (σ − σ L )/ξ − (σ − σ L ). There-
toric strains, (ρc , ρs ) are the intrinsic compressibility fore, the translation of the centre L is given by the
parameters during virgin compression and rebound, formula:
(ζ v , ηv ) are volumetric structure degradation parame-
ters and (θ q , ζ q , ηq ) are deviatoric structure degradation
parameters. Parameter θ q is used to increase the rate of
structure degradation by shearing, since shear-induced
This rule ensures that the characteristic surfaces do
structure degradation is usually dominant (compared
not intersect even for finite increments of the material
to volumetric structure degradation).
state (Kavvadas & Amorosi, 2000). The factor µ̇ is
In Equation 9 the component (α − α∗ ) represents
determined from the “consistency condition”, i.e. a
the magnitude of structure and α∗ is given by Equa-
requirement that during plastic deformation the stress
tion 8. Equation 9 ensures that the α − α∗ is constantly
decreasing down to the value of α − α∗ = 0, when all point remains on the PYE (f˙ = 0).
structure has been lost (i.e. SSE≡ISE), shifting to a
Cam – Clay type of isotropic hardening. Moreover, 4.3 Flow rule
by appropriately selecting the structure degradation
parameters, various rates of destructuring can be The plastic strain increment is determined by a non-
simulated, including collapse-type behaviour. associated incrementally linear plastic flow rule:
4.2.2 Kinematic hardening rules

The kinematic hardening rules describe the evolu-
tion of the structure-induced anisotropy during plastic
where the scalar ( ˙ f ) and the plastic potential ten-
straining by controlling the position of the centres
K and L of the SSE and PYE in the stress space. sor (Pf ) give the magnitude and direction of the
The MSS-2 model describes anisotropy by the primary plastic strain increment, (σ̇) is the corresponding effec-
(bK ≡ sK /σ K ) and secondary (bL ≡ sL /σ L ) anisotropy tive stress increment, (Hf ) is a “plastic modulus” as
tensors and uses the kinematic hardening rules of the described in a following section, and Qf ≡ ∂f /∂σ is
original MSS model (Kavvadas & Amorosi, 2000). the gradient of the PYE.
The plastic gradient (Pf ) has the following isotropic
Motion of the centre K of the SSE and deviatoric components:
For material states inside the SSE:
30
where ρs , A and B are hyper-elastic constants.
6 PLASTIC MODULUS
For material states on the SSE, the plastic modulus,

Hf , is determined from the “consistency condition”,
which ensures that the stress point remains on the SSE
during plastic loading and results in equations 19, 20
and 21:
Figure 4. Failure envelopes (FE) and phase parameter (ψσ ).
where λ1 is a positive constant, h(σ) is a conical fail-

ure envelope (FE) given by Equation 16 and shown in
Figure 4, and λ2 is a positive constant.
For material states on the PYE but inside the SSE,

the plastic modulus is determined by an interpolation
mapping rule, which ensures a smooth and continuous
transition of Hf as PYE approaches SSE:
where kc and ξ are scalar and tensorial constants
respectively, related to the critical slope (M ), kd is
a positive constant and kp = kc − kd ψσ is a variable
(similar to the parameter used by Wood et al 1994). where Hfo is the plastic modulus at point M " (where
Improving the Kavvadas & Belokas (2001) model, −−→
the MSS-2 flow rule employs the phase parameter, ψσ vector OM intersects the SSE in Fig. 3), δ is the nor-
(Fig. 2), to distinguish dilatant from contractant states malized length of MM (M is the current state), δo is
and introduces a dependence on primary anisotropy the value of the parameter δ upon yield initiation and
for deviatoric strains. γ = γ 1 + γ 2 (α/α ∗ − 1) is a variable with γ 1 and γ 2
Constants λ1 and λ2 control the direction of the positive constants.
plastic strain increment for states on the SSE. The func- Equation 22 interpolates between the values:
tion hc serves the purpose of the critical state line in Hfo = ∞ (upon initiation of yielding) and: Hfo = Hfo
terms of stresses, while the function hp allows for stress (when the stress state reaches the SSE). Point M
paths that intersect this critical state line during dilatant is calculated by the condition F(λσ; σ K , α) = 0 and
behaviour. parameter δ is defined by Equation 23 (see also
Kavvadas & Amorosi 2000). Parameter δo is reset to
the value of current δ each time yielding is re-initiated.
5 ELASTICITY
Thus, δ/δo = 1 upon initiation of yielding, δ/δo < 1 at
any later stage, and δ = 0 when the material state lies
A hyper-elastic formulation which results from an
on the SSE.
elastic strain potential V e has been incorporated:
A differentiation from previous model versions is that

variable γ includes the effect of structure on the plastic
hardening modulus.
31
ψσ , has been employed, which controls dilatancy –
contractancy and depends on the magnitude of struc-
ture, d) the deviatoric component of the plastic flow
depends on structure anisotropy and e) the plastic
hardening modulus depends on the magnitude of
bonding.
Compared to the original MSS model (Kavvadas &
Amorosi, 2000) the major advances include the incor-
poration of: a) rotated distorted ellipsoids for the
bounding and the yield surfaces, b) a different damage-
type mechanism to model structure degradation and
c) the Intrinsic Strength Envelope as a reference
Figure 5. Influence volumetric degradation parameters on envelope.
compressibility.
REFERENCES
Been K and Jefferies MG. 1985. A state parameter for sands.
Géotechnique. 35(2):99–112.
Belokas G. 2008 Modelling of the Mechanical Behaviour
of Structured and Anisotropic Soil Materials. Ph.D The-
sis. National Technical University of Athens. pp695 (in
Greek).
Belokas G and Kavvadas M. An intrinsic compressibility
framework for clayey soils. Geotechnical and Geological
Engineering, under review.
Kavvadas M. 1983. A constitutive model for clays based
on non-associated anisotropic elasto-plasticity. Proc. of
the 2nd Int. Conf. on Constitutive Laws for Engineering
Figure 6. Influence of bonding on undrained shear response. Materials, in Tucson. p. 263–270.
Kavvadas M. 1998. Hard Soils – Soft Rocks: Modelling the
soil behaviour – Selection of soil parameters, General
Report. Proc. 2nd Int. Symp. on the Geotechnics of Hard
7 EXAMPLE SIMULATIONS Soils – Soft Rock, in Napoli. p. 1441–1482.
Kavvadas M and Amorosi A. 2000. A constitutive model for
Figure 5 shows the structure degradation during radial structured soils. Géotechnique. 50(1): 263–273.
compression and Figure 6 shows structure degrada- Kavvadas MJ and Belokas G. 2001. An anisotropic elasto-
tion for various degrees of bonding (α/α∗ ) during and plastic constitutive model for natural soils. Proc. 10th Int.
undrained shear. Conf. on Computer Methods andAdvances in Geomechan-
ics (IACMAG), in Tucson, Arizona,. p. 335–340.
Leroueil S and Vaughan PR. 1990. The general and congru-
ent effects of structure in natural soils and weak rocks.
8 CONCLUSIONS Géotechnique. 40(3):467–488.
Lewin P.I. and Burland J.B. 1970. Stress-probe experiments
The formulation of an anisotropic bounding sur- on saturated normally consolidated clay. Géotechnique.
face plasticity constitutive Model for Structured Soils 20(1):38–56.
Roscoe K.H., Schofield A.N. and Thurairajah A. 1963.Yield-
(MSS-2) has been presented. It has been based on the ing of clays in states wetter than critical. Geotechnique.
Kavvadas & Belokas (2001) the major advances being 13(3):211–240.
the following: a) a hyperelastic formulation has been Wood DM, Belkheir K and Liu DF. 1994. Strain softening
employed, b) the Intrinsic Compression Curves are and state parameter for sand modelling. Technical Note.
linear in the lnv-lnσ plane, c) the phase parameter, Géotechnique. 44(2):335–339.
32
An examination of strain space versus stress space for the formulation of

elastoplastic constitutive models
K.C. Ellison & K. Soga

University of Cambridge, Cambridge, UK
B. Simpson
Arup Geotechnics, London, UK
ABSTRACT: Several elastoplastic soil models have been proposed over the years that are formulated in
strain space rather than stress space due to certain analytical and computational advantages. One such model,
BRICK (Simpson 1992), has been continuously utilized and developed for industrial applications within Arup
Geotechnics for more than two decades. This paper aims to describe the advantages and difficulties associated
with strain space modeling. In addition, it will show how recent advances in modeling the effects of stress
history, stiffness anisotropy, strength anisotropy and time-dependence in conventional stress space models can
be transferred to the BRICK model.
1 INTRODUCTION 2 BRIEF OVERVIEW OF THE BRICK MODEL
Conventional elastoplastic critical state based con- A simple plane-strain version of the BRICK model
stitutive models for soil behavior are formulated was described by Simpson (1992). The key features of
primarily in stress space with one or more yield this model will be reviewed briefly here.
surfaces defined in terms of principal stresses. Alter- In the original BRICK model, the current strain
natively, several constitutive models for geomaterials state was defined in a three-dimensional coor-
have been proposed over the years that are formu- dinate system with one volumetric strain axis,
lated entirely in strain space due to certain analytical εvol = (εx + εy )/2, and two shear strain axes, εx − εy
and computational advantages over the conventional (pure shear) and γxy (simple shear). In addition, sev-
models (e.g. Yoder & Iwan (1981) and Iwan & Chel- eral ‘bricks’ that represented a portion of the material
vakumar (1988)). While most of these have failed were attached to the strain point by ‘strings’of different
to gain traction outside of academic realms, the lengths in this strain space. When a string became taut,
BRICK model (Simpson 1992) has been continu- the corresponding brick would move toward the strain
ously utilized and developed for industrial applica- point as demonstrated by the sequence in Figure 1 and
tions within Arup Geotechnics for more than two plastic strain would develop along that path. The total
decades. plastic strain increment {ε̇p } was determined by sum-
Since the initial formulation of the BRICK model, ming the contributions of each brick while the elastic
many advances have been described in the literature strain increment was determined simply via the rela-
to capture specific features of soil behavior in stress tionship {ε̇e } = {ε̇} − {ε̇p }. Thus, when all the strings
space models; however, there is little in the litera- were loose, the soil behavior was fully elastic and when
ture to describe how these advances may be applied all the strings were taut and lined up behind the strain
in strain space. For example, recent constitutive mod- point in the orientation of the strain increment, the soil
els can explicitly consider the effects of stress history, behavior was perfectly plastic.
creep, rate-dependence, stiffness anisotropy, strength In addition to the BRICK effect described above, it
anisotropy and other factors in ways that the original was assumed that the soil’s capacity for elastic strains
single-surface models could not. During the ongoing increased as the mean effective stress increased. This
development of the BRICK model, several of these was achieved by discounting some of the volumetric
advances have been modified for use in strain space. plastic strains indicated by brick movements so that
This paper aims to describe the advantages and diffi- changes in volumetric strain followed the appropriate
culties of modeling in strain space as well as to outline normal consolidation or swelling line in εvol -ln s space
how recent advances can be applied to BRICK-type (when all the bricks were aligned) where s was the
models. mean in-plane effective stress. Increases in volumetric
33
Figure 1. Example sequence of strain point and brick movement during 1D consolidation (a→b) and undrained extension
(b→d).
strain were also associated with an increase in the by the accumulation of strain. Thus, strain space
capacity for shear strain. Therefore, shear plastic strain models may be better suited to explain the underly-
reductions were applied such that shear failure was ing mechanisms that govern the constitutive behavior
achieved at a constant stress ratio t/s = f ({Lb },{Rb }) of soil.
where t was the shear stress and {Lb } and {Rb } were the It has long been recognized that void ratio is a criti-
string lengths and material proportions, respectively, cal parameter to soil behavior. In stress space models,
for each brick ‘b’. the void ratio (or volumetric strain) is needed to iso-
The elastic stiffness in the BRICK model is late the critical state line. Moreover, plastic strains are
pressure-dependent, i.e. K e = s /ι where K e is the needed to compute the hardening of yield surfaces in
elastic bulk modulus and ι is a user-defined elastic stress space models and it seems more appropriate to
stiffness parameter. However, it is also presumed that derive these plastic strains directly from an increment
the strength and stiffness of the material will increase of strain rather than an increment of stress. Therefore,
with overconsolidation as measured by the distance strain space is simply a more consistent basis for a
of the current strain point from the normal consoli- constitutive model.
dation line in εvol -ln s space. Increased stiffness was As another example, consider the phenomenon of
introduced by decreasing the parameter ι from its user- stiffness anisotropy which arises from the preferred
specified value. This would also lead to an increase orientations of particles and particle contacts that make
in the failure stress ratio t/s which was modified by up a soil’s fabric. The accumulation of large strains
adjusting the string lengths {Lb }. might change this fabric and corresponding anisotropy
The three-dimensional (3D) BRICK model cur- even if the initial and final stress states are the same.
rently used by Arup is based on the same principles as Stress space models account for this by allowing the
the plane-strain version described by Simpson (1992) yield surface to expand, translate or rotate throughout
and summarized above.The major difference is that the the stress history; however, a more realistic descrip-
strain point and bricks are defined in a six-dimensional tion should examine strain history since this is a better
space comprised of one volumetric strain and five measure of the change in fabric.
shear strains. The details of this formulation have Furthermore, if a specimen is subjected to an abrupt
recently been described by Ellison (2009) and Clarke change in stress path, plastic strains will initially con-
(2009). tinue to develop in the direction of its recent strain
history (Atkinson et. al. 1990). If the initial stress path
is small, then its recent stress history can be ‘forgotten’
3 ADVANTAGES OF STRAIN SPACE
after a period of creep (Clayton & Heymann 2001).
However, more significant stress histories cannot be
3.1 Philosophical advantages
completely forgotten due to creep (Gasparre 2005).
As stated in the frontispiece of Professor John This observation is best described by an examination of
Burland’s PhD thesis (Burland 1967) and reiterated strain: small strains will result in a small change to the
during Brian Simpson’s Rankine Lecture (Simpson soil fabric that can be overcome by subsequent creep
1992): ‘Stress is a philosophical concept - deformation strains whereas larger strains may result in a signifi-
is the physical reality’. This quote encapsulates one of cant change to the soil fabric that cannot be overcome
the most compelling reasons that an examination of by creep. In contrast, the influence of stress changes
strain rather than stress might be more appropriate to on the soil fabric will be harder to gauge since this will
describe the evolution of soil behavior. While these vary with the current stress state.
two measures are inextricably linked, changes in soil It also makes intuitive sense to model both creep and
behavior are ultimately caused by micromechanical the related phenomenon of stress relaxation (whereby
changes in soil fabric that are reflected at the mesoscale stress decreases over time while strain remains
34
constant) in strain space since these phenomena are
thought to arise from the gradual rearrangement of par-
ticles due to bond failures at the molecular level. This
rearrangement would best be expressed directly by the
development of plastic strains rather than indirectly
through the propensity of a yield surface to reposition
itself in stress space.
3.2 Practical advantages

There are also some practical reasons to formulate con-
stitutive models in strain space. The primary reason
is that strain space models are naturally compatible
with the finite element method whereby an incre- Figure 2. Sample discretization of the S-shaped curve.
ment of strain is shuttled into the constitutive model
and an increment of stress must be produced, i.e. rate of stiffness degradation compared with conven-
{σ̇} = f({ε̇}, {k}) where {k} is comprised of the state tional flow rules. However, a finite number of bricks
variables (e.g. Chelvakumar & Iwan (1988)). On the must be used and, therefore, the stiffness curve will
other hand, stress space formulations must be rear- not necessarily be smooth. For practical purposes, 10
ranged to make strain the independent variable during bricks are normally used, but in principle the number
a numerical time step, often at the expense of computa- is not limited.
tional efficiency. This complication is apparent in the The greatest challenge in working with BRICK is
calculation of the scalar multiplier which is used in the determination of the in situ state variables. At the
conjunction with the plastic potential ∂P/∂σ to obtain start of a BRICK simulation, the bricks and strain
the plastic strain increment. point must be located along the volumetric strain axis
A second practical advantage of strain space model- or else anisotropic strength will be predicted. There-
ing is that it eliminates the need to make assumptions fore, a best estimate of the brick locations at the in
about the intersection of yield surfaces. Multi-surface situ state is usually obtained by modeling a portion of
formulations in stress space with bilinear constitutive geological history considered to have influenced the
laws require that no two surfaces intersect or else the current state. For London clay, an arbitrary reference
uniqueness of solution could be destroyed (Puzrin & state of p = 2 kPa is typically used where p is the
Houlsby 2001). However, Yoder (1981) demonstrated mean effective stress. This is inconvenient because the
that no such problem exists for yield surfaces in strain geological strain history is not always well-known and
space. Moreover, in multisurface stress space models, there is no guarantee that its simulation will result in
it can be difficult to determine the proportion of a the same horizontal pressures observed in the field.
strain increment that takes place before a new surface By comparison, conventional stress space models are
is encountered. This is usually remedied via a com- attractive because their simulations may begin at the in
putationally expensive iteration scheme that would be situ state. As a result, the initial size, shape and posi-
unnecessary in a strain space formulation. tion of their yield surfaces may be selected to agree
Finally, Yoder (1981) asserts that stress space mod- with probed stress paths from laboratory tests. How-
els can become unstable for special cases (such as ever, this advantage will disappear once a yield surface
perfect plasticity) unless they are explicitly consid- is engaged.
ered within a computer code. On the other hand, strain The ability of the BRICK model to replicate
space models can seamlessly unify the cases of hard- undrained deformations has been well-established
ening, softening and ideally plastic behavior (Iwan & (e.g. Simpson (1992) and Pillai (1996)); however,
Chelvakumar 1988). comparison of the measured and simulated tests of
high quality undrained triaxial tests on London Clay in
Figure 3 demonstrates that the model is less successful
4 DIFFICULTIES ASSOCIATED WITH THE at capturing the effective stress paths. The orientations
ORIGINAL 3D BRICK MODEL of the undrained stress paths are incorrect in this sim-
ulation partly because the model does not account for
Despite the advantages described above, strain space elastic anisotropy. Moreover, the BRICK effect occa-
models exhibit a few difficulties relative to conven- sionally results in unrealistic ‘kinks’ or rapid changes
tional models that have prevented them from gaining in the stress path direction. These kinks occur when a
widespread popularity. complex recent strain history has changed the location
Rather than defining yield surfaces and a flow rule, of the bricks on short strings relative to the cur-
a user of the BRICK model must explicitly select rent strain point without significantly affecting those
the size and weight of many kinematic surfaces by attached to the longer strings. The resulting orienta-
defining points along an S-shaped curve comprised of tion of the plastic strain increment during a monotonic
material proportions and string lengths as shown in loading sequence can change rapidly once the bricks
Figure 2. This affords the user greater control over the on the longer strings become engaged.
35
Figure 3. Simulations of undrained triaxial compression and drained triaxial extension tests on samples from a depths of
11 m in Unit B2(c) of London Clay using the original 3D BRICK model (lab data from Gasparre (2005)).
As also shown by Figure 3, the model tends to under- and the orientation of a vector connecting the current
predict dilation during drained simulations. It can be stress point to its conjugate point on a larger surface.
seen from Figure 1d that continued shear straining will This ensures that yield surfaces will only intersect
lead to additional volumetric strain until the bricks line tangentially at conjugate points.
up parallel to the shear strain axis. In this manner, the Stress history is considered by the BRICK model
model computes some plastic dilation; however, there in a similar manner except that there may be multi-
is no explicit flow rule. ple active surfaces (i.e. taut strings) at a given time.
It will be shown in a subsequent section that the In fact, as mentioned in a previous section, the for-
incorporation of stiffness anisotropy in the BRICK mulation is even simpler in strain space since there is
model significantly improves the predictions of both no need to force subsequent yield surfaces to intersect
dilation and effective stress paths. However, this is tangentially.
much more difficult to accomplish in BRICK-type
models than in conventional models. In stress space,
it is a relatively trivial task to incorporate stiffness 5.2 Creep, ageing and rate effects
anisotropy by substituting the isotropic elastic stiff- Many studies have highlighted the roles of creep, age-
ness matrix with an anisotropic one. However, in ing and rate effects on soil behavior. It has been shown
the BRICK model, strength and elastic stiffness are that creep and ageing can cause the elastic region
closely intertwined. Therefore, one cannot simply to recenter itself about the current stress state (e.g.
introduce an anisotropic elastic stiffness matrix with- Clayton & Heymann (2001) and Gasparre (2005)).
out inducing an equivalent and undesired anisotropy A related phenomenon known as isotach behavior
of strength. describes how changes in the strain rate applied to
Lastly, unlike the conventional critical state based some soils can cause a jump between different isotach
models, BRICK does not necessarily approach a criti- stress-strain curves (e.g. Suklje (1969)).
cal state line in εvol -ln p space. As a result, its primary The most popular methods to capture creep and
applications are currently limited to stiff clays that isotach behavior due to viscoplastic time-dependent
undergo strain localization before such a line would effects in geomaterials are the nonstationary flow
ever be reached. surface (NSFS) and overstress theories (e.g. Perzyna
(1966) and Naghdi & Murch (1963)). The NSFS the-
5 INCORPORATING SPECIFIC FEATURES IN ory utilizes a variant of the classical elastoplastic yield
BRICK-TYPE MODELS surface that is a function of strain rate. The over-
stress theory postulates that a dynamic yield surface
5.1 Stress history effects exists beyond the static yield surface that depends upon
the strain rate and that these surfaces will gradually
Many constitutive models have employed multiple converge as the strain rate reduces to zero.
kinematically-translating yield surfaces to describe the The principles of the overstress theory readily lend
influence of stress history on the anisotropic hardening themselves for incorporation into BRICK-type mod-
of geomaterials (e.g. Dafalias & Herrmann (1982) and els. Rather than employing rate-dependent dynamic
Stallebrass & Taylor (1997)). It is straightforward to yield surfaces, Clarke (2009) has employed rate-
convert this type of formulation to strain space and this dependent string lengths that gradually converge
is the only extraordinary feature explicitly considered upon their reference values as the stress/strain rate
by the original BRICK model. decreases. Clarke’s strain rate dependent string lengths
In conventional multi-surface models, one or more are determined by the following equation:
surfaces are usually nested within an isotropically-
expanding bounding surface. The largest yield surface
engaged at any time is the active yield surface and the
translation of this surface is a function of normality
36
where β is a material constant and the superscripts shear failure. Therefore, the following equation is
‘tar’ and ‘ref ’ refer to ‘target’ and ‘reference’ values, used:
respectively. However, to avoid a large jump in string
length due to a sudden change in strain rate, a damping
function is introduced:
It is postulated that the same type of string length

correction could be applied to introduce strength
anisotropy; however, such a formulation has yet to be
where α is a material constant, t is time and the incorporated into the model.
superscript ‘n’ refers to the current step number.
Implications of the Clarke (2009) formulation are:
(i) Creep occurs due to the accumulation of plas- 5.4 Stiffness anisotropy
tic strains as the strings contract, (ii) stiffer behavior Numerous laboratory studies have revealed that the
is predicted after creep since all of the strings will small strain linear elastic behavior of soil is not
initially become loose when the strain rate increases isotropic and that the extent of anisotropy is relatively
again, and (iii) isotach behavior can be replicated. constant up to intermediate strains (e.g. Gasparre
(2005) and Yimsiri et al. (2009)). The conventional
isotropic elastic stiffness matrix can easily be substi-
5.3 Strength anisotropy and Lode angle effects tuted with an anisotropic one in stress space models.
It has long been observed that the size of the elastic It is not so simple to incorporate stiffness anisotropy
region and the residual stress ratio achieved by labora- in BRICK-type models for two reasons: 1) elastic
tory specimens are stress path dependent. Thus, many stiffness anisotropy would result in an equivalent
formulations have been proposed to alter the shape of anisotropy of strength; and 2) the BRICK model must
the yield surface in the π-plane with respect to the Lode simulate the entire geologic history over which the
angle θ (e.g. Matsuoka & Nakai (1977)). Such formu- anisotropy has likely developed. An appropriate treat-
lations can easily be incorporated into BRICK-type ment of stiffness anisotropy should account for the
models by modifying the string lengths to be strain evolution of anisotropy with strain without altering the
path dependent. residual strength condition.
For example, the Modified Drucker Prager surface The problem of strength anisotropy arising from
may be adopted, i.e.: elastic stiffness anisotropy can be mitigated by pass-
ing a modified strain increment {ε̇mod } through the
constitutive model, i.e.,
where µ is a constant that controls the shape of the where [Daniso ] and [Diso ] are the anisotropic and
yield surface in the π-plane, q* and r* are invari- isotropic stiffness matrices, respectively.
ants derived from the deviatoric stress tensor and α A matrix [M ] can be defined to convert between
is related to the effective friction angle and controls the real strain increment and the modified strain
the slope of the yield surface in the meridian plane. increment, i.e.,
This criterion is readily employed within models where
the failure surface is imposed as a discrete boundary
around pre-failure behaviour. However, it is not read-
ily adapted to the BRICK model in which stiffness and
strength are interrelated. In this formulation both the elastic and plastic strains
The 3D BRICK model used by Arup employs a will be anisotropic. As a result, there is no guaran-
revised version of Equation 3 rewritten using equiv- tee that the perfectly plastic behavior at critical state
alent strain terms to express the ratio of the polar will be volume-preserving. Thus, in order to achieve a
distance from the volumetric strain axis to the max- constant volume condition, the model must evolve to
imum polar distance (corresponding to triaxial com- become isotropic at critical state (i.e. [M ] → [I ] where
pression), therefore: [I ] is the identity matrix).
Ellison (2009) presents one possible formulation for
stiffness anisotropy in BRICK using the above frame-
work. In this model, direction-dependent anisotropy
increases with the amount of shear strain developed
since a reference value and decreases with the devel-
q
where γbr and γb are invariants of the deviatoric strain opment of shear stress such that the model becomes
along string ‘b’. However, the full string length correc- isotropic near the residual stress ratio. This formula-
tion need not be applied unless the string has reached tion significantly improves the predictions of effective
37
Figure 4. Simulations of undrained triaxial compression and drained triaxial extension tests on samples from a depth of 11 m
in Unit B2(c) of London Clay using the 3D BRICK model with stiffness anisotropy (lab data from Gasparre (2005)).
stress paths during undrained tests and dilation during Ellison, K.C. 2009. Constitutive modeling of London Clay.
drained tests as demonstrated by the simulations in First Year PhD Report. University of Cambridge.
Figure 4. Gasparre, A. 2005. Advanced laboratory characterization of
London Clay. PhD Thesis. Imperial College London.
Iwan, W. D. & Chelvakumar, K. 1988. Strain-space con-
stitutive model for clay soils. J. Eng. Mech., 114(9):
6 CONCLUSIONS 1454–1472.
Matsuoka, H. & Nakai, T. 1977. Stress-strain relationship of
In recent decades, researchers have identified several soil based on the “SMP”. Proc. Specialty Session 9, IX
shortcomings of conventional constitutive models to ICSMFE, Tokyo, 153–162.
capture certain aspects of soil behavior, such as stress Naghdi, P. M. and Murch, S. A. 1963. On the mechanical
history effects, strength anisotropy, time-dependence behavior of viscoelastic/plastic solids. J. Appl. Meterorol.,
and stiffness anisotropy. While most methods to 30, 321–328.
address these features have been introduced to stress Perzyna, P. 1966. Fundamental problems in viscoplasticity.
Adv. App. Mech., 9, 244–377.
space models, similar advances have been applied to Pillai (Kanapathipillai), A. K. 1996. Review of the BRICK
the BRICK model in strain space. model of soil behaviour. MSc dissertation, Imperial Col-
This paper has reviewed the advantages and diffi- lege, London.
culties of using strain as the independent variable in Puzrin, A. M. & Houlsby, G. T. 2001. On the non-intersection
a constitutive model. Moreover, it has described how dilemma in multiple surface plasticity. Géotechnique,
several advances in the modeling of specific features 51(4): 369–372.
in stress space can be modified for implementation in Simpson, B. 1992. Retaining structures: displacement and
BRICK-type models. design. Géotechnique, 42(4): 539–576.
Stallebrass, S. E. & Taylor, R. N. 1997. The development
and evaluation of a constitutive model for the prediction
of ground movements in overconsolidated clay. Géotech-
REFERENCES nique, 47(2): 235–254.
Suklje, L. 1969. Rheological aspects of soil mechanics, Wiley
Atkinson, J. H., Richardson, D. & Stallebrass, S.E. 1990.
Interscience, London.
Effect of stress history on the stiffness of overconsolidated
Yimsiri, S. & Soga, K. 2009. The anisotropy of two natural
soil. Géotechnique, 40, 531–540.
stiff clays. Submitted to Géotechnique.
Burland, J. B. 1967. Deformation of soft clay. PhD Thesis.
Yoder, P.J. 1981. A strain-space plasticity theory and numer-
University of Cambridge.
ical implementation. PhD Thesis. California Institute of
Clarke, S. D. 2009. Enhancement of the BRICK constitutive
Technology.
model to incorporate viscous soil behavior. PhD Thesis.
Yoder, P. J. & Iwan, W. D. 1981. On the formulation of strain-
University of Sheffield.
space plasticity with multiple loading surfaces. J. Appl.
Clayton, C. R. I. & Heymann, G. 2001. Stiffness of geo-
Mech., 48(4):773–778.
materials at very small strains. Géotechnique. 51(3):
245–255.
38
Anisotropic small strain stiffness within the multilaminate framework
B. Schädlich & H.F. Schweiger

Computational Geotechnics Group, Institute for Soil Mechanics and Foundation Engineering, Graz University of
Technology, Graz, Austria
ABSTRACT: Using the spectral decomposition of the global compliance matrix, a novel approach to modelling
anisotropic elasticity within the multilaminate framework is presented. The new approach is implemented into
a soil model which accounts for degradation of small strain stiffness with increasing shear strain and stress
dependency of stiffness. The model is calibrated by back-analysis of element tests on London Clay and applied
in a Finite Element calculation to evaluate the influence of anisotropic small strain stiffness on deformations
connected with tunnel excavation.
1 INTRODUCTION The transformation matrix Ti contains the derivatives

of the local stress components with respect to the
The high initial elastic stiffness of soils at very small global axes, represented by the direction cosine of the
strains (<10−6 ) and its degradation with accumulation unit vector nTi = (ni,1 , ni,2 , ni,3 ) normal to the plane i and
of strain is well known since the early 1970ies. With of two unit vectors within the plane, sTi = (si,1 , si,2 , si,3 )
increasing progress in both laboratory testing and soil and tTi = (ti,1 , ti,2 , ti,3 ). Vectors ni , si and ti must form
modelling, taking that effect into account in geotech- an orthogonal system of local axes, such that ni ·si = 0,
nical engineering has become more and more common ni ·ti = 0 and ti ·si = 0.
practice within the last decade.
Still, in most practical cases soil is assumed to
behave isotropically at very small strains, although
laboratory tests on natural soils indicate strongly cross-
anisotropic behaviour (Gasparre 2005, Fioravente
2000). The local elastic strains εi,loc are calculated as
In the following study an approach to model inher-
ently cross-anisotropic elastic material and degrada-
tion of small strain stiffness is presented. The model
is applied to a tunnelling problem in order to evaluate In the case of isotropic linear elastic material, Cloc is
the effect of initial anisotropy on deformations within equal for all planes. For non-linear elasticity (small
the tunnel and at the ground surface. strain stiffness), Cloc depends on the strain history of
each plane and therefore differs from plane to plane,
resulting in anisotropic global behaviour.
Global strains are obtained by back-transformation
2 MULTILAMINATE FRAMEWORK and summation of all local strains:
Multilaminate material models are based on the con-
cept that the material behaviour can be formulated on
a distinct number of local planes with varying orien-
tation. The stress – strain state varies from plane to The factor of 3 in front of the summation can be derived
plane, resulting in loading induced anisotropy within from the principle of virtual work (Bažant & Prat
an intrinsically isotropic material. The global response 1988). The weight factors wi depend on the chosen
of the material to a prescribed load is obtained by integration rule. In this study an integration rule based
summation of the contributions of all planes. on 2 × 33 planes (Bažant & Oh 1985) is used.
The local stress vector σi,loc is obtained by project- While local stresses are a projection of the global
ing the global stress vector σgl with the transformation stress state (static constraint), local strains are in gen-
matrix Ti on plane i. eral not the projection of global strains (kinematic
constraint) in multilaminate models but they are in
so-called microplane models (Bažant & Prat 1988).
39
3 ANISOTROPIC SMALL STRAIN STIFFNESS Using the idempotent matrices E1 . . . E4 , which are
MODEL defined by the 4 eigenvectors of Cgl , the global stress
vector can be split up into its spectral components or
3.1 Concept stress modes σgl,1 . . . σgl,4 .
In previous multilaminate-type soil models it was pos-
tulated, that the local stress state could be represented
by 3 components, whose directions coincided with the
direction of the vectors n, s and t (Scharinger et al.
2008). That assumption results in a 3 × 3 local elas-
tic compliance matrix Cloc and a 3 × 6 transformation
matrix T. For elasticity it was further assumed, that on
local level normal strains are only caused by normal
stresses and tangential strains are only caused by tan-
gential stresses. Therefore, Cloc was a diagonal matrix
with elements outside the main diagonal equal to 0.
For anisotropic material, the aforementioned
assumptions can no longer be maintained. Isotropic
compression of an anisotropic material results in shear
strains on all planes which are not parallel to the global
axes, although only normal stresses are obtained on
these planes from Equation 1. Anisotropic material
behaviour can therefore not be modelled by using a
diagonal local compliance matrix.
The spectral decomposition of the global stress
vector offers the possibility to obtain local compli-
ance matrices directly. Cross-anisotropic material with
a vertical axis of symmetry is considered further
on, although the method is also applicable to fully
anisotropic material. Only the step-by-step procedure
will be demonstrated in this paper. For details on
the theoretical background see Theocaris & Sokolis
(2000) and Cusatis et al. (2008).
The global compliance matrix Cgl of a cross
anisotropic elastic material is fully defined by 5 param-
eters: two elastic moduli Ev and Eh , one independent
shear modulus Gvh , and two Poisson’s ratios, νvh and
νhh . If written in Kelvin notation, Cgl possesses four
eigenvalues, λ1 . . . λ4 .
Local stress modes σi,loc,m on plane i are obtained by

projecting each global stress mode separately using
transformation matrices according to Equation 2.
Local strain modes εi,loc,m are calculated by multiply-

ing each local stress mode σi,loc,m separately with the
corresponding eigenvalue λm (Equ. 11). The sum of all
local strain modes yields the local strain vector εi,loc ,
and the sum of all local stress modes equals the local
stress vector σi,loc (Equ. 12).
40
For this split, the transformation matrix of plane i can
be written as
Back transformation to global level and summation

of the local strains follows the same procedure as for
isotropic material (Equ. 4).
It should be noted, that in this procedure the local
stresses depend on global material parameters, while
local strains are obtained by multiplying the local
stress modes with the scalar eigenvalues. That differs
quite substantially from multilaminate models, where
local stresses only depend on plane orientation, and all
elastic material properties are described by the local
compliance matrix Ci,loc .
However, local compliance matrices can also be Ci,loc depends on the plane orientation and con-
derived directly from the global compliance matrix. tains non-zero off-diagonal elements for general non-
Combining Equations 7 and 10–12 yields. isotropic material. In the case of isotropic elastic
material, Cgl has only two unique eigenvalues (Equ.
18), yielding a diagonal local compliance matrix
(Equ. 19).
Comparing Equation 13 with Equation 1 and 3,

Equation 14 is found.
As the matrices Em are of the order 6 × 6, Equation

13 only has a unique solution if both Ti and Ci,loc also Both local stresses and local strains are projections of
are 6 × 6 matrices. In that case Equation 13 can be the corresponding global quantities. Therefore, both
transformed to the static and the kinematic constraint are fulfilled.
3.2 Stress dependency of stiffness

Laboratory test data indicate, that in cross-anisotropic
natural soils Eh depends on σh and Ev depends on σv
(Kuwano & Jardine 2002). However, using such an
That means, that the local stress vector is split up into approach in boundary value problems is prone to cause
six components. For the present study local stress and numerical problems at stress free boundaries because
strain components on plane i are defined as high stresses in one direction and stresses close to
zero in the other numerically induce extreme ratios
of anisotropy in the material, which causes the global
stiffness matrix to become almost singular.
In order to avoid such problems, in this study global
stiffness parameters Ev , Eh , and Gvh are assumed to
depend on the mean effective stress p only.
with σi,n,vol . . . volumetric normal stress, σi,n,dev . . .

deviatoric normal stress, τi,s1 and τi,t1 . . . tangential
stresses in direction of s and t resulting from global
axial stresses, τi,s2 and τi,t2 . . . tangential stresses in With that approach a reference local compliance
direction of s and t resulting from global shear stresses. matrix Cloc,ref can be established, calculated from
41
Table 1. Elastic soil properties of London Clay.
Parameter Isotropic Anisotropic Unit
Eur,ref 13000 13000 kPa

pref 100 100 kPa
m 1.0 1.0 –
νur 0.2 0.2 –
Ev0,ref 48960 30000 kPa
Eh0,ref 48960 78000 kPa
Gvh0,ref 20400 20400 kPa
νhh 0.2 0.02 –
νvh 0.2 −0.16 –
γ1 0.0025 0.0025 %
γ2 0.03 0.03 %
Figure 1. Degradation of anisotropic small strain stiffness.
global stiffness parameters Eh,ref , Ev,ref , Gvh,ref at the

reference pressure pref . The local compliance matrix
at the current stress level is then obtained according to
Equation (21).
3.3 Degradation of stiffness

Experimental data from laboratory tests show a
S-shaped degradation of the initial stiffness with accu-
mulated shear strain (e.g. Gasparre 2005). Various
functions describing the degradation of small strain
stiffness can be found in the literature, involving
trigonometric, exponential and logarithmic functions Figure 2. Degradation of equivalent shear modulus.
(e.g. Jardine et al. 1986, Benz 2007).
For the anisotropic multilaminate model, it
is assumed that the initially anisotropic material
behaviour, the reader is referred to Schweiger et al.
approaches isotropy with increasing accumulated local
(2009) for details on the plasticity part of the model.
shear strain γ. Stiffness degradation of Eh follows
Equation 22, for the other anisotropic parameters
equivalent equations apply. At shear strains larger than
γ2 , the material at local level is isotropic, described by 4 ELEMENT TESTS ON LONDON CLAY
the elastic modulus Eur and Poisson’s ratio νur .
As the development of local shear strains differs 4.1 Small strain parameters
from plane to plane, also local stiffness parameters
vary over the planes, thus resulting in a smooth transi- Recent experimental data on anisotropic small strain
tion from small to large strain behaviour on global stiffness of London Clay have been published by
level. The local shear strains γ1 and γ2 have to be Gasparre (2005). The samples were retrieved from
determined by back-analysis of laboratory tests. the site of the Heathrow Terminal 5 and tested
using bender element aided triaxial tests. For a depth
of 22.6 m, values of Ev0 = 110 Mpa, Eh0 = 285 Mpa,
Gvh0 = 75 Mpa, νvh = 0.02 and νhh = −0.16 are
reported. Reference values at 100 kPa are listed in
Table 1. Back analysis of the equivalent shear modu-
lus Geq in undrained triaxial compression tests yielded
local shear strains γ1 = 0.0025% and γ2 = 0.03% for
fitting experimental data (Fig. 2).
However, the same initial shear stiffness Geq and
degradation curve could also be obtained assuming
3.4 Plastic strains
isotropic material behaviour (Fig. 2). In that case,
Once the large strain region is reached locally, the Geq would equal the isotropic shear modulus G.
model can also account for strain hardening plastic- Setting ν = 0.2, isotropic values can be derived as
ity. As this study is focused on elastic small strain G0,ref = 20.4 Mpa and E0,ref = 48.96 Mpa.
42
Figure 3. FE-model and boundary conditions.
4.2 Elastic large strain parameters

The inclination of the unloading/reloading line in Figure 4. Vertical displacement point A.
isotropic compression of natural samples is reported
as κ = 0.029 (Gasparre 2005). With V = 1 + e = 2.12
(specific volume at reference pressure pref ), νur = 0.2
and pref = 100 kPa, the unloading reloading stiffness
can be obtained as Eur,ref = 13000 kPa according to
Equation 24.
5 INFLUENCE OF ANISOTROPIC SMALL

STRAIN STIFFNESS
5.1 FE-model and boundary conditions

The soil model described above is utilized to investi-
gate the influence of small strain stiffness anisotropy Figure 5. Vertical displacement point B.
on tunnel induced surface settlements and displace-
ments at the tunnel cross section. The model developed
by Scharinger et al. (2008) has been extended to
account for anisotropic small strain stiffness. Although
the model can also account for plasticity in the large
strain range, only elastic strains are considered in this
study.
The tunnel centre is situated at 30.5 m depth, diam-
eter of the circular tunnel is 4.75 m. Soil layering is
simplified to only one soil layer (Fig. 3).
The calculations are performed with the Finite
Element code PLAXIS2D V9.0, using triangular
15-noded elements. Three different sets of soil prop-
erties are considered:
• Set 1: no small strain stiffness
• Set 2: isotropic small strain stiffness (Table 1) Figure 6. Horizontal displacement point C.
• Set 3: anisotropic small strain stiffness (Table 1)
5.2 Results
The following boundary conditions are applied in
all calculations: The displacements due to tunnel excavation are com-
pared at the following points:
• Ground water table 5 m below ground surface
• K0 = 1.5 (constant over depth) • Point A – at ground surface above tunnel
• Hydrostatic pore water pressure • Point B – at tunnel crown
• Point C – at tunnel bench
For simplicity drained conditions are assumed in
this study. Starting from the initial stress state, the The development of vertical and horizontal dis-
nodal forces of the tunnel boundary are subsequently placements with reduction of initial nodal forces
reduced from 100% to 0. within the tunnel is summarized in Figures 4–6.
43
Table 2. Displacements [mm] at 40% stress relaxation. needs to be investigated with more sets of parameters
and also in different boundary value problems.
set 1 set 2 set 3 set 3 / set 2
Point A – uv −2.1 −1.2 −1.4 120% REFERENCES

Point B – uv −10.8 −7.0 −7.6 109%
Point C – uh 16.7 10.4 8.6 83% Addenbrooke, T.I., Potts, D.M. & Puzrin, A.M. 1997. The
influence of pre-failure soil stiffness on the numerical
analysis of tunnel construction. Geotechnique 47 (3):
The curves for point B and C become parallel at 693–712.
40–50% excavation, indicating the complete loss of Bažant, Z.P. & Oh, B.H. 1986. Efficient Numerical Integra-
small strain stiffness within the soil volume relevant to tion on the Surface of a Sphere. Zeitschrift für angewandte
these points. For point A slightly different inclinations Mathematik und Mechanik 66: 37–49.
can be found even at full relaxation. The ratio Bažant, Z.P. & Prat, P.C. 1988. Microplane Model for
Brittle-Plastic Material: I. Theory. Journal of Engineering
Mechanics 114(10): 1672–1688.
Benz, T. 2007. Small Strain Stiffness of Soils and its
Numerical Consequences. Ph.D. Thesis. Mitteilung 55 des
at 40% stress relaxation varies from 83% to 120% Instituts für Geotechnik, Universität Stuttgart.
(Table 2). Although fu,aniso follows the ratio of isotropic Cusatis, G., Beghini, A. & Bažant, Z.P. 2008. Spectral Stiff-
ness Microplane Model for Quasibrittle Laminates – Part
vs. anisotropic stiffness (with uv being governed by
I: Theory. Journal of Applied Mechanics 75(2): (021009)
1/Ev and uh by 1/Eh ), the stiffness ratio is consid- 1–9.
erably higher than fu,aniso (Ev0 / Eiso = 61%; Eh0 / Fioravante, F. 2000. Anisotropy of small strain stiffness of
Eiso = 160%). Ticino and Kenya sands from seismic wave propagation
measured in triaxial testing. Soils and Foundations 40(4):
129–142.
6 CONCLUSION Gasparre, A. 2005. Advanced laboratory characterisation of
London Clay. PhD thesis, Imperial College, London.
A new approach for modelling anisotropic, stress - Jardine, R.J., Potts, D.M., Fourie, A.B. & Burland, J.B. 1986.
Studies of the influence of non-linear stress-strain char-
dependent small strain stiffness within the multil-
acteristics in soil-structure interaction. Geotechnique 36
aminate framework has been developed. The stress (3): 377–396.
dependency of stiffness currently implemented in Kuwano, R. & Jardine, R.J. 2002. On the applicability of
the model does not fully agree with experimentally cross-anisotropic elasticity to granular materials at very
observed soil behaviour and requires further investi- small strains. Geotechnique 52 (10): 727–749.
gation. Regarding the influence of anisotropy in the Scharinger, F., Schweiger, H.F. & Pande, G.N. 2008. On a
small strain range on practical boundary value prob- multilaminate model for soil incorporating small strain
lems, the study must be seen as preliminary. Only stiffness. International Journal for Numerical and Ana-
one set of anisotropic parameters was used, and only lytical Methods in Geomechanics 33(2): 215–243.
Schweiger, H.F., Wiltafsky, C., Scharinger, F. & Galavi, V.
elastic deformations in a tunnelling problem were con-
2009. A multilaminate framework for modelling induced
sidered. However, in the case investigated isotropic and inherent anisotropy of soils. Geotechnique 59 (2):
and anisotropic small strain stiffness result in simi- 87–101.
lar displacements, if both sets of parameters fit the Theocaris, P.S. & Sokolis, D.P. 2000. Spectral decomposi-
degradation curve of the equivalent shear modulus tion of the compliance fourth-rank tensor for orthotropic
Geq . Whether this is a general trend or just coincidence materials. Archive of Applied Mechanics 70(4): 289–306.
44
Application of discontinuity layout optimization to problems involving

non-associative friction
A.F. Babiker, C.C. Smith & M. Gilbert

Department of Civil and Structural Engineering, University of Sheffield, UK
ABSTRACT: According to plasticity theory, when a frictional material is sheared, dilatancy will accompany
sliding deformation, according to the ‘associated flow rule’, or ‘saw-tooth’ idealization of friction. However,
although real granular soil materials will generally be observed to dilate when sheared, the amount of dilatancy
measured experimentally will generally be much less than assumed by classical plasticity theory, potentially
leading to non-conservative estimates of safety being obtained. For many types of geotechnical problems the
degree of non-conservatism involved will be small, but in the case of highly confined problems the error can be
more significant. A variety of means of addressing this issue have previously been proposed and in this paper
the scope for using an iterative procedure developed for application to masonry structures in conjunction with
the Discontinuity Layout Optimization (DLO) numerical limit analysis procedure is examined.
1 INTRODUCTION not be significant, differences can become large when

the degree of confinement is high.
1.1 Background For a given problem there is therefore a need to:
Limit analysis provides a powerful means of assess- 1. Identify a measure of the range of feasible non-
ing the stability of a broad range of geotechnical associative solutions.
engineering problems. One benefit of the method is 2. For each non-associative solution, determine a real-
that the bound theorems of plasticity theory can be istic failure mechanism (which does not involve
used to ensure that solutions have definite status, i.e. excessive dilation).
that the solution obtained provides either an over-
For many problems the above may simply mean
or under-estimate of the true collapse load, or ‘load
showing that the associative solution is not signifi-
factor’.
cantly affected by the flow rule used.
However, to enable the formal theorems of plasticity
theory to be applied, associated plastic flow has to be
assumed. i.e. flow occurs in a direction normal to the
failure surface, according to the ‘normality rule’. This
means that the angle of dilation ψ must be taken to be 1.2 Previous work
equal to the angle of shearing resistance φ. In practice In previous work by Gilbert et al. (2006), a numerical
this assumption leads to satisfactory analysis predic- method was proposed to allow the stability of large-
tions in the case of purely cohesive materials, where scale masonry block structures with non-associative
φ is effectively zero, but predictions are less realistic frictional joints to be evaluated. The iterative method
in the case of frictional materials. Whilst experiments proposed involved the use of a succession of modi-
have demonstrated that dilation does occur in frictional fied failure surfaces for each interface (masonry joint)
soil materials, this is very much less than predicted by in the problem, with each failure surface ascribed
the normality rule and in some cases the dilation can a fictitious cohesion to limit the shear force in line
be zero, i.e. involving shearing of soil at a constant vol- with the actual Mohr Coulomb failure envelope being
ume (Cox 1963). Thus a frictional soil should ideally used. The fictitious cohesion was calculated accord-
be modelled assuming 0 < ψ < φ. (Note that in this ing to the magnitude of the normal force computed
paper φ is used in place of φ for sake of simplicity.) in previous iterations. The iterative procedure required
It has been known for many decades that severe solution of a succession of simple linear programming
difficulties can arise when calculating limit loads in (LP) problems, and proceeded until a converged solu-
the presence of non-associative friction, and that a tion, not involving dilation, was obtained. The method
wide range of possible solutions to such problems appears to be capable of evaluating realistic non-
exist (Drucker 1954). Although indications are that for associative solutions, with computed collapse loads
many practical geotechnical problems the differences always lower, or equal to, those obtained assuming
between associative and non-associative solutions will associative friction.
45
It is of interest to investigate whether the same basic
method can be applied to soil plasticity problems. One
obvious difference the problems is that whereas the
locations of potential planes of weakness are well-
defined in the case of masonry block problems (i.e.
along the joints), the continuum nature of the material
means that this is not true in the case of soil problems.
However, the recently developed discontinuity layout
optimization (DLO) procedure would appear to pro-
vide a means of addressing this issue, as will now be
described.
1.3 Discontinuity Layout Optimization (DLO) Figure 1. Original and modified failure surfaces which pro-
DLO is a novel analysis procedure which is applicable vide the same limit on the shear stress (for constant normal
to a wide range of limit analysis problems (Smith & stress σn ).
Gilbert 2007). It is capable of determining the critical
layout of discontinuities in a body of soil at collapse. Referring to Figure 1, consider a point X lying
Using DLO an upper-bound limit load for a cohesive- on the Mohr-Coulomb failure surface (indicated by
frictional soil with an associated flow rule can be the solid line), where the normal stress is given by
calculated by first discretizing the soil body under σ = σn , and the shear stress τ = c + σn tan φ. The asso-
consideration using closely spaced nodes, and then ciated flow rule clearly requires ψ = φ (i.e. flow in the
inter-connecting each node to every other node with direction of the solid arrow), whereas the requited non-
potential slip-lines. LP can then be used to identify the associated flow, with ψ = 0, will be in the direction
critical subset of slip-lines that form at collapse. indicated by the dashed arrow.
The problem can be posed in a kinematic form, In order to ensure ψ = 0, a fictitious failure sur-
where the LP variables represent the displacements face can be constructed, represented by the dashed line
along the slip-lines, and the objective function is to (ĉ = c + σn tan φ). This still correctly limits the shear
minimize the energy dissipated at collapse; alterna- stress at X provided the normal stress remains con-
tively the dual (equilibrium) form can be posed, which stant. In principle the failure surfaces corresponding
requires that constraints are imposed on the shear to all points within a soil body would be replaced in
(T ) and normal (N ) forces along the discontinuities. this way and the problem solved again. There is how-
Whichever formulation (i.e. kinematic or equilibrium) ever no guarantee that the normal stresses will remain
is used, duality theory means that results from the constant after re-solving, and hence a number of iter-
other, dual, formulation can also be obtained. The solu- ations may be required before a converged solution is
tion also provides a load factor λ which is a multiplier reached.
on a specified live load or loads required to generate Assuming a converged solution can be obtained, this
the identified collapse mechanism. solution will satisfy both the original failure surface at
all points, and will ensure that flow is non-associative
as required.
1.4 Aim of this paper Failure typically involves transformation of a soil
body into discrete blocks of soil, separated by disconti-
The aim of this paper is to present results from a pre- nuities, and hence in DLO the requisite failure criterion
liminary study into the viability of using an approach is checked along potential slip-lines, rather than at spe-
similar to that originally described by Gilbert et al. cific points within a soil body. The approach described
(2006) in conjunction with DLO, and to identify the previously therefore needs instead to be formulated
likely future direction of the research. To achieve this, in terms of shear and normal stress resultants, T and
three simple example problems which are amenable to N respectively. The material failure criterion becomes
hand calculation are considered in the paper, in order |T | = C + N tan φ, and the modified failure condition
to allow a good understanding of the method to be (assuming zero dilation) becomes:
built up.
2 PROPOSED ALGORITHM
2.2 Iterative procedure for finding a
2.1 Underlying principles non-associative solution
For sake of simplicity a cohesive-frictional mate- The iterative solution procedure involves a number of
rial with strength parameters c and φ, and dilation steps:
angle ψ = 0 will be considered throughout this paper.
However, the same basic method can also applied to 1. Assume initial modified shear strength parameters
problems with non-zero dilation angle. cˆ0 , and φˆ0 of arbitrary value for all discontinuities,
46
where Ĉ0 = ĉ0 l and l is the length of a slip-line.
Solve the resulting DLO problem. The initial nor-
mal Nk,0 and shear Tk,0 forces can be extracted from
the solution for each discontinuity k, together with
the load factor λ0 .
2. At the next iteration i the modified shear strength
parameter Ĉk,i for discontinuity k can be computed
using the normal force from the previous iteration
Nk,i−1 as follows:
where φ is the actual material angle of friction.

(In this paper only non-associative cases with
zero-dilation are considered. Thus φ̂k,i = 0.)
3. Solve the associative DLO problem using parame-
Figure 2. Sliding block example: rigid block sliding on a
ters Ĉk,i and φ̂k,i , to obtain new values of Nk,i , Tk,i , single interface (c = 0, φ = 20◦ ).
and λi
4. If i > 1 and | λi − λi−1 | / λi < tolerance, and vio-
lation of the real failure surface is not detected at Table 1. Sliding block example: final results.
any discontinuity k, then the algorithm stops.
5. If there is no convergence the process is repeated ĉ0 φ̂0 Iterations Computed λ Analytical λ Diff%
from step 2 until convergence is reached.
2 0 2 0.364 0.364 0.0
This algorithm was implemented using a modified 0.1 0 2 0.364 0.364 0.0
version of a MATLAB formulation of the DLO code
presented by Gilbert et al. (2009), in conjunction with
the Mosek LP solver.
The problem is statically determinate and the live
load P = λ.1 required to move the block can be
3 NUMERICAL EXAMPLES calculated from Equation (3) to be 0.364.
3.1 Introduction
The following simple examples are designed to provide As is evident from Table 1, the correct solution can
insights into the proposed numerical procedure. The be obtained using the proposed algorithm in 2 itera-
first two examples are statically determinate, which tions. While this is a trivial case, the actual kinematic
should ensure that the computed load factor is not mechanism of collapse is now modelled with no dila-
dependent on the flow rule. tion along the sliding surface, in contrast to when the
For sake of simplicity the blocks used in all exam- associative friction is used.
ples are all weightless, square and of unit dimensions.
Where a dead or live load is applied it is applied as a
unit edge load. Since the computed solution is a mul- 3.3 Sliding wedge example
tiple of 1 it can be interpreted as either a limit load
P = λ.1 or as a limit load factor λ. In the problem shown in Figure 3, a unit square block is
For each example tabulated results are provided restrained at the base and right hand edge and subject
for specified initial conditions. Two calculated solu- to a normal live load P = λ.1 along the left edge. The
tions are shown for each problem, one obtained using block is only allowed to shear along the diagonal. As
a high starting value for the cohesion ĉ0 and the failure can only occur along the diagonal, when c = 1
other obtained using a low starting value for ĉ0 . In and φ = 20◦ the procedure would be expected to con-
all cases the solution converged to the same value verge to the associative solution of 3.145. This value
regardless of the starting values. The last column of can be determined from equation (4), derived from
each table of results provides a comparison between application of the sine rule to the polygon of forces
the known closed form solution and the computed shown in Figures 4a and 4b.
solution.
3.2 Sliding block example

Figure 5 shows how the failure surface for the
In this example an applied unit normal load (N = 1) is diagonal slip-line is modified at each iteration until
applied to the top edge of a square block as shown in convergence is reached. (The final iterations have been
Figure 2. omitted for sake of clarity.)
47
Figure 3. Sliding wedge example (c = 1, φ = 20◦ ).
Figure 6. Sliding wedge example: progress of iterative

solution procedure (using high and low starting cohesion
values).
Table 2. Sliding wedge example: final results.

Figure 4. Sliding wedge example: equilibrium of forces.
cˆ0 ψ̂o Iterations Computed λ Analytical λ Diff%
0.2 0 11 3.145 3.145 0.0

10 0 12 3.145 3.145 0.0
is sought. Slip may occur along top and bottom inter-

faces and/or along the single diagonal slip-line shown.
Assuming associative friction, slip must occur along
the top and bottom interfaces and also the diagonal
(Figure 7a). When non-associative friction is involved,
with ψ = 0, failure no longer has to involve slip along
the diagonal (Figure 7b).
In the latter case the mobilized angle of friction
φmob on the diagonal is unknown; however it is limited
as follows: −φ < φmob < φ. The associative limit load
Figure 5. Sliding wedge example: iterative modification of
failure envelope lines. and the bounds on the non-associative limit load can
thus be calculated by examining the force polygons
shown in Figure 8. The resultant force R acting across
It is clear from Figure 6 that the solution improves the diagonal is orientated at an angle β = 45 − φmob to
with each iteration, and converges regardless of the vertical (taking shear as clockwise positive).
whether low or high starting cohesion values are used. From considerations of overall horizontal equilib-
Table 2 shows the results from two analyses with differ- rium (Figure 8a):
ent starting values of initial cohesion ĉ0 . The solution
is the same for both cases (to within the required
tolerance). Consideration of the equilibrium of block B (Fig-
ure 8b) gives:
3.4 Wedged block example
In this problem a weightless block is wedged between
a fully rough rigid surface at the base and a smooth Finally combination of equations 5 and 6 gives:
rigid surface at its top, as shown in Figure 7. The force
P = λ.1 required to cause the block to slide against
a unit restraining force applied to the right hand side
48
Figure 7. Wedged block example: a block wedged at its base
and top (c = 0, φ = 30◦ o). Figure 9. Wedged block example: progress of iterative
solution procedure (using high and low starting cohesion
values).
Figure 8. Wedged block example: equilibrium of forces.
Table 3. Wedged block example: final results.
ĉ0 φ̂o Iterations Computed λ Analytical λ Diff%
0.2 0 18 1.155 1.155 0.0

10 0 18 1.155 1.155 0.0
The associative solution is given when φmob = φ and Figure 10. Wedged block example: iterative modification
results in a closed form solution of 3.155. The mini- of failure envelopes for (a) base and (b) diagonal slip-lines.
mum non-associative solution that can be obtained for
the limit load occurs when the value of the mobilised
shear angle φmob is (−30) degrees, the collapse load calculated closed form analytical solutions, and the
then being 1.155. associated predicted failure mechanisms involved no
The same problem was solved using the proposed dilation. Further tests have indicated that convergence
algorithm and the results presented in Table 3. It can is achieved in the case of statically determinate prob-
be seen that the final converged solution is equivalent lems regardless of the starting values used; however
to the lowest possible non-associative solution and the this remains to be proved rigorously.
failure mechanism involves the block sliding along its Further investigation of the third, statically inde-
base and upper surface only with no dilation. terminate problem indicates that convergence of the
The solution converges after 18 iterations, as shown iterative procedure is sensitive to the values of T and
in Figure 9. N on the diagonal slip-surface generated by the LP
The iterative modification of the failure envelope solver at each iteration. In this problem the LP solver
for the case with ĉ0 = 0.2 is shown in Figure 10. It is free to return a range of values for T and N .
can be seen that the solution ‘spirals’ in on the final To generate the results presented in Figures 9 and
minimum non-associative solution. 10, the solver must return the lowest possible value of
T at each iteration. To achieve this result in all cases
it would be necessary to adopt an alternative approach
4 DISCUSSION to force T to its minimum value. An approach under
consideration is to include a small additional negative
For the simple problems examined, the calculated friction angle in the failure surface at each iteration,
non-associative limit loads were in agreement with as described in Gilbert et al. (2006).
49
Future aims for this research project include: ACKNOWLEDGEMENTS
• Investigation and development of robust heuristics
The first author acknowledges the support of EPSRC
to guarantee convergence for all problems involving
(DTA studentship).
non-associative friction.
• Modelling of more challenging and practical prob-
lems (e.g. foundation footing problems). REFERENCES
• Development of techniques to allow the range of
non-associative collapse loads to be established for Cox, A. D. (1963). The use of non-associated flow rule in soil
any given problem. plasticity. R.A.R.D.E Report B 2/63.
• Development of a simple test to allow the sensitivity Drucker, D. C. (1954). Coulomb friction, plasticity, and limit
to dilation angle to be checked for a given problem. loads. J. Appl. Mech., ASME 21(1), 71–74.
Gilbert, M., H. M. Ahmed, & C. Casapulla (2006). Limit
analysis of masonry block structures with non-associative
frictional joints using linear programming. Computers and
5 CONCLUSIONS Structures 84(3), 873–887.
Gilbert, M., C. Smith, I. Haslam, & T. Pritchard (2009).
An iterative procedure originally developed for appli- Plastic limit analysis using discontinuity layout optimiza-
cation to masonry structures has been used in con- tion (DLO). In 17th UK Conference on Computational
junction with the Discontinuity Layout Optimization Mechanics (ACME-UK), 6–8 April 2009, Nottingham.
(DLO) numerical limit analysis procedure to gen- ACME-UK.
erate solutions to soil plasticity problems involving Smith, C. C. & M. Gilbert (2007). Application of discontinu-
non-associative friction has been outlined. ity layout optimization to plane plasticity problems. Pro-
The procedure is shown to be capable of predict- ceedings of the Royal Society A: Mathematical, Physical
and Engineering Sciences 463(2086), 2461–2484.
ing the lowest non-associative collapse load for three
simple example problems. However further work is
required to demonstrate robustness of the algorithm.
Once this has been undertaken the procedure will be
applied to larger, more practical, problems of interest
to practicing geotechnical engineers.
50
Associated plasticity for nonassociated frictional materials
Kristian Krabbenhoft, Andrei V. Lyamin & Scott W. Sloan

Centre for Geotechnical and Materials Modelling, University of Newcastle, NSW, Australia
ABSTRACT: A new computational plasticity scheme for nonassociated frictional materials is presented. While
general, it relies solely on well established concepts of associated plasticity. The new scheme is applied to
some common boundary value problems for which the consequences of nonassociated flow rules in terms of
localization of deformations are highlighted.
1 INTRODUCTION an associated flow rule it is here possible, by appro-

priate construction of the hardening law, to achieve
The application of plasticity theory to frictional mate- a response that with non-hardening models only is
rials such as sand, clay, rock, and concrete introduces attainable by means of a nonassociated flow rule. Cam
a number of complications, the most prominent of clay is an example of such a model. Although such
which relates to the flow rule. As is well known, the models appear to circumvent the need for a nonas-
flow rule associated with a relevant yield criterion, sociated flow rule, many of the same problems that
for example Mohr-Coulomb, predicts excessive plas- plague perfectly plastic nonassociated models are in
tic dilation. Consequently, a nonassociated flow rule fact encountered again. Indeed, the nonassociativity
must be used. Although seemingly straightforward, the of the flow rule is essentially transferred to the hard-
introduction of a nonassociated flow rule gives rise to a ening law, resulting in much the same complications
number of complications that manifest themselves par- (Hjiaj et al., 2005a; Krabbenhoft, 2009).
ticularly in the numerical solution of boundary value These facts motivate a closer look at the physical
problems. These complications can be divided into two origins of nonassociated flow rules and the numerical
categories. methods used to solve problems of frictional plastic-
Firstly, from a mathematical point of view, nonas- ity. In the following, inspired by the micromechanical
sociated flow rules often lead to a situation where origins of friction and its modeling in terms of plastic-
the boundary value problem, at some characteristic ity theory, a new approach to computational plasticity
stress state, goes from being elliptic to being hyper- for frictional (and generally nonassociated) materials
bolic (Bigoni and Hueckel, 1991). Physically, this loss is presented. The new scheme is then applied to some
of ellipticity indicates an instability where a homoge- common boundary value problems that also highlight
neous mode of deformation gives way to a localized the consequences of nonassociated flow rules in terms
deformation pattern defined by one or more shear of localization of deformations.
bands (Rice, 1976). Such localized modes of defor-
mation give rise to a number of complications related
to mesh dependence, internal length scales, etc.
Secondly, and more seriously, it has frequently been 2 FRICTION AND PLASTICITY
reported that numerical solutions to boundary value
problems involving nonassociated constitutive mod- The strength of cohesive-frictional materials can in
els are much more difficult to obtain than in the case many cases be described via a yield criterion of the
where the flow rule is associated (Manzari and Nour, type
2000; Carter et al., 2005; Clausen and Krabbenhoft,
2008; Loukidis and Salgado, 2009). Even recogniz-
ing the above complications with shear banding, mesh
dependence, and so on, it is somewhat surprising that
there should be any significant problems obtaining where F is the yield function, p and q are appropriate
solutions, be they physically relevant or not. measures of mean and deviatoric stress respectively,
As an alternative to elastic-perfectly plastic models M is a friction parameter and k is a measure of the
where a nonassociated flow rule usually is a neces- cohesion of the material. In the following, compressive
sity, hardening models are often considered. Using stresses are taken as being positive.
51
p p
Let εv and εs be the plastic volumetric and shear
strains conjugate to p and q respectively. The associ-
ated flow rule then reads:
where λ̇ ≥ 0 is a plastic multiplier and a superposed dot

denotes rate quantities. The associated flow rule thus
p p
predicts a dilation given by d = −ε̇v /ε̇s = M . For most
materials, this dilation is excessive. Consequently, one Figure 1. Illustration of the microscopic origins of friction
defines a flow potential: as plastic shearing of asperities. Coulomb friction implies that
the ultimate strength of the assembly is proportional to the
macroscopically applied pressure p. The plastic shearing is
assumed to be of the ductile, purely cohesive kind. For brittle
materials such as sand grains, this assumption is justified by
such that the plastic strain rates are given by
the very high stress level at the scale of the asperities which
effectively renders the otherwise brittle material ductile.
p p
As such, a constant dilation d = −ε̇v /ε̇s = N ≤ M may
be accounted for.
The internal dissipation for a material with yield
function F and flow potential G is given by
Figure 2. Explicit evaluation of apparent cohesion: original

and approximate yield functions. The area indicated by (a)
is non-permissible according to the original yield function
but permissible according to the approximate yield function.
This expression for the internal dissipation reveals sev- Similarly, the area indicated by (b) is non-permissible accord-
eral interesting features. Firstly, for N = M , i.e. for an ing to the approximate yield function but is within the elastic
associated flow rule, the internal dissipation is propor- domain according to the original yield function.
tional to the cohesion k. Thus, for a purely frictional p p
material (k = 0), the internal dissipation is zero which desired result, namely a dilation d = −ε̇v /ε̇s = N . In
clearly is physically problematic. Secondly, for N < M the solution of boundary value problems, the apparent
the dissipation is proportional to an apparent cohe- cohesion is of course not known a priori as it is directly
sion that comprises two terms: the internal cohesion proportional to the pressure that is to be determined
k and an contribution (M − N )p that stems from the as part of the solution. However, assuming that such
prescribed nonassociativity. The interpretation of the problems are solved incrementally via a sequence of
latter term as an apparent cohesion is consistent with pseudo-time steps, some parts of the yield function
the viewpoint that friction results from mechanical may, in principle, be evaluated implicitly while other
interaction of microscopic asperities on the surfaces parts may be evaluated explicitly. Assume that the state
of the solids in contact (Bowden and Tabor, 1973). at time tn is known. The yield condition imposed at tn+1
With the stresses at the scale of the asperities being may then be approximated as
much greater than the elastic limit of the material,
it is primarily plastic deformations at the microscale
that govern the macroscopically observed frictional
resistance (see Figure 1). This interpretation motivates where
rewriting the yield function (1) as
Again, this produces the desired result that the associ-

where ated flow rule predicts that the dilation at time tn+1
is equal to N . However, the explicit evaluation of
the apparent cohesion means that the original yield
function may be exceeded for the new stress state at
is the apparent, pressure dependent, cohesion. tn+1 . Similarly, the approximation may imply plas-
Suppose now that the apparent cohesion, k̂, is tic yielding for stress states that would otherwise be
known. The associated flow rule then produces the deemed purely elastic (see Figure 2). However, for
52
small enough increments, i.e. for tn+1 − tn → 0, it methods may require slight modification to cater for
would be expected that the error introduced by the initial stress states that do not satisfy F ≤ 0, but oth-
explicit evaluation of the apparent cohesion would tend erwise all operations would be identical to those of
to vanish. Numerical experiments largely confirm this standard associated elastoplasticity.
supposition.
3 EXAMPLES
2.1 Variational formulation
With the explicit evaluation of the apparent cohe- In the following, the new computational scheme is
sion detailed above, the governing equations that tested and the effects of nonassociativity in general
define each time increment are essentially those of are examined. In all cases, a simple linear elastic-
standard associated plasticity. As such, a variational perfectly plastic model is used. The yield function is
formulation is straightforward. of the Drucker-Prager type:
Following Krabbenhoft et al. (2007b, 2009) and
assuming linear elasticity, the relevant time-discrete
problem can be written as:
We further consider a nonassociated flow rule defined

e e
where Cpp and Cqq are the elastic constants and by the plastic potential G:
kn = k + (M − N )pn . The strain increments εv and
εs are assumed given.
The optimality conditions associated with this opti-
Assuming associated flow and plane strain conditions,
mization problem are given by
the parameters M and k are given in terms of the Mohr-
Coulomb friction angle, φ, and cohesion, c, by:
Analogously, the dilation angle, ψ, may be related to

where λn+1 is a Lagrange multiplier (the plastic the Drucker-Prager parameter N by:
multiplier).
Following Krabbenhoft et al. (2007b), the local
problem (10) may be extended to the entire domain.
After appropriate finite element discretization, the
final discrete optimization problem to be solved in Following the discussion in the previous section, the
each time step is given by flow rule is imposed by evaluating the apparent cohe-
sion, k̂, explicitly. In this way, only the yield potential,
F, in the form given above, is actually needed.
3.1 Biaxial test

The first problem considered is the biaxial test
where σ are the stresses, C e is the elastic compli- sketched in Figure 3(a). Assuming associated flow, a
ance matrix, BT is an equilibrium matrix, and Fj∗ are homogeneous state of stress, and neglecting the imper-
yield functions to be enforced at the Nσ stress points. fection indicated in the figure, it is straightforward to
This problem may be solved using either general show that the ultimate compressive stress is given by
interior-point solvers (Krabbenhoft and Damkilde,
2003; Krabbenhoft et al., 2007b) or by means of more
specialized formulations. For yield surfaces of the
Drucker-Prager type, second-order cone programming
is particularly suited while Mohr-Coulomb type con- Next, assuming a nonassociated flow rule but still a
straints can be handled efficiently using semidefinite homogeneous state of stress, the ultimate compressive
programming formulations (see Krabbenhoft et al., stress can be determined as
2007a, 2008, for details).
Alternatively, more conventional methods of com-
putational elastoplasticity are also applicable. These
53
Figure 4. Load-displacement curves for biaxial test using
200 equal size displacement increments.
cases, a mesh involving some 25,000 displacement

Figure 3. Biaxial test: problem setup (a) and localized
degrees-of-freedom was used. Further details are as
deformation solution (b).
follows: friction angle φ = 35◦ (M = 0.94), dilation
where angle ψ = 5◦ (M = 0.15), cohesion c = 1 (k = 1.35),
bulk modulus K = 1, shear modulus G = 0.5. All
problems were solved using the second-order cone
programming solver MOSEK (Andersen et al., 2003).
Figure 4 shows the load-displacement curves for the
three cases discussed above: associated flow, nonas-
sociated flow with homogeneous deformation, and
This solution is derived by requiring that the out-of- nonassociated flow with localized deformation. The
p
plane plastic strain rate, ε̇22 = λ̇∂G/∂σ22 , be equal to response in the latter case is computed using the afore-
zero at the ultimate limit state. mentioned finite element mesh.To trigger localization,
Finally, assuming nonassociated flow and consid- an imperfection in the form of an element with a
ering a localized state of deformation where the block friction angle equal to φimp = 20◦ is introduced as indi-
is traversed by a single band of intense deformation cated in Figure 3(b). The deformations at incipient
[see Figure 3(b)], the ultimate compressive stress can collapse resulting from this computation are shown
be determined as Figure 5. As seen, the collapse mechanism is highly
localized with a single band traversing the sample,
largely following the diagonals of the mesh at an angle
of 53.1◦ with the horizontal.
It is interesting to compare this response with the
where predictions of a bifurcation analysis. Following Leroy
and Ortiz (1989), such an analysis is carried by initially
assuming that the stresses and the deformations are
homogeneous. The vertical stress is increased incre-
mentally and the eigenvalues of the acoustic tensor
are gauged at each time instant. Assuming plane strain
conditions, the acoustic tensor (or matrix) is given by
This solution is derived by observing that the two
normal strains lying in the plane of the shear band
approach zero as the thickness of the band approaches
zero (Krabbenhoft et al., 2004). The inclination of with Dep being the usual elastoplastic tangent matrix
the shear band corresponding to the above solution and
is given by
where θ is the shear band inclination angle shown

The finite element analyses are carried out using in Figure 3(b). A zero eigenvalue of the acoustic
a mesh consisting of quadratic displacement trian- tensor indicates that a switch from a homogeneous
gles, arranged as indicated in Figure 3(a). In all to a localized state of deformation is possible. This
54
Figure 6. Nγ problem: setup and finite element mesh (1,655
elements, 6,904 displacement degrees-of-freedom).
Figure 5. Biaxial test: deformations at incipient collapse.
first occurs as indicated in Figure 4 at an incli-

nation angle of θ = θcr = 51.9◦ . As seen from the depends on the friction angle. The determination of
figure, this bifurcation point is in good agreement with this bearing capacity factor has been the subject of a
the bifurcation point of the finite element analysis rather large number of investigations. These are sum-
(the bifurcation here occurs slightly earlier due to marized by Hjiaj et al. (2005b) who also, based on
the imperfection). Furthermore, the critical angle of computational upper and lower bound limit analysis,
the bifurcation analysis is in good agreement with the propose the following closed-form expression for the
one observed in the finite element analysis, namely bearing capacity factor for a rough footing:
θ 53◦ . This angle, however, is slightly different
from the angle θloc = 45◦ + 12 φloc = 60.0◦ for which
the lowest ultimate limit load is attained. Hence, the
computed ultimate limit load is slightly higher.
Regarding the post-bifurcation response, the load- Following the approach in the previous section of
displacement curve shown in Figure 4 appears to be deriving effective material parameters for problems
somewhat oscillatory. We view these oscillations as of nonassociated plasticity, effective bearing capacity
being the result of that there is more than one solution factors (Nγ )hom and (Nγ )loc can be defined. These are
– after the bifurcation point, there is in fact an infi- computed by replacing φ in the above equation with
nite number of possible localized deformations modes, the effective friction angles, given in terms of Mhom
each with different load-displacement characteristics. and Mloc by (17) and (19) respectively.
Interestingly, even though they may well be numerical Regarding the finite element analysis, the problem
artifacts, similar oscillations are frequently observed is solved using two different types of finite ele-
in biaxial tests on sand (Desrues and Vigiani, 2004; ments: a standard quadratic displacement triangle and
Gajo et al., 2004). a mixed linear stress-quadratic displacement element
Finally, it should be noted that the violation of first proposed by Borges et al. (1996) for limit anal-
the yield criterion that is possible with the new solu- ysis and subsequently extended to elastoplasticity by
tion approach in all cases is moderate and appears Krabbenhoft et al. (2007b). Further details are as fol-
to be systematically reduced as the magnitude of the lows: friction angle φ = 40◦ (M = 0.94), dilation angle
displacement increment is reduced. ψ = 10◦ (M = 0.15), cohesion c = k = 0, unit weight
γ = 16, bulk modulus K = 50 × 103 , shear modulus
3.2 Nγ problem G = 25 × 103 .
The results of the analysis are shown in Figure 7. For
The second example concerns the load-deformation the associated case, the displacement finite element
behaviour of a centrally loaded strip footing on a purely overestimates the bearing capacity while the mixed
frictional soil (see Figure 6). The bearing capacity of element provides a slight underestimate. Qualitatively,
such a footing is usually expressed as the results of the nonassociated analysis follow those of
the previous example. The ultimate strength decreases,
though not to the lowest level theoretically possible.
Also, the load-displacement behaviour is slightly oscil-
where Vc is the collapse load, B is the width of the latory – despite the yield function being satisfied to
footing, and Nγ is the bearing capacity factor which within a rather tight tolerance.
55
Gajo, A., Bigoni, D., and Wood, D. M. (2004). Multi-
ple shear band development and related instabilities in
granular materials. Journal of the Mechanics and Physics
of Solids, 52:2683–2724.
Hjiaj, M., Huang, W., Krabbenhoft, K., and Sloan, S. W.
(2005a). Formulation of non-standard dissipative behav-
ior of geomaterials. Journal of Engineering Mathematics,
52:147–165.
Hjiaj, M., Lyamin, A. V., and Sloan, S. W. (2005b). Numeri-
cal limit analysis solutions for the bearing capacity factor
Nγ . International Journal of Solids and Structures, 42:
1681–1704.
Krabbenhoft, K. (2009). A variational principle of elasto-
plasticity and its application to the modeling of frictional
materials. International Journal of Solids and Strcutures,
46:464–479.
Krabbenhoft, K. and Damkilde, L. (2003). A general nonlin-
ear optimization algorithm for lower bound limit anal-
Figure 7. Load-displacement curves for Nγ problem with ysis. International Journal for Numerical Methods in
φ = 40◦ and ψ = 10◦ using 100 equal size displacement Engineering, 56:165–184.
increments. Krabbenhoft, K., Lyamin, A. V., Hjiaj, M., and Sloan, S. W.
(2004). Limit analysis of materials with nonassociated
flow rules. In Proc. Eccomas 2004, pages 1–21.
4 CONCLUSIONS Krabbenhoft, K., Lyamin, A. V., Hjiaj, M., and Sloan, S. W.
(2005). A new discontinuous upper bound limit anal-
A new computational plasticity scheme for nonasso- ysis formulation. International Journal for Numerical
ciated frictional materials has been presented. While Methods in Engineering, 63:1069–1088.
general, it relies solely on well established concepts of Krabbenhoft, K., Lyamin,A.V., and Sloan, S. W. (2007a). For-
associated plasticity. As such, a number of equally well mulation and solution of some plasticity problems as conic
established numerical methods are directly applicable. programs. International Jounal of Solids and Structures,
Moreover, techniques that have proved very powerful 44:1533–1549.
Krabbenhoft, K., Lyamin, A. V., and Sloan, S. W. (2008).
for limit analysis can be easily extended to general Three-dimensional Mohr-Coulomb limit analysis using
elastoplastic problems. These include the incorpo- semidefinite programming. Communications in Numer-
ration of discontinuous stress and velocity fields ical Methods in Engineering, 24:1107–1119.
(Krabbenhoft et al., 2005; Lyamin et al., 2005a,b) as Krabbenhoft, K., Lyamin, A. V., Sloan, S. W., and Wrig-
well as application of efficient and robust methods gers, P. (2007b). An interior-point method for elastoplas-
of nonlinear optimization such as the one used in the ticity. International Journal for Numerical Methods in
present study. Engineering, 69:592–626.
Leroy, Y. and Ortiz, M. (1989). Finite element analysis
of strain localization in frictional materials. Interna-
REFERENCES tional Journal for Numerical and Analytical Methods in
Geomechanics, 13:53–74.
Andersen, E. D., Roos, C., and Terlaky, T. (2003). On Loukidis, D. and Salgado, R. (2009). Bearing capacity of
implementing a primal-dual interior–point method for strip and circular footings in sand using finite elements.
conic quadratic optimization. Mathematical Program- Computers and Geotechnics, 36:871–879.
ming, 95:249–277. Lyamin,A., Krabbenhoft, K.,Abbo,A., and Sloan, S. (2005a).
Bigoni, D. and Hueckel, T. (1991). Uniqueness and local- General approach to modelling discontinuities in limit
ization – I. Associative and non-associative elastoplas- analysis. In Barla, G. and Barla, M., editors, Proceedings
ticity. International Journal of Solids and Structures, of IACMAG, Turin.
2:197–213. Lyamin, A. V., Sloan, S. W., Krabbenhoft, K., and Hjiaj,
Borges, L. A., Zouain, N., and Huespe, A. E. (1996). A non- M. (2005b). Lower bound limit analysis with adaptive
linear optimization procedure for limit analysis. European remeshing. International Journal for Numerical Methods
Journal of Mechanics, A/Solids, 15(3):487–512. in Engineering, 63:1961–1974.
Bowden, F. P. and Tabor, D. (1973). Friction. An Introduction Manzari, M. T. and Nour, M. A. (2000). Significance of soil
to Tribology. Anchor Press/Doubleday. dilatancy in slope stability analysis. Journal of Geotech-
Carter, J. P., Poon, M. S. B., and Airey, D. W. (2005). Numer- nical and Geoenvironnemental Engineering, 126:75–80.
ical and semi-analytical techniques for footings subjected Rice, J. R. (1976). The localization of plastic deformation.
to combined loading. In Proc IACMAG 11, Turin, pages In Koiter, W., editor, Theoretical and Applied Mechanics,
163–176. pages 239–264. North-Holland.
Clausen, J. and Krabbenhoft, K. (2008). Existence and Simo, J. C. (1998). Numerical analysis and simulation in plas-
uniqueness of solutions in nonassociated mohr-coulomb ticity. In Ciarlet, P. G. and Lions, J. L., editors, Handbook
elastoplasticity. In Proc. WCCM VIII, Venice. of Numerical Analysis, pages 179–499. Elsevier.
Desrues, J. and Viggiani, G. (2004). Strain localization in Souza de Neto, E. A., Peric, D., and Owen, D. J. R. (2009).
sand: an overview of the experimental results obtained Computational Methods for Plasticity: Theory and Appli-
in Grenoble using stereophotogrammetry. International cations. Elsevier.
Journal for Numerical and Analytical Methods in Geome-
chanics, 28:279–321.
56
Comparison of methods for calculation of settlements of soft clay
H.P. Jostad
Norwegian Geotechnical Institute (NGI), Oslo, Norway
Norwegian University of Science and Technology (NTNU), Trondheim, Norway
S.A. Degago
Norwegian University of Science and Technology (NTNU), Trondheim, Norway
ABSTRACT: Calculation of long-term settlements of soft clay generally consists of many uncertainties. By
studying back-calculated field cases, from literatures, it is therefore generally very difficult to compare the
performance of different calculation tools due to varying interpretations and assumptions in the governing input
parameters. Therefore, as part of a series of creep workshops called CREBS, the participants were invited to
analyse a set of hypothetical cases using their material models and computer program. The cases involved a
30 m thick homogeneous normally consolidated soft clay layer underlying a 10 m thick sand layer subjected to a
surface stress of either 50 or 90 kPa. Six groups submitted their contributions to this exercise. This paper presents
the main results from this exercise, compare the background of the different material models and discuss the
reasons for the characteristic differences in the obtained results.
1 INTRODUCTION 2 HYPOTHETICAL CASES
The first CREBS (CREep Behaviour of Soft clay) All cases consist of a 30 m thick soft clay layer below a
workshop was held in January 2006, at NGI in Oslo. 10 m thick sand layer as shown in Figure 1. A surface
One of the conclusions from this workshop was that load of 50 kPa (light structure) or 90 kPa (heavy struc-
even for material models based on the same frame- ture) is distributed over an area that is large compared
work it is very difficult to compare the differences to the thickness of the soft clay layer (1D condition),
in assumptions and input data since all models use except for one case. Settlements below the clay layer
somewhat different expressions. Hence, it was recom- are neglected. The ground water table (GWT) is at the
mended to establish a common set of definitions and to top of the sand layer. For the bottom boundary two
systematically compare existing calculation tools used extreme assumptions are considered, either a perfectly
in long-term settlement analyses of soft soils. drained or an impervious surface. The following cases
The participants at the second CREBS workshop, were analysed (however, only some of the results are
held in September 2007 in Pisa (Italy), were invited presented here):
to analyse a set of well defined hypothetical cases
by various calculation tools. The main purpose was 1. Normally consolidated (NC) behaviour where the
to compare variations in interpreted input data and pre-consolidation pressure is assumed to be equal
obtained calculation results and not a competition to the in-situ effective vertical stress.
in predicting the most correct results. When study- 2. Normally consolidated behaviour with an appar-
ing published back-calculations of field cases (see ent pre-consolidation pressure corresponding to
for instance Leroueil 2006) large differences may be a constant over-consolidation ratio of OCR =
obtained due to uncertainties in material properties, 1.4.
in situ pore pressure distribution, drainage conditions 3. A time history where the soil profile is pre-loaded
and earlier load histories. in a period of 25 years before increasing the load.
The results from the analyses of the hypothetical 4. The clay layer is divided into two sub-layers with
cases were briefly presented at the third CREBS work- significantly different permeabilities.
shop held in July 2009, in Gothenburg (Sweden). This 5. A load is applied on a strip foundation with limited
paper gives a more detailed presentation and evalua- width of 20 m, that gives a decreasing excess stress
tion of some of the most characteristic results from distribution with depth and induce some effect of
this exercise. shear mobilisation.
57
Figure 1. The hypothetical cases.
2.1 Soil conditions Figure 2. Stress, strain and time relationships obtained from
a standard IL-test.
2.1.1 Sand layer
In order to make it easier to compare the results,
the main properties of the sand layer were directly
given: Constrained modulus, M = 10 MPa, submerged
unit weight of soil, γ = 10 kN/m3 and permeability,
3 PARTICIPANTS
k = 1 m/year.
The following participants have analyzed the given
cases:
2.1.2 Clay layer
The soft soil layer consists of a homogeneous, nor- • Dr. Martino Leoni and Professor Pieter Vermeer
mally consolidated, fully water saturated, plastic from the University of Stuttgart. They used the com-
marine clay with approximately the same age (10,000 puter program Plaxis (www.plaxis.nl) with the Soft
years). This means that the characteristic mechanical Soil Creep (SSC) (Vermeer & Neher 1999) and the
behaviour found at one depth is assumed to be valid user defined Anisotropic Creep model (Leoni et al.
for the entire depth of the layer. 2008)
The constitutive behaviour of the clay is found from • Dr. Zhen-Yu Yin and Dr. Minna Karstunen from
a standard oedometer test with incremental loading Ecole Centrale de Nantes and University of Strath-
(IL). The results from the different load increments clyde. They used Plaxis with the user defined visco-
are shown in Figure 2. The figure shows the verti- plastic EVP-SCLAY1S model (Yin & Karstunen
cal strain increment εv = δv /ho , where δv is the 2008)
vertical displacement at the top of the sample dur- • Dr. David Nash from the University of Bristol.
ing the actual load increment and ho = 20 mm is the He used the computer program BRISCON with an
initial sample height. Most of the load increments isotache based model (Nash & Ryde 2001)
have a period of about 1440 minutes (1 day). How- • Mats Olsson and Professor Claes Alén from
ever, for load increment (180 − 280 kPa) the vertical Chalmers University of Technology. They used
stress of 280 kPa was kept for a period of about 5.5 the GeoSuite Settlement Program (www.novapoint.
days (8000 minutes). Figure 2 shows the accumulated com) with the Chalmers model Claesson
vertical strain (εv = δv /ho ) after 24 hours for all load (2003)
increments. • Per-Evert Bengtsson and Rolf Larsson from the
Swedish Geotechnical Institute (SGI). They used
The initial effective vertical stress σvo was pur-
posely not provided for the actual test. The reason the settlement program EMBANKCO with a model
is that because the actual over-consolidation ratio rather similar to the Chalmers model (Bengtsson &

OCR = σvc / σvo was specified. The results should Larsson 1997)
therefore be representative for the given OCR and not • Professor Hans Petter Jostad from NGI and Nor-
affected by the interpretation of the effective vertical wegian University of Science and Technology
(NTNU). He used the GeoSuite Settlement Pro-
pre-consolidation pressure σvc .
For models based on void ratio, the initial void ratio gram with the Krykon material model (Svanø et al.
eo is 1.17. 1991)
58
4 BRIEF DESCRIPTIONS OF MODELS USED in the finite difference program Embankco (Bengts-
son & Larsson 1997).
In order to systematically compare the different mod- The 24 hr reference strain, due to stress
els used, their behaviours in uniaxial vertical strain changes, is given by the constrained modulus,
condition (1D) are briefly presented.The vertical strain Mt = Moc for σ v < σ vc ; Mt = ML for σ vc < σ v < σ vL ;
rate is then decomposed into a component due to effec- and Mt = ML + m · (σ v − σ vL ) for σ vL < σ v ; where
tive vertical stress changes (a reference strain) and a Moc = 50 · σ vc ; ML = 3140 + m · (σ vc − 192 kPa);
component only due to time (creep): σ vL − σ vc = 78 kPa; and m = 16.5. The time resistance
is given as: R = Ro + r · t, where Ro = r · to , to = 24 hr;
r is 261 at σ vc and increases asymptotically to an
infinite value at OCR = 1.25. In the NC-regime, r
increases from 261 to about 365 at 10 % vertical strain.
where Mt is an effective stress dependent tangential
constrained or oedometer modulus; R = Ro + r · t is
Janbu’s time resistance (Janbu 1969); Ro is the initial 4.4 The Briscon model
time resistance; and r is the time resistance num-
ber. Lines with constant R-values in vertical effective The material model used in the finite difference pro-
stress–strain space are called isotaches (Šuklje 1957). gram BRISCON (Nash & Ryde 2001) is also based on
The values for the model parameters are presented the isotache concept.
based on the results submitted by the participants. The 24 hr reference strain, due to stress changes, is
given by: Mt = mr · σ v for σ v < σ vc , and Mt = m · σ v
for σ v > σ vc , with m = 13.7 and mr = 8 · m = 110.
4.1 The Krykon model The time resistance is given as R = Ro + r · teq ,
The Krykon model, implemented in the GeoSuite Set- where Ro = 0.95 year. The equivalent time teq is the
tlement program, is based on Janbu’s time resistance time required to obtain the increase in creep strain from
concept. A detailed description of the model can be the reference time line (RTL) to the current strain at
found in Svanø et al. (1991). the actual effective vertical stress. The RTL in a εv −
The reference strain after 24 hr εo due to effective ln(σ v ) plot is a line that goes through the strain at σ vc
stress changes, is given by a stress dependent tangen- by a slope defined by the modified compression index
tial constrained modulus: Mt = Moc for σv < σvc
and λ∗ = 1/m. The teq is then calculated from the equivalent
Mt = m · (σv − σvr

) for σv > σvc

(where σvr is reference creep strain εeq :
stress), with Moc = 5 · m · σ vc ; m = 16 and σ vr = 0 in
this case. The time resistance R is given as function of
the vertical strain ε:
In the over-consolidated (OC) regime the initial

where Ro = 0.8 year is the time resistance at the ref- εeq,o is defined by the OCR (σ vc /σ v ) and the modulus
erence strain εo , and r is an effective vertical stress numbers m and mr :
dependent time resistance number varying linearly
between r = 1125 at σ vo and r = 300 at σ vc . In the
NC-regime r is taken to be constant equal to 300.
4.2 The Chalmers model The time resistance numbers for the load steps in
the IL test were found to be: r = 4000, 3000, 2000,
The Chalmers model, implemented in the GeoSuite
1200, 700, 280, 400, 360 and 410. In addition, it is
Settlement program, is also based on Janbu’s time
realistic to assume that some of the creep observed in
resistance concept. A detailed description of the model
the OC-regime is due to sample disturbance and that
can be found in Claesson (2003).
the in-situ creep rate for a 10,000 year old clay is negli-
The 24 hr reference strain, due to stress changes,
gible. The pre-consolidation pressure was found to be
in this case is given by an initial stress depen-
σ vc = 152 kPa. Based on this, three different variations
dent tangential constrained modulus, Mt = 12 MPa +
of the r-value were considered: r = 343 (constant); r
0.5 MPa · (z − 10 m) for σ v < σ vc and Mt = 13.5 ·
as function of σ v /σ vc based on measured results; and a
σ vo for σ v > σ vc . The time resistance is given as:
case where r is gradually increased in the OC-regime
R = Ro + r · t, where Ro = r · to ; to = 24 hr and r is
to an unlimited value at OCR = 1.4 (only constant r is
varying linearly between r = 10, 000 at σ vo and r = 300
presented here).
at σ vc . In the NC-regime r = 300.
4.3 The Embankco model 4.5 The Soft Soil and Anisotropic Creep models
The Embankco material model is very similar to the The Soft Soil Creep model (Vermeer & Neher 1999)
Chalmers model; however, the model is implemented in Plaxis, is similar to the BRISCON model; however,
59
extended to a full 3D stress condition using the frame- is taken to be 3.6 at the NC reference line. The cor-
work of the modified Cam-Clay model. The time resis- responding r-value at the NC-line is then 412. The
tance is defined as R = Ro + r · teq , where Ro = r · to , modulus numbers were taken equal to m = 13.7 and
to = 24 hrs, r = 1/µ∗ = 333, and µ∗ is the modified mr = 13.5 · 13.7 = 185.
creep index used as input in Plaxis.
The volumetric creep strain is then related to the
expansion of the ellipse in the effective mean stress 4.7 Discussions
(p ) – deviatoric stress (q) space controlled by the mod-
All models may give approximately the same 24 hr ref-
ified compression index, λ∗ = 1/m = 1/13.7 = 0.073.
erence strain and the time dependent strain. The actual
This means that the equivalent time teq and the corre-
results are therefore dependent on how the participants
sponding creep strain is governed by the expansions
interpreted the IL test.
of the ellipse compared to the ellipse given by the cur-
The modulus number used in the NC-regime did not
rent stress state (p and q). This gives the following
differ much since almost all participants based it on
expression for the time resistance:
the slope of the εv versus log(σ v ) plot at large effec-
tive vertical stresses. The interpreted creep strain in
the NC-regime also did not differ much. Most of the
participants found the time resistance factor r from
the last part of the 8000 minutes creep phase at an
effective vertical stress above the pre-consolidation
where σ vy is the updated apparent pre-consolidation pressure. However, the parameters (N and µ) used in
stress due to creep. The elastic effective stress depen- EVP-SCLAY1S were selected in order to fit the strain
dent constrained modulus is given as: for all load steps in the IL test. Consequently the model
underestimated the strain during the 8000 minutes long
creep phase at 280 kPa.
The time resistance number in the NC-regime varies
between 261 (Embankco) and 412 (EVP-SCLAY1S) at
where νur = 0.2 is the unloading/reloading Poisson’s the 24 hr reference time. Furthermore, for the N -value
ratio; Ko is the actual effective horizontal/vertical used in EVP-SCLAY1S, the r-value increases with
stress ratio; and mr = 7.85 · m = 108. increasing strain under a constant effective vertical
For the Anisotropic Creep Model (Leoni et al. stress.
2008), a rotated ellipse based on a fabric tensor The largest differences are found in the mod-
(Wheeler et al. 2003) is used. However, ACM gives elling of the creep in the OC-regime. In Embankco
the same results as the SSC model except for Case 5. the r-value asymptotically increases with increasing
OCR to infinitely at OCR = 1.25. In the Chalmers
modelr-value increases with increasing OCR to a very
4.6 The EVP SCLAY-1S model large value (r = 10,000) at OCR = 1.4. In Krykon the
In differ to the other models that are based on the corresponding value at OCR = 1.4, is r = 1125. In
isotache concept, the anisotropic elasto-viscoplastic SSC/ACM and Briscon the r-value is independent of
model EVP SCLAY-1S (Yin & Karstunen 2008) is OCR. Instead it is the equivalent time teq that increases
based on the overstress theory (Perzyna 1966) and a with OCR (see Eq. 3) which gives an increase in the ini-
rotated Cam-Clay surface as in ACM. In this case the tial time resistance Ro = r · (to + teq ). However, based
time resistance is give by a somewhat more complex on the IL test all the r-values used in all the analyses
expression: may be considered as reasonable.
The modulus used in the OC-regime depends on
whether it was based on the initial loading from the
in-situ effective vertical stress to the pre-consolidation
stresses (which may underestimate the stiffness due
where µ is the fluidity parameter; N is the strain-rate to sample disturbance), taken from the unloading
coefficient relating to the strain-rate effect on shear sequence at the end of the test (starting from a large
strength and pre-consolidation stress; df 1 ∼ 0.7 for 1D effective vertical stress) or based on in-house experi-
condition; and OCRs is the ratio between the size of the ences. For instance in-house experiences were used for
ellipse given by the current stress state (dynamic load- the Chalmers model to extrapolate to larger effective
ing surface) and the size of an inner expanding ellipse stresses.
(static yield surface).The expansion of the inner ellipse Figure 3 shows the stress dependent 24 hr refer-
is controlled by the accumulated creep strain as for the ence strain and the time dependent strain at the top
SSC/ACM. From the above expression it is seen that and bottom (with drainage boundary) of the clay layer
the time resistance Ro at the 24 hr reference strain (OCR = 1.4) given by the different models. The curves
is controlled by OCRs – value and that the creep rate are established based on the reported input parameters.
vanish when OCRs = 1. From these plots it is clear that the calculated
By fitting the IL test the constants were found to settlements for Case 1 to 5 will be smallest by
be µ = 5 · 10−16 (1/year) and N = 13.77. The OCRs EVP-SCLAY1S and largest by SSC and Briscon. The
60
Figure 4. Calculated settlements versus time for Case 1 with
open bottom boundary.
Figure 5. Calculated settlements vs. time for Case 2 with

open bottom boundary and q = 50 kPa.
Figure 6. Excess pore pressure histories for Case 2 with

q = 50 kPa and open bottom boundary.
Figure 3. Effective vertical stress-strain-time relationships
derived from the different models at the top and the bottom r = 261. The settlements due to creep for r = 300 is
with open boundary of the clay layer. The time dependent about 1 m after 100 years.
strains are for an excess load of 50 kPa. Figure 5 also shows the calculated settlements for
Case 2 with q = 50 kPa and open bottom boundary.
results obtained by SSC and Briscon should be very The results demonstrate the effect of different assump-
similar. tions of the constitutive behaviour in the OC-regime.
Chalmers and EVP-SCLAY1S models give the small-
est settlements (0.3 m) due to the smallest creep in
5 RESULTS the OC-regime, while Krykon, Briscon and SSC give
roughly the same settlements (an average of 0.66 m).
Only the main results that demonstrate some charac- Embankco gives the largest settlements (0.89 m) and it
teristic differences are selected and presented here. seems that the solution has become somewhat unstable
Figure 4 shows the calculated settlements versus after about 60 years.
time for Case 1. The plot demonstrates the effect of From Figure 6, which shows the corresponding
different time resistance numbers used, where EVP- excess pore pressure in the middle of the soft clay layer
SCLAY1S gives the smallest settlements (0.95 m) due versus time, it is seen that the primary consolidation
to a time resistance number r of more than 400, while phase is finalized after 100 years for all cases (except
Embancko gives the largest settlements (1.62 m) with for Embankco). The differences in calculated pore
61
Furthermore, it would have been of large benefit
if all creep models have used a common set of main
input parameters. It would then be easier to understand
differences in obtained results by comparison of input
parameters and to establishing a common data base for
creep behaviour of soft clays.
ACKNOWLEDGEMENTS
Figure 7. Calculated strain profiles after 100 years for The participants of the exercise are greatly acknowl-
Case 2 with q = 50 kPa and open bottom boundary. edged for their valuable contribution and allowing the
authors for publishing the results.
pressure during the first 20 years are significant, e.g.
after 10 years the excess pore pressure varies between
7 kPa or 14% (Chalmers) and 37 kPa or 74% (Krykon). REFERENCES
It is also seen that the excess pore pressure for SSC and
Briscon initially becomes larger than 50 kPa. This indi- Bengtsson, P-E. & Larsson, R. 1997. Calculation of settle-
cates that the initial creep strain rate in these models ments for embankments on fine-grade soils. Calculation
is unrealistically large at this depth. The creep rate is of course of settlements with time. In User’s guide for
Embankco programme version 1.02. Swedish Geotechni-
also too large in Krykon, however, in GeoSuite Settle-
cal Institute, Linköping.
ment the creep strain is not allowed to reduce the initial Claesson, P. 2003. Long term settlements in soft clays. Ph.D.
effective vertical stresses. The excess pore pressure for thesis, Chalmers University of Technology, Sweden.2
EVP-SCLAY1S and Chalmers decreases very rapidly Janbu, N. 1996. The resistance concept applied to defor-
due to stiffer behaviour and less creep. mations of soils. Proc. 7th Int. Conf. Soil Mech. Found.
Figure 7 shows the calculated total strain profiles Engng, Mexico. 1: 191–196.
after 100 years, also for Case 2 with q = 50 kPa and Leoni M., Karstunen M. & Vermeer P.A. 2008. Anisotropic
open bottom boundary. From this plot it is seen that creep model for soft soils. Géotechnique 58(3): 215–226.
the strain at the top and bottom of the clay layer gen- Leroueil, S. 2006. Šuklje Memorial Lecture: The isotache
approach. Where are we 50 years after its development by
erally agree with the strains-time relationship given
Professor Šuklje? 13th Danube-European Conf. Geotech.
in Figure 3. The differences in the obtained results Engng, Ljubljana, Slovenia. 2: 55–88.
are directly results of the differences in the inter- Nash, D.F.T. & Ryde, S.J. 2001. Modelling the consolida-
preted input data. However, the strain at the top of tion of compressible soils subject to creep around vertical
the clay layer is slightly larger for Embankco and drains. Géotechnique 51(4): 257–273.
slightly smaller for EVP-SCLAY1S than found from Perzyna, P. 1966. Fundamental problems in viscoplasticity.
the curves in Figure 3. Advanced Applied Mechanics 9: 244–377.
Šuklje, L. 1957. The analysis of the consolidation process
by the isotaches method. Proc. 4th Int. Conf. Soil Mech.
Found. Engng., London. 1: 200–206.
6 CONCLUSIONS
Svanø, G., Christensen, S., and Nordal, S. 1991. A soil model
for consolidation and creep. Proc.10th Int. Conf. Soil
The main conclusion from this study is that the dif- Mech. Found. Engng, Florence, Italy. 1:269–272.
ferences in calculated settlements for the set of well Vermeer, P. A. & Neher, H. P. 1999. A soft soil model that
defined idealized hypothetical cases are rather large. accounts for creep. In R.B.J. Brinkgreve (ed.), Proc. Int.
The main reason for the differences is uncertainties Symp. Beyond 2000 in Comput. Geotech.: 10 Years of
and assumptions in the creep behaviour for stress Plaxis International: 249–261. Rotterdam: Balkema.
conditions below the initial pre-consolidation stress Wheeler, S.J., Näätänen A., Karstunen, M. & Lojander, M.
(OC-regime). The differences could have been even 2003. An anisotropic elastoplastic model for soft clays.
Canadian Geotechnical Journal 40(2): 403–418.
larger for a real case where uncertainties related to the
Yin, Z.Y. & Karstunen, M. 2008.Influence of anisotropy,
OCR profile are generally significant. destructuration and viscosity on the behaviour of an
It is therefore recommended to continue the focus embankment on soft clay. In: Singh, D. N. (ed.): Proc.
on the constitutive behaviour in the OC-regime, to 12th Int. Assoc. Comput. Methods Advances Geomech.
find suitable testing procedures and interpretation (IACMAG), Goa, India: 4728–4735.
techniques that can account for sample disturbance.
62
Effect of yield surface shape on the simulated elasto-plastic response of

cohesive soils
A.G. Papadimitriou
Department of Civil Engineering, University of Thessaly, Volos, Greece
A.D. Vranna
Department of Civil Engineering, Aristotle University of Thessaloniki, Thessaloniki, Greece
Y.F. Dafalias
University of California, Davis, USA
National Technical University of Athens, Greece
M.T. Manzari
George Washington University, USA
ABSTRACT: This paper presents a discussion on the effect of the selected yield surface shape on the simulated
drained and undrained response of cohesive soils. The discussion is made possible by sequentially implementing
two different yield surface shapes to a recently proposed elastoplastic critical state (reference) model, SANICLAY,
as alternatives to its own yield surface. The reference yield surface shape of SANICLAY is a distorted and rotated
ellipse, and the two studied alternatives have the shape of a distorted lemniscate and a distorted ellipsoid. For
each constitutive model variant, the remaining equations of SANICLAY were retained unaltered. It is shown,
that the use of the distorted lemniscate may lead to selectively more accurate simulations in comparison to the
reference, but to a less accurate overall response with the same number of model constants. On the contrary,
the use of the distorted ellipsoid provides an overall slightly enhanced simulative ability, but this at the cost of
one extra model constant.
1 INTRODUCTION 2 THE REFERENCE CONSTITUTIVE MODEL
There are many examples in the literature of com- SANICLAY is a simple anisotropic plasticity model
plicated constitutive models that do offer accurate proposed by Dafalias et al (2006) that provides rela-
simulation of soil behavior. In parallel, simpler con- tively satisfactory simulation of the rate-independent
stitutive models (e.g. employing elasto-plasticity with behavior of both normally consolidated and over-
an elastic region being defined by a convex yield sur- consolidated sensitive clays which do not exhibit
face in stress space) are still being used, especially in destructuration during loading. An extension to
boundary value problems, maybe at the cost of reduced include destructuration was recently published by
accuracy. This is due to the fact that the popularity Taiebat et al. (2009). SANICLAY builds on a
of a constitutive model is governed by the balance of modification of the associative flow rule isotropic
offered accuracy and simplicity in the equations and Modified Cam Clay (MCC) model (Burland 1965),
the calibration process. at the expense of merely three extra constants.
To this extent, “simple” models are also continu- In fact, its reference model (Dafalias 1986) con-
ously being proposed in the literature. The response of stitutes the simplest possible energetic extension
such models is governed by the shape of the adopted of the MCC model from isotropic to anisotropic
yield surface and on how it evolves with loading. This response.
paper explores the effect of the adopted yield surface SANICLAY is characterized by a non associative
shape on the simulated response of cohesive soils, by flow rule which is introduced by adopting a yield
sequentially implementing two different yield surface surface different than the plastic potential surface.
shapes to a recently proposed “simple” critical state Besides the isotropic hardening of the yield sur-
(reference) model, SANICLAY (Dafalias et al 2006), face, both surfaces evolve according to a combined
as alternatives to its own yield surface shape. kinematic and distortional hardening rule.
63
For simplicity, the formulation of the SANICLAY stress-ratio N . This is exactly the property that allows
(and its model variants) is presented here in the tri- for the undrained softening in compression to take
axial space, in terms of effective stress quantities place after Ko consolidation, and the very reason for
p = (σa + 2σr )/3, q = (σa − σr ) and strain quantities introducing an f different than g.
εv = (εa + 2εr ), εq = 2(εa − εr )/3, where subscripts a
and r denote the axial and the radial directions, 2.2 Rate evolution equations for the internal
respectively. variables
For the po variable the classical evolution law of critical
2.1 The SANICLAY model surfaces state soil mechanics is postulated as:
The plastic potential surface has the shape of a rotated
and distorted ellipse and is analytically described by
(see Fig.1):
where ein is the initial value of the void ratio e and
λ, κ are the slopes of the normal compression and
rebound lines, respectively, in the e−lnp space, while
where M is the critical stress ratio, α is a non- L is the loading index (in Macauley brackets yielding
dimensional variable which introduces anisotropy in < L > = L for L > 0 and < L > = 0 for L ≤ 0), which
the plastic potential and pα is the value of p at q = pα, is related to the partial derivatives of the f = 0 function
so that Eq. (1) is satisfied for ( p, q ) values at yield. in terms of p, q,po , and β, to the loading increments ṗ
As deduced from Figure 1, M = Mc when the and q̇ and the forms of po (Eq. 3) and β (see below),
stress ratio η = q/p > α and M = Me when η = q/p < α, based on standards methods on plasticity.
where Mc and Me are SANICLAY constants. Clearly, The rate evolution equation for α is described as:
one must have |α| < M for real-valued p, q in Eq. (1).
The SANICLAY yield surface is expressed simi-
larly to the plastic potential by:
where po , β and N respectively substitute for pα , α and where C is a SANICLAY constant, that also controls
M in Equation 1. the rate evolution equation for the β variable (Eq. 5):
In particular, β is the rotational hardening variable
of the yield surface that introduces anisotropy the same
way that variable α does in the plastic potential. N
is a SANICLAY constant similar in nature to M , but
taken the same in compression and extension for sim-
plicity. Clearly one must have |β| < N for real-valued Note that the ∂g/∂p term in Eqs (3) through (5) intro-
p, q in Eq. (2). Notice that q on f = 0 at the stress- duces the volumetric plastic strain rate and diminishes
ratio M is not the peak q stress, the latter occurring at the evolution of all surfaces at the critical state. Fur-
thermore, observe the operation of “attractors” αb and
βb in Eqs (4) and (5), that enforce the aforementioned
conditions of |a| < M and |β| < N for real valued p, q
in Eqs (1) and (2), respectively.
2.3 Reference model performance

Further details on the SANICLAY model may be found
in Dafalias et al (2006). Note here, that this model has
8 constants requiring calibration, i.e. the 5 constants of
the MCC model (Mc ,Me , κ, λ, ν) and 3 extra constants
(N ,C,x), as listed in Table 1.
As presented in detail in Dafalias et al (2006), SAN-
ICLAY leads to rather accurate simulations for low
OCR values with a possible underprediction of the
undrained shear strength (in compression, following
Ko -consolidation). Its response for high OCR values is
similar to that produced by the MCC model, i.e. it leads
to overprediction of yield stresses and hence stiffer
overall response in the early loading stages. In concept,
these problematic simulations may be corrected by
Figure 1. SANICLAY model surfaces in triaxial stress changing the yield surface shape, namely by “widen-
space. ing” it at large p/po values (leading to higher |q| values
64
at yield) and by “narrowing” it at small p/po values for q shows that in order to retain real valued p, q the
(leading to lower |q| values at yield). This exercise is following must hold:
what that gave birth to the work hereby presented.
3 ALTERNATIVE YIELD SURFACE SHAPE:

DISTORTED LEMNISCATE
In order to enforce the foregoing condition, the
absolute value of “attractor” βb in Eq. 5 must be set
3.1 Presentation of distorted lemniscate
equal to the term in the right hand side of Eq. (7).
The first alternative yield surface shape studied herein For comparison purposes, all remaining equations of
is the one introduced by Pestana & Whittle (1999) in SANICLAY remain unaltered.
their MIT-S1 model, which has the shape of a distorted Figure 2 shows an example comparison of the Dis-
lemniscate and is described by: torted Lemniscate shape (for m = 1.4, n = 0.77) to that
of the SANICLAY (for N = 0.91), as well as data
for undrained triaxial compression and extension tests
on isotropically (CIU) and Ko -consolidated (CKo U)
samples of Lower Cromer Till (LCT, see Table 1).
Observe in Fig. 2a, that Eq. (6) qualitatively offers
the aforementioned necessary changes in yield sur-
face shape for CI tests. In parallel, Fig. 2b shows that
where m, n are constants of the new yield surface Eq. (6) does the same for CKo tests, with a possible
shape that substitute for SANICLAY constant N. In exception of extension tests.
particular, the m defines the aperture of the yield sur- 3.2 Comparison with the reference model
face (in terms of η = q/p values) in the vicinity of
p = 0, whereas the n controls the yield surface width In order to ascertain whether the Distorted Lemniscate
(in terms of q values) from the q = pβ(p/po )n curve. yield surface shape has the potential to offer enhanced
Note that the higher the values of m and n the wider accuracy, pertinent simulations for LCT clay were
the yield surface shape becomes, always located along repeated with the new SANICLAY model variant and
the q = pβ(p/po )n curve. Observe that solving Eq. (6) compared to those of the original model. In all cases,
model constants were given the values of Table 1, with
the exception of constant N that no longer exists in
Table 1. SANICLAY constants and their values for LCT.
the SANICLAY variant model and is replaced by (m,
Constant Description Value n) = (1.4, 0.77).
Figures 3 and 4 present the LCT data and the
Mc Value of η at critical state in TC 1.18 simulations with the use of the SANICLAY and the
Me Value of η at critical state in TE 0.86 Distorted Lemniscate SANICLAY variant model, for
κ Compressibility of OC Clay 0.009 CIU and CKo U tests, respectively. Note that Figures 3
λ Compressibility of NC Clay 0.063 and 4 present the comparisons of simulations to data
ν Elastic Poisson’s ratio 0.2 only for OCR = 1, 2 and 7. The comparisons for the
N Shape of yield surface 0.91
remaining OCR values (see Fig. 2) lead to qualitatively
x Saturation limit of anisotropy 1.56
C Rate of evolution of anisotropy 16 similar conclusions, but are not included in the forego-
ing figures for reasons of clarity of the presentation.
Figure 2. Comparison of LCT data for undrained triaxial tests and shapes of the Distorted Lemniscate and the SANICLAY
yield surface, after: a) isotropic consolidation (CIU tests) and b) Ko-consolidation (CKo U tests).
65
Figure 3. Comparison of data to simulations with the use of the SANICLAY and the Distorted Lemniscate model variant for
CIU tests on LCT and various OCR values, in terms of: a) effective stress path, b) stress-strain response.
Figure 4. Comparison of data to simulations with the use of the SANICLAY and the Distorted Lemniscate model variant for
CKo U tests on LCT and various OCR values, in terms of: a) effective stress path, b) stress-strain response.
As inferred by the comparisons of the yield surface (2002), set on a different context. It has the shape of a
shapes in Figure 2, the new model offers better sim- distorted ellipsoid and is described by:
ulations for high OCR values. Nevertheless, for low
OCR values, although the new model offers better sim-
ulation of the peak strengths, the post-peak response
is qualitatively erroneous, since it characterized by
excessive strain softening for CIU tests and minimal
for CKo U tests, exactly opposite to what is shown by
the data and predicted by the original model. An exten- where the factor X is the novel addition to Eq.(2) in
sive parametric analysis of the simulated response terms of the additional constants z and n, for the new
offered by the new model shows that by appropriately yield surface shape.
choosing the set of (m, n) values one may selectively Notice that when z = 0 and/or n = 1, Eq. (8)
attain ameliorated simulations for the OCR range in becomes Eq. (2), which describes the reference yield
question. Yet, there is no unique set of (m, n) values surface. Hence, constant N for the distorted ellipsoid
that may offer ameliorated simulations as compared to plays the same role as for the SANICLAY yield sur-
the SANICLAY model for all OCR values. face shape, i.e. it dictates the overall width of the yield
surface (in terms of η = q/p values) around the η = β
line. This is the reason why the constant N of the
4 ALTERNATIVE YIELD SURFACE SHAPE: Distorted Ellipsoid behaves similarly with the SANI-
DISTORTED ELLIPSOID CLAY constant N . Regarding the other two constants
the following may be stated:
4.1 Presentation of distorted ellipsoid
– For values of n < 1, the ellipsoid distorts as com-
The second alternative yield surface shape studied pared to the SANICLAY ellipse, and “widens” at
herein is inspired by the work of Collins & Hilder large p/po values (>0.5) and “narrows” at small
66
Figure 5. Comparison of LCT data for undrained triaxial tests and shapes of the Distorted Ellipsoid and the SANICLAY
yield surface, after: a) isotropic consolidation (CIU tests) and b) Ko-consolidation (CKo U tests).
p/po values (<0.5), as required for ameliorated the CKo tests (for which β = 0.77 at the end of Ko
simulations, based on the discussion in section 2. consolidation).
For any given value of n < 1, an increase in the value
of constant z (to values z > n) furthers the distortion
4.2 Comparison with the reference model
in this qualitatively accurate manner.
– For values of n > 1, the distortion of the ellipsoid Again, in order to ascertain whether the Distorted
is qualitatively opposite to what is desired, since Ellipsoid yield surface shape proposed here has the
it “widens” at small p/po values (< 0.5) and “nar- potential to offer enhanced accuracy, pertinent sim-
rows” at large p/po values (>0.5), and will not be ulations for LCT clay were repeated with the new
discussed further. SANICLAY model variant and compared to those of
– It may be easily shown, that for any given set of the original model. In all cases, model constants were
β and N , there is an infinite set of (z, n) values given the values of Table 1, with the exception of con-
describing practically the same distorted ellipsoid stant N = 0.88 for the distorted ellipsoid model variant,
yield surface shape (differences in terms of q of less for which (z , n) = (0.88, 0.8) also holds.
than 2%). Figures 6 and 7 present the LCT data and the simula-
tions with the use of the SANICLAY and the Distorted
Based on all the above, the distorted ellipsoid yield
Ellipsoid variant model, for CIU and CKo U tests,
surface shape described in Eq. (8) is a potentially good
respectively, and for OCR=1, 2 and 7. As inferred
alternative to the SANICLAY ellipse, with merely two
by the comparisons, the use of the Distorted Ellip-
extra constants (N , z), since n may be considered prac-
soid as a yield surface shape provides apparently
tically fixed (e.g. n = 0.8 < 1). Again, solving Eq. (8)
improved simulations on CIU tests for both compres-
for q (and for n < 1) shows that in order to retain real
sion and extension tests and all OCR values. After
valued p, q the following must hold:
Ko -consolidation though, the new yield surface shape
shows an enhanced response only for compression
tests and this for high OCR values. On the whole, it is
deduced that the use of the Distorted Ellipsoid model
In order to enforce the foregoing condition, the abso-
variant leads to a slightly more accurate simulative
lute value of “attractor” βb in Eq. (5) must be set equal
ability, with the cost of only one extra variable. This
to the term in the right hand side of Eq. (9). For compar-
benefit is also underlined by drained triaxial test sim-
ison purposes, all remaining equations of SANICLAY
ulations on the same LCT clay (not shown here due to
remain unaltered.
paper length limitations), that shows enhanced accu-
Figure 5 shows an example comparison of the Dis-
racy after both isotropic and Ko -consolidation and for
torted Ellipsoid shape (for N = 0.88, z = 0.88, n = 0.8)
all OCR values.
to that of the SANICLAY (for N = 0.91), as well as
data for undrained triaxial compression and extension
tests on CIU and CKo U samples of LCT. Observe in
5 CONCLUSIONS
Fig. 5, that Eq. (8) qualitatively offers the aforemen-
tioned necessary changes in yield surface shape for
Based on this study, the following may be stated:
both CI and CKo tests. Nevertheless, the difference
between the two sets of yield surface shapes is not 1. Comparing yield surface shapes with undrained
large, since the N values are approximately equal for test data in the stress space offers the potential for
both sets and the value of z may not be increased too assessing the relative benefits of each candidate
much, since it would lead to a violation of Eq. (9) for shape, at least for tests that do not induce large
67
Figure 6. Comparison of data to simulations with the use of the SANICLAY and the Distorted Ellipsoid model variant for
CIU tests on LCT and various OCR values, in terms of: a) effective stress path, b) stress-strain response.
Figure 7. Comparison of data to simulations with the use of the SANICLAY and the Distorted Ellipsoid model variant for
CKo U tests on LCT and various OCR values, in terms of: a) effective stress path, b) stress-strain response.
yield surface rotations (e.g. CKo U extension test Collins, I. F., Hilder, T. 2002. A theoretical framework for
for OCR = 1). constructing elastic/plastic constitutive models of triaxial
2. Compared to the ellipse of the SANICLAY model, tests. International Journal for Numerical and Analytical
the Distorted Lemniscate is a flexible yield surface Methods in Geomechanics 26: 1313–1347.
Dafalias, Y. F. 1986. An anisotropic critical state soil plastic-
shape with merely one extra constant, that may lead ity model. Mechanics Research Communications 13(6):
to selectively accurate cohesive soil response sim- 341–347.
ulations. Nevertheless, it fails to provide enhanced Dafalias, Y. F., Manzari, M. T. & Papadimitriou, A. G. 2006.
accuracy for all loading histories (CIU and CKo U SANICLAY: simple anisotropic clay plasticity model.
tests) with the same set of values of the model International Journal for Numerical and Analytical Meth-
constants, at least when the remaining SANICLAY ods in Geomechanics, 30(12): 1231–1257.
equations remain unaltered. Pestana, J. M., Whittle, A. J. 1999. Formulation of a uni-
3. The use of the Distorted Ellipsoid hereby proposed, fied constitutive model for clays and sands. Interna-
provides a slightly improved response for all OCR tional Journal for Numerical and Analytical Methods in
Geomechanics, 23(12): 1215–1243.
values, loading histories and drainage conditions, Taiebat, M., Dafalias, Y. F., Peek R. 2009. A destructuration
as compared to the SANICLAY model. Never- theory and its application to SANICLAY model. Interna-
theless, its use requires one extra constant, and tional Journal for Numerical and Analytical Methods in
therefore it is up to the end user to decide on its Geomechanics, DOI: 10.1002/nag.841.
selection over SANICLAY.
REFERENCES
Burland, J. B. 1965. The yielding and dilation of clay.
Geotechnique 15 (2): 211–214.
68
Impact of input data on soil model calibration using Genetic Algorithms
D. Taborda & A. Pedro

Department of Civil & Environmental Engineering, Imperial College London, London, UK
P.A.L.F. Coelho
Department of Civil Engineering, University of Coimbra, Coimbra, Portugal
D. Antunes
Department of Informatics Engineering, University of Coimbra, Coimbra, Portugal
ABSTRACT: Soil model parameters can be obtained from experimental data using genetic algorithms-based
software. This paper highlights the efficiency of model calibration using this type of tool for the Modified
Cam-Clay model. The influence of the type and amount of information provided on the efficiency of the opti-
mization process and the quality of the solution obtained is also discussed. It is shown that including additional
experiments in the input may be ineffective, unless these reflect specific aspects of soil behaviour.
1 INTRODUCTION which represent mathematically the analysed problem,

generating a set of solutions. In geotechnical engi-
Genetic Algorithms (GAs) have been used in the past neering, for example, these can include stress-paths
to carry out the subjective and complex task of soil and stress-strain curves, the evaluation of settlements
model calibration. However, it is not often clear how due to the excavation of a tunnel, bending moments
many and which type of experiments are necessary to in a particular section of a retaining wall, etc. The
obtain a good quality solution. These aspects are dis- obtained solutions are then compared with the refer-
cussed herein, considering the effectiveness of genetic ence input data (laboratory testing results, instrumen-
algorithms-based software in obtaining the parameters tation values, etc.) and their quality is measured by a
required for the Modified Cam-Clay model (MCC) “fitness function”. This expression can either return
from pseudo-experimental sets of results representing the “error”, in which case the lower the value, the
a progressively larger number of different loading con- better the individual, or the “merit”, meaning that to
ditions. The results presented provide evidence about larger magnitudes correspond better solutions. Once
the viability of reaching a better solution and the more the whole population has been tested, if the best indi-
efficient way to do it. vidual satisfies certain quality criteria, then a solution
has been found.
However, the procedure described above is not very
2 GENETIC ALGORITHMS AS A
different from an automated trial-and-error procedure
CALIBRATION TOOL
and, consequently, it only performs well when very
large populations are employed (or, alternatively, when
2.1 Basic concepts
the search area is small). As a result, GAs employ a
The computational implementation of GAs requires set of techniques designed to improve substantially the
several processes to be adequately defined. Probably efficiency of the optimization algorithm. The first of
the most fundamental one concerns the mathematical these operations is termed “selection” and its purpose
representation of a solution of the analysed problem. is to decide on which individuals are discarded and
In the particular case where constitutive models are which are transferred to the next generation. Although
employed to reproduce a certain known behaviour, many methods exist, it is generally accepted that the
a vector containing the required parameters is the best should remain in the population. The second pro-
most natural choice. Each of the values is termed cedure, known as “crossing”, consists of determining
“gene” and the vector of genes composes a candi- which characteristics (i.e. genes) of the parent indi-
date solution (“individual”). The optimisation process viduals (“ancestors”) are transmitted to the offspring
is initiated by randomly generating a certain number (“descendants”). Finally, the third genetic operator,
of individuals (“population”) from within a previously termed “mutation”, introduces small changes to these
defined search area. Sequentially, each of the indi- genes when the offspring is being determined. Once
viduals is tested in a series of numerical procedures, the individuals composing the new generation are
69
known, the previously described assessment process all of the input data given to genUC consisted of
is restarted. Naturally, due to the evolutionary aspect pseudo-experimental results. These were obtained by
of the GAs, it is expected that the best individual of simulating triaxial compression and oedometer tests
a given generation corresponds to a better solution using constitutive equations similar to those imple-
than the best individual of the previous population. mented in the code. As a result, genUC should be able
In the present case, this means that the ability of a to fit perfectly the MCC to the chosen curves while
given constitutive model to reproduce a certain mate- the proximity to the known solution (i.e. the parame-
rial behaviour increases with the number of analysed ters selected to generate the input data) can be used to
generations. evaluate the global performance of the method.
The biggest advantage of this method is probably In terms of particular aspects of the algorithm, a
the fact that many of the processes described above population of 64 individuals; a variable probability
are not, and cannot, be rigidly defined. As a result, alternate crossing process and a continuous muta-
there is considerable freedom when determining how tion scheme (with +/− 5% change) have been chosen
the “crossing” or “mutation” should be carried out, based on previous studies (Azeiteiro et al. 2009).
allowing adjustments to be made during the optimisa-
tion procedure. However, this also implies that many
3.2 Input data
conclusions drawn from research studies might be
problem dependent, thus requiring all aspects of the A total of 7 pseudo-experimental tests, which are char-
method to be separately evaluated. acterised in Table 1 in terms of loading and drainage
conditions, stress level and OCR, were generated using
2.2 Developed software and previous studies the parameters listed in Table 2. Furthermore, all of the
The computer code genUC has been developed as a triaxial compression tests were assumed to start from
general platform for the application of GAs to any an isotropic stress state (K0 = 1.0).
given problem. In its current version, however, only The undrained triaxial compression tests were
the constitutive equations of the Modified Cam-Clay defined by the stress path in p’-J plane and the stress-
(MCC) model with a constant Poisson’s ratio (Potts & strain curve Ed -J, where p’ is the mean effective stress
Zdravković 1999) are implemented, since research so and J and Ed are the second invariants of the stress
far has focussed on the computational aspects of the and strain tensors, respectively (Potts & Zdravković
optimisation method, rather than on the merits of any 1999). The drained triaxial compression tests were
specific model (Azeiteiro 2008, Azeiteiro et al. 2009). characterised by the stress-strain curve Ed -J and by the
In particular, the code uses the MCC model to simulate Ed − εvol relation, where εvol is the volumetric strain.
the behaviour of soil under three distinct types of lab- Finally, the one-dimensional consolidation tests were
oratory tests: undrained triaxial compression; drained described by the obtained e-p’ curve, where e is the
triaxial compression and one-dimensional consolida- void ratio. All of the curves mentioned above were
tion (oedometer).Therefore, data corresponding to any discretised into 25 points.
combination of these tests (number and type) can be In terms of search area, Azeiteiro (2008) showed
given as an input of the calibration procedure. that, for a problem with a well-defined solution, such
In terms of features related to the GAs, the current as the one being analysed, the limits chosen did not
version uses a constant-size population of individuals influence the quality of the results, only affecting
composed of real-format genes.There are several types
of crossing operators available, both deterministic Table 1. Characteristics of the pseudo-experimental tests.
and probabilistic. Similarly, discrete and continuous
mutation schemes with different maximum allowed p’ init. p’ end
changes per gene (1% up to 10%) have been imple- Test Type (kPa) OCR (kPa)
mented. These aspects have been studied and the
obtained results can be found in Taborda et al. (2008) UD.1 UTXC 100 6.0 Critical state
UD.2 UTXC 100 1.5 Critical state
and Azeiteiro et al. (2009). Other existing capabili-
UD.3 UTXC 200 1.0 Critical state
ties, such as the fitness functions and age control are UD.4 UTXC 500 1.2 Critical state
described in Azeiteiro (2008). DR.1 DTXC 100 1.5 Critical state
DR.2 DTXC 200 1.0 Critical state
OE.1 1-DC 100 1.5 1000
3 COMPUTATIONAL STUDIES
UTXC – Undrained triaxial compression test
3.1 General aspects and methodology DTXC – Drained triaxial compression test
This computational study aims at investigating and 1-DC – one-dimensional consolidation test
characterising the effect of the input data on the effec-
Table 2. Employed Modified Cam-Clay parameters.
tiveness of GAs, when this optimisation method is
employed to estimate the parameters of the MCC κ λ ν1 ϕcs µ
model.
To eliminate the impact of the limitations of the 0.05 0.20 3.0 32.0◦ 0.25
chosen model when reproducing real soil behaviour,
70
Table 3. Limits of the searched area.
Limit κ λ ν1 ϕcs µ
Min. 0.02 0.10 2.50 25.0 0.15

Max. 0.08 0.30 3.50 42.0 0.35
slightly the performance of the algorithm. In the

present case, the adopted limits, which are defined in
Table 3, are rather wide and are thought to influence
equally all of the studied cases.
3.3 Study based only on undrained triaxial Figure 1. Evolution of the normalised error with the size of
compression tests (UD) the sample.
For this particular study, tests UD.1; UD.2; UD.3 and
UD.4 were used as input data. The computer code was
set to stop at 4 different marks (a normalised mea-
sure of the error determined internally by genUC’s
fitness function, which, in practice, evaluates the dis-
tance between the input and the calculated curves): 1,
10, 100 and 200. Naturally, the greater the mark used
as stopping condition, the larger is the expected error
when comparing the returned solution with the original
set of values. This error was first quantified for each
of the parameters:
Figure 2. Mark and error of best individuals (case UD).

allowing, subsequently, the determination of the global
error for each solution: (instead of the more economic 100) in order to ensure
that the average result was significant.
For each of the 4 cases, Figure 2 shows the respec-
tive mark and error of the best 200 individuals. It is
interesting to note that the envelope of the points tends
It is acknowledged that a perhaps more representa- to a horizontal line at lower marks, indicating that
tive measure of the error could be obtained by dividing the adoption of a more stringent stopping criterion
the result of Equation 2 by the number of parameters might not reduce the expected error of the solution.
(Nparam ). However, since Nparam is constant for all the Furthermore, it can also be seen that the triangular
analyses performed, this normalisation was thought to shape formed by the points corresponding to high
be unnecessary and the global norm of the error (Eq. 2) marks collapses into an almost rectangular shape at
was deemed appropriate. It is apparent that the mag- lower marks. In fact, when determining the relative
nitude of this error does not truly reflect the quality of frequency of the normalised errors (i.e. the error of
the solution. Furthermore, for this error to be signif- an individual divided by the maximum error of the
icant for a given analysis, it is necessary to perform entire sample obtained using an identical stopping con-
each calculation various times, thus limiting the nat- dition), it becomes clear that the distribution tends to
ural influence of all the random processes involved have almost a uniform shape for lower marks (Fig-
in the generation of the initial population. Since this ure 3). This observation is also corroborated by the
number is likely to be a function of several factors, statistical measures of dispersion, which indicate that
it is difficult to determine its magnitude in advance. the quartiles of the samples are very close to 25%, 50%
Therefore, each calculation was repeated 200 times and 75% – the values of a uniform distribution – for
and the average error of the best individuals was deter- marks below 10. In fact, it is interesting to note that
mined for samples of various sizes (limited to 200, the first quartile is close to 25% for the three analyses,
naturally). The results, normalised by the error (Eq. 2) while the third quartile suffers the largest changes –
of the full sample, are shown in Figure 1, where it is from 50.3% for 200 up to 71.7% for 10 – and is only
evident that above 100 calculations, a maximum devi- approximately 75% for the two lowest marks. This sug-
ation of +/− 5% can be guaranteed. Nevertheless, as gests that the points located below a certain error level
these conclusions are likely to be problem dependent, are distributed uniformly while the same is not true for
it was decided to conduct 200 calculations per analysis the individuals with larger errors. If it is considered that
71
Figure 3. Relative frequency of the normalised error (case Figure 4. Mark and error of best individuals (case UD-DR).
UD).
Table 4. Dispersion of the normalised error (case UD). error to be below 10% as it is for it to belong to the
interval 20–30%.
Case Q1 (25%) Median Q3 (75%)
3.4 Introduction of drained triaxial compression
Uniform 25.0 50.0 75.0
1 27.5 45.4 80.7 tests (UD-DR)
10 23.0 48.3 71.7 This study used as input data the tests labelled as
100 27.8 40.1 52.0 UD.1, UD.4; DR.1 and DR.2 in order to evaluate if the
200 25.9 36.5 50.3
trend verified in the previous set of analyses, which
indicated an inability to further reduce the error by
decreasing the mark used as stop condition, remained
Table 5. Dispersion of the normalised error with partial
samples (case UD). valid for this combination of tests. With that objective,
200 calculations for each of the 3 chosen marks (10,
Case N Q1 (25%) Median Q3 (75%) 100 and 200) were performed and the obtained best
individuals are presented in Figure 4 (mark 1 was dis-
Uniform – 25.0 50.0 75.0 regarded since for case UD there was little difference
1 200 27.5 45.4 80.7 when compared to mark 10).
10 187 22.5 47.2 70.2 Although the maximum relative error registered for
100 179 41.7 55.9 75.1 the analysis using 200 as a stopping condition is sig-
200 144 46.5 68.0 82.6 nificantly reduced by the introduction of the drained
triaxial compression test, the results show that the
N is the new size of the sample
envelope tends to the same limit of about 30% at lower
marks. Furthermore, the corresponding relative fre-
for mark 1 the distribution is uniform, then the position quency of errors, as indicated by Figure 5, show a
of this division line must coincide with the asymptotic progression towards a uniform distribution, justify-
limit defined by the maximum relative error registered ing the approximately rectangular shape formed by
when the mark was set at this level. To test this hypoth- the points obtained when the mark is set to 10. This
esis, the quartiles were recalculated for the 4 samples observation is further confirmed by the values of the
using only the individuals with error below 28.95% quartiles of the different samples, listed in Table 6.
(the maximum error registered for mark 1). The results It is interesting to note that, as in the previous case,
are listed in Table 5 and confirm that, with the excep- the first quartile is close to 25% for the three analy-
tion of the value of the first quartile for larger marks, ses, while the third quartile suffers the largest changes
the distributions of the errors below this level are more (from 56.5% for 200 up to 71.9% for 10) and is only
uniform than those obtained when the complete sam- close to uniform for the lowest mark. Therefore, to
ples are considered. Therefore, it can be speculated further test the hypothesis raised before, the quar-
that if such a division line (i.e. a value below which a tiles were recalculated using only the individuals with
given sample is distributed uniformly) is found, then an error below 29.2% (the maximum error registered
its position may correspond to a limit error which can when the employed mark was 10). The results indicate
never be eliminated when using the current input data. more uniformly distributed errors (Table 7), agreeing
In conclusion, it can be stated that by only supply- with the previous case.
ing results of undrained triaxial compression tests the In conclusion, the analysis of the collected data indi-
error in the obtained solution tends to be defined by a cates that the introduction of the pseudo-experimental
uniform distribution, with a maximum magnitude of results of a drained test did not increase the ability of
approximately 30% (5% on average for each param- the MCC model to determine with accuracy the right
eter). Therefore, as an example, it is as likely for the parameters.
72
Figure 5. Relative frequency of the normalised error (case Figure 7. Relative frequency of the normalised error (case
UD-DR). UD-DR-OE).
Table 6. Dispersion of the normalised error (case UD-DR). Table 8. Dispersion of the normalised error (case UD-DR-
OE).
Uniform 25.0 50.0 75.0
10 23.2 49.3 71.9 Uniform 25.0 50.0 75.0
100 22.7 45.2 63.4 10 40.0 55.8 76.3
200 23.2 39.9 56.5 100 57.3 74.2 81.9
200 38.8 57.7 71.5
Table 7. Dispersion of the normalised error with partial the effect of the introduction of a test characterising a
samples (case UD-DR).
different type of loading condition (one-dimensional
Case N Q1 (25%) Median Q3 (75%)
compression instead of triaxial shearing). The results
of the 200 calculations performed for each of the
Uniform – 25.0 50.0 75.0 values of the stopping criterion (mark 10, 100 and
10 200 23.2 49.3 71.9 200) are illustrated in Figure 6. It is clear from the
100 189 23.7 48.5 68.4 graph that there is an apparent relation between the
200 179 30.1 50.2 72.8 mark and the registered errors, while the asymptotic
limit identified for the previous cases cannot be dis-
N is the new size of the sample cerned. Consequently, with this input data, and unlike
the situations analysed so far, a more stringent stop-
ping criterion does lead to better results. Furthermore,
Figure 7, which illustrates the relative frequency of
the normalised errors, shows that the distributions
cannot be considered uniform. This observation is fur-
ther confirmed by the values of the quartiles listed
in Table 8 and agrees with the hypothesis proposed
before. In fact, if there is no value of error below which
the distribution is uniform, then it must mean that the
limit error using this input data is very close to 0 and
that the employed tests are able to uniquely define the
parameters of the MCC model. This result may seem
trivial as the constitutive equations were initially for-
mulated to reproduce the observed behaviour of clay
under these two tests. However, since for other mod-
Figure 6. Mark and error of the best individuals (case els this may not be as clear, the results of the present
UD-DR-OE). study suggest that the GAs may be used as part of a
process to identify the most adequate set of tests for
3.5 Introduction of oedometer tests the efficient calibration of a given constitutive model.
(UD-DR-OE)
3.6 Comparative analysis
For this group of analyses, the data from the tests
named UD.1; UD.4; DR.1 and OE.1 were given to The error envelopes of the three groups of analyses
genUC as input with the objective of investigating are shown in Figure 8. It is clear that results obtained
73
in the analyses. Moreover, it is interesting to note that,
for all the cases where the mark was set at 10, the
program registered about 13.4 iterations per second.
Consequently, the fact that the average time spent per
calculation increases with the number of different tests
modelled is related to a larger number of iterations
being required, rather than to a higher computational
cost per iteration.
4 CONCLUSIONS
The effectiveness of genetic algorithms-based soft-

ware in soil model calibration is clearly demonstrated
Figure 8. Error envelopes of the three groups of analyses. by the results presented. It is also shown that differ-
ent types of loading conditions need to be included in
the input if a very accurate solution is pursued. For
the case studied, which considers the evaluation of
parameters for MCC model, adding experimental data
corresponding to drained triaxial compression tests
to the reference set of data including only undrained
triaxial tests merely increases the efficiency of the
optimization procedure. However, the quality of this
solution cannot be improved beyond a certain limit,
irrespective of the number of iterations performed.
A more accurate solution can only be obtained if the
results of a 1-D compression test are included in the
input.
The conclusions are consistent with the obvious fact
that some experiments provide a better representation
of particular aspects of soil behaviour than others,
Figure 9. Calculation times for the different groups of
analyses. which is particularly evident in the case of MCC
model. However, this conclusion is likely to retain its
validity for more complex constitutive relations, where
by only performing triaxial shearing tests, drained or the material’s response under certain loading condi-
undrained, cannot guarantee the determination of the tions is hard to relate to specific model parameters
correct set of parameters. In fact, adding a drained which often have no clear physical meaning. Further-
triaxial compression test to the reference undrained more, the fact that the enlargement of the experimental
shearing data only decreased the error for large marks, database for the calibration of soil models may not be
maintaining the same asymptotic value of about 30% effective does not necessarily apply to practical prob-
which cannot be eliminated. This can be justified lems, where the crucial issue of soil variability must
by the fact that triaxial shearing behaviour depends be addressed through this route.
more on a combination of parameters rather than on
their individual magnitudes. However, the introduction
REFERENCES
of the results of a one-dimensional compression test
greatly improved the quality of the solutions, leading Azeiteiro, R.N. 2008. Application of genetic algorithms to the
to much smaller errors and indicating that, with this calibration of constitutive models for soils. MSc Thesis,
combination of tests, there is a direct relation between University of Coimbra, Portugal (in Portuguese).
the mark and the error. Azeiteiro, R.N., Coelho, P.A.L.F., Taborda, D., Pedro, A. &
In terms of performance, the calculation times were Antunes, D. 2009. Computational study of the perfor-
used to evaluate the effectiveness of the different cases mance of a genetic algorithms-based software. Proc. of
analysed. In Figure 9, the average error is plotted the 1st Int. Symp. on Computational Geomechanics –
COMGEO I, Juan-les-Pins, 29 April – 1 May 2009.
against the calculation time, while the labels over the
Potts, D.M. & Zdravković, L. 1999. Finite element analy-
points indicate the mark employed as stopping cri- sis in geotechnical engineering: theory. London: Thomas
terion. It can be seen that, for the case where only Telford.
undrained triaxial compressions tests were used, there Taborda, D., Coelho, P.A.L.F., Antunes, D. & Pedro, A.
is a large increase in computational cost without any 2008. Genetic algorithms as a calibration method for
reduction of the error when the mark is decreased from constitutive models. Proc. of the 11th National Confer-
10 to 1. This is undoubtedly related to the existence of ence in Geotechnics – XI CNG, Coimbra, Portugal (in
the asymptotic limit on the maximum error obtained Portuguese).
74
Influence of destructuration of soft clay on time-dependant settlements
D.F.T. Nash
Department of Civil Engineering, University of Bristol, Bristol, UK
ABSTRACT: The time-dependant settlement of soft soils following application of surface loading may be
modelled using elastic visco-plastic constitutive models to describe the soil behaviour. For applied loadings
that increase the stresses to around the in-situ yield stress, the predicted behaviour is strongly influenced by
the associated breakdown of clay structure and the way in which this is modelled. The paper describes some
predictions for a hypothetical case prediction exercise recently organised by the Norwegian Geotechnical Institute
(NGI), to compare different calculation methods used in settlement analyses of soft soil around the world. The
results presented here were made using a one-dimensional coupled consolidation analysis implemented in the
spreadsheet-based software Briscon. Parameters were obtained from an oedometer test and the results were
extrapolated over the full soil profile. Various plausible assumptions about the shape of the isotaches around
yield were explored, and it is shown that the predicted long-term settlement may vary by factors of two or more
depending on the assumptions made.
1 INTRODUCTION
1.1 Background
The prediction of long-term settlements of embank-
ments foundations and fills on soft clays requires a
good understanding of the time dependant behaviour
of the clay. While a complex three dimensional analy-
sis may sometimes be necessary, it will often be more
appropriate to undertake a one-dimensional analysis
of the centerline conditions. The study presented here
explores the implications of making different assump-
tions about the one-dimensional creep behaviour of the
clay as it is stressed towards and beyond the yield stress
as structure of the soft clay is gradually damaged.
1.2 Origins of isotache models

One-dimensional compression under constant effec-
tive stresses was first studied by Taylor & Mer-
chant (1940), who showed that secondary compression
movements decrease logarithmically with time. Taylor
(1948) stated that creep occurs during primary consol-
idation as well as subsequently, and following Taylor’s
ideas, Šuklje (1957) and notably Bjerrum (1967) pre-
sented diagrams showing a system of approximately Figure 1. Principles of settlement calculation for soft clays
parallel e vs. logσ curves (see Figure 1) that describe (Bjerrum 1967).
secondary compression behaviour. In this widely used
diagram, the lines indicate void ratio after constant
time for delayed compression. Bjerrum introduced the labelled to indicate equal times after the application of
terms instant and delayed compression to describe the the loading, but Bjerrum (1972) noted that they also
behaviour of the soil skeleton in the absence of pore represent lines of constant creep rate or isotaches. Sub-
pressure effects, and argued that delayed compression sequently Lerouiel (2006) has shown that the isotaches
(or creep) occurs during the whole consolidation pro- in one-dimensional compression are part of strain rate
cess. The parallel e vs. logσ lines in Figure 1 are dependant limit state surfaces.
75
Figure 2. Observed and predicted settlements of buildings
in Drammen, Norway. (Bjerrum 1967).
Bjerrum (1967) presented interesting observations

of settlement of buildings founded over soft clay in
Drammen (see Figure 2) and showed that the mag-
nitude and rate of settlement was strongly influenced
by the degree of loading in relation to the yield stress
or pre-consolidation pressure p/(pc − p0 ). He sug-
gested that if this ratio is less than 50% then the
creep settlement within the lifetime of the buildings
would be small. Bjerrum also predicted the long- Figure 3. Isotaches observed on reloading (dashed) merging
term settlements and suggested that eventually all the with those for normal consolidation.
time-settlement curves would be parallel as shown in
Figure 2. This behaviour follows from the shape of the spacing of the isotaches as well as the initial state of
equal time lines shown in Figure 1. the clay. In practice the clay in-situ is probably under-
going creep at a very slow rate as indicated by point B*
1.3 Shape of isotaches around yield in Figure 3, but the state of the soil must be adjusted
to take account of sample disturbance.
In Figure 1 the isotaches are parallel and their position
is not affected by the development of clay structure nor
by physical over-consolidation, although the spacing 1.4 Objectives of this study
of the isotaches decreases with elapsed time. When fur-
In predicting creep settlements resulting from limited
ther load is applied and the yield stress is approached,
loading, the choice of isotache model is obviously fun-
the state path crosses the isotaches, and the increas-
damental. The two models shown diagrammatically in
ing creep rate is associated with structural breakdown.
Figure 4 were applied in a benchmark study proposed
This is illustrated in Figure 3 (developed from Bjer-
by NGI (NGI 2007, Jostad & Degago 2010), and as
rum’s plot), from which it may be deduced that a soil
expected the choice of model significantly influences
element that has crept from state A when it was initially
the short and long-term behaviour. Combined consoli-
deposited 10000 years ago to state B today would cur-
dation and creep analyses have been undertaken using
rently be undergoing creep at a very slow rate that
the elasto-visco-plastic procedure Briscon using these
can be evaluated assuming that the parallel isotaches
isotache models; a brief outline of the procedure is
apply.
given below. The results were first presented at an
However in laboratory IL oedometer tests on high
informal workshop on creep on soft soils (CREBS III)
quality samples that are reloaded to the in-situ state
held at Chalmers University in 2009.
B, creep is observed at rates that are initially compa-
rable with those in the normally-consolidated region
but the creep rate falls off rapidly. Thus as the clay is
2 BRIEF DESCRIPTION OF ANALYSIS
loaded towards the yield stress, the creep behaviour
is not consistent with the parallel isotaches, but rather
2.1 Brief description of BRISCON
with isotaches that are curved, as shown by the dashed
isotaches in Figure 3. This implies that creep occurs Briscon is an Excel™ spreadsheet-based program for
well inside the limit state surface with the gradually analysis of large strain one-dimensional consolidation
increasing isotache spacing associated with structural in multi-layered soils exhibiting creep with vertical
breakdown. and radial flow. The program was originally developed
The long-term strains resulting from a modest for use in back-analysis of the centreline settlement
increase in stress are clearly strongly influenced by the behaviour of an embankment on soft clay near Bristol,
76
(Nash & Ryde 2001, Nash 2001). It is written in Visual
Basic for Applications (VBA), and data is read from
and written to a series of spreadsheets. The ground
profile may include alternating drained sand and con-
solidating clay layers with varying hydraulic boundary
conditions. Permeability (linked to void ratio) and
geometry may be updated throughout the analysis,
with allowance made for gradual submergence of the
fill placed at the surface.
In Briscon the consolidation equation with vertical
flow only is expressed thus:
where u is the total pore pressure, σ is the total verti-

cal stress, kz is hydraulic conductivity, mv is the elastic
compressibility and the last term ∂εtp /∂t expresses the
natural strain creep rate. Thus creep is modelled con-
tinuously throughout the consolidation. The equation
is formulated in natural strain, but void ratio is an input
parameter. A standard finite difference procedure with
a graded grid is used to represent the problem.
2.2 EVP models used in this study

In the elastic visco-plastic (EVP) models that are
implemented in Briscon the plastic strain rate ε̇tp is
uniquely related to the soil state given by void ratio e
and the current vertical effective stress σ . One such
model was developed from the work of Yin & Graham
(1996) who used the λ − κ model from critical state
soil mechanics to define the elastic-plastic behaviour
of the soil skeleton. The normal consolidation line Figure 4. Isotaches for a) model 3 and b) model 3-d.
(NCL) is replaced by a reference time line (RTL)
with slope λ on which the creep strain rate is known,
that is used to define the complete set of isotaches; the yield stress σy has been taken as equal to the pre-
equally spaced isotaches indicate a logarithmic change consolidation pressure pc determined from the 24-hour
in strain rate, with the spacing controlled by parameter IL oedometer test, and defines the intersection of the
ψ. There is no lower limit to creep. The parameter ψ, elastic reloading line and the RTL on a plot of void
has a similar meaning to the coefficient of secondary ratio against logσ . In Figure 4 the stress is normalised
compression Cα except that it is independent of start- by the initial vertical effective stress. The creep rate
ing time; at very large times Cα = ψ. ln(10). The ratio at the initial effective stress is determined from the
ψ/λ is equal to Cα /Cc in the normally consolidated appropriate isotache.
region (ie on the RTL). Details of this model (denoted
model 3 here) are given by Nash & Ryde (2001).
The two creep parameters that define the creep
tp 3 BENCH MARK PROBLEM
behaviour are thus ψ and ε̇0 on the RTL. It is conve-
nient but not essential to choose the creep strain rate 3.1 Outline of problem
on the RTL equal to that after 24 hours in incremental
load oedometer tests. In the work of Yin and Graham The stratigraphy proposed by NGI (NGI 2007,
it was assumed that the isotaches are linear on a plot of Jostad & Degago 2010) is shown in Figure 5. Five
void ratio vs effective stress. Other models have been cases were proposed for analysis which differ primar-
developed (Nash 2001) that allow for a curved RTL ily in respect of the OCR assumed for the soft clay, the
and associated isotaches. magnitude of the load increment applied, and whether
Two EVP models have been used in the analyses the base of the clay layer is an open or closed bound-
reported here as follows: Model 3 (see Figure 4a) with ary. In this paper only one case will be considered in
parallel isotaches defined by λ, κ and ψ and Model which a uniform loading of 90 kPa is applied over a
3-d (see Figure 4b) similar to model 3 but with a bilin- clay deposit with OCR = 1.4 as shown in Figure 5. It
ear RTL and curved isotaches around yield to account may be seen that the final stress straddles the initial
for the effects of gradual destructuration. In this study yield stress pc .
77
Table 1. Model parameters derived from oedometer test
no. 693.
e0 1.17
Cc /(1 + e0 ), Cc , λ 0.167, 0.362, 0.157
Cs /(1 + e0 ), Cs , κ 0.021, 0.045, 0.020
yield stress σy 152 kPa
Cα /(1 + e0 ), /(1 + e0 ), 0.0067, 0.0029, 0.0063, 0.04
, /λ
Strain rate on RTL 2 × 10−4 %/min
Figure 5. Geotechnical profile for problem analysed.
Figure 7. Observed and predicted stress strain behaviour

for oedometer test no. 693.
obtained and normalized. The chosen fit lines are

shown in Figure 6a and the parameters are given in
Table 1. The steps involved were as follows:
1. From the plot of strain vs stress (log) – Plot the best
fit 24 hour NCL; determine Cc /(1 + e0 ); Choose
best fit average slope for initial reload and swelling
lines - Cs /(1 + e0 );
Calculate Cc and Cs , λ and κ values;
Determine yield stress σy from intersection of
Figure 6a. 24-hour data from test no. 693 with fitted lines.
initial reload line with NCL.
2. From time-settlement increment data – Assess
slope of final stages of each strain vs log(time)
plot (Figure 8a) to obtain Cα /(1 + e0 ) for each load
increment; calculate values of ψ;
Select representative value of ψ for the RTL (24
hour-NCL) and calculate ψ/λ;
Calculate strain rate on 24-hour NCL.
3. From increment data for creep parameter ψ -;
Plot ψ/(1 + e0 ) against normalised stress level σ /σy
(Figure 6b) and fit sin function.
3.3 Initial modelling of oedometer test

The data were used in a simulation of oedometer test
no. 693 using Briscon with models 3 and 3-d. The
Figure 6b. Creep parameter from test no. 693 plotted results are shown in Figures 7 and 8; the stress strain
against normalized stress with fitted sin function. plot shows good agreement between the experimental
and predicted data. The small differences around yield
reflect the different shape of the isotaches – agreement
3.2 The IL oedometer test
appears to be satisfactory. The time settlement plots
The significant data for the clay was a single IL (using model 3-d1) also show fair agreement with the
oedometer test (no. 693) from which parameters were experimental data (Figure 8).
78
Figure 9. Comparison between 24-hour NCL in test no. 693
and assumed RTL.
Figure 8. a) Observed settlement and b) settlement pre-

dicted with model 3-d1 vs time for loading increments in
test no. 693.
3.4 Model parameters for the full-scale

problem
The parameters in Table 1 were then applied directly
in the analyses of the full scale problem. Since test no.
693 was assumed to be representative of a clay stra-
tum with a specified OCR = 1.4, the most significant Figure 10. Predicted settlement vs time for a) no load and
additional assumption was the location of the refer- b) 90 kPa.
ence time line (RTL). The experimental data from test
no. 693 showed that strains of more than 2.5% had 4 ANALYSES OF FULL SCALE PROBLEM
occurred during reloading to around 108 kPa (corre-
sponding to an OCR of 1.4). Much of this strain was Two cases are reported here; others are described by
assumed to arise from recompression of the sample Nash (2008). In the first no loading was applied to
after disturbance during sampling and test prepara- the ground surface and the ground was allowed to set-
tion, and in applying this data to the full scale problem tle under its own weight. Secondly a surface loading
it was decided to disregard this. It was assumed that of 90 kPa was applied to the ground surface over a
the slope of the reload and RTL lines would remain period of 0.01 years and then the ground was allowed
unchanged (λ and κ) and that they would intersect at to consolidate and settle for 500 years. Analyses were
the yield stress σy . The yield stress profile has been undertaken with creep models 3, 3-d1 and 3-d2, and
taken equal to the initial effective stress multiplied for comparison without creep (model 2 using the same
by the given OCR. The in-situ RTL is located at a NCL).
smaller strain than the experimental NCL as shown in Firstly, analyses were carried out with no loading
Figure 9. applied, so as to examine the inherent creep of the
With model 3 the isotaches are at constant spac- ground. The time-settlement behaviour is summarised
ing (Figure 4a), and the creep rate at the initial state in Figure 10a. It shows that the 500 year settlements of
is controlled by the OCR. With model 3-d the sinu- the unloaded ground were 0.36 and 0.15 m with creep
soidal variation of creep parameter shown in Figure 6b models 3 and 3-d1 respectively. With creep model 3-d2
was used, leading to a varying spacing of isotaches settlement was negligible as expected. This illustrates
(Figure 4b). The creep rate at the initial state is con- how the model used in the over-consolidated region
trolled both by OCR and by the equivalent age of the strongly influences the calculated creep.
deposit. Two ages were used: 24 hours (point B in Fig- When 90 kPa load was applied there were still sig-
ure 3) and 10000 years (point B*); these models are nificant differences between the final settlements pre-
denoted 3-d1 and 3-d2 respectively. dicted by the various models as shown in Figure 10b.
79
is broken down around the yield stress. Without this,
using a simple isotache model can result in significant
over-prediction of creep settlements. The models used
here appear to be conceptually similar to that described
by Claesson (2003).
As is also discussed by Jostad & Degago (2010),
the main limitations of isotache models such as these
reflect our lack of knowledge of the creep behaviour
below and around the yield stress, at very large times
and after unloading. It is uncertain how the creep
behaviour of laboratory samples should be interpreted
to allow for sample disturbance.
REFERENCES
Bjerrum, L. 1967. Engineering geology of Norwegian
normally-consolidated clays. Seventh Rankine Lecture.
Geotechnique 17(2): 81–118.
Bjerrum, L. 1972. Embankments on soft ground. State of
the Art report. Proc. ASCE Spec. Conf. on Performance of
Earth and Earth-supported structures, Purdue. 1: 1–54.
Claesson, P. 2003. Long term settlements of clays. PhD thesis,
Chalmers University of Technology. Götheborg, Sweden.
Jostad, H.P. & Degago, S.A. 2010. Comparison of methods
for calculation of settlements of soft clay. 7th Eur. Conf.
Numerical Methods in Geotech. Eng. Trondheim, Norway.
Leroueil, S. 2006. Šuklje Memorial Lecture: The isotache
Figure 11. Predicted variation of soil state under 90 kPa approach. Where are we 50 years after its development by
load for a) model 3 and b) model 3-d2. Professor Šuklje? 13th Danube-European Conf. Geotech.
Engng, Ljubljana, Slovenia. 2: 55–88.
Nash, D.F.T. & Ryde, S.J. 2001. Modelling the consolida-
Without creep (model 2), the long-term settlement tion of compressible soils subject to creep around vertical
was 0.26m, whereas using creep models 3, 3-d1 and drains. Géotechnique 51(4): 257–273.
3-d2, the values increased to 1.13, 0.58 and 0.41 m Nash, D.F.T. 2001. Modelling the effects of surcharge to
reduce long term settlement of reclamations over soft
respectively. Again the choice of model has a profound clays: a numerical case study. Soils and Foundations 41
influence on the long-term settlements. Figure 11 (5): 1–13.
shows how the state of the clay at the top, centre and Nash, D.F.T. 2008. Comparison of methods for calculation
base of the clay stratum varies during consolidation. of long term settlements of soft clay – hypothetical case
The data are normalized by the initial effective stresses prediction exercise. Report.on predictions using Briscon.
and are superimposed on isotaches spaced logarithmi- University of Bristol (unpublished).
cally. The state of the soil in-situ lay below that on Norwegian Geotechnical Institute 2007. Comparison of
the RTL at all times since even near to the top of the methods for calculation of long term settlements of soft
clay, the rates of strain in the field did not reach those clay – a hypothetical case prediction exercise. Report
20061075-1.
obtained after 24 hours in the oedometer. Comparison Šuklje, L. 1957. The analysis of the consolidation process by
of the two figures shows how modelling the changing the isotache method. Proc. 4th Int. Conf. on Soil Mech.
isotache spacing around yield reduces the strains in the and Found. Engng., London. 1: 200–206.
lower levels of clay where the stresses never reaches Taylor, D.W. & Merchant, W. 1940. A theory of clay con-
the yield stress (see Figure 5). solidation accounting for secondary compression. Journ.
Math. Phys. 19(3): 167–185.
Taylor, D.W. 1948. Fundamentals of Soil Mechanics. Chap-
5 CONCLUSIONS man and Hall, London and Wiley, New York.
Yin, J-H. & Graham, J. 1996. Elastic visco-plastic modelling
of one-dimensional consolidation. Géotechnique 46(3):
This study has shown the significance of modelling the 515–527.
gradually changing creep behaviour as clay structure
80

Numge 2010

Uploaded by

Copyright:

Available Formats

Numge 2010

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Numge 2010

Uploaded by

Copyright:

Available Formats

cover_NUMGE2010_D.

qxd 29-04-2010 12:44 Pagina 1

Numerical Methods in Geotechnical Engineering

Numerical Methods in Geotechnical

No claim to original U.S. Government works

International Standard Book Number-13: 978-0-203-84236-2 (eBook - PDF)

A non-associated creep model for structured anisotropic clay (n-SAC) 3

Computer codes and algorithms

Discontinuum and particulate modelling

A Coupled Eulerian-Lagrangian approach to solve geotechnical problems involving

Flow and consolidation

Unsaturated soil mechanics

Reliability and probability analysis

Dynamic problems and Geohazards

Slopes and cuts

Effect of updated geometry in analyses of progressive failure 497

Embankments, shallow foundations, and settlements

Deep excavations and retaining walls

Tunnels and caverns

Ground improvement modelling

Offshore geotechnical engineering

Numerical methods and Eurocode

Author index 953

Thomas Benz and Steinar Nordal

Scientific Committee (ERTC 7)

Local Organizing Committee at Geotechnical Division NTNU, Norway

Reviewers not being member of ERTC 7 or the local organizing committee

Conference secretariat NTNU Videre

A non-associated creep model for structured anisotropic clay (n-SAC)

1 BACKGROUND big consequence of wrongly predicting the horizontal

where p = mean stress; σ d = deviatoric stress vec-

ζi , MC , αK0NC and ηK0NC are internal model param-

Isotropic hardening of the intrinsic reference stress,

Figure 1. concept of destructuration – effect on time resis-

The mobilization under virgin loading in oedome-

where λi and κ are alternative input parameters to

where rsi = the intrinsic time resistance number and

0.15 0.55 160 16 100 800

where Dn+1 is the mean stress dependent isotropic

Figure 2. Result of test 1 in p − q space.

Typically the scheme will converge in a few iterations.

The 21 state variables, including the effective stresses,

peak or at 6% shear strain) are 0.51 and 0.70 respec-

6 CONCLUSIONS AND RECOMMENDATIONS

This paper describes a constitutive model which

A state dependent constitutive model for sand-structure interfaces

ABSTRACT: Numerical simulations of geotechnical engineering problems depend strongly on predictive

1 INTRODUCTION later introduced by Gajo & Wood (1999). Hereto-

where εn and εt are normal and tangential strains with

where Rn plays the role of dilatancy, d, in the proposed

where Kn and Kt are respectively the interface normal

Table 1. Model parameters used in model predictions shown

Kt 0 (MPa) 1.57 5.0 3.0

*: experiments reported by Mortara et al. (2007).

In equation (16), e0 and λ are model parameters.

3.3 Dilatancy function and plastic hardening

where αb and αd are respectively the bounding and

where h0 is a model parameter and αin is the initial

Mortara et al. (2007) published the results of a