Constitutive response of idealized granular media under the principal stress axes rotation

International Journal of Plasticity 24 (2008) 1967–1989

Contents lists available at ScienceDirect

International Journal of Plasticity journal homepage: www.elsevier.com/locate/ijplas

Review

Constitutive response of idealized granular media under the principal stress axes rotation Seiichiro Tsutsumi a,*, Kenji Kaneko b a

Department of Marine System Engineering, Kyushu University, Motooka 744 (W3-616), Nishi-ku, Fukuoka 819-0395, Japan Department of Environmental and Civil Engineering, Hachinohe Institute of Technology Obiraki 88-1, Myo, Hachinohe 031-8501, Japan

b

a r t i c l e

i n f o

Article history: Received 29 February 2008 Received in final revised form 9 May 2008 Available online 7 July 2008

Keywords: Granular element method Non-coaxiality Two-scale analysis B. Anisotropic material B. Elastoplastic material

a b s t r a c t This paper is dedicated to the understanding of the phenomena, which give rise to anisotropy and non-coaxiality in granular materials. In achieving three-dimensional numerical simulation under static condition of granular media, granular element method (GEM) is adopted in this study. The method has been incorporated into the so-called mathematical homogenization theory for quasistatic equilibrium problems, which enables us to obtain the macroscopic/phenomenological inelastic deformation response of a representative volume element (RVE). To examine the anisotropic macroscopic deformation properties of the assumed RVE, which is solved by granular element method (GEM), a series of numerical experiments involving the pure rotation of the principal stress axes are carried out, and its results are discussed in relation to induced anisotropy and non-coaxiality. Ó 2008 Elsevier Ltd. All rights reserved.

Contents 1. 2.

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Two-scale numerical method for granular media based on mathematical homogenization. . . . . . . . 2.1. Two-scale modeling of granular media . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.1. Variational formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.2. Two-scale boundary value problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.3. Strong form and general numerical algorithm of the two-scale boundary value problem 2.2. 3D granular element method for micro-scale problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.1. Constitutive relation for particle contacts. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

* Corresponding author. Tel./fax: +81 92 802 3464. E-mail address: [email protected] (S. Tsutsumi). 0749-6419/$ - see front matter Ó 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijplas.2008.05.001

2062 2063 2063 2063 2065 2066 2067 2067

1968

3.

4.

5.

S. Tsutsumi, K. Kaneko / International Journal of Plasticity 24 (2008) 1967–1989

2.2.2. Equilibrium equation of a particle. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.3. Treatment of periodic boundary condition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.4. Definition of the macroscopic stress . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.5. Numerical algorithm for micro-scale simulation by GEM . . . . . . . . . . . . . . . . . . . . . . On the description of anisotropic response . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1. Conventional plasticity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2. Extension of conventional plasticity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Mechanical responses . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1. Granular media and stress paths . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2. Procedures of principal stress axes rotation test. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3. Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Concluding remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2068 2068 2071 2072 2072 2073 2073 2074 2074 2076 2077 2082 2082 2082

1. Introduction The mechanical behavior of granular materials under non-proportional loading condition is a problem of great interest, since the stress path often deviates severely from proportional loading in many real situations, e.g. tidal waves, earthquakes, wheel rotations, footing penetration and plastic instability phenomena inducing shear bands and/or diffuse modes. Also, natural geomaterials are often deposited in horizontal layers and then subjected to anisotropic stresses leading to preferred orientation of the particles. As a consequence, most natural sand deposits possess some anisotropic structures, which cause variation in deformation strength characteristics as the loading direction changes. The following results have been found experimentally for proportional and non-proportional loading behavior of granular materials: (i) The largest stiffness and strength occur for loading perpendicular to the bedding plane and the lowest values are observed along the bedding plane (cf. Oda, 1972, 1993; Matsuoka and Nakai, 1974; Arthur et al., 1980, 1981). (ii) Not only the magnitude but also the direction of inelastic stretching is dependent on the direction of the stress increment, as revealed by stress probe tests (Ishihara and Towhata, 1983; Miura et al., 1986a,b; Pradel et al., 1990; Gutierrez et al., 1991). (iii) Inelastic deformation is induced by the rotation of the principal stress axes, even if all of the values of the principal stresses and invariants of stress are kept constant (Broms and Casbarian, 1965; Arthur et al., 1980, 1981; Ishihara and Towhata, 1983; Symes et al., 1984; Miura et al., 1986a,b; Gutierrez et al., 1991; Joer et al., 1998). These studies show that the deformation strength characteristics highly depend on loading condition and anisotropic structures of materials. However, the experimental tests have difficulties to derive absolute conclusions on behave of the restrictions of the test system and the variation of specimen. In an attempt to overcome these experimental problems, inelastic mechanical responses of idealized granular media have been investigated numerically, since any intermediate data can be used for different tests and any stress or strain control can be applied flexibly (Bardet, 1994; Kishino, 1989, 2002; Wren and Borja, 1997; Kuhn, 1999; Kishino et al., 2001; Kaneko, 2001; Kaneko et al., 2003; Tsutsumi et al., 2005; Nicot and Darve, 2007a). From these micromechanical investigations for the ideal granular assemblies, it has been revealed that a regular plastic flow rule (independence of the inelastic stretching direction with respect to the tangential stress rate) exists only in particular (axisymmetric) loading conditions. For more general loading conditions, the inelastic stretching direction depends on the stress rate component tangential to the yield surface. Among them, Bardet (1994) have conducted a numerical simulation on the stress probe test by using the distinct element model (DEM), originally proposed by Cundall and Strack (1979), to examine the macroscopic elastoplastic mechanical behavior of a granular assembly. It was shown that DEM might be suited mainly to dynamic problems, not to static, since it is difficult to control arbitrary some stress components under dynamic motions. On the other hands, Kishino (1989) developed the granular element method (GEM), which enables us to obtain the elastoplastic mechanical response of granular media under quasi-static equilibrium

S. Tsutsumi, K. Kaneko / International Journal of Plasticity 24 (2008) 1967–1989

1969

states. Furthermore, a method of global-local (two-scale) analysis method has been developed for quasi-static equilibrium problems by incorporating GEM into the mathematical homogenization theory (Kaneko, 2001; Kaneko et al., 2003). The mathematical theory of homogenization for heterogeneous media with periodic microstructures enables us to realize the two-scale modeling, which consistently encompasses both micro- and macro-scales together with variational statements (Terada and Kikuchi, 2001). Due to this consistency, the nonlinear mechanical behavior is easily incorporated into the numerical analysis. The two-scale numerical method combines these macro- and micro-scale problems to solve them simultaneously and it enables structural analysis that faithfully reflects the deformation characteristics of a microstructure allocated to the micro-scales. In the two-scale method of granular medium, the macro-scale problem solved by the finite element method (FEM), whereas the micro-scale problem by GEM (Kaneko, 2001; Kishino et al., 2001; Kishino, 2002; Kaneko, 2001; Kaneko et al., 2003; Tsutsumi et al., 2005). The role of the constitutive equation in FEM is performed by micro-scale analysis incorporating the GEM with periodic boundary control. In other words, according to the mathematical homogenization theory, micro-scale analysis provides a macroscopic constitutive relationship. In this study, paying attention to the above mentioned feature of the two-scale method based on mathematical homogenization, the three-dimensional GEM is adopted in order to examine their macroscopic mechanical response to pure rotation of the principal stress axes, in which all of the values of the principal stresses and invariants of stress are kept constant (Symes et al., 1984; Miura et al., 1986a,b; Gutierrez et al., 1991). A series of numerical experiments involving the pure rotation of the principal stress axes are carried out, and its results are discussed in relation to induced anisotropy and non-coaxiality. Throughout this paper, the mechanical responses under static loading conditions are considered. 2. Two-scale numerical method for granular media based on mathematical homogenization This section gives a formulation of three-dimensional two-scale modeling for granular media based on mathematical homogenization theory (Lions, 1981; Allaire, 1992) in order to explain that the micro-scale problem plays a role of a constitutive law for the macro-scale problem. The formulation is nearly the same as for two dimensions (Kaneko et al., 2003). 2.1. Two-scale modeling of granular media 2.1.1. Variational formulation Let us consider the quasi-static boundary value problem of a granular body shown in Fig. 1. The body is regarded as an assembly of periodic arrangements of micro-structural elements (RVE), which are composed of a random distribution of elastic spheres and voids. We assume that the size of the RVE is small enough to the overall structure and is represented by a normalized parameter e.

Fig. 1. A body composed of grains.

1970

S. Tsutsumi, K. Kaneko / International Journal of Plasticity 24 (2008) 1967–1989

Let Xe be an open domain of the granular body with smooth boundary oXe, and divided into two parts as follows:

Xe ¼ XeC [ C e

ð1Þ

where Ce is the totality of the internal surfaces of spheres associated with contact and friction. We also define partial open domain XeC by excluding Ce from Xe. Then, this problem is fully described by the equilibrium problem in XeC with oXe and friction and contact conditions on Ce. The equilibrium equation for the stress re(x) in XeC and boundary conditions on ouXe  oXe and o Xe  oXe are respectively given by t

e

div re ðxÞ þ b ðxÞ ¼ 0 ue ðxÞ ¼ 0 on ou Xe and re  n ¼ t on ot Xe

ð2Þ ð3Þ

where ue(x) is the displacement vector, be (x) is the body force, t is the traction vector specified on otXe with the outward unit normal vector n. Here and in the subsequent sections, we indicate the dependency on the microscopic heterogeneities by a superscript e on each variable. We also assume that individual particles and voids reveal elasticity along with the following constitutive equation:

re ðxÞ ¼ De ðxÞ : ee ðxÞ

ð4Þ

where ee(x) is an infinitesimal strain and De(x) the elasticity tensor, which is symmetric and positive definite. The strain is related to the displacement by the following relationship as usual:

ee ðxÞ ¼ rS ue ðxÞ

ð5Þ

S

where r is a gradient operator which produces a symmetric second-order tensor. In order to provide the complete set of governing equations, we define the friction and contact conditions on Ce. The outward unit normal nC and the unit tangential vectors t Cu and t Ch are defined so that (nC, t Cu , t Ch ) can be a set of base vectors for the right-hand local coordinate system (see Fig. 1c). Then, the displacement vector ue and the stress vector Te on Ce are decomposed respectively into their normal and tangential components as follows:

n oT ue ¼ uCn ; uCu ; uCh ; n oT T e ¼ T Cn ; T Cu ; T Ch ;



 uCn ¼ ue  nC ; uCu ¼ ue  t Cu ; uCh ¼ ue  t Ch   T Cn ¼ T e  nC ; T Cu ¼ T e  t Cu ; T Ch ¼ T e  t Ch

ð6Þ ð7Þ

Using these components, we have the contact condition of Kuhn–Tucker form as

T en P 0; ½½uen  P 0; T en ½½uen  ¼ 0 on C e ;

ð8Þ

where ½½uen  represents the relative normal displacement. The Coulomb’s friction law is introduced as

lT en P

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi T 2u þ T 2h on C e

ð9Þ

where l is the friction coefficient. The quasi-static boundary value problem for a granular material, whose local structure is composed of discrete spheres, is completely represented by the Eqs. (2)–(9). Then, the problem is governed by the following variational inequality (Kikuchi and Oden, 1988):

aðue ; we  ue Þ þ jðue ; we Þ  jðue ; ue Þ P lðwe  ue Þ; Z rwe : D : rue dx aðue ; we Þ ¼ lðwe Þ ¼

XeC

Z ot Xe

jðue ; we Þ ¼

t  we ds þ

Z Ce

Z ot Xe

e

b  we dx

ljT en ðue Þjj½½wet j ds;

8we ;

ð10Þ

S. Tsutsumi, K. Kaneko / International Journal of Plasticity 24 (2008) 1967–1989

1971

where we is the admissible displacement vector and ½½wet  is the tangential component of the relative displacement. 2.1.2. Two-scale boundary value problem We introduce two distinct scales; micro- and macro- scales, x and y, the latter of which is related to the former as y = x/e. Due to the introduction of these spatial scales, the domain XeC is divided into X measured by the macroscopic variable x and YC = Y/C measured by the microscopic one y as follows (see Fig. 2):

XeC ¼ X  Y C ¼ fðx; yÞjx 2 X  R3 ; y ¼ x=e 2 Y C  R3 g;

ð11Þ

where Y is the microscopic domain and C is the contact region. With this space decomposition, all the field variables with superscript e are re-defined as functions of two scale variables, x and y, e.g., re (x) = re (x,y). These variables are periodic with respect to y, i.e. Y-periodic. It has been demonstrated in Terada and Kikuchi (2001) that the theory of two-scale convergence of Allaire (1992) can be utilized in the derivation of the two-scale boundary value problem for a heterogeneous solid. In this particular situation, the similar argument would hold for the variational inequality (10) and the following formula can be obtained as a limit of the appropriate convergence study:

Z

rx ðw0  u0 Þ : hDi : rx ðu0 Þ dx þ

X

Z

Z

hry ðw1  u1 Þ : Di : rx ðu0 Þ dx þ

X

Z

rx ðw0  u0 Þ

X

: hD : ry u1 i dx þ hry ðw1  u1 Þ : D X Z Z 1 : ry u i dx þ ljT n ðu1 Þjj½½w1t j ds  ljT n ðu1 Þjj½½u1t j ds C C Z Z 0 0 0  t  ðw  u Þ ds þ hbi  ðw  u0 Þdx ot X

P 0;

ou X

8w0 and 8w1 ;

Fig. 2. Decomposition to micro- and macro-space.

ð12Þ

1972

S. Tsutsumi, K. Kaneko / International Journal of Plasticity 24 (2008) 1967–1989

where hi indicates the volume average over the unit cell domain Y, and rx and ry indicate the gradient with respect to the macro- and micro-scales, respectively. Here, u0(x) and w0(x) are macroscopic displacement and its variation and are independent of the microscopic heterogeneities. On the other hand, the microscopic displacement (trial function) u1(x, y) and its variation (test function) w1(x, y) are Y-periodic. In inequality (12), the test function w1 can be chosen as w1 = u1, while w0 can be chosen as 0 w = u0 + av0 where v0 is an arbitrary function and a is an arbitrary real number. Then, inequality (12) yields the following equality:

Z

rx w0 : hD : ðrx u0 þ ry u1 Þi dx þ

X

Z

hbi  w0 dx þ

ou X

Z

t  w0 ds;

8w0 ;

ð13Þ

ot X

On the other hand, in inequality (12), w0 can be set to w0 = u0 and w1 to w1 = u1 + av1. Here, v1 is an arbitrary function. Then the variational inequality for u1(x, y) in a unit cell (RVE) yields

Z

ry w1 : D : ry u1 dx þ

YC

Z P

Z

C





ljT n ðu1 Þj½½w1t  ds  

ry w1 : D dy : rx u0 ;

Z

C





ljT n ðu1 Þj½½u1t  ds

8w1

ð14Þ

X

where the macroscopic deformation rxu0 plays the role of constant excitation. The average behavior of the overall structure can be characterized by the macro-scale virtual work Eq. (13) under the influence of the microstructure whose mechanical behavior is characterized by the micro-scale variational inequality (14). These micro- and macroscopic simultaneous equations govern the two-scale boundary value problem for a granular material. 2.1.3. Strong form and general numerical algorithm of the two-scale boundary value problem The strong forms of the microscopic problem (14) can be identified with

div r0 ¼ 0 in Y c ½½u1n 

ð15Þ T n ½½u1n 

P 0;  T n P 0; qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2  lT n P T u þ T 2h on C

¼ 0 on C

ð16Þ ð17Þ

along with the constitutive equation in a micro-scale

r0 ðx; yÞ ¼ D : ðrx u0 þ ry u1 Þ ¼ D : ðE þ ry u1 Þ;

ð18Þ

where we have defined the macroscopic (average) strain E as

E ¼ rsx u0

ð19Þ

The actual displacement u(x, y) is expressed with the macroscopic (average) displacement vector u0 and the Y-periodic microscopic one u1 as

uðx; yÞ ¼ u0 ðxÞ þ u1 ðx; yÞ ¼ E  y þ u1 ðx; yÞ

ð20Þ

On the other hand, Eq. (13) is equivalent to the local form of the macroscopic boundary value problem as follows:

div R þ B ¼ 0 in X 0

u ðxÞ ¼ 0 on ou X and R  n ¼ t on ot X

ð21Þ ð22Þ

Here, R and B are the macroscopic stress and the body force, respectively, and are given by the volume averages of the corresponding microscopic variables over a unit cell as

RðxÞ ¼ hr0 ðx; yÞi

ð23Þ

BðxÞ ¼ hbðx; yÞi

ð24Þ

In our two-scale modeling, the boundary value problem is governed by the equations in both scales. In particular, the Y-periodic microscopic displacement u1 is a solution under the constraint conditions

S. Tsutsumi, K. Kaneko / International Journal of Plasticity 24 (2008) 1967–1989

1973

associated with contact and friction. Thus, we note that the governing equations in both scales are completely coupled, but can be solved independently by a sequential solution scheme in each loading step. A general numerical algorithm used by FEM to solve such two-scale nonlinear boundary value problems is proposed in Terada and Kikuchi (2001). In the algorithm, both the macroscopic virtual work equation and the microscopic one can be solved by the continuum-based FEM regardless of the mechanical nature of the microstructure. However, it is generally a difficult task to solve (14) by using the FEM in the case of granular type media because of their discrete nature with friction and contact. Therefore, we identify the physical model assumed in the above formulation, namely, the microstructure composed of elastic particles and voids, with that composed of rigid spheres, spring elements and frictional devices. This type of physical models is widely used in the single-scale analysis of granular media; e.g., the distinct element method (DEM; Cundall and Strack, 1979). In our numerical method, we employ the granular element method (GEM; Kishino, 1989) for the micro-scale problem and work out the detailed formulation in next section. GEM solves the quasi-static equilibrium equation unlike DEM. Therefore, it is equivalent to solving microscopic quasi-static equilibrium Eq. (15) with the constraint conditions (16) and (17). In our method, the microscopic problem is solved by GEM, while the macroscopic one is by a continuum-based finite element method. We show the concept of two-scale numerical method for each loading step in Fig. 3. The micro-scale problem for RVE is given the macroscopic strain DE by the macro-scale problem, and returns the macroscopic stress R after solving the micro-scale problem using DE. Therefore, in our two-scale numerical method, we can say that the micro-scale problem plays a role of a constitutive equation for the macroscopic boundary value problem. In this paper, we pay attention to this feature of the two-scale method based on mathematical homogenization, and use three-dimensional micro-scale analysis by GEM to estimate the macroscopic constitutive response under rotation of principal stress axis. 2.2. 3D granular element method for micro-scale problems We review the formulation of the granular element method with reference to Kishino (1989), which would provide the equivalent solution for microscopic problem. After providing the constitutive equation with constraint conditions and the equilibrium equation of a single particle, we drive the equilibrium equation with the periodic boundary conditions for unit cells (granular assembly). For our notational convenience to describe the GEM, all the tensor quantities are represented in vector or matrix forms expressed by {} or []. 2.2.1. Constitutive relation for particle contacts We first introduce the constitutive equation and the constraint conditions for particle contacts. Let rigid circular spheres i and j be in contact with each other on contact point Cij. In this section, we

Macro-scale problem of equivalent overall structure Eq. (13); FEM macroscopic strain E

macroscopic stress Σ

localization

homogenization Micro-scale problem of microstructure (RVE) Eq. (14); GEM

relative displacement [[uij ]]

contact force Tij

local constitutive law and friction-contact conditions Fig. 3. Concepts of global–local computations.

1974

S. Tsutsumi, K. Kaneko / International Journal of Plasticity 24 (2008) 1967–1989

express the outward unit normal nc = nij and the unit tangential vectors t Cu ¼ t iju and t Ch ¼ t ijh . The unit normal vector of particle i toward particle j on Cij is given by

nij ¼

yj  yi ; jyj  yi j

ð25Þ

where yi and yj are the position vectors of the center points of these particles. The kinetic variable on Cij is the contact force fT ij g ¼ fT ijn ; T iju ; T ijh gT , whereas the kinematic variable associated with {Tij} in the constitutive equation is the relative displacement f½½uij g ¼ f½½uijn ; ½½uiju ; ½½uijh gT between two particles. In the GEM, the elastic characteristics of a continuum are replaced by those of the spring devices at contact points. We employ the normal and tangential spring devices on contact point Cij as shown in Fig. 3(a). Then, we have the following simple linear relationship between {Tij} and {[[uij]]}:

fT ij g ¼ ½Sf½½uij g;

ð26Þ

where [S] is the elastic modulus matrix given as

2

sn

6 ½S ¼ 4 0 0

0

0

st 0

3

7 0 5: st

Here, sn and st are the elastic constants in the normal and tangential directions, respectively. In order to express the slip between particles, we set up the frictional contact device as illustrated in Fig. 3a. Then the contact condition on Cij is expressed as

½½uijn  ¼ 0:

ð27Þ

When this condition does not hold, we set zero to the contact force vector. Also, the friction condition between particles is written as

T ijn tan / P

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ij ij2 T ij2 u þ T h on C ;

ð28Þ

where / is the friction angle on Cij. If the tangential component of the contact force violates this condition, we correct it compulsively to the limit value. The contact behavior in GEM is summarized in Fig. 3b. 2.2.2. Equilibrium equation of a particle Let us consider the quasi-static equilibrium of a single particle i whose radii is given by ri. The translational equilibrium for particle i is given by a X fF ij g ¼ f0g;

ð29Þ

j¼1

where {F ij} is the translational force exerted on i by any particle j, which is in contact with i on a contact point C ij, and a is the total number of particles in contact with i. The rotational equilibrium of particle i takes the form a X fM ij g ¼ f0g;

ð30Þ

j¼1

where {Mij} is the moment of a contact force exerted by j. The equilibrium Eqs. (29) and (30) can be combined as a n o X ff ij g ¼ f0g; ff ij g ¼ F ij1 ; F ij2 ; F ij3 ; M ij23 =ri ; M ij31 =ri ; M ij12 =ri ;

ð31Þ

j¼1

where {fij} is the generalized force vector for a single contact point Cij. The generalized force vector {fij} defined in Cartesian axis y1, y2, y3, are related to the contact force vector {Tij} through the coordinate transformation matrix [Rij] via

S. Tsutsumi, K. Kaneko / International Journal of Plasticity 24 (2008) 1967–1989

ff ij g ¼ ½Rij fT ij g:

1975

ð32Þ

Here, the components of [Rij] are given by

2

sin u cos h cos u cos h

6 sin u sin h 6 6 6 cos u ij ½R  ¼ 6 6 0 6 6 4 0

cos u sin h  sin u  sin h cos h

0

3

 sin h

7 7 7 7 0 7:  cos u cos h 7 7 7  cos u sin h 5 cos h

sin u

0

Using Eq. (32), we can rewrite the equilibrium Eq. (31) for particle i as follows: a X ð½Rij fT ij gÞ ¼ f0g:

ð33Þ

j¼1

Since the motion of each particle is constrained from friction and contact on contact point Cij, the microscopic analysis is always nonlinear. Thus, we replace the quasi-static equilibrium Eq. (33) by the incremental form as a X ð½Rij fDT ij gÞ ¼ f0g:

ð34Þ

j¼1

The kinematic variables of the particle motion are the translational displacements, ui1 , ui2 and ui3 , and the rotation xi23 , xi31 and xi12 of particle i. Putting them into a generalized element displacement vector such that fui g ¼ fui1 ; ui2 ; ui3 ; ri xi23 ; ri xi31 ; ri xi12 gT , we shall use their increments to be consistent with incremental form Eq. (34). Then, the relationship between the incremental relative displacement {[[Duij]]} on Cij, and the generalized element displacements {Dui} and {Duj} are given as follows:

^ ij T fDuj g; f½½Duij g ¼ ½Rij T fDui g þ ½R

ð35Þ

where

2

sin u cos h cos u cos h  sin h

6 sin u sin h 6 6 6 b ij  ¼ 6 cos u ½R 60 6 6 40 0

3

 cos h

7 7 7 7 0 7: cos u cos h 7 7 7 cos u sin h 5

0

 sin u

cos u sin h  sin u sin h

cos h

Finally, using the linear constitutive law of Eqs. (26) and (35), we can rewrite the incremental form of the equilibrium Eq. (34) for particle i as follows: a a X X ii ^ij fDuj gÞ ¼ f0g; ð½Rij fDT ij gÞ ¼ ½k fDui g  ð½k j¼1

ð36Þ

j¼1

where ii

½k  ¼

a X

^ij  ¼ ½Rij ½S½R ^ ij T : ð½Rij ½S½Rij T Þ; ½k

j¼1

2.2.3. Treatment of periodic boundary condition We here introduce the periodic boundary condition to the GEM formulation. Fig. 4 schematically shows the periodicity of a unit cell with respect to adjacent cells. Here, due to the periodicity, particle i in one of the adjacent cells is identified with particle i in this unit cell. We call this particle i the master particle whereas i the slave particle. The microscopic displacement vectors of these particles are respectively expressed as follows:

1976

S. Tsutsumi, K. Kaneko / International Journal of Plasticity 24 (2008) 1967–1989

Fig. 4. Illustration of the constitutive relation for particle contacts.

fui g ¼ ½Efyi g þ fui1 g; 



ð37Þ



fui g ¼ ½Efyi g þ fui1 g:

ð38Þ i1

Here, [E]: = E is the macroscopic (average) strain matrix, and {ui1} and fu g are the Y-periodic displacement vectors of i and i, respectively. Then the periodic boundary condition is expressed as 

fui1 g ¼ fui1 g:

ð39Þ

Subtracting Eq. (37) from Eq. (38) with Eq. (39), we have 

fui1 g ¼ fui g þ ½Eðfyi g  fyi gÞ:

ð40Þ

 slave particles, the last term in the leftIf master particle i is in contact with a particles containing a hand side of the equilibrium Eq. (36) for master particle i is divided into the parts about the master particles and slave ones and Eq. (36) yields ii

½k fDui g 

aa a X X ^ij fDuj gÞ  ^ij fDuj gÞ ¼ f0g; ð½k ð½k

ð41Þ

 j¼1þaa

j¼1

in which the last term in the left-hand side is due to the periodicity constraint for particle i with respect to slave particles. Here, [kii] is rewritten as following expression: ii

½k  ¼

aa a X X   ð½Rij ½S½Rij T Þ þ ð½Rij ½S½Rij T Þ:  j¼1þaa

j¼1

Then, by substituting the incremental forms of Eq. (40) into the last term in the left-hand side in Eq. (41) with some algebraic manipulations, we arrive at the following equilibrium equation for a single particle i with the periodicity constraint associated with slave particles: ii

fDHi g ¼ ½k fDui g 

aa a X X ^ij fDuj gÞ  ^ij fDuj gÞ ¼ ½K i fDuij g; ð½k ð½k

ð42Þ

 j¼1þaa

j¼1

where the components of the element stiffness matrix [Ki] and the displacement vector {Duij} is given by ii ^i1      ½k ^ij      ½k ^ia ; ½K i  ¼ ½½k   ½k n oT fDuij g ¼ fDui gT ; fDu1 gT ; . . . ; fDuj gT ; . . . ; fDua gT :

And, here, we have defined the vector of external forces as

fDHi g ¼ ½DE

a X  j¼1þaa





^ij ðfyj g  fyj gÞ ¼ ½DE ½k

a X



jj

^ij fl gÞ; ð½k

 j¼1þaa

ð43Þ

S. Tsutsumi, K. Kaneko / International Journal of Plasticity 24 (2008) 1967–1989

1977

where jj



fl g ¼ fyj g  fyj g is the characteristic length vector, which provides the length of RVE. Finally, if there are n particles in the unit cell, we define the overall external force vector and the overall generalized velocity vector for all the particles as

n oT fDHg ¼ fDH1 gT ; . . . ; fDHi gT ; . . . ; fDHj gT ; . . . ; fDHn gT ; n oT fDug ¼ fDu1 gT ; . . . ; fDui gT ; . . . ; fDuj gT ; . . . ; fDun gT : Then, by piling up the equilibrium Eq. (41) according to the order of {Du} for all the particles within a unit cell except for the slave particles, which enjoy the periodicity constraints, we have the rate form of the equilibrium equation for the unit cell as

fDHg ¼ ½KfDug;

ð44Þ

which should be solved for the overall generalized velocity vector {Du}, where [K] is the overall stiffness matrix, which consists of the element stiffness [Ki] of all the particles with reference to the order of {Du}. The components of [K] is given by

2

11

^1i     ½k ^1j     ½k .. .. .. .. . . . . ^ii     ½k ^ij     ½k .. .. .. .. . . . . ^ji     ½kjj     ½k .. .. .. ... . . . ^ni     ½k ^nj     ½k

½k  6 . 6 . 6 . 6 ^i1  6 ½k 6 6 . ½K ¼ 6 6 .. 6 6 ½k ^j1  6 6 6 .. 4 . ^n1  ½k

^1n  3    ½k .. .. 7 7 . . 7 7 in 7 ^    ½k  7 .. .. 7 7: . . 7 7 ^jn  7    ½k 7 7 .. .. 7 . . 5 

nn

½k 

^ij  should be a zero-matrix if particles i and j do not conWe here note that the non-diagonal block ½k tact with each other. 2.2.4. Definition of the macroscopic stress We consider a unit cell, which is composed of rigid particles; see Fig. 5a. In the analysis by the GEM, the macroscopic (average) stress of this granular assembly is represented by the following relationship:

R ¼ hr0 i ¼

1 jYj

Z Y

r0 dy ¼

1 X X ij  T  y ij ; jYj i 

ð45Þ

j



where T ij is the contact force from an external particle to an internal particle on the external boundary  oYC of the unit cell, yij is the position vector of a contact point between the internal and the external particles, and jYj is the volume of the unit cell; see Fig. 5b. This definition of average stress is well known in the area of micromechanics for granular materials (see, e.g.; Oda and Iwashita, 1999). The macroscopic stress is always symmetric since all the particles within a unit cell satisfy the equilibrium condition for the moment. In this study, we simulate the deformation behavior of granular media under the pure rotation of the principal stress axes. Therefore, the micro-scale simulation in this study should be simulated by the stress control. We show the relationship between the incremental macroscopic stress and the incremental macroscopic strain as follows:

DR ¼ DG DE;

ð46Þ

1978

S. Tsutsumi, K. Kaneko / International Journal of Plasticity 24 (2008) 1967–1989

Fig. 5. Master and slave particles for periodic motion.

where

DG ¼

1 X X ijT ij ij y K l : jYj i  j

The incremental macroscopic strain for the given stress increment can be calculated by this equation. 2.2.5. Numerical algorithm for micro-scale simulation by GEM GEM enables us to control the stress accurately, because it solves the quasi-static equilibrium equation directly. The analysis flow of stress controlled GEM in each incremental loading step is presented as follows: (1) Calculate DR for given stress increment DR by Eq. (46) and update E. (2) With the value of DE, calculate {DH(DE)} by Eq. (43). (3) Solve the quasi-static stiffness equation of the granular assembly (44) for the increment of the overall generalized displacement vector {Du}. (4) With the value of {Du}, compute the contact force increment between each particle, and update all the coordinate transformation matrix [Rij] and all the contact force vectors {Tij}. (5) Check whether or not the contact condition of Eq. (27), the friction condition (28) are satisfied. (6) According to the results of (5), correct the contact force {Tij} and check the equilibrium condition (31) of all the particles and the stress convergence with corrected {Tij} and [Rij]. (7) If this condition is not satisfied, go to (1) using the residual stress and forces instead of DR and {DH(DE)}. If the condition is satisfied, output the macroscopic strain E and go to next loading step.

3. On the description of anisotropic response The phenomenological description of the anisotropic responses of geomaterials has been studied up to the present (e.g. Sekiguchi and Ohta, 1977; Hashiguchi, 1977, 2001; Hashiguchi and Chen, 1998; Nemat-Nasser and Zhang, 2002; Hashiguchi and Tsutsumi, 2001, 2003; Tsutsumi and Hashiguchi, 2005; Nicot and Darve, 2005; Dean, 2005; Zhu et al., 2006; Nicot and Darve, 2007a,b,c; Anandarajah, 2008). In this section the essential roles in the phenomenological description of the anisotropy and the non-coaxiality of granular media are reviewed briefly (Tsutsumi and Hashiguchi, 2005).

S. Tsutsumi, K. Kaneko / International Journal of Plasticity 24 (2008) 1967–1989

1979

Denoting the current configuration of the material particle by x and the current velocity by v, the velocity gradient is described as L = ov/ox in which the stretching and the continuum spin are defined as D (=(L + LT)/2) and W(=(L  LT)/2), respectively. 3.1. Conventional plasticity In traditional elastoplasticity of the rate-form the stretching D is additively decomposed into elastic stretching De and plastic stretching Dp, i.e.

D ¼ De þ Dp ;

ð47Þ

where De is given by 

De ¼ E1 r;

ð48Þ

where r is the Cauchy stress. Here, by avoiding the influence of rigid-body rotation on the constitutive relation, the following Jaumann rate is adopted: 



T T WT þ TW:

ð49Þ

Here, T denotes the arbitrary second-order tensor with the dimension of the stress, and ( ) denotes the material-time derivative. Hereafter, the rate form of Hooke’s law is adopted for the elastically isotropic materials. The plastic stretching Dp based on the consistency condition and the plastic potential is generally described as

Dp ¼

  1 of  og r ; p tr or or M

ð50Þ

where f and g are the yield and the plastic potential function, respectively. Mp is the plastic modulus. These are functions of stress and internal variables in general. The plastic constitutive Eq. (50) has the following properties: (1) The direction of plastic stretching Dp is given to the out-ward normal of the plastic potential function g (og/or) and independent of the stress rate. The magnitude of the plastic stretching is independent of the stress-rate component tangential to the yield surface, called the tangential-stress rate. (2) An isotropic function of the stress tensor obeying the orthogonal transformation f(QrQT) = Qf(r)QT (Q: orthogonal tensor) can be described by a function of principal stresses or invariants of stress. Then, if the yield function f is assumed to be given by an isotropic function of the stress tensor, plastic stretching is not induced during the principal stress axes rotation process with constant values of the principal stresses. (3) The coaxiality, i.e. the coincidence of the principal directions of plastic stretching with those of stress, holds as far as g is an isotropic function of the stress tensor. In fact, however, inelastic stretching is induced by the tangential stress-rate and its direction depends on the stress rate (cf. Broms and Casbarian, 1965; Poorooshasb et al., 1966; Lewin and Burland, 1970; Arthur et al., 1980, 1981; Ishihara and Towhata, 1983; Symes et al., 1984; Miura et al., 1986a,b; Pradel et al., 1990; Gutierrez et al., 1991; Joer et al., 1998; Tsutsumi and Hashiguchi, 2005). 3.2. Extension of conventional plasticity In order to describe aforementioned general non-proportional loading behavior of geomaterials, the conventional elastoplasticity model has been extended up to the present (Rudnicki and Rice, 1975; Yatomi et al., 1989; Hashiguchi and Tsutsumi, 2001, 2003, 2007; Tsutsumi and Hashiguchi, 2005). In this study, the formation of the unconventional model considering the so-called tangential-stress rate (or vertex) effect is viewed.

1980

S. Tsutsumi, K. Kaneko / International Journal of Plasticity 24 (2008) 1967–1989

Now, the stretching D be additively decomposed into the elastic stretching De and the inelastic stretching Di, i.e.

D ¼ De þ Di :

ð51Þ i

p

Further, let the inelastic stretching D be additively decomposed into plastic stretching D and tangential stretching Dt, i.e.

Di ¼ Dp þ Dt ;

ð52Þ p

t

on the premise that D and D are induced by the stress-rate component normal and tangential, respectively, to the yield/loading surface. Then, the inelastic stretching Di possesses the following properties: (a) Not only the magnitude but also the direction of the inelastic stretching Di exhibits a dependence on the stress rate. (b) The non-coaxiality between inelastic stretching Di and stress r is caused by the existence of Dt, even if an isotropic plastic potential surface of stress is adopted. (c) The inelastic stretching induced during the rotation of the principal stress axes can be described by the incorporation of the tangent effect, even if an isotropic yield function of stress is used. Here, consider the constitutive model with the anisotropic plastic potential surface, which exhibits non-coaxiality and in which the direction of the inelastic stretching depends only on the state variable but is independent of the rate variable. One should note the following facts: (i) Plasticity model cannot describe the dependence of the inelastic stretching on the tangential stress-rate as long as the tangent effect is not incorporated. That is to say, the non-coaxial model without the tangent effect is incapable of predicting appropriately the stress probe behavior. (ii) Further, the model cannot describe the inelastic stretching during the rotation of the principal stress axes as long as an anisotropic yield surface and/or the tangent effect is not incorporated. Inversely, inelastic stretching during the rotation of the principal stress axes is predicted even by the coaxial model with an isotropic plastic potential surface, only if an anisotropic yield surface is adopted. The anisotropic and non-coaxial response of geomaterials would be predicted by the tangential effect and/or the anisotropy in the yield condition. The anisotropy of soils can be described concisely by the concept of the rotation of the yield surface around the origin of the stress space, called the rotational hardening by Hashiguchi (1977), which was originally proposed by Sekiguchi and Ohta (1977). On the other hand, the anisotropy of metals is described by the translation and the transformation of the yield surface (e.g. Barlat et al., 2003a,b, 2007; Han et al., 2003; Kyriakides et al., 2004).

4. Mechanical responses 4.1. Granular media and stress paths Fig. 6 The body considered is a periodic assembly of micro-structural elements (RVE), in which the size of the RVE is small enough relative to the overall structure. In the two-scale method, the microscopic problem is treated numerically by means of GEM. The model under consideration, depicted in Fig. 7, has randomly created 1073 spheres of between 0.03 and 0.07 mm in diameter. The normal and tangential spring constants for the granular interaction are assumed to 100 and 70 (kN/m), respectively, while the friction angle is set to tan15°. After achieving the initial equilibrium state under confining pressure up to p0 = 100 kPa, three different stress paths are applied under the stress-controlled condition. Fig. 8 shows the initial equilibrium (isotropic) stress state p0 and the following stress paths A, B and C considered in this study. In the stress path A,

S. Tsutsumi, K. Kaneko / International Journal of Plasticity 24 (2008) 1967–1989

Fig. 6. Unit cell composed of particle.

Fig. 7. Granular element model randomly created by 1073 spheres of between 0.03 and 0.07 mm in diameter.

|| σ*|| / p proportional loading paths A and B 0.8

0.6

A-3

0.4

A-2, B-2, C -2

0.2

A-1, B-1

path C 0

pC

p0

Fig. 8. The stress paths and rotational path points.

p

1981

1982

S. Tsutsumi, K. Kaneko / International Journal of Plasticity 24 (2008) 1967–1989

the stresses are controlled from the isotropic stress state p0 so that the triaxial condition drxx = dryy = drzz/2 can be fulfilled under constant mean stress p0. The stresses are controlled in the stress path B so as to fulfill the condition drxx = 0, dryy = drzz under constant mean stress p0. On the other hands, in the stress path C, the stresses are first unloaded from the stress state p0 to the mean stress pc, and then triaxial compression is applied under the condition drxx = dryy = 0 and drzz 6¼ 0. 4.2. Procedures of principal stress axes rotation test In the experiments using the follow cylindrical specimen the four stress components, i.e. the axial stress rzz, the radial stress rxx, the peripheral stress ryy and the torsional shear stress ryz, can be applied independently under the coordinate system (x, y, z). These four stress components are described by the effective pressure p, the magnitude k k of the deviatoric stress r*, the Lode’s angle hr and the rotation angle a of the principal stress axes as

rffiffiffi 2  kr k sin hr ; 3 rffiffiffi    ryy 2  sin hr cos hr cos 2a pffiffiffi kr k ; ¼ p þ 3 3 rzz 3 rffiffiffi 2 cos hr sin 2a pffiffiffi ; ryz ¼ kr k 3 3

rxx ¼ p þ

2 3

ð53Þ ð54Þ ð55Þ

where

pffiffiffi tr r3 1 1 hr  sin 6 3 kr k3

!

2r2  r1  r3 ¼ tan1 pffiffiffi : 3ðr1  r3 Þ

ð56Þ

Eqs. (53)–(56) are derived by orthogonal transformation of the principal stresses to a coordinate system rotated around the intermediate principal stress axis (rxx = r2). The b-value is related to Lode’s angle hr as follows:

b

 r2  r3 1 pffiffiffi 3 tan hr þ 1 : ¼ r1  r3 2

ð57Þ

In the tests adopted for later examination, however, the stresses are controlled to fulfill rxx = r2 = (r 1 + r3)/2 = const., leading to b = 0.5, i.e. hr = 0 °. Then, Mohr’s stress circle for the pure principal stress axes rotation is depicted as shown in Fig. 9a. The stresses rxx, ryy, rzz and ryz are described by the principal stresses and a as follows:

rxx ¼ r2 ¼ p;  ryy r1 þ r3 r1  r3 ¼ cos 2a; 2 2 rzz r  r3 ryz ¼ 1 sin 2a: 2

ð58Þ ð59Þ ð60Þ

The state of stress induced in the tests can also be represented in the (X, Y) stress plane (Fig. 9b), where

X¼

ryy  rzz 2

;

Y ¼ ryz :

ð61Þ

In the (X, Y) plane, the length of the stress vector is equal to the radius of Mohr’s stress circle and makes twice the angle a, i.e.

tan 2a ¼

Y : X

ð62Þ

Strain increments induced by the imposed stress increments in the (X, Y) plane are also depicted taking (ezz  eyy)/2 and ezy for the X- and the Y-axis, respectively.

S. Tsutsumi, K. Kaneko / International Journal of Plasticity 24 (2008) 1967–1989

1983

τ σ yz 2α

σ3

σ zz

σ xx σ2

σ yy 2α

σ1

σ

0

σ zy

B-2b

Y

Pure principal stress rotation path

σ

B-2c

2α

A-2, B -2, C -2

0

X B-2a

B-2d Fig. 9. (a) Mohr’s circle of stress and (b) the principal stress rotation in (X, Y) plane.

4.3. Numerical results To examine the mechanical responses to the non-proportional loading conditions, a series of tests is conducted from the particular loading points shown in Fig. 9. In the non-proportional loading calculation with a pure rotation of the principal stress axes, the magnitude of the deviatoric stress kr*k, the mean stress p and the b-value are all kept constant. Fig. 10 shows the stress-strain relationships of the proportional loading path A. From the loading points A-1, A-2 and A-3 shown in Figs. 9 and 10, a series of tests on a continuous principal stress axes rotation of 2a = 0–360° is conducted. Fig. 11 shows the strain path induced during the imposed stress path. The numerical result on A-1 exhibits a circular locus, while only a slight deviation from the origin can be found on A-2. On the other hand, the result on A-3 brings about a distinct inelastic strain at the end of a stress cycle as has been observed in the experimental results (e.g. Miura et al., 1986a,b). Therefore, it can be said that the inelastic strain during the principal stress axes rotation is not induced for the stresses close to isotropic state, but is induced with the increase in magnitude of the deviatoric stress leading to the induced anisotropy. These results are supported not only by experimental evidences but also from the viewpoint of the phenomenological elastoplasticity, where these properties are discussed in relation to the non-coaxiality of the material response induced by both the anisotropic yield function and the tangential (or vertex) effect (Tsutsumi and Hashiguchi, 2005; Tsutsumi et al., 2005).

1984

S. Tsutsumi, K. Kaneko / International Journal of Plasticity 24 (2008) 1967–1989

1

1

path A

||σ *||/ p

0.8

0.8

proportional loading

0.6

0.6

A-3

||σ *||/ p

εv

A-2

0.4

0.4

εv

A-1

0.2 0

0.2 0

−0.2 0

0.4

0.8

1.2

||ε ∗|| (% )

−0.2 1.6

Fig. 10. Stress–strain relationship under the triaxial stress and constant pressure conditions together with the rotational path points A-1, A-2 and A-3.

0.4

0.2

A-1

ε zy (% )

A-2 0

A-3 −0.2

−0.4 −0.6

−0.4

−0.2

0

0.2

(ε zz − ε yy) /2 (%) Fig. 11. Strain responses in principal stress axes rotation from stress states A-1, A-2 and A-3.

Fig. 12 shows the stress–strain relationships under the proportional loading path B, and from the loading points B-1 and B-2, principal axes rotation tests are conducted. Fig. 13 shows the obtained strain path and we can find that the numerical result on B-1 exhibits a circular locus, whilst the result from B-2 brings about an inelastic strain at the end of a stress cycle. The stress-strain relationships from the first unloaded stress state pc to the followed triaxial compression condition is shown in Fig. 14. Fig. 15 shows the strain responses to the principal axes rotation tests from the loading point C-2 together with those obtained for A-2 and B-2. As can be seen from these figures, the numerical

S. Tsutsumi, K. Kaneko / International Journal of Plasticity 24 (2008) 1967–1989

1985

1

1

path B 0.8

0.8 0.6

proportional loading

|| σ *||/ p

B-2

0.4

εv

0.4

B-1

0.2

0.6

0.2

0

0

−0.2 0

0.4

0.8

−0.2 1.6

1.2

||ε ∗|| (% ) Fig. 12. Stress–strain relationship under the plane stress and constant pressure conditions together with the rotational path points B-1 and B-2.

0.2

0.1

B-1

ε zy (% )

0

B-2 −0.1

−0.2 −0.3

−0.2

−0.1

0

0.1

(ε zz − ε yy ) /2 (% ) Fig. 13. Strain responses in principal stress axes rotation from stress states B-1 and B-2.

result for the stress path C-2 coincides with that of A-2, both of which was controlled so as to fulfill the triaxial condition. Thus, it can be concluded that the strain responses during a principal stress axes rotation are not influenced by the change of mean stress history, and is uniquely given by the particular state of stress. Also, we can say that the deviation of the result of B-2 calculated under the plane stress condition from those of A-2 and C-2 with the triaxial condition is induced by the difference of loading condition and the resulting induced anisotropy. Fig. 16 shows the stress-strain relationships

1986

S. Tsutsumi, K. Kaneko / International Journal of Plasticity 24 (2008) 1967–1989

1

1

path C 0.8

0.8 proportional loading

0.6

0.6

εv

C-2

|| σ*|| / p 0.4

0.4 0.2

0.2

0

0 −0.2

0

0.4

0.8

1.2

−0.2 1.6

|| ε ∗|| ( %) Fig. 14. Stress–strain relationship under the triaxial stress condition together with the rotational path point C-2.

0.2

0.1

ε zy ( %)

A-2, C-2

0

−0.1 B-2

−0.2 −0.3

−0.2

−0.1

(ε zz − ε yy ) / 2

0

0.1

( %)

Fig. 15. Strain responses in principal stress axes rotation from stress states A-2, B-2 and C-2.

under the proportional loading path B, while their loading directions represented by b-value are quite different as shown in Fig. 9b, and principal axes rotation tests are conducted from the loading points B2a, b, c, and d in Figs. 9b and 15. As can be seen from the results in Fig. 17, the obtained strain responses for B-2a (or B-2) and B-2c are almost symmetry with identical size and shape. On the other hands, a slight difference is observed between B-2b and B-2d at the end point of strain responses. Following the results of Section 3, this should be induced by the anisotropic material property as has been investigated by numerous researchers (cf. Oda, 1972; Matsuoka and Nakai, 1974; Sekiguchi and Ohta, 1977; Oda et al., 1985; Oda, 1993; Houlsby and Puzrin, 2000; Hashiguchi, 2001; Nemat-Nasser and Zhang, 2002; Dean, 2005; Tashman et al., 2005; Sansour et al., 2006). Although, numerical tests under a large numbers of stress states and loading histories have to be done to derive absolute conclusions, from a theoretical point of view, the phenomenological description of the deformation behavior of

S. Tsutsumi, K. Kaneko / International Journal of Plasticity 24 (2008) 1967–1989

1

1987

1

path B 0.8

0.8 B-2c

0.6

|| σ*||/ p

0.6 B-2b

B-2a

0.4

0.4

εv

B-2d

0.2

0.2

0

0

−0.2 0

0.4

0.8

1.2

−0.2 1.6

||ε ∗|| (% ) Fig. 16. Stress–strain relationships under the triaxial stress condition of path B together with the rotational path points B-2a, b, c and b.

0.4

B -2d 0.2

ε zy (% )

0

B -2a (B-2)

B -2c

−0.2

B -2b −0.4 −0.4

−0.2

0

(ε zz − ε yy) /2

0.2

0.4

(% )

Fig. 17. Strain responses in principal stress axes rotation from stress states B-2a, b, c and b.

elastoplastic materials to non-proportional loading conditions has been examined by one of anthers (Tsutsumi and Hashiguchi, 2005). They adopt an elastoplastic constitutive mode based on subloading surface model with both tangent effect and anisotropy of yield surface (Hashiguchi and Tsutsumi, 2001, 2003) and concluded that both tangential stress-rate effect and anisotropy of yield condition should be introduced to describe an inelastic response under a principal stress axes rotations.

1988

S. Tsutsumi, K. Kaneko / International Journal of Plasticity 24 (2008) 1967–1989

5. Concluding remarks The macroscopic constitutive response of granular media to the principal stress axes rotation has been examined numerically under quasi-static equilibrium situation. Also, the two-scale analysis adopting GEM has been applied for a 3D unit cell with periodic boundary condition. The mechanical response of the granular element model to the non-proportional loading exhibits the following characteristics. (1) Inelastic strain from the pure principal stress rotation is not induced for the stresses close to the isotropic state, but is induced with increase in magnitude of the deviatoric stress and corresponding induced anisotropy. (2) The strain response during a principal stress axes rotation is not influenced by the mean-stress history, and is uniquely governed by the current state of stress and loading condition. These results are supported not only by experimental evidences but also from the viewpoint of the phenomenological elastoplasticity with both the anisotropic yield function and the tangential effect. However, the numerical tests under a large numbers of stress states have to be done to derive absolute conclusions for the general loading conditions. Acknowledgements The authors would like to express our sincere gratitude to Professor K. Hashiguchi, Department of Civil and Environmental Engineering, College of Technology, Daiichi University, Japan. He provided us with valuable advices and discussions on this study. References Anandarajah, A., 2008. Multi-mechanism anisotropic model for granular materials, I. Int. J. Plasticity 24 (5), 804–846. Allaire, G., 1992. Homogenization and two-scale convergence. SIAM J. Math. Anal. 23, 1482–1518. Arthur, J.R.F., Chua, K.S., Dunstan, T., Rodoriguez, J.I., 1980. Principal stress rotation: a missing parameter. Proc. ASCE 106 (GT 4), 419–433. Arthur, J.R.F., Bekenstein, S., Germaine, J.T., Ladd, C.C., 1981. Stress path tests with controlled rotation of principal stress directions. Lab. Shear Strength Soil ASTM, STP 740, 516–540. Bardet, J.P., 1994. Numerical simulations of the incremental responses of idealized granular materials. Int. J. Plasticity 10 (8), 879–908. Barlat, F., Duarte, J.M.F., Gracio, J.J., Lopes, A.B., Rauch, E.F., 2003a. Plastic flow for non-monotonic loading conditions of an aluminum alloy sheet sample. Int. J. Plasticity 19, 1215–1244. Barlat, F., Brem, J.C., Yoon, J.W., Chung, K., Dick, R.E., Lege, D.J., Pourboghrat, F., Choi, S.-H., Chu, E., 2003b. Plane stress yield function for aluminum alloy sheets - Part 1: theory. Int. J. Plasticity 19, 1297–1319. Barlat, F., Yoon, J.W., Cazacu, O., 2007. On linear transformations of stress tensors for the description of plastic anisotropy. Int. J. Plasticity 23 (5), 876–896. Broms, B.B., Casbarian, A.O., 1965. Effects of rotation of the principal stress axes and of the intermediate principal stress on the shear strength. In: Proc. 6th ICSMFE, Montreal, vol. 1, pp. 179–183. Cundall, P.A., Strack, O.D.L., 1979. A discrete numerical model for granular assemblies. Geotechnique 29, 47–65. Dean, E.T.R., 2005. Patterns, fabric, anisotropy, and soil elasto-plasticity. Int. J. Plasticity 21 (3), 5851–5861. Gutierrez, M., Ishihara, K., Towhata, I., 1991. Flow theory for sand during rotation of principal stress direction. Soils Found. 31 (4), 121–132. Han, C.-S., Chung, K., Wagoner, R.H., Oh, S.-I., 2003. A multiplicative finite elasto-plastic formulation with anisotropic yield functions. Int. J. Plasticity 19, 197–211. Hashiguchi, K., 1977. An expression of anisotropy in plastic constitutive equations of soils. In: Proc. 9th ICSMFE, Specialty Session 9, Constitutive Equations of Soils, Tokyo, JSSMFE, pp. 302–305. Hashiguchi, K., 2001. Description of inherent/induced anisotropy of soils: rotational hardening rule with objectivity. Soils Found. 41 (6), 139–145. Hashiguchi, K., Chen, Z.-P., 1998. Elastoplastic constitutive equations of soils with the subloading surface and the rotational hardening. Int. J. Numer. Anal. Meth. Geomech. 22, 197–227. Hashiguchi, K., Tsutsumi, S., 2001. Elastoplastic constitutive equation with tangential stress rate effect. Int. J. Plasticity 17 (1), 117–145. Hashiguchi, K., Tsutsumi, S., 2003. Shear band formation analysis in soils by the subloading surface model with tangential stress rate effect. Int. J. Plasticity 19, 1651–1677. Hashiguchi, K., Tsutsumi, S., 2007. Gradient plasticity with the tangential-subloading surface model and the prediction of shearband thickness of granular materials. Int. J. Plasticity 23, 767–797. Houlsby, G.T., Puzrin, A.M., 2000. A thermomechanical framework for constitutive models for rate-independent dissipative materials. Int. J. Plasticity 16 (9), 1017–1047. Ishihara, K., Towhata, I., 1983. Sand response to cyclic rotation of principal stress directions as induced by wave loads. Soils Found. 23 (4), 11–26. Joer, H.A., Lanier, J., Fahey, M., 1998. Deformation of granular materials due to rotation of principal axes. Geotechnique 45, 605– 619.

S. Tsutsumi, K. Kaneko / International Journal of Plasticity 24 (2008) 1967–1989

1989

Kaneko, K., 2001. Doctor Thesis, Tohoku University. Kaneko, K., Terada, K., Kyoya, T., Kishino, Y., 2003. Global-local analysis of granular media in quasi-static equilibrium. Int. J. Solids Struct. 40, 4043–4069. Kikuchi, N., Oden, J.T., 1988. Contact Problem in Elasticity: A Study of Variational Inequalities and Finite Element Methods. SIAM, Philadelphia. Kishino, Y., 1989. Computer Analysis of Dissipation Mechanism in Granular Media. Powders and Grains. A.A. Balkema, Rotterdam. pp. 323-330. Kishino, Y., Akaizawa, H., Kaneko, K., 2001. On the Plastic Flow of Granular Materials. Powders and Grains. A.A. Balkema, Rotterdam. pp. 199-202. Kishino, Y., 2002. The incremental nonlinearity observed in numerical tests of granular media. In: CD-ROM Proc. 15th ASCE Eng. Mech. Conf. Kuhn, M.R., 1999. Structured deformation in granular materials. Mech. Mater. 31, 407–429. Kyriakides, S., Corona, E., Miller, J.E., 2004. Effect of yield surface evolution on bending induced cross sectional deformation of thin-walled sections. Int. J. Plasticity 20, 607–618. Lions, J.L., 1981. Some Methods in the Mathematical Analysis of Systems and their Control. Science Press, Beijing, China. Lewin, P.I., Burland, J.B., 1970. Stress prove experiments on saturated normally consolidated clay. Geotechnique 20, 38–56. Matsuoka, H., Nakai, T., 1974. Stress-deformation and strength characteristics of soil under three different principal stresses. Proc. JSCE 232, 59–70. Miura, K., Miura, S., Toki, S., 1986a. Deformation behavior of anisotropic dense sand under principal stress axes rotation. Soils Found. 26 (1), 36–52. Miura, K., Toki, S., Miura, S., 1986b. Deformation prediction for anisotropic sand during the rotation of principal stress axes. Soils Found. 26 (3), 42–56. Nemat-Nasser, S., Zhang, J., 2002. Constitutive relations for cohesionless frictional granular materials. Int. J. Plasticity 18, 531– 547. Nicot, F., Darve, F.RNVO Group, 2005. A multi-scale approach to granular materials. Mech. Mater. 37, 980–1006. Nicot, F., Darve, F., 2007a. Basic features of plastic strains: From micro-mechanics to incrementally nonlinear models. Int. J. Plasticity 23, 1555–1588. Nicot, F., Darve, F., 2007b. A micro-mechanical investigation of bifurcation in granular materials. Int. J. Solids Struct. 44, 6630– 6652. Nicot, F., Darve, F., 2007c. Micro-mechanical bases of some salient constitutive features of granular materials. Int. J. Solids Struct. 44, 7420–7443. Oda, M., 1972. Initial fabrics and their relations to mechanical properties of granular material. Soils Found. 12 (1), 17–36. Oda, M., 1993. Inherent and induced anisotropy in plasticity theory of granular soils. Mech. Mater. 16, 35–45. Oda, M., Nemat-Nasser, S., Konish, J., 1985. Stress-induced anisotropy in granular masses. Soils Found. 25 (3), 85–97. Oda, M., Iwashita, K. (Eds.), 1999. Mechanics of Granular Materials: An Introduction. Balkema, A.A., Rotterdam. Poorooshasb, H.B., Holubec, I., Sherbourne, A.N., 1966. Yielding and flow of sand in triaxial compression. Canadian Geotech. J. 3, 178–190. Pradel, D., Ishihara, K., Gutierrez, M., 1990. Yielding and flow of sand under principal stress axes rotation. Soils Found. 30 (1), 87–99. Rudnicki, J.W., Rice, J.R., 1975. Conditions for localization of deformation in pressure-sensitive dilatant materials. J. Mech. Phys. Solids 23, 371–394. Sansour, C., Karšaj, I., Soric´, J., 2006. A formulation of anisotropic continuum elastoplasticity at finite strains. Part I: Modelling. Int. J. Plasticity 22 (12), 2346–2365. Sekiguchi, H., Ohta, H., 1977. Induced anisotropy and time dependency in clays. In: Proc. Spec. Session 9, 9th ICSMFE, Constitutive Equations of Soils, Tokyo, pp. 229–239. Symes, M.J., Gens, A., Hight, D.W., 1984. Undrained anisotropy and principal stress rotation in saturated sand. Geotechnique 34 (1), 11–27. Tashman, L., Masad, E., Little, D., Zbib, H., 2005. A microstructure-based viscoplastic model for asphalt concrete. Int. J. Plasticity 21 (9), 1659–1685. Terada, K., Kikuchi, N., 2001. A class of general algorithms for nonlinear multi-scale analyses of heterogeneous media. Comp. Meth. Appl. Mech. Eng. 190, 5427–5464. Tsutsumi, S., Kaneko, K., Toyosada, M., Hashiguchi, K., Kishino, Y., 2005. Non-coaxial constitutive response of idealized 3D granular assemblies to rotation of principal stress axes. J. Appl. Mech., JSCE 8, 565–571. Tsutsumi, S., Hashiguchi, K., 2005. General non-proportional loading behaviour of soils. Int. J. Plasticity 21 (10), 1941–1969. Wren, J.R., Borja, R.I., 1997. Micromechanics of granular media, Part II: Overall tangential moduli and localization model for periodic assemblies of circular disks. Comput. Meth. Appl. Mech. Eng. 141, 221–246. Yatomi, C., Yashima, A., Iizuka, A., Sano, I., 1989. General theory of shear bands formation by a non-coaxial Cam-clay model. Soils Found. 29, 41–53. Zhu, H., Mehrabadi, M.M., Massoudi, M., 2006. Three-dimensional constitutive relations for granular materials based on the dilatant double shearing mechanism and the concept of fabric, Int. J. Plasticity 22 (5), 826–857.