A three-invariant cap model with isotropic–kinematic hardening rule and associated plasticity for granular materials

Available online at www.sciencedirect.com

International Journal of Solids and Structures 45 (2008) 631–656 www.elsevier.com/locate/ijsolstr

A three-invariant cap model with isotropic–kinematic hardening rule and associated plasticity for granular materials H. DorMohammadi, A.R. Khoei

*

Center of Excellence in Structural and Earthquake Engineering, Department of Civil Engineering, Sharif University of Technology, P.O. Box 11365-9313, Tehran, Iran Received 16 October 2006; received in revised form 24 August 2007 Available online 8 September 2007

Abstract In this paper, a three-invariant cap model is developed for the isotropic–kinematic hardening and associated plasticity of granular materials. The model is based on the concepts of elasticity and plasticity theories together with an associated flow rule and a work hardening law for plastic deformations of granulars. The hardening rule is defined by its decomposition into the isotropic and kinematic material functions. The constitutive elasto-plastic matrix and its components are derived by using the definition of yield surface, material functions and non-linear elastic behavior, as function of hardening parameters. The model assessment and procedure for determination of material parameters are described. Finally, the applicability of proposed plasticity model is demonstrated in numerical simulation of several triaxial and confining pressure tests on different granular materials, including: wheat, rape, synthetic granulate and sand. Ó 2007 Elsevier Ltd. All rights reserved. Keywords: Three-invariant plasticity; Cap-plasticity surface; Isotropic–kinematic hardening; Granular materials

1. Introduction Modeling the behavior of granular materials under various loading conditions poses several important and challenging problems. It is a characteristic of granular materials that the relative motion between the grains leads to irrecoverable deformation, and the use of plasticity theory is necessary as a theoretical framework for the mechanical behavior of granular media. The characterization of the stress–strain and failure behavior of granular materials is so complex due to particulate nature. It can be particularly observed in triaxial tests that starting from different confining pressure and initial void ratio, such as: hardening, densification, dilation, or combinations of these behaviors. Thus, it is of great interest to be able to accurately predict the mechanical behavior of granular materials by an appropriate constitutive model. This study focuses on the frictional effects of granular materials based on a combined isotropic–kinematic hardening rule and associated plasticity.

*

Corresponding author. Tel.: +98 21 6600 5818; fax: +98 21 6601 4828. E-mail address: [email protected] (A.R. Khoei).

0020-7683/$ - see front matter Ó 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.ijsolstr.2007.08.019

632

H. DorMohammadi, A.R. Khoei / International Journal of Solids and Structures 45 (2008) 631–656

There are a number of plasticity models available in the literature for constitutive relations of granular materials, such as: sand, concrete and rock. A review of various plasticity models for frictional and granular materials can be seen in Desai and Siriwardane (1984), Chen and Baladi (1985) and Lewis and Schrefler (1998). Basically, there are two different approaches in constitutive modeling of granular materials; ‘the micromechanical model’ and ‘the macromechanical model’. The micromechanical models are considered as ideal approaches, and cannot be directly applied to analyze general boundary value problems (Jenkins and Strack, 1993; McDowell and Bolton, 1998; Chang and Hicher, 2005). Therefore, the continuum approaches which consist of implementation of macroscopic constitutive models based on plasticity theories have been proposed for granular materials. The crucial factor that governs the behavior of granular soils is stress-dilatancy, the soil’s ability to increase in volume due to shear stresses and geometrical effects. A plasticity theory was developed by Dorris and Nemat-Nasser (1982) for finite deformation of granular materials, which accounts for the true stress triaxiality, pressure sensitivity and dilatancy. They captured the effect of stress triaxiality by including the third deviatoric stress invariant in the yield function and the flow potential. A simple constitutive model was developed by Wan and Guo (1998) within the plasticity theory by introducing a modified stress-dilatancy law. They modeled the dependencies of granular soil behavior on void ratio and stress by a modified stress-dilatancy law. By using a void ratio-dependent factor that measures the deviation of the current void ratio from the critical one, Rowe’s stress dilatancy equation was modified. This modification was related to a new energy dissipation equation for a granular assembly which sustain the kinematical constraints under the action of stresses. The experimental results of Watson and Wert (1993) and Brown and Abou-Chedid (1994) demonstrated that a ‘two-mechanism-model’, such as: Drucker–Prager or Mohr–Coulomb and elliptical cap models can be utilized to construct the suitable phenomenological constitutive models which capture the major features of the response of geological and frictional materials. These models consist of two yield surfaces; a ‘distortion surface’ and a ‘consolidation’ or ‘cap’ surface, which has an elliptical shape. The distortion surface controls the ultimate shear strength of material and the cap-surface captures the hardening behavior of material under compression. Gudehus (1996) presented a constitutive model by incorporation of density and pressure dependencies for granular materials within hypoplasticity formulation. A cone-cap model based on a density-dependent Drucker–Prager yield surface and a non-centered ellipse was developed by Brandt and Nilsson (1999). A double-surface plasticity model was proposed by Lewis and Khoei (2001) for the non-linear behavior of porous materials in the concept of the generalized plasticity formulation. The non-linear behavior of frictional and granular materials is adequately described by double-surface plasticity models. However, it suffers from a serious deficiency when the stress-point reaches in the intersection of these two different yield functions. In the flow theory of plasticity, the transition from an elastic state to an elasto-plastic state appears more or less abruptly. For double-surface plasticity models, it is very difficult to define the location of yield surface and special treatment must be performed to avoid numerical difficulties in the intersection of these two surfaces (Khoei et al., 2004). Lade and Kim (1995) proposed a single hardening plasticity model for frictional materials, which does not have such a drawback. Their model was based on a single, isotropic yield surface shaped as an asymmetric teardrop with the pointed apex at the origin of the principal stress space. The yield surface was expressed in terms of stress invariants to describe the locus at which the total plastic work is constant. The total plastic work due to shear and volumetric strains was served as the hardening parameter, and used to define the location and shape of the yield surface. The non-associated flow rule was derived from a potential function which describes a three-dimensional surface shaped with an asymmetric cross-section. A non-associated plasticity theory was developed by Krenk (2000) and Krenk and Ahadi (2000) for granular materials based on the concept of a characteristic stress state of vanishing incremental dilation. The theory was presented based on a common format for the yield and flow potential surfaces, representing the surfaces in terms of stress invariants and a single shape function. The flow potential surface was determined by an approximate friction hypothesis, and the plastic work hardening was introduced in a linear invariant form that permits dilation before the ultimate state, by including the work associated with shape change in addition to the traditional contribution from volume change. Bigoni and Piccolroaz (2004) proposed a yield/damage function for modeling the inelastic behavior of a broad class of pressure-sensitive, frictional, ductile and brittle-cohesive materials. The yield function was based on a single, convex and smooth surface in stress space approaching as limit situations well-known criteria and the extreme limits of convexity in

H. DorMohammadi, A.R. Khoei / International Journal of Solids and Structures 45 (2008) 631–656

633

the deviatoric plane. A three-invariant isotropic/kinematic hardening cap plasticity model was presented by Foster et al. (2005) for geomaterials based on an associated flow rule. Their yield surface was an exponential shear failure function at its core with triangularity in p-plane that captures the pressure-dependence of the shear strength of geomaterials. Recently, a single cap plasticity with an isotropic hardening rule and associated plasticity was developed by Khoei and Azami (2005), which could generate some common elliptical yield surface and double-surface plasticity models, as special cases. In the present study, the cap-plasticity model developed by Khoei and Azami (2005), Khoei and DorMohammadi (2007) and Khoei et al. (2007) is extended to three-invariant cap plasticity with isotropic–kinematic hardening rule and associated plasticity for granular materials. The paper focuses on different aspects of material behavior, including: isotropic and kinematic hardening and dilation. A generalized framework for the cap plasticity is presented based on three invariants of stress states, which describes the isotropic and kinematic hardening behavior of material. The constitutive elasto-plastic matrix and its components are extracted. The model assessment and parameter determination are described. Finally, the applicability of the model is demonstrated in several numerical examples. 2. Three-invariant isotropic–kinematic cap plasticity In order to present a constitutive plasticity model for granular media, several requirements must be taken into account (Collins, 2005). Firstly, the soils and granular materials are multiple phase materials, consisting of solid grains together with voids. It is therefore necessary that the model considers the effects of porosity and voids ratio during loading conditions. Secondly, the elasto-plastic behavior depends on the current stress state and the pressure history, particularly the pre-consolidation pressure, which is the fundamental hardening parameter in granular materials. Thus, the material elasto-plasticity matrix must be defined as a function of hardening parameter. Thirdly, granular materials undergo significant plastic volume changes, and present the dilation and contraction depending on the current voids ratio and pressure. The dilation can be caused due to the fact that individual grains must ride up over each other in order to accommodate shear strains. Thus, the model must capture both the dilation and compaction behavior of material. Basically, there are various procedures in the literature for the derivation of plasticity models. Some researchers used the concept of characteristic stress state by using surfaces that can be characterized by their ‘traces’ in the deviatoric plane and in planes containing a ‘radial’ centre-line (Krenk, 2000). Other researchers applied the shear failure function which captures the pressure-dependence of the shear strength (Foster et al., 2005). The current research is based on the second procedure to derive the constitutive plasticity model following the works of senior author given in Khoei et al. (2003) and Khoei and Bakhshiani (2004). In these references, an endochronic plasticity theory was developed based on coupling between deviatoric and hydrostatic behavior for pressure-sensitive materials. Although the concept of yield surface has not been explicitly assumed in endochronic theory, it was presented by Khoei and Bakhshiani (2004) that the single-cap plasticity can be derived as a special case of their endochronic model. 2.1. The yield function The constitutive model is developed here for both isotropic and kinematic hardening behaviors based on three invariants of stress states, J1, J2D and J3D. As the first step in derivation of plasticity model, a shear-failure function is defined as  2 fd F f ¼ fd2  J 21 ð1Þ fh where J1 is the first invariant of stress tensor, and fh and fd are the hydrostatic and deviatoric material functions, respectively, and will be defined later as functions of hardening parameters. The function Ff captures the pressure-dependence of the shear strength of material, which increases with more compressive mean stresses, as shown in Fig. 1. It will be shown later that different combinations of fh and fd lead to different shapes of the yield surface. The initial yield surface f0 is offset from the pressure-dependence function Ff, as shown in Fig. 2.

634

H. DorMohammadi, A.R. Khoei / International Journal of Solids and Structures 45 (2008) 631–656

Ff

Ff = f d2

− J1

J1 = − f h Fig. 1. The shear-failure function Ff.

Ff

Ff

f0

− J1 Fig. 2. The shear-failure function Ff and the initial yield surface f0.

Considering the shear-failure function Ff, the yield function can be written in (J1, J2D) stress space as 2 F 1 ¼ J 2D  F f ¼ 0 3

ð2Þ

 2 2 fd F 1 ¼ J 2D þ J 21  fd2 ¼ 0 3 fh

ð3Þ

or

where pffiffiffiffiffiffiffiffi J2D is the second invariant of deviatoric stress tensor. In Fig. 3, the yield function F1 is presented in ð J 2D ; J 1 Þ stress space. In above relation, we define the coefficients /h and /d in order to indicate the effect of material functions fh and fd in the yield function. The yield function (3) can be therefore rewritten as  2 2 /d fd 2 F 2 ¼ J 2D þ J 21  ð/d fd Þ ¼ 0 ð4Þ 3 /h fh Finally, applying the effect of third invariant of deviatoric stress tensor J3D into Eq. (4), i.e. the triangularity of deviatoric trace along the hydrostatic axis, the yield function can be written in its final representation as  2 2 /d fd 2=3 2 F ðr; a; jÞ ¼ wJ 3D þ J 2D þ J 21  ð/d fd Þ ¼ 0 ð5Þ 3 /h fh

H. DorMohammadi, A.R. Khoei / International Journal of Solids and Structures 45 (2008) 631–656

635

J 2D

F1

− J1 pffiffiffiffiffiffiffiffi Fig. 3. The cap-surface function F1 in ð J 2D ; J 1 Þ stress space.

where w is a constant, and a and j are the kinematic and isotropic hardening parameters, respectively. It must be mentioned that the material functions fh and fd can be decomposed into two parts, the isotropic and kinematic parts, which control the shape of yield surface (5). Fig. 4 presents the 3D representation of yield surface (5) for the isotropic and kinematic hardening behavior of material. The isotropic and kinematic hardening parameters j and a evolve with plastic deformation. The evolution of j is related to the mean stress, and more directly to the volumetric plastic strain epv , i.e. j ¼ epv , while the evolution of a is related to the deviatoric plastic strain ep. As can be expected, the kinematic hardening parameter a = {a1 a2 a3}T can be decomposed in two directions J1 and J2D in meridian plane, which contains two parts as follows     a ¼ a1 exp a2 ððep ÞT : ep Þa3 m þ a4 ððep ÞT : ep Þa3 ep ð6Þ where the first term controls the movement of yield surface in J1-axis and the second term controls the movement of yield surface in perpendicular direction to J1-axis. In relation (6), m is the unit vector, defined as m = {1 1 1}T, and a1, a2, a3 and a4 are the material parameters. The three components of kinematic hardening

Fig. 4. Trace of 3D representation of yield surface (5) in principal stress space for the isotropic and kinematic hardening behavior of material.

636

H. DorMohammadi, A.R. Khoei / International Journal of Solids and Structures 45 (2008) 631–656

parameter a, i.e. a1, a2 and a3, determine the values of the yield surface movements in the directions of principal stresses r1, r2 and r3, respectively. 2.2. The isotropic and kinematic material functions As mentioned earlier, the material functions fh and fd control the size and movement of the yield surface, and are functions of hardening parameters. Then, these functions can be decomposed into two parts, i.e. the isotropic and kinematic parts, as fh ¼ fhisotropic þ fhkinematic

ð7Þ

fd ¼ fdisotropic þ fdkinematic

ð8Þ

The isotropic part of material functions fh and fd are the exponential increasing functions of the isotropic hardening parameter j ¼ epv , defined as fhisotropic ¼ ðb1 þ b2 expðb3 epv ÞÞdðepv Þ

ð9Þ

c2 expðc3 epv ÞÞdðepv Þ

ð10Þ

fdisotropic ¼ ðc1 þ

where b1, b2, b3, c1, c2 and c3 are the material parameters and dðepv Þ is defined as  1 if epv 6¼ 0 dðepv Þ ¼ 0 if epv ¼ 0

ð11Þ

In order to determine the kinematic parts of material functions (7) and (8), consider two different stress spaces ri and Ni; the former is located in the center of the yield surface before kinematic hardening and the later is placed in the center of the yield surface after kinematic hardening (Fig. 5). The distance of centers of two coordinate systems can be defined by ai, in which the relationship between the principal stresses in two stress spaces is defined as Ni ¼ ri þ ai

ð12Þ

Ξ3 Ξ2 σ3 Ξ1

α σ2 σ1

Fig. 5. The stress spaces before kinematic hardening ri and after kinematic hardening Ni.

H. DorMohammadi, A.R. Khoei / International Journal of Solids and Structures 45 (2008) 631–656

637

Defining the parameters J1a and J2Da similar to the definitions of the invariants of stress and deviatoric stress tensors as J 1a ¼ a1 þ a2 þ a3 i 1h J 2Da ¼ ða1  a2 Þ2 þ ða2  a3 Þ2 þ ða3  a1 Þ2 6

ð13Þ ð14Þ

Considering the definition of the principal stresses in two stress spaces, defined by Eq. (12), the three-invariant single plasticity (5) can be written in new stress space as  2 2 /d fd 2=3 J 2I  /2d fd2 ¼ 0 ð15Þ wJ IIID þ J IID þ 3 /h fh where JI is the first invariant of stress tensor and JIID and JIIID are the second and third invariants of deviatoric stress tensors in the second stress space, defined as J I ¼ J 1 þ J 1a J IID ¼ J 2D þ J 2Da þ J ra where J ra

ð16Þ ð17Þ

!   pffiffiffi 2 pffiffiffiffiffiffiffiffi 1 1 3 sin xða2  a3 Þ ¼ pffiffiffi J 2D cos x a1  a2  a3 þ 2 2 2 3

ð18Þ

where

! pffiffiffi 1 3 3 J 3D 1 x ¼ cos 3=2 3 2 J 2D

0 6 x 6 60

ð19Þ

Substituting relations (16) and (17) into the yield surface (15) in the absence of third invariants of deviatoric stress with respect to Eq. (5), results in  2 /d fd 2 ðJ 2Da þ J ra Þ ð20Þ ¼ 3 ð2J 1a J 1 þ J 21a Þ /h fh According to Eq. (20), the kinematic parts of material functions fh and fd can be therefore written as 1 1=2 ð2J 1a J 1 þ J 21a Þ /h  1=2 1 2 ¼  ðJ 2Da þ J ra Þ /d 3

fhkinematic ¼

ð21Þ

fdkinematic

ð22Þ

Finally, the material functions fh and fd can be defined by substituting Eqs. (9), (10), (21) and (22) into (7) and (8) as 1 1=2 ð2J 1a J 1 þ J 21a Þ /h  1=2 1 2 fd ¼ ðc1 þ c2 expðc3 epv ÞÞdðepv Þ þ  ðJ 2Da þ J ra Þ /d 3

fh ¼ ðb1 þ b2 expðb3 epv ÞÞdðepv Þ þ

ð23Þ ð24Þ

3. Model assessment In order to demonstrate the performance of the proposed plasticity model in prediction of granular material behavior, the experimental tests must be performed to determine and calibrate the parameters of material functions fh and fd, defined by (23) and (24), in the yield surface (5). These two material functions control the size and movement of the yield surface, and are decomposed into the isotropic and kinematic parts, given by

638

H. DorMohammadi, A.R. Khoei / International Journal of Solids and Structures 45 (2008) 631–656

(7) and (8), as functions of the hardening parameters j and a, or directly the plastic volumetric strain epv and the deviatoric plastic strain ep. It must be noted that the kinematic hardening parameter a indicates the movement of the yield surface in the direction of J1-axis and perpendicular direction to J1-axis. It is worth mentioning that different values of material functions fh and fd result in different aspects of the yield surface (5). Consider that the first two terms of Eq. (5) are zero, it leads to  2 /d fd J 21  /2d fd2 ¼ 0 ð25Þ /h fh The above equation generally yields to three roots, the points of intersection of yield surface with J1-axis, i.e. J1 = ± /h fh and one more from fd = 0.0, which has been defined in Eq. (24). If fd = 0.0 does not lead to any value for J1, the yield surface of Eq. (5) yields to two roots for J1, i.e. J1 = ± /hfh, which results in an elliptical shape in meridian plane. The 3D representation of the yield surface (5) is shown in Fig. 6 in principal stress space for different values of isotropic hardening, where the intersection point of yield surface with J1-axis are

σZ

20 10 0 -10 -20 -20

σX

-20 -10

-10 0

0 10

σY

10 20

20

Fig. 6. Trace of 3D elliptical yield function in principal stress space for the isotropic hardening behavior.

H. DorMohammadi, A.R. Khoei / International Journal of Solids and Structures 45 (2008) 631–656

639

/hfh and +/hfh. This representation clearly shows how the yield surface grows with densification, eventually becoming independent of the hydrostatic stress J1 at full dense material, where the von-Mises yield surface is generated. This yield surface is very similar to the elliptical yield functions developed by authors for porous metals based on an extension of von-Mises’s concept (Doraivelu et al., 1984; Oliver et al., 1996). If fd = 0.0 leads to the value of J1 between /hfh and +/hfh, the cone-cap yield surface will be produced from Eq. (5) based on the intersection points of J1 =  /hfh and the value obtained from fd = 0.0 (Fig. 3). The 3D representation of Fig. 3 is shown in Fig. 7 for different values of isotropic hardening. As can be observed from this figure the yield surface grows with densification and reduces to the Drucker–Prager yield function for full dense bodies. This yield surface is very similar to the double-surface cap models, i.e. a combination of Mohr–Coulomb or Drucker–Prager and elliptical surfaces, which has been extensively used by authors to demonstrate the behavior of granular and powder materials (Brandt and Nilsson, 1999; Lewis and Khoei, 2001; Khoei et al., 2004). It is worth mentioning that as the double-surface plasticity consists of two different yield functions, special treatment should be made to avoid numerical difficulties in the intersection of these two surfaces (Khoei et al., 2004), however – the single yield surface (5) does not have such a drawback. Fig. 8 presents the effect of parameter w in the yield surface (5) that causes triangularity of deviatoric trace along the

σZ

25 20 15 10 5 0 -5 -10 -10

σX-5

0

5

10

15

20

25 25

20

15

10

5

0

-5

-10

σY

Fig. 7. Trace of 3D cone-cap yield function in principal stress space for the isotropic hardening behavior.

640

H. DorMohammadi, A.R. Khoei / International Journal of Solids and Structures 45 (2008) 631–656

σZ

25 20 15 10 5 0 -5 -10 -10

σX

-10 -5

0

0

-5

σY

5

5 10

15

15 20

10

20 25 25

Fig. 8. 3D irregular hexagonal pyramid of the Mohr–Coulomb and cone-cap yield function in principal stress space for w = 0.1.

hydrostatic axis. This yield surface is similar to the irregular hexagonal pyramid of the Mohr–Coulomb and cone-cap yield surface employed by researchers for description of soil and geomaterial behavior (Lade and Kim, 1995). It must be noted that the single-surface plasticity (5) is not only for the isotropic compression part of the triaxial tests, and is capable to model the complete triaxial and confining pressure tests. 4. Parameter determination The important issue in prediction of granular material behavior is the identification of parameters of the proposed plasticity model. The calibration procedure for three-invariant isotropic–kinematic cap plasticity is carried out based on a series of isostatic and triaxial tests. An organized approach to determine model parameters is to utilize an optimization routine. Mathematically, an objective function and a search strategy are necessary for the optimization. The objective function, which represents the constitutive model, captures the material behavior and can be used in a simultaneous optimization against a series of experimental data. The simplest search strategy is based on the direct search approach, which is proved to be reliable and its relative simplicity make it quite easy to program into the code.

H. DorMohammadi, A.R. Khoei / International Journal of Solids and Structures 45 (2008) 631–656

641

For the proposed constitutive model, the total number of 10 material constants need to be determined for the material functions fh and fd. The parameters of isotropic parts of fh and fd (i.e. b1, b2, b3 and c1, c2, c3) are firstly evaluated using the confining pressure test, where the values of J2D and J3D in the yield surface (5) are zero. The parameters of the kinematic parts of fh and fd (i.e. a1, a2, a3 and a4) are then estimated performing the LSM method on the data obtained by a series of triaxial tests. These parameters control different aspects of predicted stress–strain curves obtained numerically by fitting the stress path to the triaxial and confining pressure tests, including: the slope, transition, expansion and contraction. a1 a2 a3 a4 b1 b2 b3 c1 c2 c3

controls the slope of ‘stress–axial strain’ and ‘radial strain–axial strain’ diagrams. controls the slope of ‘radial strain–axial strain’ diagram after dilation. controls the transition of ‘radial strain–axial strain’ diagram before dilation. controls the transition of ‘stress–axial strain’ and ‘radial strain–axial strain’ diagrams, and the slope of ‘stress–axial strain’ curve. controls the transition, expansion and contraction of ‘radial strain–axial strain’ diagram, and the slope, expansion and contraction of ‘stress–axial strain’ curve. controls the transition, expansion and contraction of ‘stress–axial strain’ and ‘radial strain–axial strain’ diagrams, and the slope of ‘radial strain–axial strain’ curve. controls the slope and transition of ‘radial strain–axial strain’ curve, and the slope of ‘stress–axial strain’ diagram. controls the slope, expansion and contraction of ‘stress–axial strain’ and ‘radial strain–axial strain’ diagrams. controls the slope, transition, expansion and contraction of ‘stress–axial strain’ diagram, and the slope of ‘radial strain–axial strain’ diagram. controls the slope of ‘stress–axial strain’ curve, and the slope and transition of ‘radial strain–axial strain’ diagram after dilation. The procedure of parameter determination is performed as follows;

Step 1: Based on the results obtained from the confining pressure test, the values of J1 are evaluated using the yield surface (5) where the values of J2D and J3D are zero. From the values of volumetric strain ev, the elastic and plastic volumetric strains, eev and epv , are estimated. The parameters b1, b2 and b3 in the isotropic part of fh are computed. The parameters c1, c2 and c3 in the isotropic part of fd are then calculated by a least square method on the data obtained from the confining pressure test. Step 2: Applying the results of triaxial tests and the isotropic parameters of fh and fd obtained from Step 1, the kinematic parameters of fh and fd, i.e. a1, a2 and a3, in the first term of relation (6) are first estimated. Parameter a4 in the second term of relation (6) is then obtained by performing the least square method on the data obtained from the triaxial tests.

5. Numerical simulation results 5.1. Triaxial tests on dry silo materials The first three examples refer to set of triaxial tests on wheat, rape and synthetic granulate materials performed by Kolymbas and Wu (1990). The device was a triaxial apparatus specially designed to test dry silo materials. The triaxial apparatus was designed in the Institute of Soil Mechanics and Rock Mechanics of the Karlsruhe University. The axial force was measured beneath the pressure chamber by a load cell with a precision of ±30 N. As the axial force was measured beneath the pressure chamber, the measurement was not influenced by the friction between the loading piston and the sealing or by the confining pressure. The lateral strain of the sample was measured directly by means of three collars which contact the sample in the upper, middle and lower parts, respectively. The axial and lateral deformations were measured during hydrostatic loading by the displacement transducer mounted between the end plates and the three collars.

642

H. DorMohammadi, A.R. Khoei / International Journal of Solids and Structures 45 (2008) 631–656

The specimens are prepared by pluviation. The specimens are cylindrical with the initial diameter of 100.0 mm and length of 100.0 mm. The axial load is exerted by moving the loading piston. The apparatus allows a maximum axial load of 100 kN. The maximum design confining pressure is 1400 kPa. 5.1.1. Compaction of wheat Based on the procedure described in preceding section, the material parameters of the yield surface calibrated for the wheat material are as follows a1 ¼ 5:0 e  11 b1 ¼ 1500:3

a2 ¼ 32:0 b2 ¼ 600:8

a3 ¼ 32:0 b3 ¼ 600:5

c1 ¼ 3000:0

c2 ¼ 500:0

c3 ¼ 83:5

a

a4 ¼ 5:0

1.8 1.6

Radial Strain (%)

1.4 1.2 1 0.8 0.6 0.4

Three-invariant isotropic-kinematic cap model

0.2 0

Experiment (Kolymbas and Wu,1990)

0

0.5

1

1.5

1

1.5

Axial Strain (%)

b

1.8 1.6

Radial Strain (%)

1.4 1.2 1 0.8 0.6 0.4 0.2 0

0

0.5

Axial Strain (%) Fig. 9. Radial strain versus axial strain in isostatic compression test for wheat; (a) a comparison between the numerical and experimental results, (b) schematically isotropic hardening behavior.

H. DorMohammadi, A.R. Khoei / International Journal of Solids and Structures 45 (2008) 631–656

643

In Fig. 9(a), the radial strain versus axial strain is presented for the hydrostatic pressure test. It shows a good agreement between the model and experimental results. In Fig. 9(b), the relevant yield surfaces corresponding to the isotropic hardening behavior has been shown schematically. This representation clearly shows how the yield surface grows isotropically by increasing the hydrostatic pressure. In this figure, the inner circle shows the initial yield surface and the outer one presents the yield surface at the specified point. In Figs. 10 and 11, the variations of volumetric strain and stress ratio (axial stress/radial stress) are depicted with the axial strain at confining pressures 50 kPa and 400 kPa, respectively, corresponding to complete triaxial tests. In Figs. 10(b) and 11(b), the relevant yield surfaces are plotted corresponding to the kinematic hardening behavior. In these figures, the isotropic hardening behavior can be clearly observed in triaxial tests due to different confining pressures. Obviously, the initial yield surface at the beginning of triaxial test with confining pressure 400 kPa is bigger than that presented at the confining pressure 50 kPa. Thus, both isotropic and kinematic hardening behavior can be observed in this example. As can be seen, the proposed model captures the behavior of wheat in complete triaxial experiment.

a

Volumetric Strain (%)

Stress Ratio

4 3 2 1

Three-invariant isotropic-kinematic cap model

0

Experiment (Kolymbas and Wu,1990)

-1 -2 -3 -4

0

5

10

15

Axial Strain (%)

b Stress Ratio

4 3 2 1

Volumetric Strain (%)

0 -1 -2 -3 -4

0

5

10

15

Axial Strain (%) Fig. 10. Stress ratio and volumetric strain versus axial strain in triaxial test for wheat (confining pressure = 50 kPa); (a) a comparison between the numerical and experimental results, (b) schematically isotropic–kinematic hardening behavior.

644

H. DorMohammadi, A.R. Khoei / International Journal of Solids and Structures 45 (2008) 631–656

a

Volumetric Strain (%)

Stress Ratio

4 3 2 1

Three-invariant isotropic-kinematic cap model

0

Experiment (Kolymbas and Wu,1990)

-1 -2 -3 -4

0

5

10

15

Axial Strain (%)

b Stress Ratio

4 3 2 1

Volumetric Strain (%)

0 -1 -2 -3 -4

0

5

10

15

Axial Strain (%) Fig. 11. Stress ratio and volumetric strain versus axial strain in triaxial test for wheat (confining pressure = 400 kPa); (a) a comparison between the numerical and experimental results, (b) schematically isotropic–kinematic hardening behavior.

5.1.2. Compaction of rape In the compaction of rape material, the material parameters calibrated for the yield surface are as follows a1 ¼ 5:0 e  4

a2 ¼ 9:0

a3 ¼ 2:0

b1 ¼ 900:0

b2 ¼ 170:5

b3 ¼ 400:0

c1 ¼ 1500:5

c2 ¼ 600:2

c3 ¼ 5:0

a4 ¼ 1:0

In Fig. 12(a), the variation of radial strain with axial strain is presented for the isostatic test. The isotropic behavior of material during the expansion of the yield surfaces is demonstrated in Fig. 12(b) by increasing the hydrostatic pressure. In Figs. 13(a) and 14(a), the variations of volumetric strain and stress ratio (axial stress/radial stress) with axial strain are depicted at confining pressures 200 kPa and 400 kPa, respectively, corresponding to complete triaxial tests. Also plotted in Figs. 13(b) and 14(b) are the relevant yield surfaces corresponding to the isotropic and kinematic hardening behavior.

H. DorMohammadi, A.R. Khoei / International Journal of Solids and Structures 45 (2008) 631–656

a

645

4.5 4

Radial Strain (%)

3.5 3 2.5 2 1.5 1

Three-invariant isotropic-kinematic cap model

0.5 0

Experiment (Kolymbas and Wu,1990)

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

3.5

4

4.5

5

Axial Strain (%)

b

4.5 4

Radial Strain (%)

3.5 3 2.5 2 1.5 1 0.5 0

0

0.5

1

1.5

2

2.5

3

Axial Strain (%) Fig. 12. Radial strain versus axial strain in isostatic compression test for rape; (a) a comparison between the numerical and experimental results, (b) schematically isotropic hardening behavior.

5.1.3. Compaction of synthetic granulate The material parameters corresponding to the yield function (4) calibrated for the synthetic granulate material are as follows a1 ¼ 0:8 e  6

a2 ¼ 15:0

a3 ¼ 14:0

b1 ¼ 100:0 c1 ¼ 30:0

b2 ¼ 70:7 c2 ¼ 910:0

b3 ¼ 150:0 c3 ¼ 0:95

a4 ¼ 3:0

In Fig. 15(a), the variation of radial strain with axial strain is presented for the isostatic test. The isotropic hardening behavior of granular material is illustrated in Fig. 15(b) due to expansion of the yield surfaces by increasing the hydrostatic pressure. Clearly, the effect of confining pressure on isotropic hardening of material can be observed in this figure. Figs. 16 and 17 correspond to the complete triaxial compression tests. The variations of the volumetric strain and stress ratio (axial stress/radial stress) are presented with the axial strain

H. DorMohammadi, A.R. Khoei / International Journal of Solids and Structures 45 (2008) 631–656

a

5

Stress Ratio

646

4 3 2 1 Three-invariantisotropic-kinematic capm odel

0

Volumetric Strain (%)

-1 Experiment(Kolymbasand Wu,1 990)

-2 -3 -4 -5 -6 -7 -8

0

5

10

Axial Strain (%)

b Stress Ratio

5 4 3 2 1 0

Volumetric Strain (%)

-1 -2 -3 -4 -5 -6 -7 -8

0

5

10

Axial Strain (%) Fig. 13. Stress ratio and volumetric strain versus axial strain in triaxial test for rape (confining pressure = 200 kPa); (a) a comparison between the numerical and experimental results, (b) schematically isotropic–kinematic hardening behavior.

in Figs. 16(a) and 17(a). Figs. 16(b) and 17(b) present the relevant yield surfaces corresponding to the isotropic and kinematic hardening behavior in complete triaxial tests. Both isotropic and kinematic hardening behavior with dilatancy can be observed in this example. 5.2. Triaxial tests on loose and dense sands The next two examples refer to set of triaxial tests on loose and dense sands performed by Krenk and Ahadi (2000). The two materials are loose Baskarp sand with an initial specific pore volume e = 0.85 (Borup and Hedegaard, 1995), and dense Lund sand with initial specific pore volume e = 0.55 (Ibsen and Jakobsen, 1996). In both tests, the initial isotropic confining pressure was p0 = 0.64 MPa. The shear modulus of the dense sand is nearly four times larger, and the elastic and elasto-plastic flexibilities are two to three times smaller. Clearly, this gives smaller strains in the dense sand. The height of the test specimens was equal to the diameter, thus preventing localized failure near the peak load (Krenk and Ahadi, 2000).

H. DorMohammadi, A.R. Khoei / International Journal of Solids and Structures 45 (2008) 631–656

a

647

Volumetric Strain (%)

Stress Ratio

4 3 2 1

Three-invariant isotropic-kinematic cap model

0

Experiment (Kolymbas and Wu,1990)

-1 -2 -3 -4

0

5

10

15

Axial Strain (%)

b Stress Ratio

4 3 2 1

Volumetric Strain (%)

0 -1 -2 -3 -4

0

5

10

15

Axial Strain (%) Fig. 14. Stress ratio and volumetric strain versus axial strain in triaxial test for rape (confining pressure = 400 kPa); (a) a comparison between the numerical and experimental results, (b) schematically isotropic–kinematic hardening behavior.

The material parameters calibrated for the loose Baskarp sand are as follows a1 ¼ 1:0 e  6 b1 ¼ 12:3

a2 ¼ 0:05 b2 ¼ 3150:0

a3 ¼ 20:0 b3 ¼ 1500:0

c1 ¼ 10:0

c2 ¼ 2650:0

c3 ¼ 2:6

a4 ¼ 0:06

The results for loose sand are shown in Figs. 18 and 19. In Fig. 18(a), the normalized stress, i.e. the axial stress–radial stress/radial stress is presented with axial strain. The relevant yield surfaces corresponding to the kinematic hardening behavior are shown in Fig. 18(b). In Fig. 19, the variation of volumetric strain with axial strain is presented together with the relevant yield surfaces corresponding to the kinematic hardening behavior. It must be noted that the cap plasticity model presented here is based on an associated flow rule, and the associated flow rule cannot address the contraction response of the loose sand once the stress-point

648

H. DorMohammadi, A.R. Khoei / International Journal of Solids and Structures 45 (2008) 631–656

a

3.5

3

Radial Strain (%)

2.5

2

1.5

1 Three-invariant isotropic-kinematic cap model

0.5

Experiment (Kolymbas and Wu,1990)

0

0

0.5

1

1.5

2

2.5

3

3.5

4

3

3.5

4

Axial Strain (%) 3.5

b

3

Radial Strain (%)

2.5

2

1.5

1

0.5

0

0

0.5

1

1.5

2

2.5

Axial Strain (%) Fig. 15. Radial strain versus axial strain in isostatic compression test for synthetic granulate; (a) a comparison between the numerical and experimental results, (b) schematically isotropic hardening behavior.

locates on the wing part of the yield surface. In this case, a non-associated flow rule is necessary. Thus, it is assumed that the stress-point locates on the cap part of the yield surface for the loose Baskarp sand. The results for the dense sand are shown in Figs. 20 and 21 with the corresponding material parameters given as follows a1 ¼ 0:9 e  6

a2 ¼ 15:0

a3 ¼ 4:0

b1 ¼ 125:0

b2 ¼ 60:1

b3 ¼ 260:0

c1 ¼ 35:0

c2 ¼ 103:0

c3 ¼ 44:9

a4 ¼ 1:0

The variation of normalized stress, i.e. the axial stress–radial stress/radial stress, is presented with axial strain in Fig. 20(a). Also plotted in Fig. 21(a) is the variation of volumetric strain with axial strain. In Figs. 20(b) and 21(b), the relevant yield surfaces corresponding to the kinematic hardening behavior are shown. It is of interest to be highlighted from the results of experimental and numerical simulation that the examples of loose and dense sands have been used here to demonstrate the capability of proposed model even in the

H. DorMohammadi, A.R. Khoei / International Journal of Solids and Structures 45 (2008) 631–656

649

Stress Ratio

a 2

1

Three-invariant isotropic-kinematic cap model

Volume Strain (%)

Experiment (Kolymbas and Wu,1990)

0

-1

0

5

10

15

Axial Strain (%)

b Stress Ratio

3

2

Volume Strain (%)

1

0

-1

0

5

10

15

Axial Strain (%) Fig. 16. Stress ratio and volumetric strain versus axial strain in triaxial test for synthetic granulate (confining pressure = 100 kPa); (a) a comparison between the numerical and experimental results, (b) Schematically isotropic–kinematic hardening behavior.

case of non-associated material, however – it must be mentioned that an associated plasticity with combined isotropic–kinematic hardening cannot be replaced by a non-associated plasticity model. In addition, it can be observed from the curves that the dilation happens in axial strain of 7% for the loose sand and axial strain of 2% for the dense sand, so it results in different material parameters for these two sands. These two examples demonstrate the capability of the proposed plasticity model to represent the combined isotropic and kinematic hardening behavior together with the dilatancy in standard triaxial tests for the loose and dense sands. 6. Conclusion In the present paper, a three-invariant cap plasticity model was developed based on an isotropic–kinematic hardening rule and associated plasticity for granular materials. Two material functions were introduced to control the size and movement of the yield surface, which are decomposed into the isotropic and kinematic parts as functions of the hardening parameters. The kinematic hardening parameter was defined to indicate

650

H. DorMohammadi, A.R. Khoei / International Journal of Solids and Structures 45 (2008) 631–656

Stress Ratio

a 2

1

Three-invariant isotropic-kinematic cap model

Volume Strain (%)

Experiment (Kolymbas and Wu,1990)

0

-1

0

5

10

15

Axial Strain (%)

b Stress Ratio

3

2

Volume Strain (%)

1

0

-1

0

5

10

15

Axial Strain (%) Fig. 17. Stress ratio and volumetric strain versus axial strain in triaxial test for synthetic granulate (confining pressure = 200 kPa); (a) a comparison between the numerical and experimental results, (b) schematically isotropic–kinematic hardening behavior.

the movement of the yield surface in the direction of J1-axis and in the direction of perpendicular to J1-axis. The constitutive elasto-plastic matrix and its components were derived by using the definition of yield surface, material functions and non-linear elastic behavior, as function of hardening parameters. The calibration procedure for three-invariant isotropic–kinematic cap plasticity was demonstrated based on a series of standard isostatic and triaxial tests. Finally, the applicability of the proposed model was illustrated in modeling of experiments on several granular materials, including: wheat, rape, synthetic granulate, loose sand and dense sand. Different aspects of material behavior, including: isotropic and kinematic hardening and dilation were investigated. The variation of the radial strain with axial strain was presented for the hydrostatic pressure test. It was shown how the yield surface grows isotropically by increasing the hydrostatic pressure. Furthermore, the variations of axial stress and volumetric strain were obtained with the axial strain at different confining pressures corresponding to the set of triaxial tests to represent the isotropic and kinematic hardening behavior

H. DorMohammadi, A.R. Khoei / International Journal of Solids and Structures 45 (2008) 631–656

a

651

90 80

Stress Difference (%)

70 60 50 40 30 20

Three-invariant isotropic-kinematic cap model

10

Experiment (Borup and Hedegaard,1995)

0

0

5

10

15

Axial Strain (%)

b

90 80

Stress Difference (%)

70 60 50 40 30 20 10 0

0

5

10

15

Axial Strain (%) Fig. 18. Normalized stress versus axial strain in triaxial test for loose Baskarp sand; (a) a comparison between the numerical and experimental results, (b) schematically isotropic–kinematic hardening behavior.

with dilatancy in granular materials. Remarkable agreements were observed between experimental and numerical results. In a later work, we will show how an extension and modification of proposed plasticity model can be effectively used to capture the softening and cyclic response of granular materials. Appendix A. Computation of material property matrix The object of the mathematical theory of plasticity is to provide a theoretical description of the relationship between stress and strain or more commonly, between increments of stress and increments of strain using the assumption that the material behaves plastically only after a certain limiting value has been exceeded. The elasto-plastic constitutive relation in its incremental form can be presented by dr = Depde, with dr denoting the incremental stress vector, de the incremental strain vector and Dep the constitutive elasto-plastic matrix. The yield surface F(r, a, j) = 0 determine the stress level at which the plastic deformation begins. The material property matrix Dep is defined as

652

H. DorMohammadi, A.R. Khoei / International Journal of Solids and Structures 45 (2008) 631–656

a

1.3 1.2 1.1 1

Volume Strain (%)

0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2

Three-invariant isotropic-kinematic cap model

0.1

Experiment (Borup and Hedegaard,1995)

0

0

5

10

15

Axial Strain (%)

b

1.3 1.2 1.1 1

Volume Strain (%)

0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0

0

5

10

15

Axial Strain (%) Fig. 19. Volumetric strain versus axial strain in triaxial test for loose Baskarp sand; (a) a comparison between the numerical and experimental results, (b) schematically isotropic–kinematic hardening behavior.

Dep ¼ De 

nT DTe De ng H þ nT D e ng

ðA:1Þ

where n = oF/or and ng = oQ/or are the normal vector to the yield and potential plastic surfaces, respectively. In present study, an associated flow rule is used so F = Q. In Eq. (A.1), H is the hardening plastic modulus defined as H¼

oF dl ol dk

ðA:2Þ

where dk is the plastic multiplier and l ¼ ðep ; epv Þ. In order to derive the constitutive elasto-plastic matrix and its components, we need to calculate De, n, ng and H in Eq. (A.1). The normal vector to the yield surface is determined by

H. DorMohammadi, A.R. Khoei / International Journal of Solids and Structures 45 (2008) 631–656

a

653

90 80

Stress Difference (%)

70 60 50 40 30 Three-invariant isotropic-kinematic cap model

20 10 0

Experiment (Ibsen and Jakobsen,1996)

0

2

4

6

8

10

8

10

Axial Strain (%)

b

90 80

Stress Difference (%)

70 60 50 40 30 20 10 0

0

2

4

6

Axial Strain (%) Fig. 20. Normalized stress versus axial strain in triaxial test for dense Lund sand; (a) a comparison between the numerical and experimental results, (b) schematically isotropic–kinematic hardening behavior.

2=3

n¼

oF oF oJ 1 oF oJ 2D oF oJ ¼ þ þ 2=3 3D or oJ 1 or oJ 2D or oJ 3D or

ðA:3Þ

where  2 oF / fd /2 f 2 2J 21 J 1a ¼ 2J 1 d þ d3 d3 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi oJ 1 /h fh /h fh ð2=3Þð2J 1a J 1 þ J 21a Þ     oF 2 fd /2 ofd ofd ¼ þ 2J 21 2 d2  2/2d fd oJ 2D 3 oJ 2D /h fh oJ 2D ! ! 2 oF ofd ofd 2 2 fd /d ¼ w þ 2J 1 2 2  2/d fd 2=3 2=3 2=3 /h fh oJ 3D oJ 3D oJ 3D

ðA:4Þ

ðA:5Þ ðA:6Þ

654

H. DorMohammadi, A.R. Khoei / International Journal of Solids and Structures 45 (2008) 631–656

a

1.5

1

Volume Strain (%)

0.5

0

-0.5

-1 Three-invariant isotropic-kinematic cap model

-1.5

Experiment (Ibsen and Jakobsen,1996)

-2

0

2

4

6

8

6

8

Axial Strain (%)

b

1.5

1

Volume Strain (%)

0.5

0

-0.5

-1

-1.5

-2

0

2

4

Axial Strain (%) Fig. 21. Volumetric strain versus axial strain in triaxial test for dense Lund sand; (a) a comparison between the numerical and experimental results, (b) schematically isotropic–kinematic hardening behavior.

where pffiffiffi pffiffiffi ofd ða1  a2 Þðcos x  3 sin xÞ 3J 3D ða1  a2 Þð sin x  3 cos xÞ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi q pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi þ ¼ oJ 2D 3 8/d J 22D ðJ 27 2 4/d 3J 2D ðJ 2Da þ J ra Þ 2Da þ J ra Þð1  4 J 3D =J 2D Þ pffiffiffi 1=3 ofd 3J 3D ða1  a2 Þð sin x  3 cos xÞ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼  ffi 2=3 8/d oJ 3D þ J Þ 1  27 J 2 =J 3 ðJ 2Da

ra

4

3D

ðA:7Þ

ðA:8Þ

2D

According to the definition of stress invariants, i.e. J1 = rii, J 2D ¼ 12 sij sij and J3D = det(sij) in which the devi2=3 atoric stress tensor defined as sij ¼ rij  13 dij rmm , the vectors oJ1/or, oJ2D/or and oJ 3D =or, defined in (A.3), can be calculated as

H. DorMohammadi, A.R. Khoei / International Journal of Solids and Structures 45 (2008) 631–656

oJ 1 T ¼ f1; 1; 0; 1g or oJ 2D T ¼ fsx ; sy ; 2sxy ; sz g or 2=3 oT oJ 3D 2 1=3 n 2 ¼ J 3D sx þ s2xy ; s2y þ s2xy ; 2sxy ðsx þ sy Þ; s2z 3 or

655

ðA:9Þ ðA:10Þ ðA:11Þ

The hardening plastic modulus H can be determined by substituting the yield surface (5) into definition (A.2) as (Khoei, 2005)     oF ofh dep oF ofd dep oF ofh depv oF ofd depv þ þ H¼ þ ðA:12Þ ofh oep dk ofd oep dk ofh oepv dk ofd oepv dk where

! 2 2 oF 2 /d fd ¼ 2J 1 ofh /2h fh3 ! 2 oF 2 /d fd ¼ 2J 1  2/2d fd ofd /2h fh2 ofh ¼ ðb2 b3 expðb3 epv ÞÞdðepv Þ oepv ofd ¼ ðc2 c3 expðc3 epv ÞÞdðepv Þ oepv   ofh 1 ðJ 1 þ J 1a Þ oa1 oa2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi q ¼  þ 2 /h ð2J J þ J 2 Þ oep oep oep 1a 1 1a pffiffiffi   pffiffiffiffiffiffiffiffi   J 2D ðcos x  3 sin xÞ oa1 oa2 ofd ða2  a1 Þ oa2 oa1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼   þ oep 3/d J 2Da þ J ra oep oep oep oep 2/d 3ðJ 2Da þ J ra Þ depv oF ¼ dk oJ 1 dep oF 1 oF  ¼ m or 3 oJ 1 dk

ðA:13Þ ðA:14Þ ðA:15Þ ðA:16Þ ðA:17Þ

ðA:18Þ ðA:19Þ ðA:20Þ

where m is the unit vector m = {1 1 1}T. References Bigoni, D., Piccolroaz, A., 2004. Yield criteria for quasibrittle and frictional materials. Int. J. Solids Struct. 41, 2855–2878. Borup, M., Hedegaard, J., 1995. Data Report 9403, Baskarp Sand No. 15, Geotechnical Engineering Group, Aalborg University, Denmark. Brandt, J., Nilsson, L., 1999. A constitutive model for compaction of granular media, with account for deformation induced anisotropy. Mech. Cohes. Frict. Mater. 4, 391–418. Brown, S., Abou-Chedid, G., 1994. Yield behavior of metal powder assemblages. J. Mech. Phys. Solids. 42, 383–398. Chang, C.S., Hicher, P.Y., 2005. An elasto-plastic model for granular materials with microstructural consideration. Int. J. Solids Struct. 42, 4258–4277. Chen, W.F., Baladi, G.Y., 1985. Soil Plasticity, Theory and Implementation. Elsevier, Amsterdam. Collins, I.F., 2005. Elastic/plastic models for soils and sands. Int. J. Mech. Sci. 47, 493–508. Desai, C.S., Siriwardane, H.J., 1984. Constitutive Laws for Engineering Materials, with Emphasis on Geologic Materials. Prentice-Hall, New Jersey. Doraivelu, S.M., Gegel, H.L., Gunasekera, J.S., Malas, J.C., Morgan, J.T., Thomas, J.F., 1984. A new yield function for compressible P/ M materials. Int. J. Mech. Sci. 26, 527–535. Dorris, J.F., Nemat-Nasser, S., 1982. A plasticity model for flow of granular materials under triaxial stress states. Int. J. Solids Struct. 18, 497–531. Foster, C.D., Regueiro, R.A., Fossum, A.F., Borja, R.I., 2005. Implicit numerical integration of a three-invariant, isotropic/kinematic hardening cap plasticity model for geomaterials. Comput. Methods Appl. Mech. Eng. 194, 5109–5138.

656

H. DorMohammadi, A.R. Khoei / International Journal of Solids and Structures 45 (2008) 631–656

Gudehus, G., 1996. A comprehensive constitutive equation for granular materials. Soils Foundations 36, 1–12. Ibsen, L.B., Jakobsen, F.R., 1996. Data Report, Lund Sand No. 0, Part 1 and Part 2, Geotechnical Engineering Group, Aalborg University, Denmark. Jenkins, J.T., Strack, O., 1993. Mean-field inelastic behavior of random array of identical spheres. Mech. Mater. 16, 25–33. Khoei, A.R., 2005. Computational Plasticity in Powder Forming Processes. Elsevier, UK. Khoei, A.R., Azami, A.R., 2005. A single cone-cap plasticity with an isotropic hardening rule for powder materials. Int. J. Mech. Sci. 47, 94–109. Khoei, A.R., Azami, A.R., Haeri, S.M., 2004. Implementation of plasticity based models in dynamic analysis of earth and rockfill dams; a comparison of Pastor-Zienkiewicz and cap models. Comput. Geotech. 31, 385–410. Khoei, A.R., Bakhshiani, A., 2004. A hypoelasto-plastic finite strain simulation of powder compaction processes with density dependent endochronic model. Int. J. Solids Struct. 41, 6081–6110. Khoei, A.R., Bakhshiani, A., Mofid, M., 2003. An implicit algorithm for hypoelasto-plastic and hypoelasto-viscoplastic endochronic in finite strain isotropic–kinematic hardening model. Int. J. Solids Struct. 40, 3393–3423. Khoei, A.R., DorMohammadi, H., 2007. A three-invariant cap plasticity with isotropic–kinematic hardening rule for powder materials: model assessment and parameter calibration. Comp. Mater. Sci., in press. Khoei, A.R., DorMohammadi, H., Azami, A.R., 2007. A three-invariant cap plasticity model with kinematic hardening rule for powder materials. J. Mater. Proc. Tech. 187, 680–684. Krenk, S., 2000. Characteristic state plasticity for granular materials, Part I: basic theory. Int. J. Solids Struct. 37, 6343–6360. Krenk, S., Ahadi, A., 2000. Characteristic state plasticity for granular materials, Part II: model calibration and results. Int. J. Solids Struct. 37, 6361–6380. Kolymbas, D., Wu, W., 1990. Recent results of triaxial tests with granular materials. Powder Tech. 60, 99–119. Lade, P.V., Kim, M.K., 1995. Single hardening constitutive model for soil, rock and concrete. Int. J. Solids Struct. 32, 1963–1978. Lewis, R.W., Khoei, A.R., 2001. A plasticity model for metal powder forming processes. Int. J. Plasticity 17, 1652–1692. Lewis, R.W., Schrefler, B.A., 1998. The Finite Element Method in the Static and Dynamic Deformation and Consolidation of Porous Media. Wiley, England. McDowell, G.R., Bolton, M.D., 1998. On the micromechanics of crushable aggregates. Geotechnique 48, 667–679. Oliver, J., Oller, S., Cante, J.C., 1996. A plasticity model for simulation of industrial powder compaction processes. Int. J. Solids Struct. 33, 3161–3178. Wan, R.G., Guo, P.J., 1998. A simple constitutive model for granular soils: modified stress-dilatancy approach. Comput. Geotech. 22, 109–133. Watson, T.J., Wert, J.A., 1993. On the development of constitutive relations for metallic powders. Metall. Trans. 24, 1993–2071.