Loading path dependence and non-linear stiffness at small strain using rate-dependent multisurface hyperplasticity model

Computers and Geotechnics 49 (2013) 100–110

Contents lists available at SciVerse ScienceDirect

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

Loading path dependence and non-linear stiffness at small strain using rate-dependent multisurface hyperplasticity model Dedi Apriadi a, Suched Likitlersuang b,⇑, Thirapong Pipatpongsa c a

Faculty of Civil and Environmental Engineering, Institut Teknologi Bandung, Bandung, Indonesia Department of Civil Engineering, Faculty of Engineering, Chulalongkorn University, Payathai Road, Pathumwan, Bangkok 10330, Thailand c Global Scientific Information and Computing Centre, Tokyo Institute of Technology, Tokyo, Japan b

a r t i c l e

i n f o

Article history: Received 31 December 2009 Received in revised form 25 September 2012 Accepted 15 November 2012 Available online 17 December 2012 Keywords: Hyperplasticity Small strain Rate-dependent Recent stress history Non-linear stiffness

a b s t r a c t Many interpretations of small-strain experiments indicate a non-linear dependence of soil stiffness on pressure. This shows that small-strain stiffness can be expressed as a power function of the mean effective stress rather than as a linear function of stress. Many cases in the field show the importance of these behaviours to load–deformation prediction in a soil-structure interaction problem. This paper presents a numerical implementation and validation of non-linear pressure-dependent stiffness in a hyperplasticity model using a strain-driven forward-Euler integration scheme. The kinematic hardening function was incorporated into a finite number of multiple-yield surfaces of Modified Cam Clay to characterise small-strain stiffness and accommodate smooth transition of stiffness corresponding to different loading conditions and stress histories. Experimental results of current state and loading history dependence in overconsolidated clay are compared to the model prediction. Ó 2012 Elsevier Ltd. All rights reserved.

1. Introduction Interpretations of experimental data show that initial soil stiffness or small-strain stiffness is a non-linear function of mean effective stress [36,21,27,37]. Stiffness is also affected by other variables, such as the void ratio, anisotropic stress state and/or the preconsolidation pressure [8,14,38,30,33]. Detailed experimental investigations into the stress–strain response of overconsolidated soil have also shown that soil characteristics are dependent on both non-linearity and the most recent loading paths [19,5,34,12]. The effect of current loading history has also been observed using a different experimental approach [18,15,32]. It has been observed that zones exist at small strains, and they can change in both shape and size as the soil is subjected to different loading histories. The importance of these two behaviours in predicting load deformation in soil-structure interaction problems has been observed in many cases in the field [17,16,7,1]. To include certain physical principles such as conservation of mass, conservation of energy and the laws of thermodynamics in soil constitutive modelling, Likitlersuang [25] and Likitlersuang and Houlsby [24] proposed the hyperplasticity kinematic hardening Modified Cam Clay (KHMCC) model based on thermodynamics principles. To characterise small-strain stiffness and accommodate smooth ⇑ Corresponding author. Tel.: +66 2 2186343; fax: +66 22157304. E-mail address: [email protected] (S. Likitlersuang). 0266-352X/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compgeo.2012.11.007

transition of stiffness corresponding to different loading conditions and stress histories, the kinematic hardening function was incorporated in both energy and yield functions into a finite number of multiple-yield surfaces. However, numerical implementation and calibration of the model were restricted to a linear stress– strain relationship with stiffness proportional to isotropic pressure [23]. Hyperplastic formulations are capable of describing a non-linear response from small to large strain level. Small-strain non-linear stiffness can be expressed as a power function of the stress level to meet the non-linear response of materials observed in experiments. The current study extends the previous work on the linear version to the non-linear version, and the difficulty lies in working out the mathematics and numerical implementation. Apriadi [2] and Apriadi et al. [3] employed a numerical implementation of the non-linear KHMCC model using a strain-driven forward-Euler integration scheme under triaxial conditions by adopting small-strain stiffness in the form of the power function of pressure [11] in an energy function. Their work also presented an approach to determine the necessary parameters, obtained from bender element tests, for regulating the small-strain stiffness characteristic of the hyperplasticity KHMCC model. Numerical model implementation was verified against an analytical solution of ideal undrained triaxial test on normally consolidated clay [31,28]. The model’s capabilities in characterising the effect of recent stress history, the dependence of effective stress path on recent stress

D. Apriadi et al. / Computers and Geotechnics 49 (2013) 100–110

history during undrained shear, and the smooth and irreversible unloading–reloading responses under triaxial stress–strain conditions were highlighted. Some comparisons with experimental data on clay soils have also been carried out, such as the response of Bangkok clay to the K0-consolidated undrained direct simple shear (CK0UDSS) test. This work performs further numerical implementation of the non-linear KHMCC model in general stress conditions and validates it with an experiment using small strain levels on clay soils. It improves an existing model by incorporating a more realistic dependency of the small-strain stiffness on the mean effective stress. Owing to the limited data available on both consolidation pressure dependent and stress path dependent stiffness, the present work focuses on the latter. However, the effect of consolidation pressure on soil stiffness can be described with a single set of parameters [23]. The model demonstrates the non-linear stiffness response to generalised stresses. Therefore, it is not only consolidation pressure but also the stress path direction that affects soil stiffness. This article aims to present the capabilities of multiple-yield surface types of the constitutive model, which is based on hyperplasticity theory. A slight modification of mathematical expressions from the previous work is made based on the experimental evidence. Further numerical implementation under the hyperplasticity framework of the non-linear dependence of soil stiffness on pressure, which is shown by many interpretations of the experimental data, is developed in this work to present an experimentally realisable implementation of the hyperplasticity constitutive model. An experiment of undrained triaxial test on Speswhite kaolin [34] was selected to calibrate model parameters. This experiment was conducted to show the dependence of non-linearity, current state and loading history on relative directions of the current and previous loading paths in overconsolidated clay. Model validation was performed with the use of the response of tangent stiffness against stresses. Further, the experimental result of unloading– reloading at small strains was compared to the model prediction, and a study of the effect of the number of yield surfaces on stress–strain response was likewise carried out. The model formulation is initially described in terms of general stress conditions rather than triaxial stress conditions. Detailed numerical implementation into three-and two-dimensional stress–strain cases and model parameter determination are then presented, followed by a brief explanation of the test programme. Finally, remarks are made on the results of the model validation. This study is expected to contribute to the refinement of load– deformation prediction in soil-structure interaction problems.

are ap = aii and aq ¼

101

qffiffiffiffiffiffiffiffiffiffiffiffiffi 2 0 0 a a , where a0ij ¼ aij  13 aii dij is a deviatoric 3 ij ij

component of the internal variable tensor. 2.2. Hyperplasticity formulation Under the hyperplasticity framework (for details, see 10), the entire constitutive model is fundamentally governed by two scalar functions: energy function and yield function or dissipation function. In a typical formulation, Gibbs free energy function is suggested to associate stresses, material parameters and material memories (internal variables) with a unique potential function in which the reference pressure for zero volumetric strain is defined at 1 atm. It is noted that the hyperplasticity framework employs total energy instead of a rate form of energy, which is commonly used in classical plasticity theory. Since Gibbs free energy is equivalent to complementary energy, the negative elastic stored energy, negative dissipation energy and positive hardening energy are combined. A non-linear version of KHMCC is created to allocate non-linear elastic moduli at a small strain [11] to the Gibbs free energy, which ^ ij Þ: is expressed in the form E ¼ Eðrij ; aij ; a

E¼ þ

p2n p e  rij aij þ  nÞð2  nÞ kð1  nÞ ! b qa b pa ^ 0ij a ^ 0ij ^ 2kk H H þ dg; 2 3

p1n kð1 r Z

1 0

ð1Þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where pe ¼ p2 þ kð1nÞ q2 is defined as an equivalent stress variable 3g for convenience. The variables k,g and n are dimensionless material constants calibrated from elastic stress–strain relations at small strain levels. Atmospheric pressure, 1 atm (approximately 100 kPa), is defined for pr as a reference pressure, while aij is the total plastic strain tensor. Integration of differential hardening energy is evaluated in terms of the internal coordinate g, which is limited to values from 0 to 1. Here, g = 0 represents the initial hardening stage (the highest hardening response), while g = 1 represents the final hardening ^ ij is the kinematic stage (zero hardening response). The variable a internal variable tensor function of g, which can be integrated to obtain aij using Eq. (2). Therefore, aij can be regarded as a definite ^ ij over the domain of g. It is integral area of the functional variable a noted that all variables with ‘‘^’’ (hat) throughout this paper are referred to as internal variable functions of g.

aij ¼

Z

1

a^ ij dg:

ð2Þ

0

2. Non-linear kinematic hardening Modified Cam Clay model 2.1. Notation The standard soil mechanics sign convention of taking compressive stresses and strains as positive is used throughout this research, and all stresses are effective stresses. The following notation is adopted: rij is the effective Cauchy stress tensor, eij is the small-strain tensor and dij is Kronecker’s delta (dij = 1 if i = j, dij = 0 if i – j, where i, j 2 {1, 2, 3}). The stress invariants are qffiffiffiffiffiffiffiffiffiffiffiffiffiffi p ¼ 13 rii ; q ¼ 32 r0ij r0ij , where r0ij ¼ rij  pdij is a deviatoric component of the effective stress tensor. The corresponding strain invariqffiffiffiffiffiffiffiffiffiffiffiffi ants are ep = eii and eq ¼ 23 e0ij e0ij , where e0ij ¼ eij  13 eii dij is a deviatoric component of the strain tensor. Similarly, generalised qffiffiffiffiffiffiffiffiffiffiffiffiffi where vp ¼ 13 vii ; vq ¼ 32 v0ij v0ij , stress invariants are

v0ij ¼ vij  vp dij is a deviatoric component of the effective generalised stress tensor. The corresponding internal variable invariants

b p; H b q are non-linear kinematic hardening functions corresponding H to isotropic and deviatoric hardening responses, respectively, as expressed in the following equations: n

1n b p ¼ kp0 pr ð1  gÞ3 ; H 2ða  1Þ n 1n 3gp 0 pr bq ¼ H ð1  gÞ3 ; 2ða  1Þ

ð3Þ ð4Þ

where p0 is a mean pre-consolidation pressure at the end of the consolidation stage and a is a material constant. In the hyperplasticity framework, total strain components are considered conjugate variables of stresses, which are derived from the Gibbs free energy function. From Eq. (1), the strains eij can be obtained via differentiation with stresses rij as shown in the following equation:

eij ¼ 

  r0ij @E 1 p ¼  1 d þ þ aij : ij @ rij 3kð1  nÞ p1n pne 2gpr1n pne r

ð5Þ

102

D. Apriadi et al. / Computers and Geotechnics 49 (2013) 100–110

However, the Gibbs free energy function expressed in Eq. (1) and the strain components derived in Eq. (5) can produce only a constant modulus (n = 0) or power function of pressure-dependent modulus (0 < n < 1). For a linear pressure-dependent modulus (n = 1), these equations will be replaced by the following equations [10]:

    Z 1 b ^2 b ^2! H p ap H q aq p p q2 1  ln E¼  rij aij þ þ dg; ð6Þ k pr 6gp 2 2 0     r0ij @E 1 p q2 dij  ln eij ¼  ¼ ð7Þ þ þ aij : @ rij k pr 6gp2 3 2gp The second derivatives of Gibbs free energy are derived in Box 1. This form is applicable to 0 6 n 6 1. Box 1: Second derivatives of Gibbs free energy (after 10).

2





@ E @ rij @ rkl

8   9 dij dkl nq2 np 1 0 0 >  n > < = þ  r d þ d r ij kl 2 2 ij kl 9 k 6gpe 1 pr 3gpe   ¼ d d nkð1nÞ kl ij pr pe > : þ1 d d  ; r0 r0 >  2g

ik jl

3

4g 2 p2e

ij

kl

b Edg;

ð9Þ

0

b ^ 2p Hpa

b ^ 2q Hqa

þ 2 . 2 r ^ ij is referred to in Eq. (10) as the derivatives of b Therefore, v E ^ ij . with respect to a

@b E

2

b qa ^ 0ij : v^ ij ¼  ^ ¼ rij  Hb p a^ p dij  H 3 @ aij

ð10Þ

The evolution rule of the kinematic internal variable functions

a^ ij is derived from the derivatives of flow potential w with respect to generalised stress variables [10]:

@w

a^_ ij ¼ ^ : @ vij

where l is the artificial viscosity coefficient in units of stress  time [6]. It is defined as 1 kPa  second for rate-dependent calculations if the SI physical unit of stress in kPa and unit of time in seconds are employed. The operator h i is called Macaulay brackets, which de 0;  < 0 fine hi ¼ . ;  P 0 However, l can actually be defined as a true viscosity coefficient to incorporate rate-dependent behaviours of soil, such as creep problems. According to Eqs. (5) and (7), total strain components are derived as functions of rij and aij. Following the flow rule, the rate of change of total strain components can be written by the following equation:

# #

Z 1" @2E @2E @w dg:   fr_ kl g ¼ fe_ ij g  ^ kl ^ kl @ v @ rij @ rkl @ rij @ a 0

where M is a frictional critical state parameter [31], which is the vab pa ^ p g is a lue that the stress ratio q/p0 attains at critical state; ^c ¼ H hardening stress which represents the size of the yield surface func^ ij are generalised stresses. When the stress is inside the tion; and v ^ < 0Þ, materials behave elastically, but they behave yield surface ðy ^ P 0Þ. plastically when the yield surface is active ðy ^ ij are defined as the changing of the The generalised stresses v free energy functional with respect to the internal kinematic variables. The Gibbs free energy functional b E can be defined as the internal energy function of g. So the integration of b E throughout the domain of g, as shown in Eq. (9), is the Gibbs free energy.

p ^ ij þ where b E ¼  p1n kð1nÞð2nÞ þ kð1nÞ  rij a

ð12Þ

ð13Þ

"

ð8Þ

p2n e

;

By using Eqs. (11) and (12), the rate form of the stress–strain relationship obtained in Eq. (13) can be expressed.

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2 v^ q ^ 2p þ ^¼ v y  ^c; M

1

2l

# # Z 1" @2E @2E ^_ kl gdg: fa fe_ ij g ¼   fr_ kl g þ ^ kl @ rij @ rkl @ rij @a 0

@2E ¼ dij dkl ¼ identity matrix ^ @ aij @ rkl

Z

^i2 hy ¼ 2l

"

In the non-linear KHMCC model, the yield function in terms of generalised stress variables is defined as follows:

E¼

w¼

*rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi +2  2 v^ 2p þ v^Mq  ^c

ð11Þ

One may associate this kind of evolution rule with the flow rule used in classical plasticity. The flow potential w is defined based on elasto-viscoplastic theory (for details, see 40) in relation to the ^, as shown in Eq. (12). yield function y

ð14Þ

The derivative of flow potential w in Eq. (11) is shown in the following equation:

@w ¼ ^ ij @v

*rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi +  2 b pa ^p g v^ 2p þ v^Mq  H

l

v^ p dij þ M32 v^ 0ij rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2ffi : ^ 2p þ v^Mq 2 v 2 3

ð15Þ

3. Numerical implementation 3.1. Incremental stress–strain response The rate form of the stress–strain relationship based on continuous hyperplasticity is described in the previous section. However, continuous hyperplasticity is suitable for monotonic loading with simple hardening functions. In order to handle complicated loading conditions, multisurface hyperplasticity is employed in numerical implementations. Non-linearity is expressed by multiple piecewise responses (see also [26]). Therefore, multiple internal variables play a role as discrete memories of materials, and the smooth transition between piece-wise responses depends on the finite number of internal variables. For linear and non-linear KHMCC models, continuous yield surface is discretised to a finite number of yield surfaces [25,24,3]. The integration operator simply turns to the summation operator without losing the general meaning. Each yield surface has its own state variables, which are generalised and hardening stresses. A finite number of hardening stresses can be considered to be multiple material memories, which are updated when the multiple-yield surfaces are active. An illustration of multiple-yield surfaces in a three-dimensional stress space is shown in Fig. 1. The rate-dependent incremental response is obtained by integrating the incremental stress–strain relation using the strain-driven forward-Euler integration scheme. This indicates that an increment of the variable based on the rate form x_ ¼ f ðxÞ can typically be written in the manner xi+1  xi = f(xi)Dt, where

103

D. Apriadi et al. / Computers and Geotechnics 49 (2013) 100–110

2. Determine plastic parameters. for m 2 {1, 2, . . . , N} if y(m,t) 6 0 then a_ ijðm;tÞ ¼ 0 @w else a_ ijðm;tÞ ¼ @ v ðm;tÞ ijðm;tÞ 3. Compute incremental variables.

"

Drkl

@2E ¼ @ rij @ rkl

#1 "

N 1X e_ ij  N m¼1

"

# !# @2E  a_ klðm;tÞ Dt @ aklðm;tÞ @ rij

4. Update state variables.

rijðtþ1Þ ¼ rijðtÞ þ Drij for m 2 {1, 2, . . . , N} Fig. 1. Non-linear KHMCC yield surfaces in three-dimensional stress space [3].

aijðm;tþ1Þ ¼ aijðm;tÞ þ a_ ijðm;tÞ Dt Dt = ti+1  ti with i + 1 representing the current step number. Eq. (14) combined with (12) is used to update stress components in any increment of strain, and can be used to update the generalised stress components using Eq. (10) after the internal variable has been updated from the evolution rule given by Eq. (11). The algorithm of the rate-dependent numerical implementation is explained in Box 2. The subscripts m and t are positive integers representing the index of the yield surfaces and the increment of specific variables, respectively. A subscript index enclosed by parentheses () is used to distinguish between the tensor expression index. For an arbitrary variable tensor , the component of the tensor in accordance with the yield surface-mth at increment-ith is represented by ij(m,t), where N is the number of yield surfaces. All variables with ‘‘^’’ are removed and notified by a subscript index (m), as well as replacing a term g with a ratio m/N and dg with a ratio 1/N in a numerical process of discretisation. Box 2: Strain-driven forward-Euler integration scheme. For given Gibbs free energy E; e_ ij ; Dt; for m 2 {1, 2, . . . , N}: 1. Initialise.

eijðtþ1Þ ¼ eijðtÞ þ e_ ij Dt; aijðtÞ ¼

eijðtÞ ; rijðtÞ ; aijðm;tÞ

N 1X aijðn;tÞ N n¼1

for m 2 {1, 2, . . . , N}

kp0 pr1n  m3 ; 1 N 2ða  1Þ n

HpðmÞ ¼

HqðmÞ ¼

 3gpn0 p1n m3 r 1 N 2ða  1Þ

2 3

vijðm;tÞ ¼ rijðtÞ  HpðmÞ akkðm;tÞ dij  HqðmÞ a0ijðm;tÞ m cðm;tÞ ¼ HpðmÞ akkðm;tÞ ; yðm;tÞ N sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 vkkðm;tÞ 3 v0ijðm;tÞ v0ijðm;tÞ  cðm;tÞ ¼ þ 2 9 M2

wðm;tÞ

hyðm;tÞ i ¼ 2l

2

2 3

vijðm;tþ1Þ ¼ rijðtþ1Þ  HpðmÞ akkðm;tþ1Þ dij  HqðmÞ aijðm;tþ1Þ 5. Go to the next incremental step.

3.2. Two-dimensional stress–strain implementation For general stress and strain, there are six independent stress and strain components, i.e., the components of the stress tensor rij, where rij = rji, and the strain tensor eij, where eij = eji. Based on Eq. (14), the incremental stress–strain relationships of the rate-dependent strain-driven forward-Euler integration scheme under general stress can be generally formulated. The components of compliance tensor Cijkl (second derivatives of Gibbs free energy with respect to stresses) is calculated from Box 1. The evolution rule of kinematic internal variable function tensor aij(m,t) is determined from Eq. (2), using Eqs. (11) and (15). In two-dimensional cases, usually four stress tensor components are involved: three normal stresses r11, r22, r33 and one shear stress r12. Then, the increment in stress–strain relationships of the rate-dependent strain-driven forward-Euler integration scheme can simply be shown by the following equation:

2

C 11

C 12

C 13

6C 6 21 6 4 C 31

C 22

C 23

C 32

C 33

C 41

C 42

C 43

8 9 8 9 9 38 a_ 1ðm;tÞ > C 14 > Dr1 > > De1 > > > > > > > > > > > > > > N < _ < = < De = 1 X a2ðm;tÞ = C 24 7 2 7 Dr2 Dt: ¼  7 > N m¼1> a_ 3ðm;tÞ > C 34 5> Dr3 > De3 > > > > > > > > > > > > > : : ; : ; ; a_ 4ðm;tÞ C 44 Dr4 De4 ð16Þ

For computer programming purposes, r1, r2, r3, r4 are defined as r11, r22, r33, r12 and e1, e2, e3, e4 are defined as e11, e22, e33, c12, where c12 = e12 + e21 is an engineering shear strain. A similar definition is also applied to v1, v2, v3, v4, which are defined as v11, v22, v33, v12, and to a1, a2, a3, a4, which are defined as a11, a22, a33, a12 + a21. The Gibbs free energy and yield function can be rewritten as shown in Box 3. The component of compliance matrix Cij in Eq. (16) can be simply defined in expanded form as shown in Box 4. The effective generalised stress tensor vi(m,t) is calculated using Eqs. (17) and (18). The kinematic internal variable function tensor ai(m,t) is calculated using Eqs. (19) and (20).

viðm;tÞ ¼ 

@E @ aiðm;tÞ

i ¼ 1; 2; . . . ; 3;

¼ riðtÞ  Hpðm;tÞ

3 X

2 3

akðm;tÞ  Hqðm;tÞ a0iðm;tÞ ;

k¼1

ð17Þ

104

D. Apriadi et al. / Computers and Geotechnics 49 (2013) 100–110

v4ðm;tÞ ¼  a_ iðm;tÞ

@E 1 ¼ r4ðtÞ  Hqðm;tÞ a04ðm;tÞ ; @ a4ðm;tÞ 3

@wðm;tÞ ¼ @ viðm;tÞ

¼

ð18Þ

yðm;tÞ

hyðm;tÞ i

!

3 X

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 3  2 v X m ¼ vpðm;tÞ þ qðm;tÞ  HpðmÞ aiðm;tÞ ; N M i¼1

where

vkðm;tÞ þ M32 v0iðm;tÞ

2 9

k¼1

l 2

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ;  P3  02  P3 2 02 v

k¼1 kðm;tÞ

2v4ðm;tÞ þ

þ 32

9

vkðm;tÞ

k¼1

vpðm;tÞ ¼

M2

i ¼ 1; 2; . . . ; 3;

a_ 4ðm;tÞ ¼

Yield function:

ð19Þ

3 1X v ; 3 i¼1 iðm;tÞ

vqðm;tÞ ¼

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  X3  3  02 v02iðm;tÞ : 2v4ðm;tÞ þ i¼1 2

@wðm;tÞ @ v4ðm;tÞ hyðm;tÞ i

6

v04ðm;tÞ

ð20Þ

Box 4: Compliance matrix in expanded form for two-dimensional implementation.

Box 3: Gibbs free energy and yield function for two-dimensional implementation.

8   n n  o npr0 nq2 1 > 1 pr 1 1 > þ 3gp  3gp2i þ 3g  nkð1nÞ r02i ; i ¼ j 2 f1; 2; 3g 2 > k 9 4g 2 p2e > pr pe e e > >   8 9 > >   < > =  np r0i þr0j > n > nkð1nÞ 0 0 nq2 1 1 pr 1 1 >

>  n n o > 0 npr > 1 pr >  3gp24  nkð1nÞ r0i r0j ; i ¼ 4 xor j ¼ 4 > > pr pe 2g 2 p2e e > >   n o > n > 1 pr > 1 :  nkð1nÞ r024 ; i ¼ j ¼ 4 p p g g 2 p2

¼

l

2

M sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi :  P3  02  P3 2 2v02 þ v v 4ðm;tÞ kðm;tÞ k¼1 k¼1 kðm;tÞ 2 þ 32 9 M2

Gibbs free energy for n – 1:

p2n

p

eðtÞ ðtÞ E ¼  p1n kð1nÞð2nÞ þ kð1nÞ  r

riðtÞ aiðtÞ

i¼1

0 þ N1

4 X

N B X BHpðm;tÞ @

P3

a i¼1 iðm;tÞ

2

Hqðm;tÞ

  a0

4ðm;tÞ 2

s2

þ

2

2 P  3 þ

a02 iðm;tÞ

i¼1

3

m¼1

1

r

C C: A

pðtÞ k



ln

 

0 N B X BHpðm;tÞ þ N1 @ m¼1

pðtÞ pr

e

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi note: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  P3 02  P3 2 1 3  pe ¼ p2 þ kð1nÞq r ; q ¼ 2r02 ; p ¼ i¼1 i i¼1 ri 4 þ 3g 3 2  i = j 2 {1, 2, 3} ) C11, C22, C33  i – j 2 {1, 2, 3} ) C12, C13, C21, C23, C31, C32  i = 4 xor j = 4 ) C14, C41, C24, C42, C34, C43  i = j = 4 ) C44

Gibbs free energy for n = 1:

E¼

e

4  X q2 riðtÞ aiðtÞ  1  6gpðtÞðtÞ  i¼1

P3

a i¼1 iðm;tÞ

2

2

Hqðm;tÞ

 a0 2

4ðm;tÞ 2

2 P  3

þ

þ

i¼1

a02 iðm;tÞ

3

1

4. Model parameter determination

C C; A

The non-linear KHMCC model has a small number of parameters, namely three dimensionless material constant parameters (g, k and n), one parameter for critical state (M) and a kinematic hardening parameter (a). These parameters are obtained through the processes of parameter calibration from the experimental results. The dimensionless material constants g, k and n are related to elastic stiffness, which can be directly determined from experimental measurements of small-strain stiffness such as in bender element tests, or determined from empirical equations for preliminary analysis, as summarised in Table 1. The critical state frictional parameter M is determined at the stage of failure

where

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 3 kð1  nÞq2ðtÞ 1X peðtÞ ¼ p2ðtÞ þ riðtÞ ; ; pðtÞ ¼ 3 i¼1 3g rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  X3  3  02 qðtÞ ¼ 2r4ðtÞ þ r02iðtÞ : i¼1 2

Table 1 Small-strain stiffness empirical equation G = AF(e)pn (after [39]). Soil type

A

F(e)

N

Void ratio e

Test method

References

Reconstituted NC kaolin

3270

ð2:973eÞ2 1þe

0.5

0.5–1.5

Resonant column

Hardin and Black [9]

Several undisturbed NC clays Several undisturbed silts and clays

3270

ð2:973eÞ2 1þe

0.5

0.5–1.7

Resonant column

Hardin and Black [9]

1726

ð2:973eÞ2 1þe

0.46–0.61

0.4–1.1

Resonant column

Kim and Novak [20]

Undisturbed NC clays Undisturbed Italian clays

90

ð7:32eÞ2 1þe

0.6

1.7–3.8

Cyclic triaxial

Kokusho et al. [22]

4400–8100

ex ðx ¼ 1:11  1:43Þ

0.4–0.58

0.6–1.8

Resonant column and bender element

Jamiolkowsky et al. [13]

105

D. Apriadi et al. / Computers and Geotechnics 49 (2013) 100–110

Fig. 2. Determination of initial stiffness G from specified stress–strain curve. Fig. 5. Prediction of tangent stiffness Gt against deviatoric stress q of the 3-SKH model compared to experimental results (after 34).

Table 2 Summary of 3-SKH model parameters of Speswhite kaolin (after 34). M

k⁄

ecs

j⁄

G (kPa)

T

S

w

0.89

0.073

1.994

0.005

1964p0:65 R0:2 o

0.25

0.08

2.5

Fig. 3. Determination of parameter a using steps (1)–(3).

Fig. 6. Series of stress points with different stress paths for overconsolidated clay.

G ¼ gpn pr1n ; n

K ¼ kp p1n : r

ð22Þ ð23Þ

Apriadi et al. [3] showed that Eq. (22), which is derived from the constitutive model as shown in Fig. 2, is identical to the experimental observation from bender element tests [36,21,27,37]: Fig. 4. Determination of parameter a using a small parametric study.

G ¼ Cpn ;

conditions. It can be approximated by the relationship in the following equation:

M¼

6 sin / : 3  sin /

ð21Þ

It can be shown that the initial stiffness G and K, after completion of isotropic consolidation, are given by:

ð24Þ

where C and n are constants. This equation can be normalised into the form suggested by Houlsby et al. [11], shown in Equation (25), to obtain the dimensionless material constants g and n. The dimensionless material constant k can be determined using the elastic relationship as presented in Eq. (26). It is assumed that during isotropic consolidation, the ratio of K/G is constant regardless of consolidation pressure. The value of Poisson’s ratio m is made to range in a certain limit and can be related to the coefficient of earth pressure at rest K0, as shown in Eq. (27).

106

D. Apriadi et al. / Computers and Geotechnics 49 (2013) 100–110

(0 6 n 6 1) in modelling of the non-linear responses of materials; therefore, the accuracy on the initial stiffness is improved from the previous version by dropping the limited assumption on specifying n = 1. According to Houlsby and Puzrin [10], the kinematic hardening parameter a is calibrated to fit the stress–strain curve for specific test data (such as triaxial undrained or drained tests), as shown in Fig. 3. This parameter characterises the first loading curve (also referred to by many as backbone curve) of the stress–strain response with a hyperbolic curve. The following is the general step-by-step procedure for determining parameter a: (1) Normalise undrained or drained stress at the vertical axis by undrained or drained strength c. (2) Multiply the strain response at the horizontal axis by the appropriate initial stiffness and divide by c. (3) Find a as the inversion slope (1/secant stiffness) of a point at 50% of the normalised stress.

Fig. 7. Parametric study of kinematic hardening parameter a.

Table 3 Summary of non-linear KHMCC model parameters of Speswhite kaolin. k g n M a

1012 467 0.65 0.89 15

Dimensionless material constants

Slope of critical state line in the q–p plane Kinematic hardening parameter

However, steps (1)–(3) are applicable if the laboratory stress– strain curve can be approached with a simple hyperbolic stress– strain curve. In many practical scenarios, kinematic hardening parameter a can be directly obtained from a small parametric study to match the stress–strain curve of the specific test, as shown in Fig. 4. 5. Model validation

Fig. 8. Position of yield surfaces at initial conditions.

 n G p ¼g ; pr pr K k 2ð1 þ mÞ ¼ ¼ ; G g 3ð1  2mÞ K0 m¼ : 1 þ K0

ð25Þ ð26Þ ð27Þ

As emphasised in Table 1, the non-linear initial stiffness with a power pressure dependency of stiffness using material dimensionless parameter of n fits the experimental data better than that of a linear pressure dependency. Therefore, the value of n shown in Eq. (25) is an input parameter determined from bender element test, resonant column test or cyclic triaxial test, not a fitting parameter to particular loading conditions and stress histories as illustrated in this study. According to Apriadi [2], comparisons on the initial stiffness with a variation of n convinces that linear dependent stiffness (n = 1) is essentially inferior to power pressure dependent stiffness

Atkinson et al. [5], Jardine [18,15], Smith et al. [23], Stallebrass and Taylor [34] and Houlsby [12] observed an additional influence on the stress–strain behaviour concerning the loading path dependence or immediate stress history of soil, described by the most recent loading. This may take the form of an extended period of rest or a sudden change in the direction of the stress path. They found that zones exist at small strains, which can change in both shape and size as the soil is subjected to different loading histories. Stallebrass [35] and Stallebrass and Taylor [34] have conducted an experimental undrained triaxial test on Speswhite kaolin and simulated the results using the Three-Surface Kinematic Hardening (3-SKH) model [4]. They plotted the results with tangent shear stiffness Gt against deviatoric stress q, as shown in Fig. 5. Table 2 shows the 3-SKH soil model parameter in this simulation. The experiment was simulated with a given series of stress points (A, B, C, D and E), as shown in Fig. 6. Four cases of consolidated stress paths are given, i.e., ACB, ACAB, ACBDB and ACBEB, which represent a sudden change in the direction of the stress path before the undrained shearing stage at point B, which has an Ro, OCR in terms of isotropic pressure, of 2.4 (p0 = 300 kPa and pm = 720 kPa). Also shown in Fig. 6 are all stress histories located within a bounding yield surface. Hereafter, a symbol p labelled to figures and equations represents an effective mean stress. In this simulation, the small-strain stiffness parameters g, k and n are obtained using stiffness relationships normalised to Equation (25) with 1 atm (100 kPa) as a reference pressure. Eq. (28) shows the small-strain stiffness relationship used in this simulation.

 0:65 G p ¼ 467 : pr pr

ð28Þ

The constant k can be determined using Eq. (26), where Poisson’s ratio m is 0.3. The kinematic hardening parameter a is determined from a small parametric study of the actual stress–strain curve, as shown in Fig. 7. Table 3 summarises the non-linear KHMCC model in the simulation. Figs. 8 and 9 show the fields of yield surfaces at initial conditions and after the four different stress paths, respectively. The overlapping of yield surfaces observed in

D. Apriadi et al. / Computers and Geotechnics 49 (2013) 100–110

107

Fig. 9. Position of yield surfaces after four different stress histories.

Fig. 10. Prediction of tangent stiffness Gt against deviatoric stress q of the nonlinear KHMCC model compared to experimental results.

Fig. 11. Simulation of stress paths of the non-linear KHMCC model after different loading histories.

Fig. 9 is permissible. The possibility of overlapping of yield surfaces has already been discussed based on the laws of thermodynamics [29]. Fig. 10 shows the simulated results reflecting the plot of normalised tangent shear stiffness Gt against the deviatoric stress. This plot clearly shows that small-strain stiffness is affected by loading path history. Prediction of the model is also in good agreement with the experimental result. However, slightly non-smooth responses appear in the curve corresponding to the paths after ACBEB as a result of numerical discretisation. Complicated stress

paths activate more yield surfaces than do other paths; therefore, discretised responses are observed, but can be reduced by increasing the number of yield surfaces. Though measured data are not available, Fig. 11 aims to represent simulations of the behaviours of undrained effective stress paths due to stress path rotations after compression and after swelling. The figure shows the resulting four undrained stress paths from the non-linear KHMCC model. It is obvious that the non-linear KHMCC model can predict loading path dependence during undrained shear.

108

D. Apriadi et al. / Computers and Geotechnics 49 (2013) 100–110

that stiffness reduction for larger stress changes is in good agreement with the experimental results. The simulations used 100 yield surfaces, except for the unloading–reloading response, which was carried out using 50 yield surfaces. This simple computation was run using Intel Pentium M 1.6 GHz. 6. Effect of number of yield surfaces

Fig. 12. Prediction of unloading–reloading response at small strain levels of the non-linear KHMCC model compared to experimental results.

This additional numerical study aims to determine the effect of the number of yield surfaces on the response of different stress paths and tries to find the optimum number of yield surfaces to obtain accurate results. The numerical study is performed using 1, 2, 10 and 100 yield surfaces, and the effect is shown in stress–strain and stress path responses in Figs. 13 and 14. The figures clearly show that the stress–strain and stress path responses are highly influenced by the number of yield surfaces. It can be concluded that 20 yield surfaces are adequate to obtain accurate results. 7. Conclusions

Finally, the model was validated to simulate un/reloading behaviour at small strains compared to the experimental test data, as shown in Fig. 12. It is noted that the calibration of parameter a in Fig. 7 is carried out only for the first loading path, but the validation in Fig. 12 also shows un/reloading path forming a hysteresis loop, which is close to the experimental result. The model prediction shows that it can model a non-recoverable strain during an unloading–reloading response. It also shows

A numerical implementation of non-linear pressure-dependent stiffness in the hyperplasticity model under general stress conditions using a strain-driven forward-Euler integration scheme was presented. Further model verification with a small-strain experiment on overconsolidated clay soil was conducted. Verification clearly shows that the multisurface hyperplasticity model with a small number of parameters can accurately simulate the key features of small-strain behaviour using a unified framework. The

Fig. 13. Effect of number of yield surfaces on the stress–strain response.

D. Apriadi et al. / Computers and Geotechnics 49 (2013) 100–110

109

Fig. 14. Effect of number of yield surfaces on the stress path response.

unloading–reloading response in the model prediction shows a non-recoverable strain and loading path dependence, which are in good agreement with the experimental observation. In summary, this study achieved its two main purposes. First, it showed that using the multiple-yield surface could explain the effect of stress history (i.e., the different stress paths give the different non-linear stiffness responses). The more yield surfaces are used, the smoother the transition of stiffness that is provided. Second, it provides the actual non-linear pressure-dependent stiffness even inside the first yield surface. This is because when the stress state is located inside the first yield surface, the multiple-yield-surface model still offers a purely elastic response. These developments could benefit cases of small strain levels, especially during cyclic loading such as traffic load, small vibrations or a small earthquake.

Acknowledgements This work was supported by AUN/SEED-Net JICA. The authors would like to express their sincere appreciation for the grant.

References [1] Addenbrooke TI, Potts DM, Puzrin AM. The influence of pre-failure soil stiffness on the numerical analysis of tunnel construction. Geotechnique 1997;47(3):693–712. [2] Apriadi D. Development and numerical implementation of thermodynamicsbased soil model. Ph.D. thesis, Chulalongkorn University; 2009.

[3] Apriadi D, Likitlersuang S, Pipatpongsa T, Ohta H. On the numerical implementation of hyperplasticity non-linear kinematic hardening modified cam clay model. IES J Part A: Civil Struct Eng 2009;2(3):187–201. [4] Atkinson JH, Stallebrass SE. A model for recent stress history and non-linearity in the stress–strain behaviour of overconsolidated soil. In: Proc. 7th int. conf. on computer methods and advances in geomechanics, Cairns 1; 1991. p. 555– 60. [5] Atkinson JH, Richardson D, Stallebrass SE. Effect of recent stress history on stiffness of overconsolidated soil. Geotechnique 1990;40:531–40. [6] Griffiths DV. Finite element analyses of walls, footings and slopes. Ph.D. thesis, University of Manchester; 1980. [7] Gunn MJ. The prediction of surface settlement profiles due to tunnelling. In predictive soil mechanics. In: Houlsby GT, Schofield AN, editors. Proc. writh. memorial symp.. Oxford, London: Thomas Telford; 1993. p. 304–16. [8] Hardin BO. The nature of stress–strain behaviour for soils. Earthquake engineering and soil dynamics. In: Proceedings of the ASCE geotechnical engineering division specialty conference, Pasadena; 1978. p. 3–90. [9] Hardin BO, Black WL. Vibration modulus of normally consolidated clays. J Soil Mech Found, ASCE 1968;94(SM2):353–69. [10] Houlsby GT, Puzrin AM. Principles of hyperplasticity: an approach to plasticity theory based on thermodynamics principles. London: Springer-Verlag London Ltd.; 2006. ISBN-10: 184628239X, ISBN-13: 9781846282393. [11] Houlsby GT, Amorosi A, Rojas E. Elastic moduli of soils dependent on pressure: a hyperelastic formulation. Géotechnique 2005;55(5):383–92. [12] Houlsby GT. A model for the variable stiffness of undrained clay. In: Proc. int. symp. on pre-failure deformation of soils, Torino, 26–29 September, Balkema, vol. 1; 1999. p. 443–50 [ISBN: 9058090760]. [13] Houlsby GT, Wroth CP. The variation of the shear modulus of a clay with pressure and overconsolidation ratio. Soils Found 1991;31(3):138–43. [14] Jamiolkowski, M, Lancellotta, R, Lo Presti, DCF. Remarks on the stiffness at small strains of six Italian clays. Shibuya S, Mitachi T, Miura S, editors. vol. 2: A.A. Balkema; 1994. p. 817–836. [15] Jardine RJ. Some observations on the kinematic nature of soil stiffness. Soils Found Jpn Soc Soil Mech Fdn Eng 1992;32(2):111–24. [16] Jardine RJ, St. John HD, Hight DW, Potts DM. Some practical applications of a non-linear ground model. Proc. 10th Eur. conf. on soil mech. and found. eng. florence 1991;1:223–8.

110

D. Apriadi et al. / Computers and Geotechnics 49 (2013) 100–110

[17] Jardine RJ, Potts DM, Fourie AB, Burland JB. Studies of the influence of nonlinear stress–strain characteristic in soil–structure interaction. Géotechnique 1986;36(2):377–96. [18] Jardine RJ. Investigations of pile-soil behaviour with special reference to the foundations offshore structures. Ph.D. thesis, University of London; 1985. [19] Jardine RJ, Symes MJ, Burland JB. The measurement of soil stiffness in triaxial apparatus. Géotechnique 1984;34(3):323–40. [20] Kim TC, Novak M. Dynamic properties of some cohesive soils of Ontario. Canadian Geotech. J. 1981;18(3):371–89. [21] Kohata Y, Tatsuoka F, Wang, Lo, Jiang GL, Hoques E, et al. Modelling the non-linear deformation properties of stiff geomaterials. Géotechnique 1997; 47(3):563–80. [22] Kokusho T, Yoshida Y, Esashi Y. Dynamic properties of soft clay for wide strain range. Soils Found. 1982;22(4):1–18. [23] Likitlersuang S, Houlsby GT. Predictions of a continuous hyperplasticity model for Bangkok clay. Geomech Geoeng 2007;2(3):147–57. [24] Likitlersuang S, Houlsby GT. Development of hyperplasticity model for soil mechanics. Int J Numer Anal Method Geomech 2006;30:229–54. [25] Likitlersuang S. A hyperplasticity model for clay behaviour: an application to Bangkok clay. D.Phil. thesis, Oxford University; 2003. [26] Mroz Z, Norris VA. Elastoplastic and viscoplastic constitutive models for soils with application to cyclic loading. In: Pande GN, Zienkiewicz OC, editors. Soil mech. cyclic transient loads. NY: John Wiley & Sons; 1982. p. 173–218 [chapter 8]. [27] Pennington DS, Nash DFT, Lings ML. Anisotropy of G0 shear stiffness in Gault clay. Géotechnique 1997;47(3):391–8. [28] Potts DM, Ganendra D. An evaluation of substepping and implicit stress point algorithm. Comput Methods Appl Mech Eng 1994;119:341–54. [29] Puzrin AM, Houlsby GT. On the non-intersection dilemma in multi-surface plasticity. Geotechnique 2001;51(4):369–72. [30] Rampello S, Silvestri F, Viggiani G. The dependence of G0 on stress state and history in cohesive soils. In: Proc. 1st int. conf. on pre-failure deformation characteristics of geomaterials, Sapporo 2; 1994. p. 1155–60.

[31] Roscoe KH, Burland JB. On the generalised stress–strain behaviour of wet clay. Engineering plasticity. Cambridge University Press; 1968. p. 535–610. [32] Smith PR, Jardine RJ, Hight DW. The yielding of Bothkennar clay. Géotechnique 1992;42(2):257–74. [33] Soga K, Nakagawa K, Mitchell JK. Measurement of stiffness degradation characteristics of clays using a torsional shear device. In: First international conference on earthquake geotechnical engineering, Tokyo, November 14–16; 1995. p. 107–12. [34] Stallebrass SE, Taylor RN. The development and evaluation of a constitutive model for the prediction of ground movements in overconsolidated clay. Géotechnique 1997;47(2):235–54. [35] Stallebrass SE. The effect of recent stress history on the deformation of overconsolidated soils. Ph.D. thesis, City University; 1990. [36] Tanizawa F, Techavorasinskun S, Yamaguchi J, Sueoka T, Goto S. Measurement of shear wave velocity of sand before liquefaction and during cyclic mobility. In: Shibuya S, Mitachi T, Miura S, editors. Proceedings of the international symposium on pre-failure deformation characteristics of geomaterials, vol. 1. Rotterdam: Balkema; 1994. p. 63–8. [37] Techavorasinskun S, Amornwithayalax T. Elastic shear modulus of Bangkok clay during undrained triaxial compression. Géotechnique 2002;52(7): 537–40. [38] Viggiani G. Small-strain stiffness of fine grained soils. Ph.D. thesis, City University, London; 1992. [39] Yimsiri S, Soga K. Anisotropy of highly-overconsolidated clay in small-and intermediate-strain levels. In: Proc. of the 14th Southeast Asian geotechnical conference, 14th SEAGC, Hong Kong; 2001. [40] Zienkiewicz OC, Cormeau IC. Viscoplasticity, plasticity and creep in elastic solids – a unified numerical solution approach. Int J Numer Methods Eng 1974;8:821–45.