A time dependent constitutive model for soils with isotach viscosity

Computers and Geotechnics 38 (2011) 809–820

Contents lists available at ScienceDirect

Computers and Geotechnics journal homepage: www.elsevier.com/locate/compgeo

A time dependent constitutive model for soils with isotach viscosity Teresa M. Bodas Freitas a,⇑, David M. Potts b, Lidija Zdravkovic b a b

Instituto Superior Técnico, Departamento de Engenharia Civil e Arquitectura, Avenida Rovisco Pais, 1049-001 Lisboa, Portugal (formerly at Imperial College London) Imperial College London, Department of Civil & Environmental Engineering, Skempton Building, London SW7 2AZ, UK

a r t i c l e

i n f o

Article history: Received 5 September 2010 Received in revised form 18 May 2011 Accepted 19 May 2011 Available online 23 June 2011 Keywords: Constitutive relations Viscoplasticity Critical state line Time dependence Creep

a b s t r a c t This paper presents a three-dimensional elastic viscoplastic constitutive model that is able to reproduce the time dependent behaviour of soils with isotach viscosity. This constitutive law is based on the overstress theory and incorporates some important features, namely: (i) a non-linear creep law with a limit for the amount of creep deformation under isotropic stress conditions; (ii) a flexible loading surface that is capable of reproducing a wide range of shapes in p0  J stress space and incorporates the Matsuoka– Nakai failure criterion in the deviatoric stress space and (iii) assumes that the viscoplastic scalar multiplier is constant on a given loading surface to ensure that critical state conditions are reached. A full description of the model, its governing equations and implementation in a finite element program are presented. The model is then used to simulate the results of laboratory tests, highlighting the model’s abilities and shortcomings in reproducing the time dependent behaviour of soft clays. Furthermore, the paper investigates the importance of accounting for creep non-linearity and the consequences of isotach viscosity on the critical state line. Ó 2011 Elsevier Ltd. All rights reserved.

1. Introduction The study of the time dependent behaviour of soils was initiated in the 1920s, prompted by the practical need to quantify the longterm settlement of soft clays. As a result, the first mathematical formulations of the creep phenomenon were developed for onedimensional stress conditions. In this respect the works of Buisman [1], Bjerrum [2] and Garlanger [3] are of special note. Subsequent laboratory studies have focused on the identification and characterisation of the soil’s time and rate dependent behaviour under general stress conditions, namely the influence of the applied strain rate on the pre-consolidation pressure and on the size of the bounding surface [4–6], and the characterisation of undrained creep and the conditions for the occurrence of undrained creep rupture [7,8]. These laboratory studies highlight the importance of considering the time dependent behaviour of soils in the design of geotechnical structures and in the derivation of the strength and deformability parameters appropriate for field conditions, to ensure that both serviceability and ultimate limit states are verified [9]. Driven by the increasing use of numerical methods in geotechnical design, various constitutive models have been proposed to simulate the time dependent behaviour of soils that follow ⇑ Corresponding author. Tel.: +351 21 8418419; mobile: +351 96 6028382; fax: +351 21 8418427. E-mail address: [email protected] (T.M. Bodas Freitas). 0266-352X/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.compgeo.2011.05.008

isotach viscosity, under general stress conditions [10]. Due to its simplicity and flexibility, most of these are elastic viscoplastic (EVP) models based on Perzyna’s overstress theory [11]. According to Perzyna’s overstress theory the size of the viscoplastic strain increment is a function of the distance between the current loading surface and the static (yield) surface. This function – referred to in the literature as the viscous nucleus or viscoplastic scalar multiplier – and thus, the size of the viscoplastic strain increment, is implicitly constant for all stress states located on a given loading surface. The individual strain components are then obtained from a plastic potential that is coincident with the current loading surface. A significant difference between available EVP overstress-based models concerns the mathematical form of the viscoplastic scalar multiplier. In some earlier works, the viscoplastic scalar multiplier is assumed to be an exponential function of overstress and the function constants are obtained by fitting laboratory test results [12,13]. An alternative is to relate the viscoplastic scalar multiplier to the observed drained creep behaviour under one-dimensional or isotropic stress conditions [14–17]. This procedure simplifies considerably the derivation of the model parameters. But when applying this latter approach some researchers [15,16,18,19] assume that the volumetric strain rate is constant on a given loading surface, as suggested by some laboratory test data [20,21]. However such an assumption implies that critical state conditions cannot be reached [22]. The model presented herein, assumes that the viscoplastic scalar multiplier is constant on a given loading surface [13,14,17].

810

T.M. Bodas Freitas et al. / Computers and Geotechnics 38 (2011) 809–820

In most EVP overstress-based constitutive models the shape of the loading surface and plastic potential are fixed and mathematically described by either the original Cam Clay yield surface ([23] in [12,13]) or the elliptical shape of the Modified Cam Clay surface ([24] in [16,17]). In this case, once the angle of shearing resistance is set, the undrained shear strength under a given shearing mode is predefined. For many soils this results in a significant divergence between predicted and available undrained shear strength. To overcome this limitation, Yin et al. [18] proposed to describe the loading surface using two functions, one for the dry and the other for the wet side of the critical state, which can be adjusted through the model parameters, to ensure that both the soil’s drained and undrained shear strength are simultaneously well reproduced. Although the approach has some merit, the existence of two functions increases computing time and significantly complicates the derivation of the model parameters and the model implementation. In the model presented herein, the loading surface and plastic potential are described by the mathematical function proposed by Lagioia et al. [25] that allows the simulation of a wide range of shapes without the drawbacks mentioned above. The model described in this paper – termed hereafter the ET model – incorporates a hyperbolic law to describe the time delayed compression of soils under isotropic stress conditions as proposed by Yin [26]. In contrast with a semi-logarithmic creep law, widely used in engineering practice and in constitutive modelling, a hyperbolic creep law incorporates a limit to the amount of volumetric drained creep (that is specified by the user) and predicts that the evolution of volumetric deformation with time, during a drained creep period, plots as a non-linear curve in semi-logarithmic space, in agreement with experimental observations made both in the laboratory [6,26] and in the field [27]. As a rule, the stress–strain response of natural clays is subjected to persistent strain rate effects up to failure, such that their behaviour can be well described by a unique stress–strain–strain rate relationship. This behaviour is termed in the literature as isotach viscosity (e.g. [28]). Concurrently the undrained shear strength, Su, is found to increase with the applied strain rate; however the effective stress failure envelope at large deformations (i.e. the angle of shearing resistance at critical state) appears to be unique and time independent [8,29–31]. The fact that the undrained strength is rate dependent, associated with creep driven pore water pressures, suggests that the critical state line (CSL) is not unique in terms of void ratio, since for a given void ratio an infinite number of failure stress states can be obtained during undrained shearing, depending on the applied strain rate. On the other hand, it may be that such tests have been terminated prior to reaching critical state conditions and the observed dependency of the undrained strength with strain rate corresponds instead to the rate dependency of the bounding surface. This raises an issue regarding the uniqueness of the CSL in both effective stress space and in terms of void ratio. Whereas the first is well corroborated by laboratory data published in the literature, there is little information regarding the latter. Sorensen [32] has carried out an extensive laboratory testing programme on natural and reconstituted London clay and has found that normally consolidated reconstituted (NC) London Clay has a unique strain rate independent CSL both in stress space and in terms of void ratio. These observations are consistent with the fact that NC reconstituted London Clay shows no persistent strain rate effects at high stress levels (i.e. if during shearing an increase of strain rate is applied, the soil will initially show an increase in the mobilised stress, however with further straining the stress–strain response will tend to a unique rate-independent response; this behaviour is referred to in the literature as TESRA – Temporary Effects of Strain Rate Acceleration, see e.g. [28]). On the other hand, intact London Clay, which has isotach viscosity up to failure, is reported to have a un-

ique CSL in effective stress space, but the uniqueness of the CSL in terms of void ratio was not established [32]. Other research that may be relevant to this topic is that on Champlain Sea Clay from Quebec [41], that shows that the CSL is temperature dependent in void ratio–mean effective stress space, but not in q  p0 space. The relevance of this work is that temperature and strain rate effects are two aspects of the soil’s viscous nature and thus the same type of viscosity is expected to be observed for the two variables. In fact, the effects of both temperature and strain rate on the soil’s stress–strain response are found to be persistent up to failure, i.e. the Champlain Sea Clay has isotach viscosity. This is consistent with the findings regarding the effect of temperature on the CSL, since for a given void ratio an infinite number of failure stress states are possible and thus the CSL is not unique in void ratio–mean effective stress space. One would expect to obtain the same overall result when considering the effect of strain rate on the CSL. In the light of the above, the ET model is used to simulate a series of undrained triaxial compression tests on normally consolidated and overconsolidated samples to investigate the predictive capabilities of the ET model under undrained shearing and the implications of the overstress theory on the CSL. 2. General model framework The constitutive model is mostly presented in terms of the individual stress and strain components, using indicial notation. However, sometimes it is convenient to separate the model response into its volumetric and deviatoric components, and stress and strain invariants are used. The stress invariants used herein are the mean effective stress, p0 , the deviatoric stress, J, and the Lode’s angle, h, given by Eqs. (1)–(3), respectively.

p0 ¼

J¼

1 0 ðr þ r0y þ r0z Þ 3 x

ð1Þ

 1=2 1 ðsij : sij Þ 2

2 pffiffiffi 1 1 63 3 h ¼ sin 4 h 3 2

ð2Þ 3 detðsij Þ ð12 ðsij : sij ÞÞ1=2

7 i3 5

ð3Þ

where sij is the deviatoric stress tensor, its components being calculated as:

sij ¼ r0ij  p0  dij

ð4Þ

where r0ij is the effective stress tensor and dij is the Kronecker’s delta. The corresponding strain invariants are the volumetric strain, evol, and the deviatoric strain, Ed, given by Eqs. (5) and (6), respectively.

ev ol ¼ ex þ ey þ ez   ev ol 2  ev ol 2 þ ey  Ed ¼ 2 ex  3 3 1=2   2 ev ol 1 1 1 þ c2xy þ c2yz þ c2zx þ ez  2 2 2 3

ð5Þ

ð6Þ

Compressive stresses and strains are positive. It is assumed that the soil deformation, associated with an effective stress increment, {Dr0 }, over a time increment, Dt, is decomposed into an elastic and a viscoplastic part:

fDeT g ¼ fDeel g þ fDev p g

ð7Þ

T.M. Bodas Freitas et al. / Computers and Geotechnics 38 (2011) 809–820

where the elastic strain increment, {Deel}, is instantaneous and thus time independent and the viscoplastic strain increment, {Devp}, is time dependent and irreversible. The elastic strain increment, {Deel}, can be determined by inverting Eq. (8).

fDr0 g ¼ ½D0   fDeel g

ð8Þ

where [D0 ] is the elastic constitutive matrix. The elastic response is assumed to be isotropic and is fully characterised by two elastic parameters, a stress dependent bulk modulus (defined by Eq. (9)) and a second elastic parameter that can be either the Poisson’s ratio, m, or the elastic shear modulus, G.

K¼

V  p0

ð9Þ

j

According to the viscoplastic theory proposed by Perzyna [11], the viscoplastic strain increment can be evaluated as:

fDev p g ¼ fe_ v p g  Dt ¼ c  hUðFÞi 



@fd @ r0

  Dt

ð10Þ

where hUðFÞi ¼ UðFÞ if F > 0

fe_ v p g is the viscoplastic strain rate tensor; c is a fluidity parameter and U(F) is a function of the overstress, F, the quantity F being defined as the normalised distance between the current dynamic loading surface, fd, and the static loading surface, fs, which defines the region of time-independent and purely elastic behaviour (i.e. a yield surface). In the present case, U is a function of the distance of the current loading surface to both a reference loading surface (characterised by a finite value of viscoplastic strain rate) and the static loading surface (or yield surface), and viscoplastic strains are predicted provided the stress state lies above the yield surface. However, the significance of the yield surface in this formulation is different to that in other critical state based elasto-plastic models, e.g. Modified Cam Clay, where it is associated with large scale yielding and the normally consolidated (NC) state. Herein NC states are located on the reference loading surface and the yield surface, being inside the reference loading surface, is the locus of a set of heavily overconsolidated soil states, with the associated degree of overconsolidation depending on the value adopted for the model parameters. This ensures that time and rate dependent behaviour is predicted for NC and OC states, provided the stress states lie above the static yield surface. Some, if not most, overstress based EVP models published in the literature [10] do not include a yield surface and in that case viscoplastic strains are always predicted. The introduction of a yield surface, and thus a region of pure elastic behaviour, implies that there is a limit for the amount of drained creep, which is consistent with the adopted hyperbolic creep law. The viscoplastic strain increment is calculated as:

fDev p g ¼ hUi 



@g @ r0



 Dt

3. The concept of equivalent time The concept of equivalent time was introduced by Yin and Graham [33] to overcome the problems that arise when the soil’s delayed compression is related to real loading time. Fig. 1 shows schematically the model framework under isotropic stress conditions in volumetric strain–mean effective stress space, evol  ln p0 . Equivalent time, te, is defined as the time that a soil element would take to creep from a reference compression line to the current stress state under constant mean effective stress [33]. The reference compression line (or reference time line), corresponds to the locus of the soil states with zero equivalent time. In practice, it corresponds to the soil’s isotropic normal compression line (NCL). Soil states located below the reference line have a positive equivalent time, between 0 and infinity. Above the reference compression line, the equivalent time varies between zero and the value t0, where t0 is a model parameter and corresponds to the real time associated with the reference compression line. Noting that the strain origin is arbitrary, the reference compression line is assumed to pass through the point [evol = 0; p0 = 1 kPa] and is then given by the following equation:

eref v ol;m ¼

hUðFÞi ¼ 0 if F 6 0

811

0

k pm ln V 1 kPa

ð12Þ

where the superscript ref denotes a quantity evaluated on the reference compression line; the subscript m denotes an isotropic stress 0 0 state; eref v ol;m is the volumetric strain at p = pm ; and k=V is a model constant. An equivalent time is related to a unique viscoplastic (creep) strain rate, with large equivalent times being associated with smaller viscoplastic strain rates. The stress–strain–equivalent time relationship (or volumetric viscoplastic strain rate) is unique and independent of the loading history. The projection of the yield surface in evol  ln p0 space is a line parallel to the reference compression line, at a vertical distance evv pol;m;Limit . This line, also termed the limit time line [33], is the locus of the soil states that are attained after an infinite drained creep period, being characterised by an equivalent time equal to infinity and a viscoplastic strain rate equal to zero. If a linear logarithmic creep law is adopted there is no limit to the amount of volumetric creep deformation and evv pol;m;Limit is equal to infinity. According to the model formulation, when a soil element is subjected to a mean effective stress increment, instantaneously, it experiences purely elastic deformations. Graphically, in evol  ln p0

ð11Þ

where hUi ¼ U if the stress state lies above the yield surface hUi ¼ 0 if the stress state lies on or below the yield surface and g is the plastic potential, in principle different from the current loading surface, fd. The value of the viscoplastic scalar multiplier, U, is subsequently derived.

Fig. 1. Schematic representation of the model framework under isotropic stress conditions.

812

T.M. Bodas Freitas et al. / Computers and Geotechnics 38 (2011) 809–820

space, the soil state moves along a line with inclination j/V passing through the current soil state. Yin and Graham [33] termed such lines as instant time lines. 4. Model description 4.1. Constitutive law for isotropic stress conditions Based on Fig. 1, the volumetric strain, evol,m, of a given soil element under a mean effective stress, p0m , can be written as: vp ev ol;m ¼ eref v ol;m þ ev ol;m

ð13Þ

where eref v ol;m is the volumetric strain on the reference compression line at the current mean effective stress p0m and evv pol;m is the viscoplastic volumetric strain. Assuming that the soil time delayed deformation is described by the non-linear logarithmic law proposed by Yin [26], evv pol;m is then:

evv pol;m ¼

w0 V

ln



1 þ Vev pw0

t 0 þt e t0

ln

v ol;m;Limit





t 0 þte t0



ð14Þ

where w0/V is a model parameter, te is the current equivalent time and t0 is the real time associated with the reference NCL, conventionally taken as 1.0 day. Based on Eqs. (13) and (14) the volumetric viscoplastic strain rate, e_ vv pol;m , can be shown to be given by Eq. (15) [34]:

3

2

e_ vv pol;m ¼

eref w0 v ol;m  ev ol;m  1þ V  t0 evv pol;m;Limit

!2

6 6V  exp 6 6w 4 0

7 eref 7 v ol;m  ev ol;m 7 ! ref ev ol;m  ev ol;m 7 5 1þ evv pol;m;Limit ð15Þ

4.2. Extension to general stress space Consider a general state (p0 , J, evol) above the yield surface as represented in Fig. 2. The stress state (p0 = p0i ; J = Ji) is on a loading surface characterised by the quantity p0m – the mean effective stress at J = 0. Based on Eq. (11) the incremental viscoplastic volumetric strain is:

Devv pol ¼

@g U  Dt  0 @p

ð16Þ

Combining Eqs. (15) and (16) and assuming that U is constant on a given loading surface [22], the value of U can be written as: 3

2

U¼

e w0  1þ V  t0

1     @g  @p0 p0 ¼

ref v ol;m  v ol;m vp v ol;m;Limit

e

e

!2

6 6V  exp 6 6w 4 0

7 ref 7 v ol;m  v ol;m !7 7 ref  5 v ol;m v ol;m þ vp v ol;m;Limit

e

1

e

e

e

Fig. 2. Model framework in general stress space.

ev ol;m ¼ ev ol þ

p0m J ¼ 0

   @g  where @p is the absolute value of the partial derivative of the 0 p0 ¼p0 m

V

ln



p0m p0i

ð18Þ

where evol is the current total volumetric strain. The geometric significance of the quantity evol,m is illustrated in Fig. 2b. The problem of determining the viscoplastic strain increment at a general state (p0 , J, evol) is reduced to evaluating U at the equivalent isotropic state (p0 = p0m , J = 0, evol,m). With reference to Fig. 2b, the condition for a stress state to be located above the yield surface can be mathematically expressed as:

eref v ol;m  ev ol;m >0 evv pol;m;Limit

1þ

ð19Þ

Eq. (20) describes the soil time dependent behaviour under general stress conditions, for stress states above the yield surface. Given its incremental form it is suitable for direct implementation in a finite element code. For stress states below the yield surface Eq. (20) can still be used by setting U equal to zero, and only elastic deformations are predicted.

e

ð17Þ

j

T ij g

fDe

1 j Dp0  Dsij þ ¼   dij þ U  2G V  p0 3

(

@g @ r0ij

)  Dt

ð20Þ

Where is given by Eq. (17) 4.3. Plastic potential and loading surface

J¼0

plastic potential function, g, in relation to the mean effective stress, p0 , evaluated at the equivalent isotropic stress state (p0 ¼ p0m , J = 0), the absolute value function being introduced to ensure that the viscoplastic scalar multiplier is always a positive quantity; eref v ol;m is the volumetric strain on the reference NCL at p0 ¼ p0m and evol,m is the volumetric strain on the current elastic line at p0 ¼ p0m , being calculated as:

In general stress space, the loading and the plastic potential surfaces are described by Eq. (21).

  K2 g ð1lÞðK 1 K 2 Þ 1 þ K2 p f or g ¼ 0   ¼0  K1 pm  ð1lÞðK 1 K 2 Þ 1 þ Kg1 0

ð21Þ

T.M. Bodas Freitas et al. / Computers and Geotechnics 38 (2011) 809–820

where p0m is the mean effective stress at J = 0, K1 and K2 are constants determined as:

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi!

l  ð1  aÞ 4a  ð1  lÞ K 1=2 ¼  1 1 2  ð1  lÞ l  ð1  aÞ2

ð22Þ

g ¼

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi J 2g J 2g;failure

ð23Þ

J2g is the stress ratio squared i.e. J2g ¼ ðJ=p0 Þ2 and J2g,failure is the stress ratio at failure squared, which is calculated using Eq. (24) assuming the Matsuoka–Nakai failure criterion in the deviatoric plane.

2 pffiffiffiffiffiffi  C  sinð3hÞ  ðJ 2g;failure Þ3=2 þ ðC  3Þ  J 2g;failure  ðC  9Þ 27 ¼0

9  M2 2M3 27

2

 M3 þ 1

6 sin /0 3  sin /0

where Dex, Dey and Dez are the strain components in the x, y and z direction, respectively. In addition, it is also of interest to know the value of specific volume (or void ratio) during the analysis procedure, as it is a common state parameter and it may be required by other aspects of the analysis (e.g. if used in combination with a void ratio dependent permeability model). The specific volume is initiated as Vi = 1 + ei, where ei is the prescribed initial void ratio. The change in the current specific volume is equal to DV = De, where De is the change in void ratio. 4.5. Model implementation

ð25Þ

with M being the gradient of thepcritical state line (CSL) in p0  q ffiffiffi stress space (deviator stress q ¼ 3  J) under triaxial compression (h = 30°), which is related to the angle of shearing resistance, /0 as follows:

M¼

ð26Þ

Eq. (21) was proposed by Lagioia et al. [25] to describe the yield surface and the plastic potential for use in constitutive modelling of soils, allowing the simulation of a wide range of shapes (e.g. the Modified Cam Clay model ellipse is recovered by setting a = 0.4 and l = 0.9). This equation was obtained by integrating a postulated general variation of dilatancy, d (the ratio of volumetric to deviatoric plastic strain increments) with stress ratio g (g = J/p0 ) and ensures that at critical state the dilatancy is zero and under isotropic stress conditions (g = 0) the dilatancy is infinity (i.e. only incremental volumetric plastic strains are predicted). Four parameters are required to fully define this function; M, l and a are model parameters while p0m defines the size of a particular surface and depends on the current stress state. For more details regarding the plastic potential and the loading surface function being used see Lagioia et al. [25].

The model described above was implemented in the Imperial College Finite Element Program, ICFEP [35]. A modified Newton– Raphson scheme, with an error controlled sub-stepping algorithm was used as the non-linear solver. To achieve this implementation it was necessary to establish the relationship between the incremental stresses and strains and to determine the contribution of the viscoplastic strains to the governing finite element equations. 4.5.1. Relationship between incremental stresses and strains This relationship is required in the error controlled substepping stress point algorithm where stress changes must be calculated from strain steps. Noting that the incremental elastic strains are related to the incremental stresses by the elastic constitutive matrix [D0 ], then Eqs. (7) and (11) can be combined to give:

fDr0 g ¼ ½D0  

4.4. Hardening parameters

ð29Þ

4.5.2. Contribution to governing equations Following a similar procedure to that described by Potts and Zdravkovic [35] the following governing finite element equation for describing the mechanical behaviour can be derived by invoking the principle of minimum potential energy. N N X X ½Ki  ðfDdgn Þi ¼ ðfDRgn Þ

ð30Þ

i¼1

where

The application of Eq. (20) to describe the behaviour of a soil element implies the knowledge of the current total volumetric strain, evol, at any instant, t. The current total volumetric strain, evol, is initialised at the start of the analysis and its value is continuously updated during the analysis process. The initial volumetric strain, einitial v ol , can be evaluated as:

0

0

k pmc j p ln  ln mc V 1 kPa V p0i

 

T @g De  hUi  Dt  @ r0

which is the required relationship.

i¼1

einitial v ol ¼

ð28Þ

ð24Þ

where

C¼

function of the applied strain rate but for simplicity, the OCR is defined in relation to the reference compression line. During the analysis procedure the change in volumetric total strain is calculated as:

Dev ol ¼ Dex þ Dey þ Dez

 is the generalised norwhere l and a are model parameters and g malised stress ratio:

813

ð27Þ

where p0i is the initial mean effective stress and p0mc is the size of the loading surface corresponding to the largest NC stress state that the soil element has experienced, and it can be determined from the initial stress state, the overconsolidation ratio (OCR) and the model parameters. The time dependent nature of the model implies that the pre-consolidation pressure or yield vertical effective stress is a

R ½K ¼ v ol ½BT ½D½B  dVol is the element stiffness matrix; {Dd}n is the vector of incremental nodal displacements; R R R fDRgn ¼ Vol ½NT fDFg  dVol þ Vol ½BT  fDrc g  dVol þ Srf ½NT fDTg dSrf is the right hand side load vector; fDrc g ¼ 12  hUi  Dt  ½D  f@@grg is the contribution from viscoplastic strains; [N] is the shape functions matrix; [B] is the shape functions derivatives matrix; {DF} is the body forces vector; {DT} is the surface tractions vector; Srf is that part of an element boundary over which the surface tractions are applied and N is the number of finite elements. This equation differs from that for elastic (or elasto-plastic) material behaviour in that the viscoplastic strains contribute to the right hand side load vector.

814

T.M. Bodas Freitas et al. / Computers and Geotechnics 38 (2011) 809–820

Table 1 Model parameters. Model Units parameter

Definition

Aspect of soil behaviour

k=V

–

Slope of the reference time line in evol  ln p0 space.

Large scale compressibility

j/V

–

Elastic response

m or G

–/ kPa

Slope of the elastic line in evol  ln p0 space. Second elastic parameter, either the Poisson’s ratio l or the elastic shear modulus G

a l

Loading surface parameter Loading surface parameter q/p0 value at CS, under triaxial compression

Loading surface in p0  J space

M

– – –

w0/V

–

Slope of the evol  ln t curve at t = t0 = 1 day Limit to the amount of volumetric viscoplastic strain

Creep law under isotropic stress conditions

evv pol;m;Limit –

the quantities k=V and j/V are the slope of the reference NCL and the elastic line in evol  ln p0 space, under isotropic stress conditions. However, noting that in the NC range the coefficient of earth pressure at rest, K0, remains constant, and at low OCR values the error introduced by such approximation is not significant, the above parameters can be estimated under one-dimensional (1D) stress conditions, from oedometer test data (either constant strain rate or incremental load type). The parameter k=V is the slope of the compression curve in evol  ln p0 space, in the NC range and the parameter j/V should be determined as the slope of an unloading–reloading loop. It is unadvisable to use the recompression data at the start of the test to determine j/V, as at low stress levels the laboratory data is often affected by the apparatus compliance and sample disturbance. Based on laboratory testing data [36,37] the parameters j/V and k=V for reconstituted Carse clay have been derived. Fig. 3 shows the model predictions of a 24 hour incremental load (IL) oedometer test obtained by means of a coupled finite element analysis (FEA) using the parameters included in Table 2, in which the combined effects of loading and consolidation are taken into account. 5.2. Second elastic parameter The elastic part of the model is characterised by a stress dependent bulk modulus, K, defined by Eq. (9) and a second elastic parameter, either the Poisson’s ratio, m, or the elastic shear modulus, G. The elastic shear modulus, G, can be estimated from the stress–strain data of undrained triaxial tests, at small strains. 5.3. Time dependent parameters w0/V and evv pol;m;Limit

Fig. 3. Oedometer tests on reconstituted Carse Clay; laboratory data and numerical predictions (Note: R = reconstituted sample).

5. Model parameters Assuming that the plastic potential coincides with the loading surface (i.e. f = g), the ET model requires 8 input parameters, which are listed in Table 1. The model parameters can be derived from common laboratory tests: ideally two consolidated undrained triaxial tests (one in compression and one in extension) and an isotropic compression test, preferably with a long drained creep period. As discussed below, for practical purposes, the isotropic compression test can be replaced by an oedometer test, the latter being more often available in engineering practice. In addition to the model parameters in Table 1, one needs to specify the initial stress conditions, the initial void ratio, e0 and the OCR.

The parameters w0/V and evv pol;m;Limit together with the parameter t0 (conventionally set as 1.0 day) define the soil’s creep behaviour under isotropic stress conditions. However, the time dependent volumetric deformations of soils under isotropic and 1D stress conditions, as observed experimentally [21] and predicted by elastic– viscoplastic models [22], are not that dissimilar and the above parameters can be derived from an oedometer test. The parameter w0/V is the slope of the evol  ln t curve, obtained during a load increment in the NC range, at a time t = t0. w0/V can be related to the coefficient of secondary compression, Cae, through Eq. (31).

w0 C ae ffi V V  ln 10

plastic strain that is reached after an infinite drained creep period under constant effective stress, and should be derived by curve fitting a long creep test. The required length of such a creep test is difficult to ascertain, as it is likely to depend on the particular soil being tested; as a rule it is suggested that evv pol;m;Limit should be derived from a creep test, at least, 10 to 100 times t0 long. Yin et al. [18] suggest that in the absence of a long creep test the value of evv pol;m;Limit could be taken as:

evv pol;m;Limit ¼

5.1. Parameters k=V and j/V The ET model is formulated in evol  ln p0 space, assuming that k=V and j/V are model parameters and soil constants. By definition

ð31Þ

The parameter evv pol;m;Limit is by definition the volumetric visco-

e0 1 þ e0

ð32Þ

where e0 is the initial void ratio. The value given by Eq. (32) corresponds to the volumetric strain required to reach the situation when voids no longer exist in the soil.

Table 2 Model parameters for Carse Clay. e0 (–)

V (–)

k/V (–)

j/V (–)

m (–)

w0/V (–)

evv pol;m;Limit (–)

a (–)

l (–)

k (m/s)

1.85

1.462

0.0842

0.0105

0.26

0.00378

0.065

0.05

0.85

5  109

T.M. Bodas Freitas et al. / Computers and Geotechnics 38 (2011) 809–820

815

Fig. 4. Drained creep tests: laboratory data (from [37]) and model predictions. (L = Laval sample; S = Sherbrooke sample).

Fig. 6. Variation of the yield stress ratio with strain rate: laboratory test data (from [36]) and model predictions (Note: for IL oedometer tests the strain rate is that observed at the end of the load increment).

Fig. 4 shows the results of two identical creep tests carried out on good quality natural samples, in which the samples were consolidated in the triaxial apparatus at 1% axial strain per day, while maintaining K0 stress conditions, up to a mean effective stress of 100 kPa, and then allowed to creep at constant effective stress [37]. Also shown in Fig. 4 are model predictions obtained by means of coupled FEA using the parameters in Table 2, that illustrate the influence of the parameter evv pol;m;Limit . If evv pol;m;Limit is set to the value given by Eq. (32) then the ET model predictions coincide with a linear logarithmic law for the time interval considered (and for time intervals of interest to most engineering problems); on the other hand the adoption of a lower value of evv pol;m;Limit (say 10% of the value given by Eq. (32)) increases the

creep law non-linearity, its influence becoming more significant for larger creep times. To model the soil creep non-linearity that is observed experimentally [6,26,27] the parameter evv pol;m;Limit needs to be given a value significantly lower than that given by Eq. (32). These results highlight the importance of having experimental creep data for periods well in excess of 1.0 day if a good prediction of the soil’s long-term behaviour (or of a geotechnical structure) is to be recovered. A linear logarithmic creep law predicts that, following excess pore water pressure dissipation, the soil compression under maintained load in evol  ln t space follows a straight line, with slope equal to w0/V. In this case the vertical spacing in evol  ln p0 space, between any two lines of constant equivalent time te,i and te,j such

Fig. 5. Consolidated undrained triaxial tests on intact Carse Clay: laboratory data (from [37]) and model predictions (L = Laval sample; S = Sherbrooke sample).

816

T.M. Bodas Freitas et al. / Computers and Geotechnics 38 (2011) 809–820

Table 3 Model parameters for HKMD. e0 (–)

V (–)

k/V (–)

j/V (–)

m (–)

w0/V (–)

evv pol;m;Limit (–)

a (–)

l (–)

1.504

1.265

0.0793

0.0025

0.3

0.00521

0.6

0.20

0.56

Fig. 7. Undrained triaxial compression test with stepwise change in strain rate combined with stress relaxation on reconstituted HKMD – laboratory test data (from [18]) and numerical simulations.

that t0 + te,j = 10  (t0 + te,i) (i.e. e_ vv pol;i ¼ 10  e_ vv pol;j ), is constant and equal to ln 10  w0/V. If a hyperbolic law is used, the slope of the compression curve decreases with time from the value w0/V at te = 0 to zero at te = infinity and so the spacing between the same two lines of constant equivalent time will decrease for larger values of te (i.e. away from the reference compression line). Examination of Eq. (20) shows that if w0/V is set equal to zero, then U = 0 and no viscoplastic strains are predicted – i.e. pure elastic behaviour. On the other hand, if w0/V is set to a very small value, say 0.001, then the response of the model is almost unaffected by changes to the applied strain rate, as the distance between compression lines corresponding to different applied strain rates is proportional to w0/V, and thus very small. In evol  ln p0 space, the soil state follows closely the reference compression line independently of the applied strain rate and, for practical purposes, the model response may be then taken as time-independent. 5.4. Parameters M, a and l The parameter M is the stress ratio q/p0 at critical state under triaxial compression (i.e. h = 30°) and can be related to the angle

of shearing resistance, /0 , through Eq. (26). The parameters a and l are determined by curve fitting the effective stress path of undrained triaxial tests, aiming to reproduce the soil’s undrained shear strength, Su. Given that the shape of the failure surface and plastic potential in the deviatoric plane is fixed, the variation of the undrained shear strength with the intermediate principal stress (i.e. with shearing mode) cannot be controlled by the user. If the parameters a and l are calibrated to recover Su in triaxial compression then the Su under other shearing modes (different h values) is set. Comparison of laboratory data and model predictions for other shearing modes (e.g. triaxial extension) will enable the user to assess the error involved. It is recommended that the parameters a and l are derived to recover the undrained shear strength under the shearing mode likely to dominate the problem being analysed (e.g. triaxial extension for an excavation problem). Fig. 5 shows the effective stress paths obtained during five undrained triaxial compression and extension tests on intact samples of Carse clay from 5.3 m and 6.2 m depth. The samples were anisotropically consolidated to the in situ effective stresses following a swelling loop to mimic the recent stress history, and were then sheared at 1% or 4.5% axial strain per day. The parameters a and l were derived by curve fitting the effective stress path of sample LCU1. Also shown in Fig. 5 are the predicted effective stress paths for the remaining samples, obtained by means of FEA that simulate the complete loading sequence and use the parameters given in Table 2. The proposed constitutive model reproduces reasonably well the peak undrained shear strength under triaxial compression and its variation with applied axial strain rate; however it is unable to reproduce the full stress–strain response, in particular the post peak strength loss due to destructuration. Because the parameters a and l were derived to match Su under triaxial compression, significant divergence is observed between measured and predicted values of Su under triaxial extension. Fig. 6 shows the variation of the yield stress ratio (YSR) with the applied axial strain rate for intact Carse Clay, determined from laboratory oedometer tests and predicted by the ET model. The model predictions were obtained by means of single element FEA using the model parameters given in Table 2. The ET model reproduces well the general trend of the variation of the YSR with strain rate, the scatter in the laboratory data being associated with the sensitivity of Carse Clay to sample disturbance during sampling and handling [38].

6. Numerical simulation of laboratory tests on hong-kong marine deposits The strengths and limitations of the ET model are further investigated in a series of FEA in which the model is employed to reproduce laboratory tests on reconstituted samples of Hong Kong Marine Deposits (HKMD) using the model parameters given in Table 3. For further details on HKMD and testing procedures see Yin et al. [18] and Zhu et al. [30]. Fig. 7 shows the laboratory test results and the numerical predictions of an undrained triaxial compression test on a NC sample, isotropically consolidated to 300 kPa, during which the applied axial strain rate is changed in a stepwise manner and combined with stress relaxation periods.

T.M. Bodas Freitas et al. / Computers and Geotechnics 38 (2011) 809–820

Fig. 8. Undrained triaxial creep tests on reconstituted HKMD – laboratory test data (from [18]) and numerical simulations. Variation of axial strain with time: (a) in logarithmic scale and (b) in natural scale.

Fig. 9. CRS undrained triaxial compression test on NC reconstituted HKMD – laboratory test data (from [18]) and numerical simulations.

817

Fig. 10. Numerical predictions of CRS undrained triaxial compression tests on NC samples: (a) effective stress path and (b) variation of deviatoric stress and excess pore water pressure during shearing.

The ET model predictions approximate well the experimental data. The ET model is capable of simulating the variation of the sustained deviator stress, q, with the applied strain rate during the shearing stages and the decay in deviator stress during the stress relaxation periods. The model is able to reproduce the permanent effects of strain rate, as it is based on a unique stress–strain–viscoplastic strain rate relationship. Fig. 8 shows the results of three undrained triaxial creep tests on reconstituted samples of HKMD in terms of the variation of the axial strain with time, both in semi-logarithmic and natural time scales. The samples were initially isotropically consolidated to a mean effective stress of 400 kPa, and then loaded instantaneously with a deviator stress, q, equal to 134, 189 and 243 kPa, respectively, and restrained from undergoing any volume change while the applied stresses were kept constant. Fig. 8a suggests that the ET model recovers reasonably well the laboratory test results for q = 134 kPa and q = 189 kPa and underestimates the time that it took the sample at q = 243 kPa to fail. However, when the data is plotted in natural scale (Fig. 8b), it becomes apparent that, after an initial jump in axial strain associated with the application of the deviator stress increment, the ET model

818

T.M. Bodas Freitas et al. / Computers and Geotechnics 38 (2011) 809–820

Table 4 Model parameters for evaluation of the significance of the CSL in EVP models. e0 (–)

V (–)

k/V (–)

j/V (–)

G (kPa)

w0/V (–)

evv pol;m;Limit (–)

a (–)

l (–)

1.5

1.2872

0.088

0.0088

1700

0.00521

0.6

0.4

0.9

Fig. 11. Numerical predictions of CRS undrained triaxial compression tests on NC and OC samples: (a) effective stress path and (b) variation of deviatoric stress and excess pore water pressure during shearing for tests sheared at 10% axial strain per day.

predicts that the strain rate (i.e. the inclination of the axial strain– time curve) decreases with time in all cases, and thus it cannot reproduce accelerating creep and undrained creep rupture, which are observed in the laboratory data on the sample under q = 243 kPa. As noted by Adachi et al. [39] constitutive models based on the overstress theory are unable to simulate undrained creep rupture (associated with an increase in creep strain rate), because in their formulation an increase in the viscoplastic (creep) strain rate leads to an increase in the soil’s undrained strength. During an undrained creep test the applied deviatoric stress remains constant throughout the test and at failure an increase in the axial strain rate would imply an increase in the soil’s undrained strength, which would inhibit failure. This has important implications in the application of this type of model in the analysis of boundary value problems in conditions close to failure [34,40]. Fig. 9 shows the stress path of a constant rate of strain (CRS) undrained triaxial compression test sheared at 1.5% axial strain per hour, on a NC reconstituted sample of HKMD isotropically consolidated to 100 kPa, and the ET model predictions using the parameters given in Table 3. In addition, Fig. 9 shows the numerical results given by the ET model, when the parameters a and l are set equal to 0.4 and 0.9, respectively, so that the loading surface and plastic surface have the elliptical shape of the MCC model, and the remaining parameters are left unchanged. The data in Fig. 9 illustrates the error that models with a fixed loading surface shape in q  p0 space may incur, as the drained and undrained soil strength cannot be controlled separately.

7. The critical state line This section investigates the predictive capabilities of EVP models during undrained shearing under triaxial compression. Fig. 10a) shows the stress paths predicted by the ET model during CRS undrained triaxial compression tests on NC samples, isotropically consolidated to 100 kPa. The data has been obtained from single

Fig. 12. State paths predicted by the ET model during undrained triaxial compression tests on NC and OC samples in evol,m  ln p0m space.

element FEA using the model parameters given in Table 4. Fig. 10b shows the variation of the deviatoric stress and excess pore water pressure during the analyses. NC samples are predicted to have contractive behaviour, with p0 decreasing as J increases monotonically up to the CSL. The ET model predicts an increase in the undrained shear strength for higher values of applied axial strain rate. However, the effective stress critical state envelope remains unchanged, as the differences in the predicted undrained shear strength result from differences in the predicted excess pore water pressure. A further series of FEA analyses was carried out to obtain numerical predictions of CRS undrained triaxial compression tests on overconsolidated samples. Fig. 11 shows the results of the numerical simulations, normalised by the initial pre-consolidation pressure, equal to 100, 150 and 400 kPa, for the samples with an OCR equal to 1.0, 1.5 and 4.0, respectively. The ET model recovers the undrained shear strength strain rate dependency of both NC and OC samples. For heavily overconsolidated samples (OCR = 4.0) the ET model predicts an overall dilative

T.M. Bodas Freitas et al. / Computers and Geotechnics 38 (2011) 809–820

behaviour, with the generation of negative pore water pressures. The behaviour of lightly overconsolidated samples (OCR = 1.5) depends on the applied strain rate; at lower strain rates the behaviour is contractive throughout while for higher strain rates dilative behaviour is predicted beyond the pseudo-elastic stress range, showing that an increase in applied strain rate is equivalent to an increase in the sample overconsolidation ratio. This is in agreement with the fact that the observed yield stress ratio (in the laboratory and in the field), and thus the size of the apparent elastic stress range, depends on the applied strain rate (increases). Fig. 12 shows the state paths followed during the above CRS undrained triaxial compression tests in evol,m  ln p0m space, where p0m is the size of the current loading surface and evol,m is the associated volumetric strain, calculated using Eq. (18). As noted previously, the distance of the equivalent isotropic state (evol,m, p0m ) to the reference compression line defines the magnitude of U (and thus of the viscoplastic strain increment). Initially, the soil states follow closely the elastic line that passes through the initial state, tending gradually to the NCL line appropriate to the applied strain rate. At critical state (at the end of shearing), soil states for samples sheared at the same applied axial strain rate are found to define a line parallel to the reference compression line, further to the right for higher values of strain rate. The ET model assumes that the CSL is unique in stress space (Figs. 10 and 11) but is rate dependent in void ratio–mean effective stress space (Fig. 12), noting that all the samples had the same initial void ratio and the relationship between void ratio and volumetric strain is direct. Inspection of Fig. 12 shows that the vertical distance between critical state lines per logarithmic cycle of strain rate is about 0.01 volumetric strain (i.e. 0.025 in void ratio in this case). This value is very small, comparable to the error encountered in void ratio laboratory measurements, suggesting that the experimental verification of the CSL strain rate dependency is pending until higher accuracy in laboratory testing is achieved. The same time dependency of the CSL in void ratio–mean effective stress space is predicted by the ET model under drained shearing [34]. 8. Concluding remarks – This paper presents the development, implementation and validation of an elastic viscoplastic model based on the overstress theory that aims to simulate isotach viscosity. The model requires the definition of eight parameters that can be easily determined from common laboratory tests, the derivation of which has been described in detail. – The model has been used to simulate a series of laboratory tests on natural and reconstituted clayey samples. Based on the results of these numerical simulations it can be concluded that the developed model reproduces well persistent strain rate effects and thus is adequate to simulate most of the features of the time and rate dependent behaviour of natural clays up to failure. – Due to the overstress formulation, and the underlying assumption of a unique stress–strain–viscoplastic strain rate relationship, the model cannot reproduce accelerating creep; this has implications on the use of overstress models in the analysis of boundary value problems close to failure, as such a feature inhibits the occurrence of undrained failure. – The paper highlights the importance of controlling independently the soil drained and undrained shear strength, and its variation with the intermediate principle stress (i.e. shearing mode) if good predictions of failure conditions are to be obtained. – With the introduction of an additional creep parameter, in relation to a linear logarithmic creep law, the hyperbolic creep law is more flexible and provides predictions of the delayed com-

819

pression of soils, in closer agreement with laboratory and field data. – The model predicts that the critical state line is unique and time independent in stress space, but time dependent in terms of void ratio. Existing laboratory data corroborates the former but is inconclusive about the latter. For the case analysed, the ET model predicts that the CSL moves to higher values of void ratio with applied strain rate, about 0.025 per logarithmic cycle of strain rate, this being within the expected scatter in experimental measurements of void ratio. Further laboratory testing is needed to clarify the adequacy of these numerical predictions. Acknowledgments The authors would like to thank Fundação Ciência e Tecnologia, Portugal for sponsoring the doctoral programme of T.M. Bodas Freitas at Imperial College, London, during which this work was developed. References [1] Buisman K. Results of long duration settlement tests. In: Proceedings 1st ICSMFE, vol. 1. Cambridge; 1936. p. 103–07. [2] Bjerrum L. Engineering geology of Norwegian normally-consolidated marine clays as related to the settlements of buildings. Geotechnique 1967;17:81–118. [3] Garlanger JE. The consolidation of soils exhibiting creep under constant effective stress. Geotechnique 1972;22(1):71–8. [4] Tavenas F, Leroueil S. Effects of stress and time on yielding of clays. In: Proceedings of the 9th ICSMFE, vol. 1; 1977. p. 319–26. [5] Leroueil S, Tavenas F, Samson L, Morin P. Pre-consolidation pressure of Champlain clays. Part II. Laboratory determination. Can Geotechn J 1983;20:803–16. [6] Leroueil S, Kabbaj M, Tavenas F, Bouchard R. Stress–strain–time rate relation for the compressibility of sensitive natural clays. Geotechnique 1985;35(2):159–80. [7] Arulanandan K, Shen CK, Young RB. Undrained creep behaviour of a coastal organic silty clay. Geotechnique 1971;21(4):359–75. [8] Vaid YP, Campanella RG. Time-dependent behaviour of undisturbed clay. ASCE J Geotech Eng Div 1977;103(7):693–709. [9] Leroueil S, Marques MES. State of art: Importance of strain rate and temperature effects in geotechnical engineering. Measuring and modelling time dependent behaviour of soils. ASCE, Geotechnical Special Publication, No. 61; 1996, p. 1–60. [10] Liigaard M, Augustesen A, Lade PV. Characterization of models for the time dependent behaviour of soils. Int J Geomech 2004;4(3):157–77. [11] Perzyna P. Constitutive equations for rate sensitive plastic materials. Quart Appl Math 1963;20(4):321–32. [12] Adachi T, Okano M. A constitutive equation for normally consolidated clay. Soils Found 1974;14(4):55–73. [13] Adachi T, Oka F. Constitutive equations for normally consolidated clay based on elasto-viscoplasticity. Soils Found 1982;22(2):57–70. [14] Kutter BL, Sathialingam N. Elasto-viscoplastic modelling of the rate dependent behaviour of clays. Geotechnique 1992;42(3):427–41. [15] Vermeer PA, Neher HP. A soft soil model that accounts for creep. Beyond 2000 in Computational Geotechnics – 10 years of Plaxis international. Rotterdam: Balkema; 1999. [16] Yin J-H, Graham J. Elastic viscoplastic modelling of the time-dependent stress– strain behaviour of soils. Can Geotech J 1999;36:736–45. [17] den Haan EJ, van den Berg P. Evaluation of creep models for soft soils under axially symmetric conditions. Report CO-710203/21, GeoDelft, The Netherlands; 2001. [18] Yin J-H, Zhu J-G, Graham J. A new elastic visco-plastic model for timedependent behaviour of normally and overconsolidated clays: theory and verification. Can Geotech J 2002;39:157–73. [19] Leoni M, Karstunen M, Vermeer PA. Anisotropic creep model for soft soils. Geotechnique 2008;58(3):215–26. [20] Tavenas F, Leroueil S, La Rochelle P, Roy M. Creep behaviour of an undisturbed lightly overconsolidated clay. Can Geotech J 1978;15:402–23. [21] Boudali M. Comportement tridimensionnel et visqueux des argiles naturelles. PhD thesis. University of Laval, Québec; 1995. [22] Bodas Freitas TM, Potts DM, Zdravkovic L. The extension of elastic-viscoplastic models to general stress space. Geotechnique, Submitted to for publication. [23] Roscoe KH, Schofield AN. Mechanical behaviour of an idealised ‘‘wet-clay’’. In: Proceedings of 2nd European conf on soil mechanics and found eng, vol. 1; 1963. p. 47–54. [24] Roscoe KH, Burland JB. On the generalised stress–strain behaviour of ‘wet’ clay. Engineering plasticity. Cambridge University Press; 1968. p. 535–609. [25] Lagioia R, Puzrin AM, Potts DM. A new versatile expression for yield and plastic potential surfaces. Comput Geotech 1996;19(3):171–91.

820

T.M. Bodas Freitas et al. / Computers and Geotechnics 38 (2011) 809–820

[26] Yin J-H. Non-linear creep of soils in oedometer tests. Geotechnique 1999;49(5):699–707. [27] Aboshi H. Long-term effect of secondary consolidation on consolidation settlement of marine clays. Advances in geotechnical engineering: The Skempton conference. London: Thomas Telford; 2004. p. 345–55. [28] Tatsuoka F. Keynote Lecture: Inelastic deformation characteristics of geomaterial. Soil stress–strain behaviour: measurement, modelling and analysis. In: Proceedings of the geotechnical symposium in Rome; 2006. [29] Sheahan TC, Ladd CC, Germaine JT. Rate dependent undrained behaviour of saturated clay. J Geotech Eng 1996;122(2):99–108. [30] Zhu JG, Yin J-H, Luk ST. Time-dependent stress–strain behavior of soft Hong Kong marine deposits. Geotechnical Testing Journal, ASTM 1999;22(2): 118–26. [31] Augustesen A, Liingaard M, Lade PV. Evaluation of the time-dependent behaviour of soils. Int J Geomech 2004;4(3):137–56. [32] Sorensen KK. Influence of viscosity and ageing on the behaviour of clays. PhD thesis. University of London; 2006. [33] Yin J-H, Graham J. Viscous-elasto-plastic modelling of one-dimensional timedependent behaviour of clays. Can Geotech J 1989;26:199–209.

[34] Bodas Freitas TM. Numerical modelling of the time dependent behaviour of clays. PhD thesis. Imperial College, University of London; 2008. [35] Potts DM, Zdravkovic L. Finite element analysis in geotechnical engineering: theory. London: Thomas Telford; 1999. [36] Nash DFT, Sills GC, Davison LR. One-dimensional consolidation testing of soft clay from Bothkennar. Geotechnique 1992;42(2):241–56. [37] Smith PR. The behaviour of natural high compressibility clay with special reference to construction on soft ground. PhD thesis. Imperial College, University of London; 1992. [38] Hight DW, Bike R, Butcher AP, Clayton CRI, Smith PR. Disturbance of the Bothkennar clay prior to laboratory testing. Geotechnique 1992;42(2): 199–217. [39] Adachi T, Oka F, Mimura M. Mathematical structure of an overstress elastoviscoplastic model for clays. Soils and Foundations 1987;27(3):31–42. [40] Losacco N. Raising of embankments on soft clay. MSc thesis. Department of Civil and Environmental Engineering. Imperial College, London; 2007. [41] Marques MES, Leroueil S, Almeida MSS. Viscous behaviour of St-Roch-del’Achigan clay, Quebec. Can Geotech J 2004;41:25–38.