A study of localized deformation pattern in granular media

Comput. Methods Appl. Mech. Engrg. 193 (2004) 2719–2743 www.elsevier.com/locate/cma

A study of localized deformation pattern in granular media Konrad N€ ubel b

a,*

, Wenxiong Huang

b

a Ed. Z€ublin AG, Albstadtweg 3, D-70567 Stuttgart, Germany School of Engineering, The University of Newcastle, NSW 2308, Australia

Received 27 March 2003; received in revised form 6 October 2003; accepted 7 October 2003

Abstract Shear localization in granular materials are studied numerically and compared with experimental results. General agreements are shown for the localization pattern in granular body behind a retaining wall in active and passive limit state and under a strip footing. Quantitative agreement is also shown for the limit load in strip footing. Numerical simulations are based on a micro-polar (Cosserat) continuum approach. An extended hypoplastic model is used for constitutive description of granular media, which includes rotational degrees of freedom for material particles and takes into account the effects of grain size and inter-granular friction. For a more realistic modeling of material response, fluctuation in state variables are introduced based on physical argument. It is shown through numerical examples that realistic localized deformation pattern can be initiated with a fluctuation in initial density. Ó 2004 Elsevier B.V. All rights reserved. Keywords: Shear localization; Granular materials; Cosserat continuum; Hypoplasticity; State fluctuation

1. Introduction The development of finite element methods and other computer based numerical schemes has reached a stage where solutions for many different geotechnical problems are possible. In these continuum mechanicsbased numerical processes, constitutive modeling of geomaterials plays a key role. A recently developed framework for constitutive modeling of granular materials, known as hypoplasticity [28,29], as an alternative approach of elasto-plasticity, has demonstrated the capability to give realistic description of behavior of granular materials with relatively simple formulations. A practical and comprehensive hypoplastic model was proposed by Gudehus [15] and Bauer [1] which contains stresses and void ratio as state variables and can capture many important features of the behavior of cohesionless granular materials for a wide range of pressure and density with only eight state-independent material parameters. However, for conventional safety assessment of geotechnical structures, the gap between service states and limit states has not yet been closed. Granular materials have a discrete nature, which is characterized by their grains possessing independent degrees of freedom in translation and rotation. In phenomenological *

Corresponding author. Fax: +49-711-7883-650. E-mail address: [email protected] (K. N€ ubel).

0045-7825/$ - see front matter Ó 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2003.10.020

2720

K. N€ubel, W. Huang / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2719–2743

modeling based on conventional (Boltzmann) continuum theory, the materials are idealized as continuum media whose particles possess translation degrees of freedom only. This assumption can lead to rather good approximations in modeling the deformation of granular materials, in which sharp strain gradient variation or discontinuity are not present. But a class of deformation pattern, namely the localized deformation pattern, can often be observed in granular bodies when limit states are approached or inhomogeneous deformation has been fully developed. In a localized deformation pattern, intensive shear deformation develops in narrow zones commonly called shear bands, which have a finite thickness depending mainly on the grain size [39]. Very sharp strain gradient variation accompanied by pronounced grain rotation and dilation can be observed across the localized zones [51]. Conventional continuum theory meets the limitation in modeling shear localization phenomenon. Due to the lack of an intrinsic length, conventional constitutive models can not scale the thickness of shear bands. Consequently, numerical solutions of the rate boundary value problems exhibit a spurious mesh sensitivity [7]. In numerical simulation of homogeneous deformations leading to a limit state, numerical procedure becomes unstable after the bifurcation point has been reached. A promising way to surpass the limitation of the conventional continuum is to introduce the rotational degrees of freedom to material particles, which leads to the theory of micro-polar (Cosserat) continuum. Along with the rotational degrees of freedom, Cosserat continuum is associated with couple stresses in addition to traditional stresses, which enables an intrinsic length to enter constitutive models and the nonlocal effects to be taken into account. The capability of Cosserat continuum models for capturing the phenomenon of shear localization in granular materials has been demonstrated, for instance, by M€ uhlhaus [33], de Borst [7], Tejchman and Wu [49], Steinmann [42], Ehlers and Volk [10], who used an elasto-plastic approach, and by Tejchman [47], Tejchman and Gudehus [48], Bauer and Huang [3], N€ ubel [35], Huang and Bauer [24], who adopted a hypoplastic approach. The present paper is focused on studying the localized deformation pattern in granular media under various boundary conditions. Experimental results for localization pattern behind a retaining wall in active and passive limit state and in granular foundation for strip footing [35] are presented and compared with numerical results. Different visualization techniques are used in these experiments. Numerical simulation is based on a Cosserat continuum approach. The mechanical behavior of the granular material is described with an extended hypoplastic model, which typically includes stresses, couple stresses and void ratio as state variables, the mean grain diameter as an intrinsic length and an additional material parameter related to the polar effect. Different formulations for Cosserat continuum extensions of the hypoplastic model for conventional continuum were proposed by Tejchman [47], Tejchman and Gudehus [48], Bauer and Huang [4], Huang et al. [26]. All of these models have been shown to be able to capture shear localization phenomenon with a width independent of element size. The formulation proposed by Huang et al. [26] is used in this study, which has the following advantages over the earlier formulations: (1) Stationary states or the so-called critical states are consistently embedded in the model for the evolution of the stress, the couple stress, and the void ratio. Under monotonic shearing, the void ratio tends to the critical void ratio independent of initial state. (2) The limit stress and the limit couple stress are related in a rational manner, which allows a physical interpretation for the additional polar parameter. (3) The thickness of the localized zones is scaled by a characteristic length which is proportional to the mean grain diameter and also related to the property of inter-granular friction. Numerical tests with the micro-polar hypoplastic model showed that a localized deformation mode is not guaranteed to emerge in the specimen with a homogeneous initial state under an uniform straining or stressing condition. This is confirmed by a recent analysis [25], which shows that the bifurcation condition

K. N€ubel, W. Huang / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2719–2743

2721

will never be met in homogeneous biaxial compression tests modeled with the micro-polar hypoplastic model, whereas the bifurcation condition is met before the peak stress state in the conventional continuum hypoplastic description. While localized deformation mode is almost always observed in experiments, this result based on Cosserat continuum approach suggests that shear localization develops as a consequence of deformation inhomogeneity and strain softening, rather than as a result of the loss of uniqueness in mathematical description of a boundary value problem. From a physical point of view, local inhomogeneity is inevitable in a granular material which has a micro-structure. On the micro-scale, void ratio varies significantly from point to point, while a macroscopic homogeneous condition is maintained. Shahinpoor [41] showed that even within a granular body packed with equal-sized spheres, the void ratio or porosity will not be uniform throughout, but rather have a geometrical distribution. It has also been shown (e.g. [36,54]), theoretically and experimentally, that the void ratio alone is not sufficient to describe the mechanical behavior of granular materials. Frequency distribution of void ratio is therefore suggested to be added as another important factor to evaluate the mechanical behavior of granular materials. In this work, fluctuations are introduced to the void ratio in numerical simulations by the use of frequency distributions. Certainly, fluctuations can also be introduced to other state variables without difficulty. However, with a frequency distribution in void ratio, fluctuation in stresses and couple stresses will be induced immediately in a deformation process and therefore it is enough to prescribe a fluctuation only in initial void ratio and realistic localization pattern can be initiated in finite element simulations.

2. A description of constitutive model and FE implementation 2.1. Kinematics for Cosserat continuum Cosserat continuum kinematics differs from conventional continuum theory by introducing additional rotational degrees of freedom to material points in the three-dimensional Euclidean space (e.g. [11,31]). _ which defines the Consequently, the motion of continuum particles are characterized by the velocity field u, so-called macro-motion, and the rate of Cosserat rotation field w_ c , which defines the so-called micro-motion [11]. The rate of deformation can be described with the velocity gradient L and the Cosserat spin tensor Wc . The general expression of the velocity gradient represents the rate of deformation resulting from the macromotion: _ L ¼ ru;

Lij ¼

ou_ i ; oxj

ð1Þ

where xj is the position vector of the current configuration. The Cosserat spin tensor, which corresponds to the micro-motion, is defined as Wc ¼
w_ c  ;

Wijc ¼

kij w_ ck ;

ð2Þ

where  represents the permutation tensor in Euclidean space. With respect to the identity:
kij
klm ¼ dil djm
dim djl ;

ð3Þ

where dij is known as the Kronecker delta, the following relation is valid: w_ c ¼
12 : Wc ;

w_ ck ¼
12
klm Wlmc :

ð4Þ c

The deformation of a Cosserat continuum can be measured with a generalized stretching tensor D and a micro-curvature rate K. The generalized stretching tensor can be expressed as a sum of the symmetric part and the skew-symmetric part of the velocity gradient and the Cosserat spin tensor:

2722

K. N€ubel, W. Huang / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2719–2743

D c ¼ D þ W
Wc ¼ L
Wc ;

Dcij ¼

ou_ i þ
kij w_ ck : oxj

ð5Þ

The micro-curvature rate is defined as the gradient of the rate of Cosserat rotation: K ¼ rw_ c ;

Kij ¼

ow_ ci : oxj

ð6Þ

Analogous to the relation between Cosserat spin tensor and the rate of Cosserat rotation, the rate of macrorotation can be determined from the macro-spin tensor: w_ ¼
12
: W;

w_ k ¼
12
klm Wlm :

ð7Þ

Note that the rate of the micro-rotation w_ c is generally not in compliance with the rate of the macrorotation w_ and thus the generalized stretching Dc is non-symmetric. In the case where the rate of microrotation equals to the rate of macro-rotation, the generalized stretching tensor Dc coincides with the standard stretching tensor D for conventional continuum. 2.2. Equilibrium In a Cosserat continuum, static quantities associated with the generalized stretching Dc and the microcurvature rate K are the Cauchy stress tensor T and the couple stress tensor M, which are defined in the current configuration of a continuum body. With existing of the couple stresses, the static equilibrium condition can be expressed by 9 oTij > > rTþb¼0 þ bi ¼ 0; = oxj ð8Þ oMij > rM
:Tþc¼0

ijk Tjk þ ci ¼ 0; > ; oxj where b denotes the body force per unit volume and c the body couple per unit volume. Note that due to the existence of couple stresses, the stress tensor is generally non-symmetric, too. The equilibrium of a Cosserat continuum body with a volume X can be described alternatively in the weak form as Z Z Z Z ðdDc : T þ dK : MÞ dV ¼ du_  t ds þ dw_ c  m ds þ ½du_  b þ dw_ c  c dV : ð9Þ oXt

X c

oXm

X

c

_ dw_ , dD and dK represent a virtual kinematic field attributed to the Cosserat continuum body Herein, du, in consideration, and oXt and oXm denote the parts of the boundary of the granular body where the surface traction t and the surface couple m are prescribed, respectively. On these boundaries the following relations are valid:
t¼Tn ti ¼ Tij nj on Xt ; ð10Þ m ¼ M  n mi ¼ Mij nj on Xm ; with n being the unit normal vector for the boundaries. 2.3. A micro-polar hypoplastic constitutive law The micro-polar extension of the hypoplastic model is developed based on the model proposed by Gudehus [15] and Bauer [1]. The model was slightly modified by von Wolffersdorff [52] and Bauer [2] to incorporate the Matsuoka and Nakai and other criteria as predefined limit conditions at the critical

K. N€ubel, W. Huang / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2719–2743

2723

states. The extended model incorporates couple stress M as an additional state variable. The historydependent behavior of the granular media is now described with three state variables, namely, the stress tensor T, the couple stress tensor M and the void ratio e. The tensor-valued functions have been extended to define an interaction in evolution of stresses and couple stresses while asymptotic behavior and three attractors [17] are retained. Mean grain diameter d50 is taken as an intrinsic length for the material description [34]. The evolution of the state variables reads:  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 9  1 > c d > ^ ^ ^ > T :¼ fb fe LðT; M; D ; K Þ þ fd NðTÞ ^a2 kDc k2 þ d^2 kKd k2 > > ^ ^ T:T  = q ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi  : ð11Þ 1 ^ M; ^ Dc ; Kd Þ þ fd Nc ðMÞ ^ ^a2 kDc k2 þ d^2 kKd k2 > Lc ðT; M :¼ d50 fb fe > > ^ :T ^ > T > ; e_ :¼ ð1 þ eÞtr D ^ :¼ T=tr T denotes the normalized stress tensor, and M ^ :¼ M=ðd50 tr TÞ the normalized couple Herein T d stress, K :¼ d50 K is the scaled rate of micro-curvature. The factors fb , fe and fd are functions of the void ratio e and the mean pressure ps :¼
tr T=3, and comprise the dependence on the mean pressure (barotropy) and density (pyknotropy):  a e
ed ; ð12Þ fd :¼ ec
ed fe :¼

 e b c

e

ð13Þ

;

and  fb :¼

ei0 ec0

b

hs 1 þ e i n ei




tr T hs

1
n 

  pffiffiffi ei0
ed0 a
1 3þa
a 3 : ec0
ed0 2

ð14Þ

The void ratios ei , ed and ec , which characterize respectively the loosest state, densest state and critical state of a granular medium on macroscopic level, decrease with increasing mean pressure. The following relation proposed by Bauer [1] for ei based on experiments and postulated for ed and ec by Gudehus [15] are used:   n  ec ed ei
tr T ¼ ¼ :¼ exp
: ð15Þ hs ec0 ed0 ei0 Herein hs , called granular hardness, and n are two material parameters served to fit the compression curves, ei0 , ec0 and ed0 are the maximum, the critical and the minimum void ratio for a stress-free state Eq. (15) can be shown graphically by the so-called granular phase diagram (Fig. 1, [16]) With respect to the pressuredependence of the maximum, minimum and critical void ratio, a pressure-dependent density index Id for granular media introduced by Cudmani and Osinov [6] is adopted, reading: ec
e Id :¼ : ð16Þ ec
ed Id ¼ 1 represents the densest state while Id ¼ 0 in the critical state. Collapsible zones with Id < 0 can also occur but are note explicitly marked in the contour plot (refer to Fig. 1). ^ a and d^ in Eq. (11) are related to the limit stress and limit couple stress at critical states. By assuming a constant value for ^ a, the Drucker–Prager limit stress condition is embedded in the model. With the following interpolation for ^ a, the Matsuoka-Nakai limit stress condition is incorporated [3]:

2724

K. N€ubel, W. Huang / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2719–2743

0

1 Id

e ei

ec

ed

ln (ps / hs) Fig. 1. Decrease of the maximum void ratio ei , the critical void ratio ec and the minimum void ratio ed with an increase of the mean pressure ps .

3 2sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi pffiffiffi ^ s k2 þ ð 6=2ÞkT ^ s k3 cosð3hÞ sin uc 4 ð8=3Þ
3kT ^ s k5; pffiffiffiffiffiffiffiffi s ^ aðhÞ :¼
kT ^ k cosð3hÞ 3
sin uc 1 þ 3=2kT

ð17Þ

^ s :¼ T ^ s
1 1 ¼ Ts
1 1 is the deviator of the normalized symmetric part of stress T ^ s , and h is the where T 3 3 tr T Lode angle defined in the principal stress space of the symmetric part of the stress tensor [23] as ^ s  T ^ s  T ^ s Þ pffiffiffi tr ðT cos 3h :¼
6 : 3=2 s s ^ :T ^  ½T

ð18Þ

The constitutive tensors L, Lc , N and Nc for micro-polar model (11) have the following representations: ) ^ : Dc þ M ^ : Kd ÞT; ^ ^ þT ^ L :¼ ^ a2 Dc þ ðT N :¼ T : ð19Þ ^ : Dc þ M ^ : Kd ÞM; ^ ^ Lc :¼ d^2 Kd þ ðT Nc :¼ 2M ^
1 1 is the normalized stress deviator. In the case of M ¼ 0 and K ¼ 0 in Eq. (11), the non^  :¼ T Herein T 3 polar formula of the hypoplastic model by Gudehus and Bauer is recovered. On the other hand, if condition kDc k ¼ 0 is fulfilled, the evolution equation for the couple stress tensor i  1 h^2 d ^ ^ d kÞM ^ : Kd þ 2fd dkK M ¼ d50 fb fe ð20Þ d K þ ðM ^ :T ^ T resembles the non-polar evolution equation for the deviatoric stress T , i.e. i  1 h  ^ : ^ : D þ 2fd ^ ^ T ¼ fb fe akDkÞT aD þ ð T ^:T ^ T

ð21Þ

Therefore, stationary flow is embedded also for the couple stress so that the response of the couple stress to a pure increase of micro-curvature is hypoplastic (i.e., the couple stress tends to a limit value asymptotically as micro-curvature keeps increasing).

K. N€ubel, W. Huang / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2719–2743

2725 

 A stationary state is reached as changes in state variables vanish simultaneously, expressed as T ¼ 0, M ¼ 0 and e_ ¼ 0, for continuing deformation. The present micro-polar hypoplastic model is characterized by the following limit conditions on the stress, the couple stress, and the void ratio at stationary states [26]:

^ 2 kMk ^ 2 kTk þ ¼ 1 and ^ a2 d^2

e ¼ ec ;

ð22Þ

which are independent of the initial state. Eq. (22) combines the limit value of normalized deviatoric stress ^ ^a can therefore be and the limit value of the normalized couple stress by means of parameters ^a and d. ^ interpreted as the inter-granular frictional resistance to grain sliding and d the inter-granular frictional resistance to grain rotation. The ratio between inter-granular frictional resistance to grain rotation to sliding ^ g :¼ dðhÞ=^ aðhÞ

ð23Þ

is assumed to be constant, independent of the Lode parameter cosð3hÞ. This ratio defines a frictional parameter related to grain shape and roughness. The original non-polar hypoplastic model contains eight parameters, namely, uc , hs , n, ed0 , ec0 , ei0 , a, b. They can be easily determined from simple index tests and element tests in laboratory by the procedure proposed by Herle and Gudehus [21]. The present micro-polar extension comprises only two additional parameters, i.e. the frictional parameter g and the mean grain diameter d50 . As the parameter g is related to the limit value of normalized couple stress, which has no proposed means to measure, a back analysis can be performed for determination of g as proposed by Huang et al. [26]. 2.4. Numerical implementation The proposed micro-polar hypoplastic constitutive model has been implemented into the commercial finite element program ABAQUS for numerical analysis. For completeness of the paper, a brief description of numerical aspects is given here. Details can be found in [23,24]. For plane strain problems, a four-node compatible displacement element with translational and rotational degrees of freedom at each node is formulated. Finite element equations for the discrete system are obtained based on the virtual work equation, where geometric non-linearity is taken into account using a Lagrangian description. Four Gauss points are used for numerical quadrature of element stiffness and for integration of constitutive equations. In order to improve the poor performance of the element at nearly incompressible limit, the selective reduced integration technique is incorporated by adopting the so-called B-bar method of Hughes [27]. The strongly non-linear finite element equations are solved incrementally with the Newton iteration method. A fictitious time scale is introduced for loading increments. Within each time increment, an implicit algorithm known as the h-method is used in combination with an automatic sub-time stepping scheme for integration of constitutive equations. The technique allows a global convergence to be achieved for the same error tolerance with sub-time increments initiated in areas of high strain gradient. In modeling a practical problem, an initial state needs to be defined. On boundaries, either surface couple or rotation needs to be prescribed in addition to a prescription of surface traction or translation displacement.

3. State fluctuation In this study, only quasi-static processes are considered. Thus a state of a granular body refers to an equilibrium state only. In the present micro-polar hypoplastic modeling, material parameters and state

2726

K. N€ubel, W. Huang / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2719–2743

variables are strictly separated. The void ratio e, the (effective) stress tensor T and the couple stress tensor M are the only variables characterizing the state of a granular material. It is experimentally evident that void ratio is an average quantity measured within a certain volume. Void ratio may fluctuate from subvolume to sub-volume within a system while macroscopic homogeneity is maintained. Void ratio fluctuation can even be observed in a granular package comprised by equal-sized spheres [41]. Similarly, stresses and couple stresses may fluctuate in a homogeneous granular system with void ratio fluctuation, while global equilibrium holds. The state fluctuation may significantly influence the global response of a granular body, it is therefore considered in this work. 3.1. Fluctuation from numerical processes Before introducing a frequency distribution to state variables, we present results from numerical modeling of biaxial compression tests with the constitutive model described in the foregoing section. Such tests are heuristic as plane strain conditions are imposed and boundary conditions are well defined. Experimental studies can be found in, for instance, Vardoulakis et al. [50], Tatsuoka et al. [44], Yoshida et al. [53], Desrues et al. [8] and Finno et al. [13]. For numerical modeling, we consider a dry granular cuboid with a height h ¼ 20 cm, a breadth b ¼ 8 cm and an unit depth. The specimen is laterally confined by a membrane with an external pressure ps ¼ 78:4 kPa and kept between two smooth opposing horizontal plates. In finite element calculations, initial homogeneous state is assumed with initial stresses being set to T22 ¼ T33 ¼ T11 ¼ 78:4 kPa. The finite element mesh is displayed in Fig. 2. We take the material parameters for Hostun sand (Table 1) determined by Herle [21] with respect to the non-polar hypoplastic model and assumed for the polar parameter. Numerical tests are performed for a loose sample with an initial pressure-adjusted density index Id ¼ 0:04, which corresponds to an initial void ratio very close to the critical void ratio, and for a dense sample with Id ¼ 0:87. It is found that, for the loose sample, homogeneous deformation is maintained up to the end of the calculation, and for the dense sample, shear localization occurs in the course of the loading process. This is reflected by the evolution curves of the deviatoric stress shown in Fig. 3. These curves are taken from an element in the finite element calculations (using the micro-polar continuum model) and compared with those obtained from integration of the corresponding non-polar hypoplastic equations

20 cm

T11 = const

u

8 cm Fig. 2. Finite element mesh for the biaxial compression.

K. N€ubel, W. Huang / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2719–2743

2727

Table 1 Material parameters for Hostun sand with respect to the micro-polar hypoplastic model uc ¼ 31° ed0 ¼ 0:61 ec0 ¼ 0:91 ei0 ¼ 1:09 d50 ¼ 0:31  10
3

hs ¼ 1000 [MPa] n ¼ 0:29 a ¼ 0:13 b ¼ 2:00 g ¼ 1:00

400

T11 - T22 [-]

non-polar

d

300

polar 200 c

100

polar and non-polar

0 0

2

4

6

ε22 [%]

8

10

Fig. 3. Calculated biaxial compression with micro-polar and non-polar hypolasticity and initially homogeneous state, dense (d) or loose (c).

under biaxial compression condition. Note that the latter (from direct integration of the non-polar model) are not structure responses. For the case with a loose initial-state, the two curves coincide completely (curve c). In contrast, for the case with a dense initial-state, the two curves branch after the peak in response to the development of shear localization in the specimen (curve d). Results from bifurcation analysis indicates that, shear bifurcation condition will not be encountered in a micro-polar hypoplastic continuum specimen under homogeneous deformation [25]. That is to say, that as long as the deformation remains completely homogeneous, localization could not initiate itself and the material response would be that of the homogeneous solution. Further examination shows that the spontaneous localization develops as a consequence of numerical error accumulation, which leads to a minute fluctuation in otherwise uniformly distributed field variables. Such numerical inhomogeneities create a random imperfection in elements. The results showed that the homogeneous deformation mode is rather sensitive to such imperfection when softening is more pronounced. 3.2. Frequency distribution of void ratio As has been demonstrated, fluctuations may play an important role for the macroscopic behavior of granular bodies. In particular, the evolution of shear localization is closely related to initial state-perturbations and therefore the micro-structure of the grain assembly. The observed width of shear bands is in the range of 5–20 grain diameters. (e.g. [39,51]). Large void ratio fluctuations have been observed at this scale (e.g. [37]). As the subdomains increase in volume and in the number of grains within a homogeneous system, the fluctuations can be expected to decrease. A state v, described by the state variables, scatters around its mean
v. The deviation of the mean depends on the size of the subdomain. Fluctuations of state variables can be described by probabilities. In the following fluctuation in void ratio is prescribed at the initial state in the following manner.

2728

K. N€ubel, W. Huang / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2719–2743

Suppose a probability density function pðeÞ is normalized as Z ei pðeÞ de ¼ 1:

ð24Þ

ed

The range of possible void ratio values is limited by the minimum void ratio ed and the maximum void ratio ei . In general the relation ei > ei holds, which means that some local areas may be in a looser state than the average loosest packing characterized by ei . Consequently void ratio values above ei can only be achieved by an increasing proportion of very loose subdomains. Obviously, the average void ratio ei ¼ ei cannot be obtained in a dry grain skeleton since the loose areas are only present because of a proportion of dense areas stabilizing the grain skeleton. The transmission of forces by force chains in the grain skeleton works out to be a complicated feature of grain assemblies and will not be discussed here. Due to the condition that void ratio fluctuations in the critical state are at a maximum, the probability density function at this state pðec Þ must correspond to an uniform distribution, so that pc ¼

1 : ei
ed

ð25Þ

The critical void ratio in this case can be obtained as Z ei pc e de ¼ 12ðei þ ed Þ: ec ¼

ð26Þ

ed

The relation is useful for determination of ei , since ec and ed of granular material can be easily determined [21]. As a reasonable interpolation between the densest state and the critical state, a probability density function proposed by Shahinpoor [41] is adopted, which reads pðeÞ ¼

k expð
keÞ : expð
ked Þ
expð
kei Þ

ð27Þ

The continuous probability density is of the negative exponential type skewed towards the population of denser cells (Figs. 4 and 5) and assumed for the two- and three-dimensional case. The constant k is dependent on the soil density and can be obtained from a given average void ratio e (in the following e refers always to the mean value hei) by the following equation:

relative frequency

0.04

0.03 Id = 0.53 (e = 0.75) 0.02

0.01

0.00 0.6

0.8

1.0 1.2 void ratio e

1.4

Fig. 4. Theoretical distribution of local void ratio e for Id ¼ 0:53.

K. N€ubel, W. Huang / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2719–2743

2729

relative frequency

0.04

0.03 Id = 0.20 (e = 0.85) 0.02

0.01

0.00 0.6

0.8

1.0 1.2 void ratio e

1.4

Fig. 5. Theoretical distribution of local void ratio e for Id ¼ 0:20.

e ¼ k
1 þ




ed expð
ked Þ
ei expð
kei Þ : expð
ked Þ
expð
kei Þ

ð28Þ

Since the probability distribution is dependent on the soil density, void ratio fluctuations are also correlated with the pressure-dependent density index Id . Examples for frequency distributions generated numerically for different densities (for Hostun sand Table 1) are shown in Figs. 4 and 5. Until now we have been concerned primarily with fluctuations without any typical length scale for a macroscopic volume element. A volume element comprises a particular number of grains forming a particular number of voids. The void ratio is smeared out over the volume element which means that an average is taken. We assume that in the critical state all void ratios appear with the same probability. Clearly, the larger the values of void ratio the larger the size of the lattice of grains forming the void. That means that a shear band must at least be able to incorporate a lattice of grains forming the maximum local void ratio ei . Thus it will be supposed in the following that the width of the shear band represents a maximum value for the fluctuation length. The random void ratios are calculated at points xn in a field. If the void ratios of two points (x1 ; x2 ) in a random field are influenced by each other they are considered to be correlated and termed as such. As there is no experimental evidence to the authorsÕ knowledge of a typical correlation length, we can assume that the random field is supposed to be uncorrelated. It is a technical simplification if the finite element mesh corresponds to the random field discretization, which we apply in the following when modeling shear localization. For the random field discretization we use a mesh size below 10d50 which is slightly smaller than the shear band width. 3.3. Evolution of scattered state variables The deviation of the void ratio from its mean value can be used as a measurement of the fluctuation: n 1X e2
e2 ; ð29Þ s2 :¼ n j¼1 j where n denotes the total number of points and j a particular point in the random field of void ratios. As we have just noted, the random field of void ratio will correspond to the finite element mesh, so that n also denotes the total number of elements and j a particular finite element.

2730

K. N€ubel, W. Huang / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2719–2743

Fluctuations in other state variables can be measured similarly. Since the stress T is a tensor, fluctuation in stress may be referred to the principal stress space. When dealing with non-symmetric stress states, it is convenient to split the stress in a symmetric- and a skew-symmetric part: Ts ¼ 12ðT þ TT Þ and

Ta ¼ 12ðT
TT Þ:

ð30Þ

We want to express fluctuations of the stress tensor normalized to the pressure level tr T. The norm of ^ ¼ T=tr T can then be expressed as (refer to [23]): T ^ 2 ¼ 1 tr 2 T ^ s k2 þ kT ^ a k2 ¼ n 2 þ n2 þ n2 : ^ s þ kT kTk 1 2 3 3

ð31Þ

Using the scalar quantities ni ¼ n1 ; n2 ; n3 , the fluctuations of the stress tensor can be defined as n 1X 2 n2ij
ni : w2i :¼ n j¼1

ð32Þ

The variation of state fluctuation in a loading process can be shown by evolution of the deviations of state variables or the configuration entropy. Consider a granular sample under plane strain compression with the boundary conditions D11 ¼ 0 and D22 ¼
1. (Fig. 6). The parameters for Hostun sand are used with initial pressure-adjusted densities Id ¼ 0:04 ðep0 ¼ 0:90Þ and Id ¼ 0:87 ðep0 ¼ 0:65Þ, respectively. Note that ep0 refers to the void ratio related to p ¼ 0. Numerical calculations are performed with and without fluctuation in initial void ratio. It can be observed that the mean void ratio decreases with increasing mean pressure ps for both homogeneous and inhomogeneous initial density states (Fig. 7). The fluctuation in void ratio (Fig. 8) decrease as the mean pressure increases. Fluctuations are also observed in stresses, though an uniform

8 cm

x2

x1 20 cm

e [-]

Fig. 6. Finite element mesh for the compression test.

0.85 0.80 0.75 0.70 0.65 0.60 0.55 0.50

ep0 = 0.90 (loose)

h1 i1

ep0 = 0.65 (dense) h2 10

100

i2 1000 ps [kPa]

10000

Fig. 7. Compression of a dense (1) and a loose (2) sample with homogeneous (h) and inhomogeneous (i) initial void ratio distribution.

K. N€ubel, W. Huang / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2719–2743

2731

0.16 loose

s [-]

0.12 0.08 dense

0.04 0.00 0

2

4

6

ε22 [%]

8

10

Fig. 8. Evolution of void ratio fluctuations s in the compression test.

wi [-]

0.004

w2/3, loose

0.002

w1, dense w1, loose

w2/3, dense

0.000 0

2

4

6

ε22 [%]

8

10

Fig. 9. Evolution of stress fluctuations wi in the compression test.

stress state is initially assigned. As we travel along the curves of deviatoric stress fluctuations w2 (Fig. 9), a rapid increase followed by a slow decrease are observed. Based on these results we are able to conclude that the granular body is not only densified but also homogenized by compression. It will be shown now, that under particular boundary conditions the irregularity in a granular body may also increase. We adopt the same biaxial compression test as in the example above (Fig. 2, dense Hostun sand, Id ¼ 0:87). Random fluctuations are attributed to initial void ratio in numerical simulations. Localized shear deformation and dilation are produced from numerical simulations. Fig. 10 shows a comparison between results from the calculations and from an experiment performed by Yoshida et al. [53] with the same initial and boundary conditions. The location of the shear band cannot be predicted a priori. Assuming the upper and lower platen to be ideally smooth, the shear band does not necessarily run through the center of the biaxial body (Fig. 10c), but may be reflected at the boundary (Fig. 10d). Some field quantities calculated for the biaxial compression are presented in Fig. 11. The evolution of some quantities and fluctuations are shown in Fig. 12. At the rather diffuse onset of localization (A), only ffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
¼ m231 þ m232 indicate the forthe volume strain rate tr D, and to lesser extent the mean couple stress m mation of shear bands in two main directions. At this stage the stress deviator T11
T22 is close to its peak. The overall void ratio indicates a slight expansion after a slight initial compression. At the peak (B) a lateral offset on the sides is not yet visible, and the increase of couple stress and volumetric strain rate indicates the

2732

K. N€ubel, W. Huang / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2719–2743

Fig. 10. (a) Localized zone in the biaxial test by Yoshida et al. [53]; (b) incremental evolution of localization in the experiment; (c) and (d) numerical simulation of the biaxial test.

evolution of several shear bands. Beyond the peak (C) the offset on both sides is marked. One shear band becomes dominant, and shearing and polarization are strongly localized in it. The behavior of the model corresponds well to an experiment by Tatsuoka et al. [44]. Interestingly, while the overall fluctuation in void ration s and stress wi increases slightly, a significant increase of the fluctuations inside the shear band is observed after the peak (Fig. 12). These results are in

K. N€ubel, W. Huang / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2719–2743

2733

Fig. 11. Calculated evolution of state during biaxial compression: (a) deformed mesh, (b) density index Id , (c) volume strain rate tr D, (d) mean couple stress m.

accordance with the experimental observations by Oda and Kazama [37], who reported strong density fluctuations inside the shear band. It is remarkable that the numerical results as well as experimental data showed that not only the density but also the inhomogeneity of a granular body can be decreased (as in one-dimensional compression) or increased (as in biaxial compression). The results also suggest that fluctuation in stresses (and similarly in couple stresses) can be induced from the fluctuation in void ratio. Therefore in the following finite element calculations, fluctuations in initial state is considered by prescribing a fluctuation in void ratio only.

2734

K. N€ubel, W. Huang / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2719–2743

Fig. 12. Calculated evolution of state during biaxial compression. Curve co shows the average over the whole sample while curve in shows the average inside governing shear band.

4. Evolution of localized deformation patterns 4.1. Visualization of shear bands in experiments Shear bands in granular material have been observed by many researchers with different kinds of techniques. One of the simplest methods to visualize shear bands is to mark the body with colored layers or other markers and observe accumulated strain through glass walls, rubber diaphragms or on the surface of the body. Roscoe [39] examined localization by using this method in the passive earth pressure case.

K. N€ubel, W. Huang / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2719–2743

2735

Another widely used technique is to observe volume changes of the granular body by means of X-rays (e.g. [14,18,30,37]). Shear band patterns can also be visualized by stereo-photogrammetric methods. (e.g. [5,9,12,13,19,32]). An advancement of this method is the so-called particle image velocimetry (PIV) which has been used, among other methods, in our experiments. PIV is a measuring system with which field quantities of a complete two-dimensional field can be monitored. In the case of deformation measurements the whole deformation field can be evaluated. In a granular body the grains serve as tracers. For evaluating a two-dimensional deformation field of a granular body, digital images are taken by means of a CCD-camera (charge couple device) from consecutive deformation states of the granular body. The digital images are transferred to a computer in order to start the evaluation of the deformation field. A so-called area of interest (AOI) is cut out of the digital image and divided into small disjointed sub-areas called interrogation cells. Each interrogation cell comprises a particular number of digital pixels. Due to exposure to light, each pixel of a black and white image represents a particular gray scale intensity forming a characteristic pattern in each interrogation cell. If the deformation of the granular body between two consecutive images is sufficiently small it can be supposed that the patterns of the interrogation cells do not substantially change their characteristics but only their location. Therefore, a local displacement vector can be determined for each interrogation cell between two consecutive digital images by means of the cross-correlation. As an example we show the visualized shear band pattern due to active wall translation. In the active case the wall translates ÔawayÕ from the granular body. Sand is poured in a box with smooth and parallel walls so that plane deformations can be realized. The side walls are out of glass to minimize friction between the sand and the boundary. In order to achieve a high density ðId ¼ 0:8Þ the specimen was prepared by raining as described in detail by Rad and Tumay [38]. Applying this preparation procedure, a relatively homogeneous specimen of Karlsruhe sand was obtained. The specimen of size l ¼ 18 cm, h ¼ 9 cm, b ¼ 15 cm is shown in Fig. 13a. A total of 30 pictures was sufficient to monitor the whole deformation. The parallel shift of the right (u ¼ 0:6 cm) wall away from the sand specimen corresponds to an active earth pressure case. Fig. 13c presents the shear localization at the end of the experiment in a narrow band which is obtained with PIV. A PIV-layer is put over the image, where the magnitude of deviatoric shear deformation ep ¼

ffi pffiffiffiffiffiffiffiffi eij eij ;

with

eij ¼ eij
13ekk dij

ð33Þ

is represented by a gray scale contour plot. Fig. 13b demonstrates the efficiency of PIV by showing the same picture as in Fig. 13c but without PIV-layer. The deformation field can be exactly evaluated by this technique. A more detailed representation of the situation at the end of the experiment is illustrated in Fig. 14. For the shear band width a typical value of 11
15d50 is measured. The width increases slightly towards the free surface. The angle between the straight localization strip and the bottom of the specimen is # ¼ 68°. (Fig. 14). 4.2. Numerical modeling In this part, numerical investigations are presented to show the ability of the applied constitutive model at reproducing the complicated kinematics of localized deformation in sand behind the retaining wall and under strip footing. At first, we perform a simulation of retaining wall under an active earth pressure by moving the the wall away from a sand body. The simulation is performed with a fine mesh (5000 quadratic elements of the same size) so that the localized zone can be covered with sufficient accuracy. It is assumed that the surface of the

2736

K. N€ubel, W. Huang / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2719–2743

Fig. 13. (a) Sand specimen at initial state and location of the interrogation window; (b) sand specimen at the end of the test (c) ep evaluated with PIV.

rigid wall is ideally smooth so that no shear stresses and couple stresses can be transmitted. The parameters calibrated for Karlsruhe sand (Table 2, [20]) are used. The initial state is generated by activating the gravity force out of a state of very low pressure. Hence, a so-called K0 -state is achieved as described in N€ ubel [35]. The initial mean density index is the same as in the experiment, viz. Id ¼ 0:8. Initial fluctuations of void ratio are applied. Like in the experiment, the evolution of shear band pattern is visualized by contour plots of ep (Fig. 15). After a translation of u ¼ 0:1 cm, i.e. u=h ¼ 1:25% the slip surface has already emerged as in the experiment. From Fig. 15, one can see that numerical and experimental results are in good agreement. Next, we study the localized deformation pattern behind a retaining wall in a passive situation. Schwing [40] performed a number of such tests and monitored the volume expansion by radiographs. For numerical simulation, we take the same initial condition as in the above calculation and change the direction of motion of the wall. There is no report on successful numerical reproduction of the test results. It is now shown that numerical simulation can capture such a complicated shear banding processes with the proposed constitutive model (Fig. 16). As in the experiment, the rigid rough wall was pushed inwards and slightly downwards. We use rough surfaces which correspond to constrained translation and rotational degrees of freedom along the wall and the bottom. Shear banding starts at the foot of the wall with a horizontal localization zone. A slightly upwardly curved shear band appears almost simultaneously with the radial set of shear bands starting at the top of the wall. In the final state, the experiment as well as the results of the computation show very clearly the kinematics due to such a wall movement in the granular body. The most activated shear band is the lower one ascending slightly curved. The other localized zones are less pronounced indicating again that different friction angles up are mobilized in the shear bands. The proposed model captures not only the kinematics of soil behavior. Agreement in quantitative measures between model response and experimental results can also be obtained, provided that the material

K. N€ubel, W. Huang / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2719–2743 15 d50

13 d50

(2)

1.0

1.0

0.8

0.8

εp / max εp

εp / max εp

(1)

0.6 0.4 0.2

2737

0.6 0.4 0.2 0.0

0.0 13

14

15 x [cm]

16

13

17

14

15 16 x [cm]

17

(1) (2)

(3) (4)

11 d50

11 d50

(4)

1.0

1.0

0.8

0.8

εp / max εp

εp / max εp

(3)

0.6 0.4 0.2

0.6 0.4 0.2 0.0

0.0 14

15

16 x [cm]

17

18

15

16

17 x [cm]

18

19

Fig. 14. Deformation pattern and evaluated thickness of the shear band in the active case in the experiment.

parameters and initial state are obtained carefully. We compare computations with an experiment of a strip footing on Silver Lighton Buzzard (SLB) sand performed by Tatsuoka et al. [43,46]. For modeling the footing experiments by Tatsuoka et al. [43,46] a high resolution mesh of 20 000 elements with a mesh width of 5d50 ¼ 3  10
3 m was used. The dimension of the discretized solution space was according to the experiment l ¼ 0:6 m, h ¼ 0:3 m. The hypoplastic parameters for the SLB sand have been determined by Herle and Tejchman [22] (see Table 3). For the initially dense SLB sand a mean void ratio of ep0 ¼ 0:55 was used, which corresponds to Id ¼ 0:86. Fluctuation of initial void ratio was attrib-

Table 2 Material parameters for Karlsruhe sand with respect to the micro-polar hypoplastic model uc ¼ 30° ed0 ¼ 0:53 ec0 ¼ 0:84 ei0 ¼ 1:00 d50 ¼ 0:40  10
3

hs ¼ 5800 [MPa] n ¼ 0:28 a ¼ 0:13 b ¼ 1:05 g ¼ 1:50

2738

K. N€ubel, W. Huang / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2719–2743

Fig. 15. Numerical calculation in comparison to the experiment, active case.

uted to the solution space and the gravity field was generated by an equilibrium loading step with removal of incurred displacement. The element row just beneath the footing serves to simulate an interface. The polar hypoplastic constitutive law is applied also for those interface elements. The only difference lies in a IF
5 variation of the material parameters uIF m. The bearing capacity coefficient c ¼ 1=3uc and d50 ¼ 1  10 2 Nb ¼ 2Pmax =cB versus the normalized settlement u=B is displayed in Figs. 17 and 18. As in the experiment a decrease of Nb with increasing breadth B is obtained. For the small foundation (B ¼ 1:2 cm) the peak load is slightly overestimated. This can be attributed to the influence of friction on the boundary layer for relative small foundations with respect to the grain size as discussed by Herle and Tejchman [22]. To increase the quality of the prediction for such small foundations (compared to the grain size) a study of interface behavior (here: foundation to soil) seems to be necessary. However, it can be expected that with increasing breadth of the foundation the influence of the boundary friction on the bottom decreases. With increasing breath of the foundations the results fit quantitatively well with the experimental ones. In the footing problem very large strains are necessary to reach a steady state. High strains are concentrated in the shear bands. Particularly in the shear bands directly at the corner of the footing the ele-

K. N€ubel, W. Huang / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2719–2743

2739

Fig. 16. Calculated and experimental evolution of volume strains when moving wall inwards and slightly downwards.

Table 3 Material parameters for Leighton Buzzard (SLB) sand with respect to the micro-polar hypoplastic model uc ¼ 29° ed0 ¼ 0:51 ec0 ¼ 0:79 ei0 ¼ 0:86 d50 ¼ 0:60  10
3

hs ¼ 300 [MPa] n ¼ 0:40 a ¼ 0:16 b ¼ 1:05 g ¼ 0:50

ments become too much distorted. Therefore all calculations have been truncated before reaching a steady state. Note, that even in the experiment a steady state is only indicated and not reached. Tatsuoka et al. [45] visualized the shear band pattern in the experiment by colored sand markers and by making the sand moist after the test (Fig. 19a). The same significant evolution of shear bands can be obtained in the numerical calculation. It is clearly visible that a non-symmetric shear band pattern is formed. The dominating shear band on the left side propagates deeper into the ground than the one on the right side. Note that not only the major slip surfaces are covered by the numerical model but also the less pronounced shear bands. Similarities between the triangular pattern beside the wedge and the less clear shear bands running down from the edge of the foundation to the curved shear band are clearly visible.

K. N€ubel, W. Huang / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2719–2743

Fig. 17. Experimental influence of the footing breadth on 2P =ðcB2 Þ by Tatsuoka et al. [43,46].

400

2P/(γ B2)

2740

B=1.2 cm 200

B=2.4 cm B=10 cm

0 0

0.2 u/B

0.4

Fig. 18. Calculated influence of the footing breadth on 2P =ðcB2 Þ.

Fig. 19. (a) Plane strain bearing capacity test by Tatsuoka et al. [43,46]: (a) experiment and (b) calculation.

K. N€ubel, W. Huang / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2719–2743

2741

5. Summary Localization of shear deformation in narrow shear bands is a fundamental phenomenon in granular materials. Therefore, when modeling the behavior of granular material shear localization has to be considered. The model must include a proper constitutive law with an intrinsic length determining the width of shear bands as well as a representation of microscopic irregularities triggering shear localization. The material behavior of non-cemented grain skeletons can be described using a hypoplastic law. As the void ratio is included as an independent variable, this incrementally non-linear constitutive model allows a strict distinguish between material parameters and state variables. Fluctuations of state variables on the grain level are a decisive property of grain skeletons. State variables must be physically justified which requires a micro-mechanical interpretation. Void ratio fluctuations can be introduced by a physically justified probability density function. The main assumption is that at the densest state fluctuations are minimal and at the critical state maximal. From a mechanical point of view, stress and density distribution are related with each other. If void ratio fluctuations are prescribed, then stress fluctuations emerge satisfying all boundary conditions. Due to fluctuations in void ratio, strain is not homogeneous even under either homogeneous stress or displacement boundary conditions. Shear bands can thus evolve spontaneously in the interior of a granular body. The strain first localizes at zones of low density. Due to dilation, positive feed-back begins to lead to shear band patterns. The onset of shear localization is always below the maximum deviatoric stress. Fluctuations may also change during computation. Shear localization emerges under certain initial and boundary conditions as spontaneous growth of fluctuations. Localization in granular bodies is related to an increase in both void ratio and stress fluctuations. In critical states fluctuations are maximal, so that localization can not evolve anymore. The density of the grain skeleton is a crucial property of the granular body. It influences strongly strength and stiffness. Due to rearrangements of the grain skeleton, not only changes of the mean density occur, but also of density fluctuations. For cyclic shear deformation with small amplitudes, the grain skeleton is compacted and homogenized. A non-polar theory without fluctuations suffices whenever fluctuations decrease, which is the case for compaction. For strain or stress paths leading to shear localization, e.g. monotonic shearing with constant pressure, fluctuations must be taken into account, and a polar theory should be used. In order to scale the width of a shear band, an intrinsic length scale must be included in the theory. This can be done by means of a Cosserat continuum where the additional degree of freedom can be associated with a grain rotation. In mathematical terms the boundary value problem then remains well-posed, even when limit states are reached. With the extended formulation of the constitutive equation including the evolution of non-symmetric stresses and couple stresses, inter-granular friction as resistance to grain sliding and rotation is taken into account by a single parameter. A linear dependency between this frictional parameter and the shear band width is achieved, enabling to determine the frictional parameter with a back analysis. In addition, the new formulation is characterized by approaching the critical void ratio when a critical state is approached. Since evolution of couple stresses can occur spontaneously out of minute fluctuations, a prescription of initial fluctuating couple stresses is not necessary. It suffices to consider void ratio fluctuation at the initial state.

Acknowledgements The authors would like to express their deep thanks to Prof. G. Gudehus of the University of Karlsruhe and to Prof. S.W. Sloan of the University of Newcastle, Australia for their support of this work.

2742

K. N€ubel, W. Huang / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2719–2743

References [1] E. Bauer, Calibration of a comprehensive hypoplastic model for granular materials, Soils Foundations 36 (1) (1996) 13–26. [2] E. Bauer, Conditions for embedding Casagrande’s critical states into hypoplasticity, Mech. Cohes.-Frict. Mater. 5 (2) (2000) 125– 148. [3] E. Bauer, W. Huang, Numerical study of polar effects in shear zones, in: Pande, Pietruszczak, Schweiger (Eds.), Numerical Models in Geomechanics––NUMOG VII, Balkema, 1999. [4] E. Bauer, W. Huang, A numerical study of polar effects in shear zones, in: J. M€ uhlhaus (Ed.), International Conference on Localization and Bifurcation, Perth, Balkema, 2001. [5] R. Butterfield, R.M. Harkness, K.Z. Andrawes, A stereo-photogrammetric method for measuring displacements fields, Geotechnique 20 (3) (1970) 308–314. [6] R. Cudmani, V.A. Osinov, The cavity expansion problem for the interpretation of cone penetration and pressuremeter tests, Can. Geotech. J. 3 (4) (2001) 622–638. [7] R. de Borst, Simulation of strain localization: a reappraisal of the Cosserat-continuum, Engrg. Comput. 8 (1991) 317–332. [8] J. Desrues, R. Chambon, M. Mokni, F. Mazeroll, Void ratio evolution inside shear bands in triaxial san specimens studied by computed tomography, Geotechnique 46 (3) (1996) 524–526. [9] J. Desrues, J. Lanier, P. Stutz, Localization of the deformation in tests on sand sample, Engrg. Fract. Mech. 21 (4) (1985) 909–921. [10] W. Ehlers, W. Volk, On theoretical and numerical methods in the theory of porous media based on polar and non-polar solid materials, Int. J. Solids Struct. 35 (1998) 4597–4616. [11] A.C. Eringen, Polar and Non-local Field Theories, fourth ed., Academic Press, 1976. [12] R.J. Finno, W.W. Harris, M.A. Mooney, G. Viggiani, Strain localization and undrained steady state of sand, J. Geotech. Engrg., ASCE (1996) 462–473. [13] R.J. Finno, W.W. Harris, M.A. Mooney, G. Viggiani, Shear bands in plane strain compression of loose sand, Geotechnique (1) (1997) 149–165. [14] B. Graf, Theoretische und experimentelle Ermittlung des Vertikaldrucks auf eingebettete Bauwerke, Number 96, Ver€ offentlichung des Instituts f€ ur Boden- und Felsmechanik der Universit€at Karlsruhe, 1984. [15] G. Gudehus, A comprehensive constitutive equation for granular materials, Soils and Foundations 36 (1) (1996) 1–12. [16] G. Gudehus, Attractors, percolation thresholds and phase limits of granular soils, in: Behringer, Jenkins (Eds.), Proceedings of Powders and Grains, Balkema, 1997, pp. 169–183. [17] G. Gudehus, Forced and spontaneous polarization in shear zones, in: M€ uhlhaus (Ed.), International Conference on Bifurcation and Localization, Perth, Balkema, 2000. [18] C. Han, I. Vardoulakis, Plane-strain compression experiments on water-saturated fine-grained sand, Geotechniques 41 (1) (1991) 49–78. [19] W.W. Harris, G. Viggiani, A.M. Michael, R.J. Finno, Use of stereophotogrammetry to analyze the development of shear bands in sand, Geotech. Test. J., ASTM 18 (4) (1995) 405–420. [20] I. Herle, Hypoplastizit€at und Granulometrie einfacher Kornger€ uste, Number 142, Ver€ offentlichung des Instituts f€ ur Boden- und Felsmechanik der Universit€at Karlsruhe, 1997. [21] I. Herle, G. Gudehus, Determination of parameters of a hypoplastic constitutive model from grain properties, Mech. Cohes.-Frict. Mater. 4 (5) (1999) 461–486. [22] I. Herle, J. Tejchman, Effects of grain size and pressure level on bearing capacity on sand, in: Oka Asaoka, Adachi (Eds.), Deformation and Progressive Failure in Geomechnics, Pergamon, Nagoya (Japan), 1997, pp. 781–786. [23] W. Huang, Hypoplastic modelling of shear localisation in granular materials, Ph.D. Thesis, Graz University of Technology, 2000. [24] W. Huang, E. Bauer, Numerical investigation of shear localization in a micro-polar hypoplastic material, Int. J. Numer. Anal. Methods Geomech. 27 (2003) 325–352. [25] W. Huang, M. Hjiaj, S.W. Sloan, Bifurcation analysis for shear localization in non-polar and micro-polar hypoplastic continua, J. Engrg. Math., submitted. [26] W. Huang, K. N€ ubel, E. Bauer, A polar hypoplastic model for granular material, Mech. Mater. 34 (2002) 563–576. [27] T.J.R. Hughes, The Finite Element Method, Prentice-Hall, Englewood Cliffs, New Jersey, 1987. [28] D. Kolymbas, A generalized hypoelastic constitutive law, in: Proceedings of 11th International Conference of Soil Mechanics and Foundation Engineering, vol. 5, Balkema, 1985, p. 2626. [29] D. Kolymbas, Introduction to Hypoplasticity, A.A. Balkema, Rotterdam, Netherlands, 2000. [30] D. Lesniewska, Z. Mr oz, Study of evolution of shear band systems in sand retained by flexible walls, Int. J. Numer. Anal. Methods Geomech. 25 (9) (1994) 909–932. [31] L.E. Malvern, Introduction to the Mechanics of a Continuous Medium, Prentice Hall, Englewood Cliffs, New Jersey, 1969. [32] A. Mooney, G. Viggiani, R.J. Finno, Undrained shear deformation in granular material, J. Engrg. Mech., ASCE (1997) 577–585. [33] H.B. M€ uhlhaus, Scherfugenanalyse bei granularem Material im Rahmen der Cosserat-Theorie, Ingenieur Arch. (56) (1986) 389– 399.

K. N€ubel, W. Huang / Comput. Methods Appl. Mech. Engrg. 193 (2004) 2719–2743

2743

[34] H.B. M€ uhlhaus, I. Vardoulakis, The thickness of shear bands in granular materials, Geotechnique 37 (3) (1987) 271–283. [35] K. N€ ubel, Numerical and experimental investigation of shear localization of granular material, Number 159, Ver€ offentlichung des Instituts f€ ur Boden- und Felsmechanik der Universit€at Karlsruhe, 2002. [36] M. Oda, Co-ordination number and its relation to shear strength of granular material, Soils Foundation 17 (2) (1977) 29–42. [37] M. Oda, H. Kazama, Microstructure of shear bands and its relation to the mechanisms of dilatency and failure of dense granular soils, Geotechnique 48 (4) (1998) 465–481. [38] N.S. Rad, M.T. Tumay, Factors affecting sand specimen preparation by raining, Geotech. Test. J., ASTM 1 (1987). [39] K.H. Roscoe, The influence of strains in soil mechanics, Geotechnique 20 (2) (1970) 129–170. [40] E. Schwing, Standsicherheit historischer St€ utzw€ande, Internal report: Institut f€ ur Boden- und Felsmechanik der Universit€at Karlsruhe, 1991. [41] M. Shahinpoor, Statistical mechanical considerations on storing bulk solids, Bulk Solids Handling 1 (1) (1981). [42] P. Steinmann, A micropolar theory of finite deformation and finite rotation multiplicative elastoplasticity, Int. J. Solids Struct. 31 (8) (1993) 1063–1084. [43] F. Tatsuoka, S. Goto, T. Tanaka, K. Tani, Y. Kimura, Particle size effects on bearing capacity of footings on granular material, in: Asaoka, Adachi, Oka (Eds.), Deformation and Progressive Failure in Geomechnics, Pergamon, Nagoya (Japan), 1997, pp. 133– 138. [44] F. Tatsuoka, S. Nakamura, C.C. Huang, K. Tani, Strength anisotropy and shear band direction in plane strain tests on sand, Soils Foundations 30 (1) (1990) 35–54. [45] F. Tatsuoka, M. Okahara, T. Tanaka, K. Tani, T. Moritomo, M.S.A. Siddiquee, Progressive failure and particle size effects in bearing capacity of footing on sand, in: Geotechnical Engineering Congress, vol. 27, ASCE, Boulder, USA, 1991, pp. 788–802. [46] F. Tatsuoka, M.S.A. Siddiquee, T. Tanaka. Link among design, model tests, theories and sand properties in bearing capacity of footing on sand, in: Proceedings of 13th ICSMFF, vol. 5, Boulder, USA, 1994, pp. 87–88. [47] J. Tejchman, Modelling of shear localization and autogeneous dynamic effects in granular bodies, Number 140, Ver€ offentlichung des Instituts f€ ur Boden- und Felsmechanik der Universit€ at Karlsruhe, 1997. [48] J. Tejchman, G. Gudehus, Shearing of a narrow granular layer with polar quantities, Int. J. Numer. Anal. Methods Geomech. 25 (1) (2001) 1–28. [49] J. Tejchman, W. Wu, Numerical study on shear band pattern in a Cosserat continuum, Acta Mech. 96 (1993) 61–74. [50] I. Vardoulakis, M. Goldscheider, G. Gudehus, Formation of shear bands in sand bodies as a bifurcation problem, Int. J. Numer. Anal. Methods Geomech. (2) (1978) 99–128. [51] I. Vardoulakis, B. Graf, Calibration on constitutive models for granular materials using data from biaxial experiments, Geotechnique 35 (3) (1985) 299–317. [52] P. von Wolffersdorff, A hypoplastic relation for granular materials with a predefined limit state surface, Mech. Cohes.-Frict. Mater. 1 (1996) 251–271. [53] T. Yoshida, F. Tatsuoka, M.S.A. Siddiquee, Y. Kamegai, Shear banding in sands observed in plane strain compression, in: Proceedings of 3rd International Workshop an Localisation and Bifurcation Theory for Soils and Rocks, Aussois, France, Balkema, Rotterdam, 1994. [54] T.L. Youd, Compaction of sands by repeated shear straining, J. Geotech. Engrg. Div., ASCE 103 (GT8) (1972) 918–922.