Predicting the behavior of unbound granular materials under repeated loads based on the compact shakedown state

Accepted Manuscript Predicting the Behavior of Unbound Granular Materials under Repeated Loads based on the Compact Shakedown State E. Badakhshan, A. Noorzad, A. Bouazza, S. Zameni PII: DOI: Reference:

S2214-3912(17)30198-8 https://doi.org/10.1016/j.trgeo.2018.05.001 TRGEO 170

To appear in:

Transportation Geotechnics

Received Date: Revised Date: Accepted Date:

29 October 2017 5 April 2018 14 May 2018

Please cite this article as: E. Badakhshan, A. Noorzad, A. Bouazza, S. Zameni, Predicting the Behavior of Unbound Granular Materials under Repeated Loads based on the Compact Shakedown State, Transportation Geotechnics (2018), doi: https://doi.org/10.1016/j.trgeo.2018.05.001

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Predicting the Behavior of Unbound Granular Materials under Repeated Loads based on the Compact Shakedown State E. Badakhshan1*, A. Noorzad1, A. Bouazza2, S. Zameni2 1. Faculty of Civil, Water and Environmental Engineering, Shahid Beheshti University, Tehran, Iran 2. Department of Civil Engineering, Infrastructure Systems, Monash University, Melbourne, Australia *Corresponding Author: Email: [email protected], http://www.sbu.ac.ir. ORCID number: 0000-0001-7462-0200, Tel: + 98(21)7393-2487, Fax: +98(21)7700-6660

Abstract The behavior of unbound granular material under repeated loads primarily depends on the nature and arrangement of constituent particles, particle size and shape. However, very few analytical models are available for predicting the behavior of unbound granular material under repeated loads. This paper presents a generalized failure criterion with a non-associated flow rule which uses from the concept of compact shakedown state for predicting the stress-strain behavior of unbound granular material. At the compact shakedown state, the response of the unbound granular sample to a stress change is reversible and behaves as an elastic material. In an extension to original CANAsand constitutive model, two concepts namely critical state and compact shakedown state play paramount roles as all of the moduli and coefficients are related to these states. In contrast to the critical state, the unbound granular material samples do not experience any plastic deformation at the compact shakedown state and behaves as a purely elastic material and the response to a stress change is reversible. For stress reversal two additional concepts as reflected plastic potential and bounding surface are employed which distinguishes between virgin loading and secondary loading. The extended model is examined by simulating different typical drained tests in conventional tri-axial tests. The applicability of the present model is evaluated through comparisons between the predicted and the measured results. Keywords: Unbound granular material; repeated loads; Critical state; Compact shakedown state. Dedicated to the Memory of Prof. H. Poorooshasb and his supervisor, the late Prof. K.H. Roscoe. 1

1. Introduction The behavior of unbound granular materials subjected to repeated loads remains not fully understood, and is currently the subject of interesting research worldwide. Advanced laboratory testing techniques are specifically developed to understand and characterize soil response to repeated loads [1]. According to Zhao et al. [1], when the applied cyclic load is above the yield limit but lower than a critical load limit, termed as shakedown limit, the structure may exhibit some initial plastic deformation; however, after a number of load cycles, the structure ceases to experience any further plastic strain and responds purely elastically to the subsequent load. In the past few decades, shakedown analyses were generally carried out based on two fundamental shakedown theorems, known as Koiter’s kinematic (upper bound) shakedown theorem and Melan’s static (lower bound) shakedown theorem respectively [36,37]. The advantage of using the shakedown theorems is that detailed stress and strain history are not required during calculation. Numerical tools as well are adapted to this particular problem, to solve complex multi-physics equation systems [38-41]. Soil response to shear loading includes a tendency to change its volume. Volume change is a result of the rearrangement of grains by slip-down (repacking into denser state) and roll-over (repacking into looser state) motions. The dilative or contractive nature of the volumetric behavior of sand is linked to the initial state of the porous medium, in terms of void ratio and stress state. For medium and dense sand monotonically sheared, a contractive phase starts before the dilative motion. The characteristic state or phase transformation state [2] is the moment which marks the limit between contractive and dilative behavior. These two competing ideas are considered equivalent [3]. The concept of critical state [4] is at the base of many constitutive models of current soil mechanics [5]. Whether the initial state of the soil was loose or dense, drained or un-drained, monotonic and cyclic stress paths join the critical state line at sufficiently large strains. After that moment, the volumetric strain rate is equal to zero while large shear strain develops. Many constitutive models may be selected for an overview of existing tools for dynamic modeling of 2

soils [42,43]. Advanced constitutive models suitable for repeated loads were also adapted from existing theories of plasticity, elastoplasticity, viscoplasticity, hypoplasticity [7-9]. Those complete effective stress models describe soil behavior on a range of stress paths. Their development is based on incremental relations between strain and effective stress. Such models predict transient and permanent deformations within the soil mass, and are then much more interesting than the linear equivalent model for modeling catastrophic displacements or collapse of soil structures [6]. A few constitutive models are suitable for modeling unbound granular soils under repeated loads. The viscoelastic equivalent linear method, which is widely used with the program SHAKE91, takes into account soil non-linearity but does not allow strain accumulation or the pore pressure increase. Among advanced elastoplastic models, the multi-surface theory [7] can reproduce cyclic behavior with either two surfaces or an infinite number of nesting surfaces which can account for them. An example of a multiple surface plasticity model is presented by Elgamal et al. [8] for unbound granular soils. Another family of models is the bounding surface theory. Khalili et al. [9] applied it to cyclic loading of sand. Using the concept of multi-mechanisms plasticity [10], the total irreversible deviatoric strain increment is induced by coupled dissipative processes: an isotropic and three deviatoric plastic mechanisms. On the negative side, the elastoplastic Hujeux model [10] does not solve all the problems of the linear equivalent method. In particular, the large strain behavior is not quantitatively well simulated because axial stain rate stabilizes instead of reaching the threshold for cyclic loadings such as repeated load. The set of parameters which fairly reproduces medium strain behavior and stiffness degradation is not compatible with a good representation at large strain, around the failure threshold. This could be a limit of the model. However, it may also be possible to obtain a set of model parameters with better adequacy to the test results, by performing a systematic parametric study or implementing advanced algorithms, e.g. artificial neural networks [11,12]. The difficulty of calibration by curve fitting of such a high number of parameters is the major problem of the multi-mechanism elastoplastic law for repeated 3

loads. Although some converged shakedown limits have been obtained using the static and kinematic shakedown theorems, most of them are calculated on the basis of an associated flow rule and are assumed that the materials behavior as elastic perfectly-plastic. However, it is well known that unbound granular materials exhibit non-associated plastic behavior [44,45]. For an elastic-perfectly-plastic model, the stress state satisfying the yield condition will become unchanged, but will cause unlimited plastic strain. In other words for perfect plasticity, model does not include strain hardening or softening effects for the unbound granular material. This paper aims to develop a new general failure criterion for unbound granular material based on consideration of the new concept of compact shakedown state. It is adopted the CANAsand criterion [13] with non-associated flow rule to quantify the mobilized strength in a frictional material. The CANAsand criterion features a clear elsatoplastic failure mechanism which extended to consider the reversal of the loading. The combination of the two approaches will lead to a new two-dimensional failure criterion featuring not only a clear elastoplastic concept for the failure of materials, but also all the parameters introduced in the criterion can be conveniently determined by conventional laboratory tests. The new failure criterion offers reasonable predictive capability on describing the stress-strain behavior for unbound granular material under repeated loads, as will be demonstrated in the subsequent sections.

2. CANAsand constitutive model The constitutive models simulate the soil behavior in a single material point, i.e. the resulting stress increment due to an applied strain increment. In this study, the analysis was performed by using an extended CANAsand constitutive model to consider dependency of dilatancy on stress level and void ratio. The extended non-associated constitutive model assumes that the state of an unbound granular element can be completely represented by the state variables e, the void ratio and parameters related to the state of effective stress acting on the sand element that are strongly depended with pressure as shear stress (τ), normal stress (σ) and analogous to lode’s angle (θ) 4

ranging from –π/6 to π/6. Hence the state of an unbound granular element is represented by the three quantities (τ, σ, θ) and the void ratio e of the element and the state space is four dimensional of (e, τ, σ, θ). Conventional tri-axial compression shear mode is employed as a reference state (shear mode). A brief introduction of this CANAsand constitutive model is helpful for the subsequent description of extended criterion. The response of unbound granular material in the proposed model is an elastic-plastic constitutive model derived from the CANAsand constitutive model, which uses a non-associated flow rule along with the concept of the state boundary surface possessing a critical and a compact shakedown state. Therefore, the CANAcand failure criterion is assumed to govern the failure stage of the unbound granular material. The CANAcand constitutive model is fairly robust in solving complex elastoplastic problems in geotechnical engineering. It was first proposed by Poorooshasb et al. [13] and ever since it has been developed and extended by Poorooshasb and Noorzad [14] that the model can accommodate quite large deformation of unbound granular material. It is an incremental elastic-plastic model with the plastic shear strain as the hardener which utilizes a non-associated flow rule. In this model, the flow rule of material described in terms of state parameters are considered as an effective stress that consists of stress invariants and void ratio of unbound granular element. Furthermore, capabilities of the main features of this model were validated against a number of laboratory results [15-18,28-32] based on the approach of integro-differential equation, which was written in finite difference form and was developed to estimate the total displacement field of the system. In the present study, an incrementally scheme is developed.

2.1. Critical state surface The concept of the critical state which also known as the ultimate state influenced by the studies of the Casagrande [35], in the three dimensional frame of reference (p, q, e), which noted the existence of a line which was termed the critical state line. In the four dimensional space the line transforms into a surface which refers to as the critical state surface. The concept of the critical 5

(ultimate) state postulates the existence, in the four dimensional state space, of a surface for which unlimited distortion of the sample may take place without any change in the state parameter. This surface is referred to as the critical state surface. There exists a surface such that for all states on this surface

   e    0    

(1)

where, parameter ε, which is not a state parameter, represents the shear deformation of the sample and is derived from the second invariant of the deviator strain tensor.

2.2. Compact shakedown state Soils behave with a hysteretic component, which means that some energy is dissipated during repeated load. The representation of such behavior in the stress-strain plane resembles a loop, which is either closed or open depending on loading conditions (Fig. 1). It is often referred to as a hysteretic loop. Under repeated loads and with stress-controlled motions of constant amplitude, different mechanical behavior can be observed depending on initial conditions and loading amplitude as shakedown, accommodation and ratchet [15]. In elastic shakedown vision the stressstrain loops stabilize on a straight line, giving perfectly linear elastic response without energy dissipation. The material response is plastic while the ultimate response is purely elastic. In plastic shakedown vision the stress-strain cycles, which evolve at the beginning of the loading, eventually stabilize, and become closed loops with no more plastic strain accumulation. The dissipative behavior of the material is maintained. It can take thousands of cycles to reach a true plastic shakedown state. The applied cyclic load is slightly less than the limit to produce collapse after the incremental accumulation of permanent strain. In previous studies it was acknowledged that such behavior exist because a small amount of energy dissipation always appears during repeated loads. Capturing the shakedown behavior for unbound granular material is done with the aid of a new concept called the compact shakedown state, which was conceived after recalling the 6

results of some experiments by Wroth [19] and using the so-called one gravity similarity laws [35]. Worth [19] conducted a series of experiments on samples of an unbound granular medium composed of one-mm diameter steel balls. The samples were tested in the simple shear apparatus and subjected to a large number of stress reversals. Worth [19] observed that as the loading proceeded, the sample became denser, and after a number of cycles the steel balls, which had been placed in the apparatus in a random pattern rearranged themselves in an absolutely regular triangular pattern that is equivalent to the most compact form of assemblage of spherical particles. This observation, which is the basis of the thoughts presented in this study, was only reported in Wroth's thesis and did not attract the attention it deserved. It worth to note that, recently Meidani et al. [20] divided the total volume of voids into two hypothetical fractions active voids volume and inactive voids volume. The active voids volume defined as the fraction of voids that can be reduced, and eventually diminished by particles rotation and displacement. In other words, the fraction of voids is reactive to the kinematic process of particle rearrangement. The other fraction, which is not reactive to the kinematic process of particle rearrangement, is named inactive voids volume. Unbound granular media are not composed of equal diameter of spherical particles and hence, they can rearrange themselves in such a way as to a more compact shakedown state as illustrated in Fig. 2. This figure shows the behavior of an unbound granular medium subjected to a simple shear loading process. It is clear that to reach to compact shakedown state, the shear strain of the unbound granular material must be equal to the γcompact. It is reasonable to assume that the void ratio after all active voids which are diminished (compact void ratio) can be correlated to the minimum void ratio of the material. To make it a useful concept in analysis the following two hypotheses are made: 1. The compact shakedown state traces a line in the Casagrande space of e vs. Ln σ parallel to critical state line. This is an essential consideration if laws of similarity which have been proposed and tested by Roscoe and Poorooshasb [34] and Scott [35] are to be followed. A brief 7

description of similarity law is presented by Poorooshasb [17]. It must be emphasized that this choice of compact shakedown is indeed dictated by certain laws which govern the similarity of behavior of samples formed at the same effective void ratio and subjected to similar loading paths. The choice is not unique but the model presented above appears to be the simplest. Fig. 3 shows the isometric representation of the critical state surface in the (τ, σ, e) as determined by Poorooshasb [14]. Also, in this paper is postulated that the compact state also lie on a surface in the state space and this surface is paralle1 to the critical state surface. The isometric view of the compact state surface is similar in shape to the critical state surface. This is a direct consequence of the first hypothesis and conforms to the law of similarity for an unbound granular material when subjected to the range of loads likely to be encountered in traffic loads. 2. At the compact shakedown state, the response of the unbound granular sample to a stress change is reversible. That is, it behaves as an elastic material. The stabilization of accumulated plastic deformations is termed as shakedown. A significant feature of shakedown are residual stresses in the material that are self-equilibrating which remains constant with time (or number of cycles). At the compact shakedown state, the interlocking between the particles is so strong that slippage does not take place. Consequently, plastic deformation, which is the result of particle slippage, does not occur. Thus, at the compact shakedown state one has ε(τ, σ, e, θ) is equal to zero, where all the parameters between the brackets are at the compact shakedown state. In view of the first hypothesis, a close tie exists between the critical state and compact shakedown state. Since all functions involved in the model are expressed in terms of η and e then the principle of similarity, would be upheld. That is, within the bounds of approximation pointed out in the paper, elements with equal η and e would exhibit similar flow behavior when subjected to similar stress path increments. Consequently, the relationship between the void ratio and σ, at the compact shakedown state may also be written in the form of Eq. (2), where, ec = eh - c and c is a constant value representing the vertical distance between the compact void ratio line and the Casagrande 8

critical void ratio line, eh is the initial void ratio and ec is void ratio in compact shakedown state as shown in Fig. 3.

e  ec  .ln 

(2)

2.3. State boundary surface All possible states that an element may assume are enclosed within a boundary called the state boundary, which an unbound granular sample fails at a point on that surface [9]. A simplified version to analyze the maximum shear stress at the state boundary surface under a constant stress condition is shown in Fig. 4. State paths end at the critical state where the loose samples contracting and the dense sample expanding. The total energy stored in this system would be equal to the energy spent in friction plus the work done in moving the block and therefore the work would be:

     v

(3)

and which upon dividing the equation by δε, it would result in τmax. Thus τδε=μσδε is situation at the critical state and μ is tanφcritical. It may be noticed that the peak value of the state path is in the proximity of the state boundary surface. Hence, at peak stress, for a given value of σ one can write the following equation:

 max   SBS   critical   (dv / d )   .tan critical   (dv / d )

(4)

where, (dv/dε) represents the rate of dilation, which in other words is the rate of change of volume with respect to the shearing strain. The component σ(dv/dε) is known as the energy correction factor. This factor is related to the behavior of dense gravel and would be explained by considering two rough blocks sliding at an angle relative to the direction of the imposed movement, which is the direction of shearing. The particles in loose sand tend to seek a more compact shape on application of the shearing stress, whereas the volume of the dense sample tends to increase, because the particles must either fracture or raise out of their position to pass by one another, thus leading to an increase in volume. It can be concluded that the increase in volume of 9

dense sand is due to the strong interlocking between the particles at the compact state. Thus by considering the situation in Fig. 4, one can obtain the angle of interlocking between the particles (δ) in a (τ, e) space, where the void ratio is denoted by ecompact at the compact state and ecritical at the critical state. Thus, by considering the situation relating to the peak stress point in Fig. 4 which occurs in compact state and is enclosed within the state boundary surface, one may write:

tan   (  critical ) / (ecritical  e)

(5)

where, φcritical is the angles of friction at the critical state. If the void ratio e of a sample is greater than the critical void ratio, then the expression within the brackets is negative and indicates a tendency of the unbound granular material to contract, and if the void ratio is smaller than the critical void ratio then the expression is positive indicating the material has a tendency to dilate. If a soil is very dense to dense, in a drained test, the associated stress path would approach the state boundary surface and then onwards, if the loading proceeds, moves on it until the critical state is reached.

3. Generalized constitutive model In a two dimensional deformation, the invariants of stress are τ and σ and those of strain are ε and ν, along with their associated increments, dτ, dσ, dε and dν. Strain increments experienced by an element have elastic and plastic components. Moreover, the strain increment is composed of two strain increments which are the shear strain dε (associated with the change of shape) and the volumetric strain dν (associated with the change of volume) and each have two components composed of elastic and plastic parts. The associated flow rule postulates the normality of the incremental plastic strain vector to the failure surface. This normality follows from the convexity of the failure surface. Accurate modelling of dilatancy requires that the incremental plastic strain vector is not normal to the failure surface. The best way to describe this mathematically is to introduce an additional family of surfaces, which is called the plastic potential. It is assumed in classical plasticity theory that the plastic potential function takes the same form as the yield 10

function; i.e., associated flow is employed. Experimental evidence from tests on several types of frictional materials has clearly indicated that the use of associated flow rules results in prediction of too large volumetric expansion [13,14]. To characterize the volume change correctly, it is necessary to a non-associated flow rule is employed. The plastic potential surfaces therefore do not coincide with the yield surfaces. In this study, the utilized plastic potential ϕ(σ,τ,e) which defines the direction of the incremental plastic strain vector in the tri-axial stress space, has two properties. First, through each point of the failure surface passes one of the surfaces of the plastic potential family and second, the incremental plastic strain vector at each point is normal to the plastic potential and not to the failure surface. Thus, the plastic strain increment is derived from the plastic potential in the form of    ( ) , where η is the stress ratio (τ/σ) and a particular form of the function referred to as the critical state formulation is  ( )   exp( /  )   0  0 , where μ is the slope of the projection of the critical state line into the (τ, σ) plane and σ0 represents the abscissa of the point of intersection of a ϕ(η)= constant contour with the σ axis. In its simplest form, the yield function is a straight ray passing through the origin of stress space (τ, σ) and is given by the simple expression of Eq. (15):

f ( ij )   ( ij ) 

 [( x   z )2  ( xz   zx )2 ]0.5   ( x   z )

(6)

where, η stated earlier as the stress ratio.

Shearing and volumetric strain increments in an unbound granular material are derived from the plastic potential using Eq. (7):

   dv p  d   d p  d 

(7)

where, ϕ is the plastic potential, dεP the plastic shear strain, dνP the plastic volumetric strain and dβ is: 11

f  f  d  h d  d     

(8)

in which, h is the loading index (or proportionality modulus during loading) as defined, dβ is an incremental quantity or constant relating strain components to the gradient of plastic potential (as the consistency parameter) and f is the yield function. The parameter h, which is denoted in Appendix A, is proposed as follow,

  aSBS h  2  (SBS   ) exp( /  ) 

(9)

where, a is function of state, η the stress ratio and ηSBS is the value of the stress ratio at the state boundary surface (SBS) corresponding to the current value of σ and e. The physical interpretation of the loading index h depends on whether the soil element is loading, unloading or reloading according to Poorooshasb [21]. The dependence of h on the void ratio is through the parameter a. Parameter a in the model formulation describes the isotropic-kinematic hardening law for the plastic moduli. It is a function of void ratio as:

a  a0

(e  ecompact ) (10)

(ebase  e)

where a0 is assumed to have a value of 0.1 [18,21] and e is the void ratio. Superimposing the two space stresses (τ, σ) and (εp, vp) on each other, the combine space cover with a set of plastic potential curves together with a set of yield loci is shown in Fig. 5. The normal to each of these curves at each point stress point would indicate the direction of unit vector demonstrating the rate of plastic flow (normality condition). The directions of the plastic strain increment vectors are also shown. The difference from the directions implied from normality is dramatic. Thus yielding at low stress ratio implies volumetric compression but the rate of volumetric compression steadily decreases as the stress ratio increases. From the point of view of soil stiffness, two types of evolution are possible upon repeated loads (Fig. 6) as cyclic hardening which stiffness of the soil increases with higher secant moduli and cyclic softening which stiffness decreases and the 12

inclination of the cycles decreases. The evolution of the secant modulus, either increase or decrease, does not necessarily correspond to true hardening or softening of the soil, respectively. The definition of hardening, with respect to the evolution of the yield surface, can be expressed as any stress path that produces plastic deformation inside the failure surface results in hardening [22]. More precisely, the stress path can push out the yield surface, i.e. hardening occurs, while the secant modulus decreases, because of the coupling of the volumetric and deviatoric shear. Anyway, for the description of repeated load tests results, the expressions cyclic softening and cyclic hardening will be used as a simple way to evoke decrease and increase of soil stiffness, respectively. In an incremental process, plastic strain increments are derived from by the following set of Eqs. (11) to (13): d  xp 

h     ( d x  d z  2 d xz )   x  x  z  xz

d  zp 

h     ( d x  d z  2 d xz )   z  x  z  xz

(12)

d  xzp 

h     ( d x  d z  2 d xz )   xz  x  z  xz

(13)

(11)

A Mohr’s circle of stress is considered with the stress components taken as invariants of the stress tensor in two dimensions. The shear stress, τ, is the radius of the Mohr’s stress circle and σ is the distance between its center and the origin of the stress space. The quantities in above equations can be derived with respect to each stress component with the plastic potential as follows.

x 

 1   x   z       x    

z 

 1  ( x   z )       z    

 xz 

 1  2    xz   xz    

(14)

(15)

(16)

13

x 

    1    z    exp   1   x    x       

z 

    1  ( x   z )    exp   1      z       

xz 

    1  2   exp     xz  xz     

(17)

(18)

  

(19)

Also from the theory of elasticity, stresses are related to strains by an elastic constitutive law. For a homogeneous isotropic elastic material, there exist only two materials constants (E, the Young’s modulus, and ν, the Poisson’s ratio). In an isotropic elastic material the matrix of coefficients relating stress to strain is symmetric and contains only two constants E and ν. In the elastic-plastic model, not only the matrix is non-symmetric but it is also, containing nine independent coefficients. The elastic strains have to be incorporated in the matrix; otherwise the matrix will become singular. Thus the relation between the strain and stress increments can be represented in a matrix form.  1 h  E   x x  d x    d      h    z  E  z x  d  xz    h   xz x 







h

x z

E  1 h  z z E  h



xz x

     d x   h z 2 xz   d z     d xz   1 h  xz 2 xz  2G   h

x 2 xz

(20)

where, ηx, ηz and ηxz are the derivative of the yield function with respect to σx, σz and τxz, and ϕx, ϕy and ϕxz are the derivative of plastic potential with regard to σx, σz and τxz, too. This matrix helps in obtaining the stress increments as a function of strain increments and will be useful in the next section where the substitution of these stress increments will occur in the incremental equations of equilibrium of an element within the unbound granular material.

4. Characteristics of the model The proposed model consists of general properties and specific model parameters that are 14

presented in Table 1. To characterize the elastoplastic response of unbound granular material, model parameters are required to calibrate the constitutive equations. Some index parameters and empirical coefficients are considered to obtain a plausible stress-strain relationship. Adopting that the angle of friction is a function of the void ratio given by

  critical  (compact  critical )

(ecritical  e) c

(21)

where e is the void ratio, φcompact and φcritical are the angles of friction at the compact shakedown state and the critical state, respectively. Parameter c is the vertical distance between e critical and e compact shakedown in the e-logσ graph or the difference in the Casagrande critical void ratio line and the compact shakedown void ratio line, which in this study it is assumed to be 0.5 in accordance with Poorooshasb [21]. The Young’s modulus is taken to be a function of the current state of stress and other index parameters, such as void ratio, identifying the current state of the sand. In the model, the relationship chosen to characterize this variation of Young’ modulus with stress and void ratio is defined as E  E0 (1  

(ebase  e) ) (e  ecompact )

(22)

where, E0 is the initial young modulus and σ is the distance between the center of Mohr’s circle of the stress space. ebase and ecompact are void ratios at the baseline and at the compact shakedown state. In the complete representation of the state boundary surface for two-dimensional stress cases, the baseline designates the loosest possible state a sample may have. The void ratios higher than the one indicated by the baseline cannot be experienced by an element. ebase and ecompact are written as ebase  ehbase   ln 

(23)

ecompact  ehcompact   ln 

(24)

The value of ehbase and ehcompact are void ratios at the baseline for σ=0 (kPa) and void ratios at the baseline for σ=0 (kPa) at compact shakedown state and λ is assumed to be very small [18,21]. The 15

form of the isotropic-kinematic hardening function, h, was derived earlier from the plastic strain curve expressed by hyperbolic stress-strain relation. ηSBS is the maximum value η which η is the stress ratio (τ/σ), the sample may attain for a given value of e. The stress ratio at the state boundary surface, ηSBS, is expressed as

SBS    (compact   )

(ecritical  e) c

(25)

where, μ is the slope of the normal projection of the critical state line and μcompact is the slope of the normal projection of the compact shakedown state line in the (τ, σ) space and are also related to ϕcompact and ϕcritical. ϕ is the angle of friction of the material. Critical void ratio (ecritical) is void ratio, at which a cohesionless soil can undergo deformation or actual flow without any volume change. It is a function of confining pressure, i.e. an increase in confining pressure yields a lower value of critical void ratio. A soil with a void ratio above the critical void ratio line for the confining pressure is subject to flow failure if undergoes sufficient stress and it is written as ecritical  eh   ln 

(26)

where, eh is considered to be void ratio at critical state for p=0 (kPa) and λ has a very low magnitude. Similarly, the other assumption used for a void ratio higher than the void ratio at the critical state and closer to the void ratio at the baseline, is expressed as Eq. (27), where b is the vertical distance between ebase and ecritical in the e−logσ graph and it may be chosen as 0.2 in accordance with Poorooshasb [21].

SBS  

(ebase  e) b

(27)

Since this is an incremental analysis, the incremental volume increase/decrease or the change in volume will be added to the void ratio of the material as e  eold  de

(28)

where, eold is the void ratio of the material at latest increment and de is the incremental change in volume. The state space is in the form of a four-dimensional space and therefore it was virtually impossible to utilize a two-dimensional graphical presentation illustrating different surfaces in 16

this model. For reducing the number of dimensions, a state variable maintains constant. Thus in a constant void ratio condition the yield and state boundary surfaces, the positions of compact shakedown and critical state surfaces (which were already decreased to bands in this state) are shown in Fig 7a. In this study it is assumed that the size of the yield surface is governed by the size of the bounding surface i.e. its size is a fraction of the size of the bounding surface. Also, the bounding surface must, necessarily, be enclosed by the state boundary surface and hence it must conform to the same curvature that the state bounding surface is subjected to. Thus, it is postulated that the bounding surfaces are similar in shape to the state boundary surface except that their relative position is controlled by the current state of the sample assuming the loading to be in its virgin phase. Also Fig. 7b shows an isometric view of the state boundary surface (SBS) and the critical and the compact shakedown state lines plotted in the state space (τ, σ, e). The portion of the SBS that lies to the left of the critical state line is straight meaning that for a constant σ, the SBS produce a straight line passing through the critical and the compact shakedown state points. The slope of the line increases as the value of σ increases. The portion of the SBS that lies to the right of the critical state line is curved. Note that the base line indicates the loosest possible state a specimen may have, i.e., void ratios higher than those indicated by the base line cannot be experienced by an element. If the void ratio of a soil is near the compact shakedown void ratio, the soil is referred to as very dense and if the soil void ratio in smaller, but close to the critical void ratio, then it is considered loose.

5. Stress reversal for material To deal with stress reversal, two additional concepts are employed. The first is the concept of reflected plastic potential and the second that of the bounding surface [23,24] which distinguishes between virgin loading (σ1 > σ3) and secondary loading (σ1 < σ3) such as stress reversal and repeated loads. The behavior of the soil within the boundary surface is guided by the properties of 17

a point on the bounding surface. During virgin loading, the bounding surface moves (expand) as the stress point traces a curve in the stress space. During a stress reversal the yield surface moves within the bounding surface and its kinematic is guided by two particular stress vectors known as the conjugate and datum stress vectors located on the π plane passing through the current stress points are known as the conjugate and datum stress points. Referring to Fig. 8a, let AB represent a portion of the plastic potential contour passing through the current stress point. Then by definition, the reflected plastic potential is the contour A′B′ which is the mirror image of AB about the line passing through the origin of the stress space. In a stress reversal process, the plastic strain increment vector is given by:

r f f ( d  d )     f f dv p  hr r ( d  d )    d  p  hr

(29)

where, the ϕr is the reflected plastic potential function. The loading function h is replaced by hr to indicate that a stress reversal is in progress. The magnitude of hr is related to the value of hc associated with the conjugate stress point on the bounding surface and datum stress which is also located on the bounding surface through the Eq. (30). The datum stress point is opposite to the conjugate stress point and is the intersection of the line joining the conjugate stress point, a point on the axis of the bounding surface with the mean stress equal to that of the conjugate stress, and the yield surface on the opposite side. hr  R n hc

(30)

where, R is defined in Fig. 8b and n is the coefficient of the interpolation function which must be determined empirically. It is observed that the material response is not particularly sensitive to mild variation of the n parameter. The loading factor hc is a function of state. To incorporate the compact shakedown state, the plastic strains are not present, the factor hc must satisfy the condition that: e  ec ,

(31)

hc  0 18

Thus, if a new parameter of e* = e − ec is defined, then h = h (τ, σ, e*), where h (τ, σ, 0) = 0 can be accepted as a rational formulation. The simplest form of the function hc according to Appendix A, is:

hc 

ae*SBS [2  2 ]0.5 (SBS   )2 

(32)

where

 

  ,    

(33)

Here ηSBS is the stress ratio τ/σ at the state boundary surface. This function clearly satisfies the condition of Eq. (31) and in addition it follows the called hyperbolic form of deformation response.

6. Programming the generalized constitutive model The flow diagram of the generalized algorithm is summarized in Fig. 9. 1- Based on the stress increment and the elastic stiffness matrix, compute the strain increment and the new total stress. 2- According to the new stress state, check whether f ≤ 0 or f > 0, if f ≤ 0, the point experiences elastic loading or unloading, then go to step 4; if f > 0, the point is loading elastoplastically, then go to step 3. 3- Compute the compact shakedown void ratio. If e < ecompact the response of the unbound granular sample to a stress change is reversible and it behaves as an elastic material. Then go to step 4. 4- For f ≤ 0, the hardening parameters remain unchanged and go to step 6. Thus, since the yield function f = f (σij), no yielding take place if [

f f   ]  0 with no associated plastic flow.  

5- For virgin loading if f > 0, the point is yielded at the start of the current increment. In other words, a material is said to undergo plastic deformation if the state of stress rate is such that it 19

increases the magnitude of the yield function f. In its simplest form, the yield function is assumed to be a function of state of stress and lately it has been assumed to be a function of state of sample i.e. τ, σ and e. If [

f f   ]  0 leading the plastic flow (deformation). Thus  

the strain can be calculated as Eq. 34, where the symbol < > stands for singularity brackets or as Macauley bracket, i.e. < δ > = δ if δ > 0 and < δ > =0 if δ ≤ 0. For secondary loading such as stress reversal in repeated loads, the ϕ and h in Eq. 34 should be replaced with the reflected plastic potential function (ϕr) and hr in each increment.

d 

d  f f h  d  d  G    (34)

d  f f dv  h  d  d  K   

6- Record and save the stress state, strain, effective plastic strain and other state variables. This increment ends. Since this is an incremental analysis, all the increments defined above are to be added to the increments at each stage of calculations.

7. Validation of the proposed constitutive model The extended constitutive model has a total of 15 constants which all the relevant model parameters in the criterion can be conveniently determined through conventional laboratory tests, such as tri-axial tests and simple shear tests. In this work, different sets of resistance data and volumetric behavior data, obtained from different published experimental studies and discrete element method work, are processed and analyzed. Since the aim of this research is to predict the behavior of unbound granular material under repeated loads, the raw data have been taken from the works that have reported the resistance data. The details of the processed datasets are summarized in Table 2. The procedures for the validation are hence presented in two separate subsections as follows, for repeated loads, respectively. In each subsection, different tri-axial and direct shear tests packages are used to demonstrate the process for materials. 20

7.1. Toyoura sand The mean size of the particle for Toyoura sand is 0.17 mm and the particle shape of this sand under microscope (40x) from Yin et al. [25] for repeated load tests is shown in Fig. 10. As presented in Table 2, the values amounts of parameters corresponding to the general framework are kept the same in all simulations. Fig. 11 shows comparisons between experiments and predictions for the drained repeated load tests on loose Toyoura sand with void ratio e= 0.845 under constant confining pressure σ= 98 kPa. Fig. 12 indicates comparisons for a dense Toyoura sand with void ratio e= 0.653 under constant confining pressure σ= 98 kPa. Good agreement, only with slight overestimation of the stress was achieved between experimental data and simulations results.

7.2. Micromechanical behavior using DEM The micromechanical behavior of unbound granular materials during repeated load using the discrete element method (DEM) was studied by Sazzad and Suzuki [27]. Oval particles were used to model the samples. Dense samples were subjected to traffic loads with constant axial strain amplitude (εy= ±0.50%). During loading, the sample height was slowly reduced vertically downward by maintaining a lateral stress constant (100 kPa) with a small strain increment of 0.0001% in each step. During unloading, the stress direction was reversed, and the sample height was slowly increased vertically upward by maintaining again the same lateral stress with the same strain increment. The soil parameters are selected based the Sazzad and Suzuki [27] study and are presented in Table 2. Due to lack of results of laboratory tests, the required parameters are determined by trial and error to fit the DEM data with the response of the proposed model. The stress–strain behavior of unbound granular materials during repeated loads from DEM simulations and extended model in current study are depicted in Fig. 13a, where, τ is the shear stress and σ is the mean stress where σ11, σ22 are the average stresses in x and y directions, respectively. 21

Fig. 13b shows the evolution of the volumetric strain with the axial strain. The positive value of the volumetric strain indicates compression, while the negative value indicates dilation. As a consequence, it is observed that the proposed model can reasonably predict the behavior of samples responses in both loading and reverse loading phases.

7.3. Toyoura sand in monotonic triaxial tests The values of constant parameters for optimal simulation of dry-deposited Toyoura sand response, as reported by Yoshimine et al. [26] are presented in Table 2. This type of sand has high angularity and is widely used in different studies. Particle shapes of this sand sample under microscope (40x) is displayed in Fig. 14 [33]. G0 and E0 are related to the elastic branch of behavior. They can be determined by the results of very small strain tests such as resonant column or bender element tests. In the absence of such data, one may estimate these parameters by constructing tangent to the very beginning parts of shear stress-shear strain and volumetric-shear strain curves in drained triaxial or simple shear tests. The critical state parameters can be obtained directly from the critical state stress ratio in triaxial compression and extension tests and the location of the critical state line in the e-σ plane (for ecritical and λ). The parameters a0, b, c and n can be determined by trial and error to fit the monotonic triaxial compression tests. It is found that these parameters are closely related to the particle constitution of sand such as gradation, maximum and minimum void ratio [21,24]. The parameters ehbase and ehcompact can be directly obtained based on the location of the void ratio for isotropic compression in the e−logσ space and it can be estimated with the conventional index tests in a geotechnical laboratory [ehcompact= (1.2−0.8)emin]. Since there is no test data available for compact shakedown state of Toyoura sand in constant stress ratio compression, parameter φcompact is assumed about 45° based on the best fitting of the isotropic compression stress-strain curve. Fig. 15 presents the results of triaxial compression test and simulations. This is clear proof of how well the model simulates the measured response, based on the effect of compact shakedown state on the failure criteria and 22

deformations. It can be concluded that, the proposed model follows well the trend of experimental test in term of stress-strain curve. The result of numerical simulation of on loose and medium dense samples having an initial void ratio of 1.05 and 0.75, respectively, is presented in Fig. 16. As al1 the simulations are under a constant normal stress of 100 kPa, at pressure, the critical void ratio can be calculated about 0.85. It is chosen to demonstrate the capabilities of the mode1 in verifying the proposed formulation. Fig. 16a shows the entire stress-strain curves for drydeposited Toyoura sand under the 2700 and 1500 cycles of stress paths, for both loose and medium dense samples, respectively. This curves show the tendency of being more elastic as more load cycles are applied. Essentially, the sand under cycles of stress paths was stabilized and elastic shakedown was reached. In other words, it reaches the compact shakedown state and behaves as elastic solid. The sand subjected to cycles of stress paths where a higher void ratio e0= 1.05 is used shows more plastic deformation. In the other words medium dense samples tended to approach the compact shakedown state in lower cycles than loose sample. The axial strain versus the number of cycle curves for Toyoura sand is demonstrated in Fig. 16b. It shows that the proposed compact shakedown state is capable of capturing the long-term behavior in permanent deformation under repeated loads. In particular, most axial strains are taking place in the first 1,500 load cycles. When the unbound granular materials are contained of higher void ratio, more load cycles were needed to reach a stable state shakedown. It may be noted that, although the medium dense sample still contracts (its void ratio decreases as the number of cycles of loading increases), nevertheless, the degree of contraction is not high in comparison with loose sample. The results show unbound granular material with a void ratio (e0= 1.05) above the critical value for the confining pressure is subject to shakedown if it undergoes sufficient stress. However, if the void ratio is lower than the critical value (ecritical= 0.85) for the given confining pressure, the soil will experience compact shakedown state under stress increase. Deposits denser than the critical void ratio could also meet shakedown limit. The tendency to contract exists even in these deposits. Dilatation (increase in volume) manifests itself only if the deposit is subjected to large 23

amplitude repeated loads. For a very dense sample which is subjected to a high amplitude repeated load at the initiation of each cycle, there is a small tendency for the sample to contract. As the magnitude of the imposed stress increases, the sample dilates and the overall effect is an increase in the void ratio of the sample. At some stage, however, it reaches the state boundary surface and fails. After this limit if the straining continues, the sample would ride on the state boundary surface and behave as a work softening material without any shakedown limit.

8. Conclusions Non-linear behavior of unbound granular material and the stress-strain relationship have been described and derived along with the constitutive equations to analyze of elastoplastic behavior of the unbound granular material under repeated loads. The formulation of proposed constitutive equations has been developed based on the compact shakedown and critical states along with the state boundary surface. The compact shakedown state is a new concept that denotes when the interlocking between the particles is in stronger mode that slippage does not take place, the response of the unbound granular sample to a stress change become reversible and behaves as an elastic material. The proposed failure criterion has been applied to the prediction of failure for a number of different sands and the obtained results indicate that the predictions by new criterion are in good agreement with the experimental data. Also, model predictions for drained monotonic and repeated load tests on Toyoura sand and DEM simulations have demonstrated that the present proposed approach is capable of modeling the behavior of unbound granular material under repeated loads the results shown the two additional concepts as reflected plastic potential and bounding surface well captured the trend of drained behavior and stress reversals under repeated loads. Also medium dense samples tended to approach the compact shakedown state in lower cycles than loose sample. In the other words, when the unbound granular materials are contained of higher void ratio, more load cycles were needed to reach a stable state shakedown. Also, if the representation of the void ratio of a sample from an unbound granular material deposit falls above 24

the critical state line, then the shakedown of the deposit, when subjected to repeated loads, is certain. If the representation indicates a point below the critical state line, it does not mean that the deposit does not meet shakedown limit. The shakedown of the deposit against repeated loads depends how far the state point is from the line presenting the compact shakedown state and the magnitude of the repeated loads.

Appendix A. Calculation of the parameters h and hc The value of parameter h depends on whether the element is loading, unloading or reloading. In a loading process, its value is determined using the plastic strain curve, which is expressed by hyperbolic stress-strain relation in the form: 

 p  a



  SBS

   

(A.1)

It may be noted that when η → ηSBS, ε → ∞ and ηSBS indicates the failure of the element at the state boundary surface. When using Eq. A.1, the strain hardening and strain softening of the medium can be modeled and in the case of the strain softening, the loaded element may have to ride on the state boundary surface till it reaches the critical state. In a constant normal stress test the following is derived from the plastic strain curve:

(   )   aSBS d p  a SBS  2 d (SBS   ) (SBS   )2

(A2)

and

d p 

aSBS d (SBS   )2

(A.3)

On the other hand, for a constant σ test, using the critical state formulation of the plastic potential:

d p  h

   1 h d  h exp( /  ) d  exp( /  )d     

By equating the right hand side of Eqs. A.3 and A.4:

25

(A.4)

aSBS h d  exp( /  )d 2 (SBS   ) 

(A.5)

Therefore the parameter h calculates as:

h

(SBS

aSBS   )2 exp( /  )

(A.6)

For stress reversal condition, the magnitude of hc is related to the conjugate stress point on the bounding surface through the following equation:

 / exp( /  )  (e  ecomp ) 1   2

 

1

 / exp( /  )  e*

(A.7)

1  [ ]2

 [

1



(A.8)

]

exp( /  ) [

 / exp( /  )  e* [

1



1



]2  [



  2

]2 (A.9)

exp( /  )]

By using the function of plastic potential  ( )   exp( /  )   0  0 , it can be concluded that: [

 / exp( /  )  e*

1



exp( /  )]2  [ [

1





  2

exp( /  )]2  e*

exp( /  )]

[2  2 ]0.5



(A.10)

Thus, from Eq.A.6, the parameter hc can be calculated as:

 aSBS  [2  2 ]0.5 hc   2    (SBS   ) 

(A.11)

References [1] Zhao J, Sloan SW, Lyamin AV, Krabbenhoft K. Bounds for shakedown of cohesive-frictional materials under moving surface loads. International Journal of Solids and Structures 2008;45(11):3290-3312. [2] Ishihara K, Tatsuoka, F, Yasuda S. Undrained deformation and liquefaction of sand under cyclic stresses. Soils and Foundations 1975;15(1):29-44. [3] Chu J. An experimental examination of the critical state and other similar concepts for granular soils. 26

Canadian Geotechnical Journal 1995;32(6):1065-1075. [4] Schofield A, Wroth P. 1968. Critical State Soil Mechanics. McGraw-Hill, London. [5] Muir Wood D. The magic of sands. Canadian Geotechnical Journal, 2007;44(11):1329-1350. [6] Muir Wood D. Deformation properties of soils for dynamic analyses. 6th SECED Conference on Seismic Design Practice into the Next Century, Oxford, England, 1998:15-24. [7] Mroz Z, Norris VA, Zienkiewicz OC. An anisotropic, critical state model for soils subject to cyclic loading. Geotechnique 1981;31(4):451-469. [8] Elgamal A, Yang ZH, Parra E. Computational modeling of cyclic mobility and post liquefaction site response. Soil Dynamics and Earthquake Engineering 2002;22(4):259-271. [9] Khalili N, Habte MA, Valliappan S. A bounding surface plasticity model for cyclic loading of granular soils. International Journal for Numerical Methods in engineering 2005;63:1939-1960. [10] Hujeux JC. Une loi de comportement pour le chargement cyclique des sols. In Génie Parasismique, Paris. 1985:287-302. [11] Obrzud RF. Numerical modeling and neural networks to identify constitutive parameters from in situ tests. EPFL, Lausanne, PhD thesis, 2009. [12] Rascol E.Cyclic properties of sand: dynamic behavior for seismic applications. EPFL, Lausanne, PhD thesis, 2009. [13] Poorooshasb HB, Holubec I, Sherbourne AN. Yielding and flow of sand in tri-axial compression. Canadian Geotechnical Journal 1966; (4):179–190. [14] Poorooshasb HB, Noorzad A. The compact state of the cohesionless granular media. International Journal of Science and Technology, Scientica Iranica 1996;3(1):1–8. [15] Roscoe KH, Poorooshasb HB. A theoretical and experimental study of strains in tri-axial tests on normally consolidated clays. Geotechnique 1963;13(1):12-38. [16] Poorooshasb HB, Ishihara K, Yang QS. Action of standing waves on seabed sands. Proceedings of the International Conference on Computational Engineering Mechanics, Beijing, China 1987:146-151. [17] Poorooshasb HB. One gravity model testing. Soils and Foundations 1995;35(3):55–59. [18] Noorzad A, Poorooshasb HB. The numerical simulation of flow of bulk solids using CANAsand constitutive model. International Journal of Civil Engineering 2005;3(3):129–139. [19] Worth CP. The behavior of soils and other granular media when subjected to shear. Ph.D. thesis, Cambridge University, U.K. 1958. [20] Meidani M, Chang CS, Deng Y. On active and inactive voids and a compression model for granular soils. Engineering Geology 2017;5(16):70–78. [21] Poorooshasb HB. Subsidence evaluation of geotextile reinforced gravel mats bridging a sinkhole. Geosynthetics International 2002;9(3):259–282. [22] Pastor M, Zienkiewicz OC, Chan AHC. Generalized plasticity and the modeling of soil behavior. International Journal for Numerical and Analytical Methods in Geomechanics 1990;14(3):151-190. [23] Dafalias YF, Manzari MT. Simple plasticity sand model accounting for fabric change effects. J. Eng. Mech. 2004;10(3):622–634. 27

[24] Dafalias YF. A.G. Papadimitriou, X.S. Li, Sand plasticity model accounting for inherent fabric anisotropy. J. Eng. Mech. 2004;11(13):1319–1333. [25] Yin ZY, Chang CS, Hicher PY. Micromechanical modeling for effect of inherent anisotropy on cyclic behavior of sand. International Journal of Solids and Structures 2010;47:1933–1951. [26] Yoshimine M, Ishihara K, Vargas W. Effects of principal stress direction and intermediate principal stress on undrained shear behavior of sand. Soils and Foundations 1998;38(3):179–188. [27] Sazzad M, Suzuki K. Micromechanical behavior of granular materials with inherent anisotropy under cyclic loading using 2D DEM. Granular Matter, 2010;12:597–605. [28] Badakhshan E, Noorzad A. Effect of footing shape and load eccentricity on behavior of geosynthetic reinforced sand bed. Geotextiles and Geomembranes 2017;45:58–67. [29] Badakhshan E, Noorzad A. Load eccentricity effects on behavior of circular footings reinforced with geogrid sheet. Journal of Rock Mechanics and Geotechnical Engineering 2015;7:691–699. [31] Badakhshan E, Noorzad A. Investigation on pullout behavior of geogrid-granular trench using CANAsand constitutive model. Journal of Rock Mechanics and Geotechnical Engineering 2017:1–15. [32] Badakhshan E, Noorzad A. A state boundary surface model for improving the dilatancy simulation of granular material in reinforced anchors. Arabian Journal of Geosciences 2017;10:1–13. [33] Qin X, Zeng H, Ming H. Influence of fabric anisotropy on seismic responses of foundations. Journal of Rock Mechanics and Geotechnical Engineering 2015;1–8. [34] Roscoe KH, Poorooshasb HB. A theoretical and experimental study of strains in tri-axial tests on normally consolidated clays. Geotechnique 1963;13(1):12–38. [34] Scott RF. Essais en centrifuge et technique et de la modelisation. Revenue Francaise de Geotecbnique, Paris 1989;48:15–34. [35] Roscoe KH, Schofield AN, Wroth CP. On yielding of soils. Geoteehnique 1958;8:22-53. [36] Melan E. Der spannungsgudstand eines henky-mises schen kontinuums bei verlandicher belastung. Sitzungberichte der Ak Wissenschaften Wie 1938;147(73):112-123. [37] Koiter WT. General theorems for elastic-plastic solids. In Progress in Solid Mechanics, edited by R. Hill. eds I. N. Sneddon. North-Holland, Amsterdam, 1960. [38] Yu HS, Hossain MZ. Lower bound shakedown analysis of layered pavements using discontinuous stress fields. Computer Methods in Applied Mechanics and Engineering 1998;167(4):209-222. [39] Werkmeister S, Dawson AR, Wellner F. Pavement design model for unbound granular materials. Journal of Transportation Engineering 2004;130(5):665-674. [40] Wang J, Yu HS. Shakedown analysis and design of layered road pavements under three-dimensional moving surface loads. Road Materials and Pavements Design 2013;14:703-722. [41] Sharp RW, Booker JR. Shakedown of pavements under moving surface loads. Journal of Transportation Engineering, ASCE 1984;110:1-14. [42] Li HX. Kinematic shakedown analysis under a general yield condition with non-associated plastic flow. International Journal of Mechanical Sciences 2010;52(1):1-12. 28

[43] Liu S, Wang J, Yu HS, Wanatowski D. Shakedown solutions for pavements with materials following associated and non-associated plastic flow rules. Computers and Geotechnics 2016;78:218-266. [44] Lade PV, Nelson RB, Ito YM. Non-associated flow and stability of granular materials. Journal of Engineering Mechanics 1987;113(9):1302-1318. [45] Lade PV, Pradel D. Instability and plastic flow of soils. I: experimental observations. Journal of Engineering Mechanics 1990;116(11):2532-2550.

29

τ

Elastic shakedown shakedown

τ

Plastic shakedown Fig. 1.. Cyclic behavior of sand under stress stress-controlled controlled motions of constant amplitude.

(a)

(b) Fig. 2. Schematic view of simple shear test on o granular material, material, a) before shearing h=2R′ and b) after shearing and reaching to compact state h′=2R′cosγcompact (where, R R′ is the radius of granular particles) particles).

Void ratio, e

Casagrande CVR line

c

Compact shakdown void ratio line (b) 0.1

1

10

Mean stress, σ, MPa

Fig. 3. The compact shakedown void ratio line.

Fig.. 4. Shear stress at the state boundary surface [21]. 7 Yeild loci

6

Plastic potentials

5 τ , εp

4

Plastic lastic strain vectors

3 2 1 0

0

5

σ, vp/2

10

15

Fig. 5. Curves of yield loci, plastic potentials and plastic strain vectors. vectors

τ

(a) Cyclic softening

τ

(b) Cyclic hardening Fig. 6. Two types of evolution of soil stiffness.

(a)

(b) Fig. 7. a) Isometric sometric view of the yield and ultimate state surfaces and b) the state boundary surface (SBS) for a granular material with critical and compact void ratios. ratios

(a)

(b) Fig. 8. a) Reflected plastic potential and b) bbounding ounding surface and the position of the datum and conjugate stress point.

Fig. 9. Flow chart for generalized algorithm.

◦ Fig. 10 10. Particle shape of Toyoura sand under microscope (40x) from Yin et al. [25].

2

data simulation

1.5

Stress Ratio, τ/σ

1

0.5

0 -2

-1.5

-1

-0.5

0

0.5

1

1.5

1

1.5

2

-0.5

-1

-1.5

Axial Strain (%) 3

data simulation

2.5

εv (%)

2

1.5

1

0.5

0 -2

-1.5

-1

-0.5

0

0.5

Axial Strain (%)

Fig. 11 Drained Drained cyclic test on loose Toyoura sand ((e0= 0.845).

2

2

data simulation

1.5

Stress Ratio, τ/σ

1

0.5

0 -2

-1.5

-1

-0.5

0

0.5

1

1.5

2

1

1.5

2

-0.5

-1

-1.5

Axial Strain (%) -1

data

-0.8

simulation -0.6 -0.4

εv (%)

-0.2 -2

-1.5

-1

-0.5

0

0.5

0 0.2 0.4 0.6 0.8 1

Axial Strain (%)

Fig. 12 Drained cyclic test on loose Toyoura sand (e0= 0.653).

1 ξ=0° data(data)

0.8

ξ=0° (simulation) simulation

0.6

Stress Ratio, τ/σ

0.4 0.2 0 -0.8

-0.4

0

0.4

0.8

-0.2 -0.4 -0.6 -0.8

(a)

-1

Axial Strain (%) -0.6 ξ=0° data(data) ξ=0° (simulation) simulation

Volumetric Strain, (%)

-0.4

-0.2

-0.9

-0.6

-0.3

0

0.3

0.6

0.9

0

0.2

(b) 0.4

Axial Strain (%)

Fig. 13 a) Stress–strain response during cyclic loading and b) evolution of volumetric strain.

Fig.. 14. Particles of Toyoura sand under microscope (40x) from Qin et al. [33].

250

Simulation

τ (kPa)

200

150

100

Experiment 50

0 0

0.4

0.8

1.2

1.6

2

γ= ε1- ε3 (%) 250

Simulation

200

τ (kPa)

150

100

Experiment

50

0 0

50

100

150

200

σ (kPa)

Fig. 15. Experimental results and model simulations at e=0.82.

100

e0=0.75

80

Loose element at compact shakedown state

e0=1.05

Shear stress, τ, kPa

60 40 20 0 0

1

2

3

4

5

6

7

-20 -40 -60 -80

Medium dense element at compact shakedown state

(a)

-100

Shear strain, ε, ×10-3 7

(b) Axial strain, ε, ×10-3

6 5 4 3 2

loose sample, void ratio= 1.05 1

Medium dense sample, void ratio= 0.75 0 0

500

1000

1500

2000

2500

3000

Number of cycles Fig. 16. a) Simulation results for complete stress-strain response and b) axial strain versus number of cycles for Toyoura sand.

Table 1 Role and physical meaning of the model parameters with typical ranges Parameters

Role and/or physical meaning

E ν μ

Young modulus (MN/m2) Poisson's ratio Slope of normal projection of critical state line in the q–p plane under triaxial tests Function of state Vertical distance between the compact void ratio line and the ultimate state line Void ratio at critical state for p=0 (kPa) Void ratios at the baseline for p=0 (kPa) Void ratios at the baseline for p=0 (kPa) at compact state Angle of dilatation of the medium (degree) Friction angle (degree) Compact state angle (degree) Critical state angle (degree) Slope of Casagrande line Vertical distance between ebase and ecritical in the e−log p graph Coefficient of the interpolation function

a0 c eh ehbase ehcompact δ φ φcompact φcritical λ b n

Typical ranges 40−150 0.1–0.3 0.8–1.9 0.1–0.9 0.4–0.6 0.7–1.2 0.9–1.2 0.3–0.6 0.0-45 35–55 38–68 25–35 0.01–0.15 0.1–0.5 0.01–0.05

Table 2 Model parameters for soils used in the present paper Parameters Soil’s name Compaction E ν μ a0 c eh ehbase ehcompact δ φ φcompact φcritical λ b n

Yin et al. [25] Toyoura sand 70 0.25 1.38 0.5 0.485 0.934 0.997 0.512 8 38 45 28 0.019 0.2 0.020

Sazzad and Suzuki [27] Oval particles dense sample 98 0.21 1.45 0.4 0.226 0.895 0.721 0.495 12 42 48.8 29.5 0.04 0.2 0.035

30

Dafalias et al. [24] Toyoura sand 80 0.25 1.25 0.1 0.424 0.834 0.875 0.451 15 44.3 49.6 31 0.02 0.2 -