Finite element modelling of rate-dependent ratcheting in granular materials

Computers and Geotechnics 38 (2011) 105–112

Contents lists available at ScienceDirect

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

Finite element modelling of rate-dependent ratcheting in granular materials A. Karrech a,b,⇑, A. Seibi c, D. Duhamel d a

CSIRO Earth Science and Resource Engineering, 26 Dick Perry Ave., Kensington, WA 6151, Australia University of Western Australia, School of Civil and Resource Engineering, 35 Stirling Hwy, Crawley Perth, WA 6009, Australia c Petroleum Institute, Mechanical Engineering Program, P.O. Box 2533, Abu Dhabi, United Arab Emirates d Universit´e Paris-Est, UR Navier, Ecole Nationale des Ponts et Chauss´ees, 6 et, 8 Avenue Blaise Pascal, Cit´e Descartes, Champs sur Marne, 77455 Marne-La-Vall´ee, Cedex 2, France b

a r t i c l e

i n f o

Article history: Received 3 January 2010 Received in revised form 2 August 2010 Accepted 5 August 2010 Available online 20 November 2010 Keywords: Viscoplasticity Ratcheting Granular materials Mixed hardening Numerical implementation

a b s t r a c t The present paper introduces a comprehensive model that is capable of describing the behaviour, under cyclic loading, of the granular materials used in railway tracks and road pavement. Its main thrust is the introduction of the ‘‘Chicago’’ law in a continuum approach to account for the ratcheting effects. It also emphasizes rate-dependency as a dissipative mechanism that acts independently or jointly with the ratcheting effect as well as the non-associated plasticity. The numerical procedure is based on the return mapping algorithm, where Newton’s method is used to calculate the nonlinear consistency parameter of the flow rule and to obtain a consistent tangent modulus. The model was applied to specific numerical examples including multi-axial and cyclic loading conditions. Ó 2010 Elsevier Ltd. All rights reserved.

1. Introduction Granular materials have received increasing interest from modern industry due to their complex mechanical behaviour. In particular, the ballasts used in railway tracks or similar agglomerates used in road pavement represent key materials for the transport industry. Understanding the behaviours of these materials, especially under cyclic loading, represents a challenging subject requiring considerable research effort. Different modelling approaches, based on discrete or continuum descriptions, have been used to study the responses of these materials, especially under cyclic loading. Discrete element methods such as molecular dynamics (Cundall and Strack [8,9], Gallas et al. [13], Roux [35], Karrech et al. [20]) and contact dynamics (Azema et al. [2], Saussine et al. [36]) are well adapted to study samples consisting of a finite number of particles. Although the discrete element methods accurately describe the local mechanisms of energy dissipation caused by ratcheting effects, a loss of contacts, and granular flow, these techniques cannot treat samples made of a large number of grains. Continuum approaches represent convenient time-saving alternatives, though they are limited in predicting local inter-granular mechanisms such as contact loss (Karrech et al. [19]). Research studies on cyclically loaded granular materials of fine granulometry, such as sands and silts, are rather well established ⇑ Corresponding author at: CSIRO Earth Science and Resource Engineering, 26 Dick Perry Ave., Kensington, WA 6151, Australia. E-mail address: [email protected] (A. Karrech). 0266-352X/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved. doi:10.1016/j.compgeo.2010.08.012

(Eekelen and Potts [12], Kovacevic and Vaughan [22], Lehane et al. [24]). Within the VErification of Liquefaction Analysis by Centrifuge Studies (VELACS; Arulanandan and Scott [1]) research project, several experimental and numerical studies have been conducted on these materials to assess their responses under repeated loading. However, these materials are different in terms of the micro-structure and water content than ballasts, which is of interest in our study. Studies on the cyclic loading of ballasts revealed that ratcheting, rate-dependency, and non-associated flow are necessary features that should be taken into account in the formulation of constitutive models (see Karrech [18] and the references therein). The experimental studies of the same spirit showed that ratcheting can be described with a logarithmic function known as the Chicago law (Ben-Naim et al. [3], Nowak et al. [31]). This empirical law was confirmed numerically using discrete element methods by Rosato and Yacoub [34], Philippe and Bideau [32], Rémond [33]. Both the experimental and the numerical results showed that granular materials exhibit a gradual decrease of the volume of voids, resulting in a nonlinear permanent deformation with respect to the number of cycles. A similar decay was observed by both Guérin et al. [14], Bodin et al. [4] when studying the settlement of granular materials of railway platforms due to repeated loading. In addition, Saussine et al. [36] developed a discrete element model based on contact dynamics that showed the same trend. More recently, Karrech et al. [20] developed a numerical algorithm based on molecular dynamics and showed that the logarithmic law is appropriate for long-term settlement prediction in granular materials. The obtained empirical law proved to be an

106

A. Karrech et al. / Computers and Geotechnics 38 (2011) 105–112

accurate description of the decaying phenomenon of cyclically loaded granular materials. The present work is based on several existing computational procedures that involve the above mechanisms of energy dissipation within the framework of continuum mechanics. The purpose of this work is to combine some of these techniques and come up with a comprehensive numerical tool that encompasses the observed logarithmic decay as well as the rate-dependency and frictional aspects for cyclically loaded ballast samples or similar materials. In this context, Lu and Wright [28] suggested a numerical model to describe the behaviour of asphalt concrete under cyclic loading, where an exponential empirical law was introduced to describe the permanent deformation. However, the cyclic effect was included in the rate-dependency mechanism rather than ratcheting. Lorefice et al. [26] introduced a different model that was intended to describe the rate-dependent behaviour of quasi-brittle materials like concrete without taking into account the ratcheting effect. Recently, Cvitanic et al. [10] formulated a theoretical and numerical framework based on non-associated plasticity with isotropic hardening. The numerical integration included the return mapping technique, which consists of projecting the trial stress on the surface of yielding at every increment (see also Key [21], Dafalias [11], Chen and Han [7], Cernocky and Krempl [6], Marques and Owen [29], Loret and Prevost [27]). This technique, which was first introduced by Wilkins [39], was improved by Simo and Taylor [37], Hofstetter et al. [15], who proposed the consistent tangent modulus. Liu [25], Kang [16], Kumar and Nukala [23] proved that return mapping coupled with a consistent tangent modulus is an effective and robust tool to integrate nonlinear constitutive equations. Although the cited models covered particular combinations of rate-dependency, ratcheting, non-associated flows, and local compaction, none of them contained a comprehensive approach to address these phenomena altogether to predict the response of granular materials. The present study is developed in light of the above mentioned contributions while taking into account the (i) non-associated flow rule, (ii) Drucker–Prager yield and plastic potential functions (instead of the frequently used Von-Mises yield function), (iii) logarithmic ratcheting effect, and (iv) robust integration method based on the return mapping algorithm. This paper is divided into three major sections. Section 2 gives a brief description of the problem under consideration, defines the constitutive continuous model, and highlights the assumptions adopted in the study. Section 3 describes the numerical procedure used to integrate the constitutive model. Section 4 is dedicated to the application of the developed model on particular examples. 2. Constitutive continuous model The constitutive model developed in this paper is based on a linear combination of kinematic and isotropic hardening, where the strain increment is defined by:

dij ¼ deij þ diij þ dkij ¼ deij þ dijv p

ð1Þ

stress, respectively. The limit of elasticity is described by a Drucker– Prager yield function written as:

^ ij Þ ¼ f ðrij  aij Þ ¼ a1^I1 þ f ðr

qffiffiffiffiffiffiffi 3^J 2 ¼ q

ð3Þ

where f is a convex function. The first and second invariants are expressed by ^J2 ¼ 12 ^Sij ^Sij and ^I1 ¼ rii  aii , respectively. The terms sij ¼ rij  pdij ; ^sij ¼ sij  aij þ 13 aii dij ; p, and dij represent the deviatoric stress, kinematic stress, hydrostatic pressure, and Kronecker identity, respectively. In this model, we express the translation or backstress aij using Prager’s hardening rule as follows:

daij ¼ cdjij

ð4Þ

In order to describe the ratcheting effect, we express the Prager parameter c in the above equation in terms of the number of loading cycles N. For a given periodic loading of frequency x, a cycle corresponds to a period of time 2px. The Chicago logarithmic law (obtained by Ben-Naim et al. [3], Nowak et al. [31], Rosato and Yacoub [34], Philippe and Bideau [32], Rémond [33]) is used to express the kinematic hardening as follows:

cðNÞ ¼

c1 A   1 þ c2 ln Ns þ 1

ð5Þ

where c1, c2, and s are material constants that can be estimated experimentally. This logarithmic effect is mainly due to the compaction of granular materials subjected to repeated loading, as explained by Guérin et al. [14], Bodin et al. [4], Saussine et al. [36], Karrech et al. [20]. The term q in Eq. (3) can be expressed in terms of the yield strength using the particular case of uniaxial loading as follows:

^ e Þ ¼ a1 r ^e ¼ q ^e þ r f ðr

ð6Þ

Therefore, it can be deduced that:

r^ e ¼

a1^I1 þ

qffiffiffiffiffiffiffi 3^J 2

1 þ a1

ð7Þ

Furthermore, a plastic potential in accordance with the Drucker–Prager equation is considered:

^ ij Þ ¼ a2^I1 þ gðr

qffiffiffiffiffiffiffi 3^J 2

ð8Þ

If the constants al and a2 are different, the material is non-associated. The flow rule in this case is defined by:

_ ijv p ¼ k_

@g @ rij

ð9Þ

_ a positive parameter that can be evaluated by selecting where k_is the appropriate viscous potential. In the present study, Perzyna’s model is adopted to describe the rate-dependency. Therefore, the consistency parameter can be expressed as follows:

wðf Þ k_ ¼ f

ð10Þ

The superscripts e and vp denote respectively the elastic and viscoplastic strains. The viscoplastic strain is decomposed into an isotropic strain, diij ¼ Mdvij p , and a kinematic strain, dkij ¼ ð1  MÞdvij p where the constant M, describes the combination of both hardening effects. It takes values of 1 when the hardening is fully isotropic and of 0 when it is fully kinematic. The elastic behaviour is described with a linear law given by:

where f is a fluidity parameter, and w is a monotonically increasing function defined by w(f) = f if f P 0 and w(f) = 0 otherwise. By defin_ and using the principle of maximum dissipation, ing f ¼ f  w1 ðkfÞ it can be seen that the Kuhn-Tucker conditions hold during the viscoplastic flow:

vp drij ¼ C ijkl ðdkl  dkl Þ

_ as will These conditions can be used to determine the factor k, be shown in the next sections. Note that unlike the model of Suiker and Borst [38], where the plastic strain evolution is expressed in terms of the number of cycles, the rate-dependency in our

ð2Þ

where Cijkl is a fourth-order tensor describing the elastic behaviour, ij and rij are second-order tensors denoting the strain and Cauchy

k_ f ¼ 0;

k_ P 0; and f 6 0

ð11Þ

107

A. Karrech et al. / Computers and Geotechnics 38 (2011) 105–112

approach requires an explicit dependency of the inelastic strain on time. This means that ratcheting and rate-dependency can act as independent mechanisms. 3. Finite element implementation The algorithm used to integrate the developed nonlinear material behaviour is similar to those adopted for rate-independent cases. A computational procedure based on the implicit return mapping scheme developed by Simo and Taylor [37] for associated plasticity and extended by Cvitanic et al. [10] for non-associated plasticity is adopted in the present work to take into account the rate-dependency of granular materials. A detailed description of the algorithm developed in this study is presented.

The return mapping starts from a converged configuration, biC n ðrn ; vn p Þ, where the state variables are known to deduce the   p new configuration, biC nþl rnþl ; vnþ1 , which satisfies the material behaviour. There are various techniques that can be used to integrate the behaviour of materials (Simo and Taylor [37], Hofstetter et al. [15]). Most of them are based on predictor–corrector schemes, which generally involve two steps. The first step uses the converged configuration Cn, Hooke’s law (2), and a deformation increment D, which can be evaluated using the incremental deformation theory (Cvitanic et al. [10]) to calculate a trial elastic stress tensor as follows:

ð13Þ

The term on the right-hand side of Eq. (13) is the correction or trial stress. The viscoplastic strain increment can be obtained from Eq. (9) as follows:

Dvij p ¼ Dk

@g @ rij

ð14Þ

3.2. Consistency factor

CðDk þ dðDkÞÞ ¼ 0

ð18Þ

Expanding the above equation in a Taylor Series to the first-order about the approximate consistency factor results in the following expression:

CðDkÞ ¼ dðDkÞ

@ CðDkÞ @ðDkÞ

ð19Þ

dðDkÞ ¼

@f @f @f @g _ f ðrnþ1 ij Þ  q  kf  Dk @ rkl C ijkl @ rkl  Dkcð1  MÞ @ rij @ rij @ CðDkÞ @ðDkÞ

ð20Þ

Eq. (20) represents the correction of the consistency factor, which can be solved using Newton’s method. The iterations consist of calculating d(Dk) and updating the consistency factor in an additive sense, such that Dknew = Dkold + d(Dk) until convergence is achieved. Once the consistency factor is known, the stresses and strains can be deduced. 3.3. Consistent tangent modulus Differentiating the incremental Eq. (13) and using Eq. (14) results in:

dDðrij Þnþ1 ¼ C ijkl dDkl  dDk

@g @2g  Dk dDðrij Þnþ1 @ rij @ rij @ rkl

! ð21Þ

Therefore, it can be deduced that:

  @g dDðrij Þnþ1 ¼ Dijkl dDkl  dDk @ rij

ð22Þ

 1 @2 g where Dijkl ¼ C 1 . The consistency condition f_ ¼ 0 ijkl þ Dk @ rij @ rkl can be applied to deduce the consistent tangent modulus:

  Bij B dDðrij Þnþ1 ¼ Dijkl  kl dDkl H

ð23Þ





Bkl ¼ Dijkl @@rf ; Bkl ¼ Dijkl @@g and H ¼ @@rf Dijkl @@g rij , rij þ ij ij p ffiffiffiffiffi 2 ^ ^ _ Mð1þa1 Þ ða2 I1 þ 3J 2 ÞHp þ Dfkt. In order to calculate ð3a1 a2 þ 3=2Þð1  MÞc þ q Dijkl, it is necessary to note that:

Dk@ 2 g 3 Dk dik djl þ dil djk dij dkl 3 ^Sij ^Skl ¼   @ rij @ rkl 2 X 2 X2 2 3 where X ¼

ð24Þ

pffiffiffiffiffiffiffi 3J2 . Hence, it can be deduced that:

C 1 ijkl þ

Dijkl ¼

!

3 Dk dik djl þ dil djk dij dkl 3 ^Sij ^Skl   2 X 2 X2 2 3

!!1 ð25Þ

It can be noted that the above tensor can be written as follows:



E1 ijkl 

Dk v ij v kl X

1

ð15Þ

Dijkl ¼

ð16Þ

After some simple mathematical identification, it can be seen that:

using Eqs. (3), (7), and (13), it can be seen that vp vp ^ Tij  C ijkl Dkl CðDkÞ ¼ f ðr  cð1  MÞDkl Þ¼0

Assume that after a given iteration of the Newton’s method, an approximation Dk to the consistency factor is obtained. Let d(Dk) be the difference between this value and the exact solution to the consistency condition. This implies that:

where

One of the main features of the computational procedure is the computation of the consistency factor, which allows the projection of the trial state onto the surface f ¼ 0, for correction. Because the evaluation of the stress–strain relationship takes a considerable amount of computing time, the efficiency of the algorithm used to integrate the material behaviour is crucial. In particular, the algorithm suggested herein uses Newton’s method to calculate the consistency parameter and deduce the equivalent viscoplastic strain, which is evaluated only once. Jointly with the consistent tangent modulus, which will be developed in the next subsection, this technique guarantees a high rate of convergence. The procedure is based on the enforcement of the consistency condition, which can be symbolically written in a scalar form as follows:

^ ij Þnþ1 Þ ¼ 0 CðDkÞ ¼ f ððr

ð17Þ

ð12Þ

If the elastic trial state does not violate the limit of elasticity, then it represents the solution of the problem. However, if the trial state is not admissible, a correction is required. The trial stress represents a starting point of the viscoplastic correction step, which leads to the following result:

 vp vp ðrij Þnþ1 ¼ ðrij Þn þ C ijkl Dkl  Dkl ¼ rTij  C ijkl Dkl

@ f @g @f @g C ijkl  Dkcð1  MÞ ¼0 @ rij @ rkl @ rij @ rij

Therefore, combining Eqs. (17) and (19) leads to:

3.1. Return mapping

rTij ¼ ðrij Þn þ C ijkl Dekl

^ Tij Þ  Dk CðDkÞ ¼ f ðr

Eq. (16) is approximated to the first-order using the Taylor Series as follows:

  d d þd d 2 3GDk Eijkl ¼ c K  G dij dkl  cKdij dkl þ 2cG ik jl il jk 3 X 2

ð26Þ

ð27Þ

108

A. Karrech et al. / Computers and Geotechnics 38 (2011) 105–112

X where c ¼ Xþ3G . Since DXk  1, therefore, a first-order development Dk results in c  1  3GXDk which means that 3GXDk c  3GXDk. Therefore,

  dik djl þ dil djk 2 Eijkl ¼ K  cG dij dkl þ 2cG 3 2

ð28Þ

when DXc  1 approaches zero, c approaches 1. Therefore, it can be seen that the inverse of Eijkl can be reduced to Cijkl in the above mentioned particular case. Applying the Sherman–Morisson formula, Eq. (26) can be rewritten as:

Dijkl ¼ Eijkl 

v ij Eijkl Eijkl v kl v ij Eijkl v kl  DXk

ð29Þ

3^s

where v ij ¼ 32 Xij . Substituting vij and Eijkl in Eq. (29) leads to the following expression:

Dijkl ¼

  ^Sij S^kl dik djl þ dil djk 2 K  cG dij dkl þ c2G þb 3 2 X X 2

cGÞ where b ¼  3ð3cG X . Since Dk

Dk X

 1, then b   ð3cGÞX

2

Dk

ð30Þ Using Eq. (30),

the terms of Eq. (23) representing the consistent operator can be deduced:

^sij @f ¼ 3Ka1 dij þ ð3cG þ bÞ @ rkl X ^skl @g Bkl ¼ Dijkl ¼ 3Ka2 dkl þ ð3cG þ bÞ X @ rij @f @g Dijkl ¼ 9bKa1 a2 þ ð3cG þ bÞ @ rij @ rkl

Bij ¼ Dijkl

ð31Þ ð32Þ ð33Þ

Therefore, 

H ¼ 9bKa1 a2 þ ð3cG þ bÞ þ ð3a1 a2 þ 3=2Þð1  MÞc  qffiffiffiffiffiffiffi Mð1 þ a1 Þ2 a2^I1 þ 3^J 2 Hp fk_ þ þ q Dt

material subroutine UMAT. In this section, we first examine the coherence of the model through simple examples and then apply it to a particular case that is relevant to granular materials subjected to repeated loading. The material parameters used for simulation are summarized as follows: Young’s Modulus: E = 30 MPa; Poisson’s ratio: m = 0.3; fluidity coefficient: f = 5 MPa s, friction angle: u = 35°, dilation angle: w = 10°, mixing hardening constant: M = 0.5, and Chicago law parameters: c1 = 20 MPa, c2 = 2, and s = 20. Isotropic hardening was described by the following stress–strain input (evp, re (MPa)) = {(0, 0.2), (0.001, 0.3), (0.005, 0.34), (0.01, 0.4), (0.05, 0.6)}. These parameters can be identified experimentally using triaxial tests at different rates of loading. The constants a1 and a2 can be related to the friction and dilation angles u and w, respectively:

a1 ¼ 

2 sinðuÞ 3ð3  sinðuÞÞ

and a2 ¼ 

2 sinðwÞ 3ð3  sinðwÞÞ

It can be noted that the model is non-associated, as a1 is different than a2. 4.1. Cylindrical sample under monotonic triaxial loading This example consists of a cylindrical sample with unit radius and length (see Fig. 1). Radial and axial boundary conditions are applied on the upper and circumferential faces, while the bottom side is constrained from moving in the axial direction. It is assumed that f = c = 0 in this particular case, in order to examine the capability of the model in predicting the frictional Drucker–Prager behaviour. We first calculate the analytical solution and then use the necessary simplifications to compare it with the numerical results. The displacement field that governs the motion takes the form:

x ¼ /ðX; tÞ ¼ X  aðtÞZez  bðtÞRer ð34Þ

It is worth noting that this work presents a general formulation of the problem that can be reduced to known particular cases: (i) If al = a2 = 0, the yield and potential functions reduce from a Drucker–Prager’s to a Von-Mises’s form, which is suitable for metallic materials. Most of the existing formulations where rate-dependency and mixed hardening are implemented can be classified within this category. For instance, Mayama et al. [30] developed a constitutive model for cyclically loaded viscoplastic stainless steel at room temperature. Kanga et al. [17] suggested a similar model where isotropic and kinematic hardening effects were considered and validated the model through an experimental investigation on stainless steel at ambient temperature. Kumar and Nukala [23] introduced a constitutive model based on multi-component forms of kinematic and isotropic hardening variables in conjunction with the Von-Mises yield criterion. (ii) If the constant M equals 1, then our hardening description reduces to the isotropic case, where material ratcheting does not play a tangible role. The formulation then reduces to the elasto-viscoplastic model of Ceia [5], where inviscid description follows Drucker–Prager’s equation. This description was used to predict the response of reinforced concrete structures subjected to dynamic loads (see also Cvitanic et al. [10] for the integration of an elasto-plastic behaviour without ratcheting). (iii) When f reduces to zero, then the model produces the ratcheting effects in rate-independent materials (see, for instance, Chen and Han [7]).

ð36Þ

Therefore, the total rate of stretches can be written as follows:

d¼

a_ b ðer  er þ eh  eh Þ  ez  ez 1a 1b

ð37Þ

The rate of plastic stretches can be expressed as: p

d ¼ k_

  @g 3 sij ¼ k_ a2 dij þ @r 2 X

4. Numerical examples The numerical procedure at hand was implemented into the general–purpose finite element code ABAQUS using the user

ð35Þ

Fig. 1. Deformable cylinder under triaxial loading.

ð38Þ

109

A. Karrech et al. / Computers and Geotechnics 38 (2011) 105–112



2 3



r_ ¼ K  G trðd  dp Þ1 þ 2Gðd  dp Þ

ð39Þ

If the response is purely elastic, then dp = 0. After the integration of Eq. (39), it can be seen that the elastic solution within the hypothesis of small perturbations reads:



 2 3   2 ¼  K  G ð2a þ bÞ  2Ga 3

rzz ¼  K  G ð2a þ bÞ  2Gb rrr

ð40Þ

0.2

Hydrostatic pressure, −p, MPa

Hence, knowing that the mapping of Eq. (36) creates no rotation, it can be deduced that the evolution of Cauchy stresses can be written as follows:

rhh ¼ rrr ; rrh ¼ rrz ¼ rhz ¼ 0



 2 Szz 3 X   2 Srr ¼  K  G ð2a þ bÞ  2Ga þ 3a2 kK þ 3kG 3 X

rzz ¼  K  G ð2a þ bÞ  2Gb þ 3a2 kK þ 3kG rrr

ð41Þ

rhh ¼ rrr ; rrh ¼ rrz ¼ rhz ¼ 0 Hence, the solution of the problem can be presented in terms of the pressure and deviatoric stress as follows:

qdev ¼ 2Gðb  aÞ þ 3kG p ¼ Kð2a þ bÞ  3ka2 K

ð42Þ

0

0.5

1

1.5

0.35 0.3 0.25 0.2 0.15 0.1 0.05

Numerical Analytical

0 0

1

2

3

4

Radial strain, − εrr, % Fig. 4. Variation of the quadratic stress with respect to the radial strain in a case of quasistatic triaxial loading.

0.2

Hydrostatic pressure, −p, MPa

Deviatoric stress, qdev, MPa

Numerical Analytical

Fig. 3. Variation of the hydrostatic pressure with respect to the axial strain.

0.35 0.3 0.25 0.2 0.15 0.1 0.05 0 −0.5

0.05

Axial strain, εzz, %

qðp Þ ¼ 3a1 ð2a þ bÞK  9a1 a2 k þ 2ðb  aÞG þ 3kG Figs. 2–5 show the analytical and numerical results obtained for this particular example. The loading scenario starts with subjecting the sample to a confining pressure so that it elongates in the axial direction and shrinks in the radial direction (because the axial surface area is smaller than the circumferential surface area). This explains why the pressure starts from zero and the sample axial and

0.1

0 −0.5

Deviatoric stress, qdev, MPa

Note that we used the first-order development ln(1 + x)  x for x  1, because the displacements are small. The above equation can be used simultaneously with the first limit of elasticity pffiffiffiffiffiffiffi a1 ðrrr þ rhh þ rzz Þ=3 þ 3J 2 ¼ q0 to evaluate the solution at the plasticity threshold. The inelastic phase involves the plastic rate of displacements. Assuming that the small deformation theory holds and combining Eqs. (38) and (37)–(39) results in the stress state:

0.15

Numerical Analytical

0

0.5

1

1.5

Axial strain, εzz, % Fig. 2. Variation of the deviatoric stress with respect to the axial strain in a case of quasistatic triaxial loading of a cylindrical sample with unit radius and length.

0.15

0.1

0.05 Numerical Analytical

0

0

1

2

3

4

Radial strain, − εrr, % Fig. 5. Variation of the hydrostatic pressure with respect to the radial strain in a case of quasistatic triaxial loading.

110

A. Karrech et al. / Computers and Geotechnics 38 (2011) 105–112

Δ t = 0.001

radial strains decrease, as can be seen in the three figures. Once the desired pressure is reached, a deviatoric stress is applied. The sample then plastifies and results in a nonlinear pattern, as depicted by Figs. 2–5. It can also be seen that a good agreement between both the analytical and the numerical methods is obtained. This demonstrates the accuracy of the numerical approach developed in this study in predicting the response of an elasto-plastic material obeying the yield and plastic functions of Drucker–Prager.

4

In this paragraph, we study the response of the same sample (see Fig. 1) but we subject it to a confining pressure of 0.07 MPa and a deviatoric stress of 0.275 MPa. We also use the same material properties and add a fluidity parameter f = 5 MPa s to take into account the rate-dependency. Unfortunately, an analytical solution turned out to be difficult to obtain because of the additional nonlinear viscous terms. Fig. 6 shows the sample’s response under different strain rates. The curves indicate that the material becomes stiffer when the loading rate increases. It is worth noting that curve (a) in the figure represents the rate-independent case, which is identical to the results shown in Fig. 2.

Δ t = 0.05

Δ t = 0.1

x 10

2 (ε ) zz

0

Strain

4.2. Cylindrical sample under monotonic triaxial loading rates

Δ t = 0.01

−3

−2

(εrr )

−4 −6 −8

0

0.2

0.4

0.6

0.8

1

Incremetal time Fig. 7. Numerical accuracy of the model in a case of typical cyclic triaxial testing.

−3

4

This example consists of subjecting a cylindrical sample to a confining pressure of 0.07 MPa and a sinusoidal deviatoric loading of the form S(t) = S0(1  cos(x0t)), where the loading amplitude is S0 = 0.09 MPa and the circular frequency is x0 = 4p s1 (see Fig. 11). The rate-dependent effect is not taken into account in this particular example. To demonstrate the accuracy of the model, different time increments were used. The shortest time increment produces results closer to the exact solution. Fig. 7 shows the strain response corresponding to the applied load, which indicates good accuracies for different time increments. The analysis also treated the variation of the radial and axial displacements as well as the volumetric change with respect to the number of cycles, as shown in Figs. 8–10. Figs. 8 and 9 show that the permanent displacements decrease with respect to the number of cycles. This result was also observed by Karrech [18], using the molecular dynamics (discrete element method) to study the ratcheting effect in granular materials subjected to cyclic

2

x 10

Radial displacement, u

r

4.3. Cylindrical sample under cyclic triaxial loading

0 −2 −4 −6 −8 −2

−0.5

0 5

x 10

Fig. 8. Cyclic radial displacement versus deviatoric stress.

−3

x 10

3

x 10

2.5

(f) (e) (d) (c) (b) (a)

2.5

Axial displacement, uz

, MPa dev

Deviatoric stress, q

−1

Applied deviatoric stress, q

5

3

−1.5

2 1.5 1

2 1.5 1 0.5 0

0.5 −0.5 0

0

0.005

0.01

0.015

0.02

0.025

Axial strain, ε

zz

Fig. 6. Effect of the rate-dependency on the response of a confined granular sample.

−1 −2

−1.5

−1

−0.5

Applied deviatoric stress, q Fig. 9. Cyclic axial displacement versus deviatoric stress.

0 5

x 10

111

A. Karrech et al. / Computers and Geotechnics 38 (2011) 105–112 −3

−3

x 10

0

Axial displacement, u(m)

Volume change, Δ V

5

0

−5

−10

−15 0

5

10

15

−0.5

−1 f = 5 Hz f = 10 Hz f = 50 Hz f = 100 Hz

−1.5

20

x 10

0

2

Incremetal time Fig. 10. Variation of the volume with respect to time, and logarithmic overall decay.

6

8

10

Fig. 12. Response of a partially confined cylindrical sample to cyclic triaxial loading under different frequencies.

loading. On the other hand, Fig. 10 shows the logarithmic compaction with respect to the number of cycles in accordance with the introduced Chicago law.

2

5. Application: partially confined viscoplastic sample under cyclic loading

0

−4

Axial displacement, u(m)

In this section, the developed model is now used to simulate the response of a partially confined cylindrical sample under cyclic loading. The sample has a radius R = 0.075 m and a length L = 0.15 m. It is subjected to a cyclic loading of the form S(t) = S0(1  cos(x0t)), where the circular frequency x0 is varied. Fig. 11 shows that the bottom surface of the sample is constrained from moving in the axial direction ez, and the lateral surface is constrained from moving in the radial direction er. However, the top surface of the sample is partially loaded with an axial stress within the region of r < 0.5R. The material properties presented in the last section are used for simulation. One of the important applications of the proposed model is the prediction of the accumulated permanent deformation in track platforms due to repeated loading. Fig. 12 shows the variation of axial displacements at the central point (r = 0, z = L), with respect

4

Number of cycles, N

x 10

−2

−4

N = 10 N = 100 N = 200 N = 300 N = 400 N = 500

−6

−8

0

0.02

0.04

0.06

0.08

Sample radius, r(m) Fig. 13. Accumulated vertical permanent deformation along the surface of the sample for different frequencies.

to the number of cycles for various loading frequencies. It can be seen that the axial displacement is decaying for all frequencies as the number of cycles increases. This leads to an increasing accumulation of permanent deformation. This result is in agreement with the phenomena observed by Guérin et al. [14], Bodin et al. [4]. The figure also shows an increase of the permanent deformation slope when the frequency increases from 1 to 100 Hz. The accumulation of permanent deformation along the surface of the sample is plotted in Fig. 13 for a frequency of 10 Hz and different numbers of cycles. 6. Conclusion

Fig. 11. Loading and boundary conditions of the flexible sample.

A comprehensive viscoplastic constitutive model describing the behaviour of non-associated ratcheting granular materials is developed and implemented using a finite element system. One of the main features of the model is the incorporation of Chicago’s law, which describes the compaction of granular materials such as the ballasts that are widely encountered in railway tracks. The model

112

A. Karrech et al. / Computers and Geotechnics 38 (2011) 105–112

was based on the Drucker–Prager yield and potential functions and involves the rate-dependent behaviour described by Perzyna’s theory. The integration of the model was performed by applying the well-known predictor–corrector algorithm used for stress state estimation coupled with Newton’s method for computing the consistency factor as well as the consistent tangent viscoplastic operator. The benchmarking of the model was performed through specific examples. Good agreements between the analytical and computational results were obtained in the case of multi-axial loading of a frictional sample. This demonstrated the capability of the model to predict one of the relevant features of granular materials. The model was then used to describe the effects of ratcheting and ratedependency on a specific field case. The logarithmic permanent deformation observed experimentally in cyclically loaded materials was produced. The model developed in this paper requires several increments within a single cycle. This aspect might be seen as a limitation in terms of computational time, but it allows the viscous terms to play an important role in ratcheting. The model can therefore be suitable to analyze local events where limited numbers of cycles are involved. Future research will take into account the dynamic aspects, as recent experimental work has proven that the frequency affects the density profiles in loaded samples. The transition from low to high loading frequencies is of particular interest, as the materials undergo phase transformation from solid-like to fluid-like behaviour due to vibration. References [1] Arulanandan K, Scott RF, editors. In: Conference proceedings on verification of numerical procedures for the analysis of soil liquefaction problems, vol. 1. Balkema, Rotterdam, The Netherlands; 1993. [2] Azema E, Radjai F, Peyroux R, Dubois F, Saussine G. Vibrational dynamics of confined granular material. Phys Rev E 2006;74:031302. [3] Ben-Naim E, Knight J, Nowak E. Slow relaxation in granular compaction. Physica D 1998;123:380–4. [4] Bodin V, Tamagny P, Sab K, Gautier P-E. Détermination expérimentale d’une ioi de tassement du ballast des voies ferrées soumises à un chargement latéral. Can Geotech J 2006;43:1028–41. [5] Ceia J. Analysis of reinforced concrete structures subjected to dynamic loads with a viscoplastic Drucker–Prager model. Appl Math Model 1998;22: 495–515. [6] Cernocky EP, Krempl E. A non-linear uniaxial integral constitutive equation incorporating rate effects creep and relaxation. Int J Non-Lin Mech 1979;14:183–203. [7] Chen W, Han D. Plasticity for structural engineers. New York: Springer-Verlag; 1987. [8] Cundall PA, Strack ODL. A computer model for simulating progressive large scale movement of blocky rock systems. In: Proceedings of the symposium of the international society of rock mechanics, vol. 1(8), Nancy, France; 1971. [9] Cundall PA, Strack ODL. A discrete numerical model for granular assemblies. Geotechnique 1979;29:47–65. [10] Cvitanic V, Viak F, Lozina Z. A finite element formulation based on nonassociated plasticity for sheet metal forming. Int J Plast 2008;24:646–87. [11] Dafalias YF. Elasto-plastic coupling within a thermodynamic strain space formulation of plasticity. Int J Non-Lin Mech 1977;12:327–37. [12] Eekelen HV, Potts DM. The behaviour of drammen clay under cyclic loading. Geotechnique 1978;28:173–96.

[13] Gallas J, Herrmann H, Sokoiowski S. Convection cells in vibrating granular media. Phys Rev Lett 1992;69(9):1371–5. [14] Guérin N, Sab K, Moucheront P. Identification expérimentale d’une loi de tassement du ballast. Can Geotech J 1999;36:523–32. [15] Hofstetter G, Simo JC, Tayior RL. A modified cap model: closest point solution algorithms. Comput Struct 1993;46(2):203–14. [16] Kang G. A viscoplastic constitutive model for ratcheting of cyclically stable materials and its finite element implementation. Mech Mater 2004;36: 213–99. [17] Kanga G, Ohnob N, Nebuc A. Constitutive modeling of strain range dependent cyclic hardening. Int J Plast 2003;19:1801–19. [18] Karrech A. Comportement des materiaux granulaires sous vibration: application au cas du ballast. PhD thesis. Paris (France): Ecole Nationale des Ponts et Chaussèes; 2007. [19] Karrech A, Duhamei D, Bonnet G, Roux JN, Chevoir F, Canou J, et al. A causality analysis of settlement versus contacts loss in vibrated granular materials. Granular Matter 2008;8(1–4):195–204. [20] Karrech A, Duhamei D, Bonnet G, Roux JN, Chevoir F, Canou J, et al. A computational procedure for the prediction of settlement in granular materials under cyclic loading. Comput Meth Appl Mech Eng 2007;197(1–4): 80–94. [21] Key W. A finite element procedure for the large deformation dynamic response of axysimmetric solids. Comput Meth Appl Mech Eng 1974;4:195–218. [22] Kovacevic N, Vaughan DMPP. The effect of the development of undrained pore pressure on the efficiency of compaction grouting. Geotechnique 2000;50: 683–8. [23] Kumar QP, Nukala V. A return mapping algorithm for cyclic viscoplastic constitutive models. Comput Meth Appl Mech Eng 2006;195:148–78. [24] Lehane BM, Jardine RJ, McCabe B. Response of a pile group in clay to first-time one-way cyclic tension loading. In: Conference on advances in geotechnical engineering. London: Thomas Telford; 2004. p. 700–9. [25] Liu GPY. Implicit consistent and continuum tangent operators in elastoplastic boundary element formulation. Comput Meth Appl Mech Eng 2001;190: 2157–79. [26] Lorefice R, Etse G, Caroi I. Viscoplastic approach for rate-dependent failure analysis of concrete joints and interfaces. Int J Solids Struct 2008;45(9): 2686–705. [27] Loret B, Prevost J. Accurate numerical solutions for Drucker–Prager elastic– plastic. Comput Meth Appl Mech Eng 1986;54:259–77. [28] Lu Y, Wright P. Numerical approach of visco-elastoplastic analysis for asphalt mixtures. Comput Struct 1998;69:139–47. [29] Marques JMMC, Owen DRJ. Some reflections on elastoplastic stress calculation in finite element analysis. Comput Struct 1984;18(6):1135–9. [30] Mayama T, Sasaki K, Ishikawa H. A constitutive model of cyclic viscoplasticity considering changes in subsequent viscoplastic deformation due to the evolution of dislocation structures. Int J Plast 2007;23:915–30. [31] Nowak ER, Knight JB, Ben-Naim E, Jaeger HM, Nagei SR. Density fluctuations in vibrated granular materials. Phys Rev E Stat Phys 1998;57(2):1971–82. [32] Philippe P, Bideau D. Numerical model for granular compaction under vertical tapping. Phys Rev E 2001;63:0513041–05130410. [33] Rèmond S. Simulation of the compaction of confined mono-sized spherical particles systems under symmetric vibration. Physica A 2003;329:127–46. [34] Rosato AD, Yacoub D. Microstructure evolution in compacted granular beds. Powder Technol 2000;109:255–61. [35] Roux JN. Elasticity, quasistatic deformation, and internal state of sphere packings. In: 17th ASCE engineering mechanics conference, Newark, DE, June 13–16 2004. [36] Saussine G, Choiet C, Gautier P, Dubois F, Bohatier C, Moreau J. Modelling ballast behaviour under dynamic loading. Part1: A 2d polygonal discrete element method approach. Comput Meth Appl Mech Eng 2005;195(19–22): 2841–59. [37] Simo JC, Taylor RL. Consistent tangent operators for rate-dependent elastoplasticity. Comput Meth Appl Mech Eng 1985;48:101–18. [38] Suiker ASJ, Borst RD. A numerical model for the cyclic deterioration of railway tracks. Int J Numer Meth Eng 2003;57:441–70. [39] Wilkins ML. Calculation of elasto-viscoplastic flow. New York: Academic Press; 1964.