# Evaluation of two homogenization techniques for modelling the elastic behaviour of granular materials

Pierre-Yves Hicher, Ching S. Chang

► To cite this version:

Pierre-Yves Hicher, Ching S. Chang. Evaluation of two homogenization techniques for modelling the elastic behaviour of granular materials. J. Eng. Mechanics, ASCE, 2005, 132 (8), pp.1184-1194. 10.1061/(ASCE)0733-9399(2005)131:11(1184) . hal-01006824

HAL Id: hal-01006824

https://hal.science/hal-01006824

Submitted on 2 Feb 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Evaluation of Two Homogenization Techniques for Modeling the Elastic Behavior of Granular Materials

P.-Y. Hicher¹ and C. S. Chang²

¹Professor, Research Institute in Civil and Mechanical Engineering, UMR CNRS 6183, Ecole Centrale Nantes and Univ. of Nantes, France.

²Professor, Dept. of Civil and Environmental Engineering, Univ. of Massachusetts, Amherst, MA 01002. E-mail: chang@ecs.umass.edu

**Abstract:** This paper discusses the capabilities of two homogenization techniques to accurately represent the elastic behavior of granular materials considered as assemblies of randomly distributed particles. The stress-strain relationship for the assembly is determined by integrating the behavior of the interparticle contacts in all orientations, using two different homogenization methods, namely the kinematic method and the static method. The numerical predictions obtained by these two homogenization techniques are compared to results obtained during experimental studies on different granular materials. Relations between elastic constants of the assembly, interparticle properties, and fabric parameters are discussed, as well as the capabilities of the models to take into account inherent and stress-induced anisotropy for different stress conditions.

**keywords:** Granular materials; Stress strain relations; Anisotropy; Constitutive models; Elasticity; Micromechanics; Homogeneity.

The determination of the elastic parameters of the granular assembly based on the characteristics of the particle contacts will be discussed. The predictions of the two models under different stress loading conditions will be compared to experimental results obtained from different granular materials. Recent developments in measuring techniques have allowed more accurate investigations for the behavior of particulate media at very small strains, giving a more comprehensive insight into their elastic properties along different stress paths.

## Introduction

A granular material can be considered as a collection of particles of different sizes and shapes. This discontinuous medium can be represented by an equivalent continuous medium whose mechanical properties should depend on geometric arrangement and contact interactions between interacting particles.

The mechanical properties of a granular assembly can be related to the properties of the microstructure by using homogenization techniques. In this paper, we seek to evaluate how two homogenization methods can be applied to the modeling of the elastic behavior of randomly packed granular materials.

A first model was developed based on the so-called kinematic hypothesis, which states that every particle displaces in accordance with a uniform deformation field (Walton 1987; Chang 1988; Chang et al. 1989; Rothenburg and Bathurst 1989). A second model was developed using a static hypothesis (Cambou et al. 1995; Chang and Gao 1996; Liao et al. 1997), which relates the average stress of the granular assembly to a mean field of particle contact forces. Using a static hypothesis, the particle displacements are considered to fluctuate around a mean displacement field that represents the best fit for the actual displacements (Liao et al. 1997). This mean displacement field can be derived by using the least-squares method, which leads to a more relaxed kinematic condition in comparison to the kinematic hypothesis.

## Micromechanical Approach

In recent years, several studies have been devoted to the elastic characteristics of granular materials within the micromechanical framework (Walton 1987; Chang 1988; Rothenburg and Bathurst 1989; Jenkins 1988; Chang and Liao 1994; Cambou et al. 1995). All these studies fundamentally presume the following two relationships.

1. Microscale contact law: This provides the relationship between force and displacement at an interparticle contact, which describes the material behavior at a microscale.

2. Link between micro and macro variables: This provides relationships between macroscale variables and micro-scale variables. There are two hypotheses of micro-macro relationships:

a. Static hypothesis provides a relationship between interparticle contact forces and the assembly stress, and

b. Kinematic hypothesis provides a relationship between interparticle contact displacements and the assembly strain.

## Interparticle Behavior

Let *N* be the total number of interparticle contact orientations in the packing. The elastic contact stiffness of two particles at the *c*th contact is defined by
<figure><img src="image_2.png" /><figcaption>Fig. 1. Local coordinate at interparticle contact</figcaption></figure>

$$ \Delta f_i^c = k_{ij}^c \Delta \delta_j^c \quad (1) $$

where $k_{ij}^c$=interparticle contact stiffness tensor. For two particles in contact, a local coordinate system can be constructed for each contact with three orthogonal base unit vectors: **n** is normal to the contact plane; **s** and **t** are tangential to the contact plane as shown in Fig. 1.

Let $k_n^c$ be the compressive contact stiffness in normal direction and $k_i^c$ the shear contact stiffness. Assuming the shear contact stiffness is the same in **s** and **t** directions and that there is no coupling effect between normal and shear directions, the contact stiffness tensor $k_{ij}^c$ can then be expressed in terms of the unit vectors **n**, **s**, and **t**, as

$$ k_{ij}^c = k_n^c n_i^c n_j^c + k_i^c (s_i^c s_j^c + t_i^c t_j^c) \quad (2) $$

For each particle contact, the corresponding auxiliary local coordinate system is related to the global coordinate system according to (see Fig. 1)

$$ \mathbf{n} = (\cos \gamma, \sin \gamma \cos \beta, \sin \gamma \sin \beta) $$

$$ \mathbf{s} = (-\sin \gamma, \cos \gamma \cos \beta, \cos \gamma \sin \beta) $$

$$ \mathbf{t} = (0, -\sin \beta, \cos \beta) \quad (3) $$

The vector **s** is located on the plane consisting of **x** and **n**. The vector **t** is perpendicular to this plane and can be obtained by the cross product of **n**×**s**. The rolling resistances between two particles are not discussed in this paper.

The value of the stiffness for two elastic spheres can be estimated from a Hertz–Mindlin formulation (Mindlin and Deresiewicz 1953). For sand grains, a revised form can be adopted (Chang et al. 1989), given by

$$ k_n = k_{n0} \left( \frac{f_n}{G_g l^2} \right)^n; \quad k_t = k_{t0} \left( \frac{f_n}{G_g l^2} \right)^n \quad (4) $$

where $G_g$=elastic modulus for the grains; $f_n$=contact force in normal direction; $l$=branch length between the two particles; and $k_{n0}$, $k_{t0}$, and $n$=material constants. For two spherical particles, the branch length is twice that of particle radius $l=2r$. Let the exponent $n=1/3$ and

$$ k_{n0} = G_g r \left( \frac{\sqrt{12}}{1 - v_g} \right)^{2/3} \quad (5) $$

the normal contact stiffness in Eq. (4) is equivalent to Hertz's contact formulation.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Schematic plot of two particles in contact</figcaption></figure>

The stress-strain relationship for an assembly can be determined from integrating the behavior of interparticle contacts in all orientations. In the integration process, a micro-macro relationship is required.

## Micro-Macro Relationship: Kinematic Hypothesis

Granular material is envisioned as a collection of particles. Under deformation, particles in the material undergo translation and rotation. In order to establish a link between the discrete particle system and its equivalent continuum system, a continuum displacement field $u_i(\mathbf{x})$ is constructed in such a way that the displacement at the centroid of the nth particle, $u_i^n$, coincides with the displacement field, i.e.,

$$ u_i(\mathbf{x}^n) = u_i^n \quad (6) $$

where $\mathbf{x}^n$=location of the nth particle.

In classic continuum mechanics, a linear displacement field is employed to describe the deformation of a representative volume element. For a granular material, the size of a representative volume element is relatively large. Thus we approximate the displacement field by

$$ u_i(\mathbf{x}) = u_i + u_{i,j}x_j \quad (7) $$

where $u_i$, $u_{i,j}$=constants for the representative volume. The second-rank tensor $u_{i,j}$ has nine components.

In the further derivation, for the sake of convenience, we define the branch vector $l_i^c$ as the vector from the centroid of particle “a” to that of particle “b” as shown in Fig. 2.

$$ l_i^c = x_i^b - x_i^a \quad (8) $$

Given the position vector of the contact point as $x_i^c$, the following may be written:

$$ x_i^a = x_i^c - r_i^{ac}; \quad x_i^b = x_i^c - r_i^{bc} \quad (9) $$

Thus $l_i^c = r_i^{ac} - r_i^{bc}$.

Between the two individual convex-shaped particles, as shown in Fig. 2, the interparticle compression $\delta_i^c$ at the contact point c can be obtained as follows:

$$ \delta_i^c = \delta_i^{ab} = u_i(\mathbf{x}^b) - u_i(\mathbf{x}^a) \quad (10) $$

Using the polynomial expression in Eq. (7), we can describe the interparticle compression using the continuum field of displacement.

$$ \delta_i^c = u_{i,j} l_j^c \quad \text{or} \quad \Delta \delta_i^c = \Delta u_{i,j} l_j^c \quad (11) $$

The representative volume V is subjected to forces on its boundary surface S. Neglecting the body force of particles, the
work done per unit volume of the discrete system can be expressed as a summation of the work done over all interparticle contacts in the unit volume.

$$ \Delta W = \frac{1}{V} \sum_{i} f_{i}^{c} \Delta \delta_{i}^{c} = \sigma_{iq} \Delta u_{q,i} \quad (12) $$

Then, by substituting the interparticle compression $\Delta\delta_i^c$ with Eq. (11), it yields

$$ \sigma_{iq} = \frac{1}{V} \sum_{i} f_{qi}^{c} \quad \text{or} \quad \Delta\sigma_{iq} = \frac{1}{V} \sum_{i} \Delta f_{qi}^{c} \quad (13) $$

### Micro-Macro Relationship: Static Hypothesis

From the static hypotheses proposed by Liao et al. (1997), the relation between the global strain and interparticle displacement in an incremental form can be written as follows (the finite strain condition is not considered):

$$ \Delta u_{j,i} = A_{ik}^{-1} \sum_{c=1}^{N} \Delta \delta_{j}^{c} l_{k}^{c} \quad (14) $$

where the branch vector $l_k^c$ is defined as the vector joining the centers of two particles, and the fabric tensor is defined as

$$ A_{ik} = \sum_{c=1}^{N} l_{i}^{c} l_{k}^{c} \quad (15) $$

Utilizing Eq. (12) and substituting Eq. (14), the mean force on the contact plane of each orientation can be obtained as

$$ f_j^c = \sigma_{ij} A_{ik}^{-1} l_k^c V \quad \text{or} \quad \Delta f_j^c = \Delta \sigma_{ij} A_{ik}^{-1} l_k^c V \quad (16) $$

### Stress-Strain Relationship

For the case of kinematic hypothesis, the stress-strain relationship can be derived using the interparticle contact law [Eq. (1)], the stress definition [Eq. (13)], and the continuum expression of $\Delta\delta_i^c$ [Eq. (11)].

$$ \Delta\sigma_{iq} = C_{iqkl}\Delta u_{k,l} \quad \text{where} \quad C_{iqkl} = \frac{1}{V} \sum_{c=1}^{N} l_{i}^{c} l_{k}^{c} k_{qk}^{c} \quad (17) $$

For the case of static hypothesis, the stress-strain relationship can be derived using Eqs. (1), (14), and (16), it yields

$$ \Delta u_{i,j} = C_{ijmp}^{-1} \Delta \sigma_{mp} \quad \text{where} \quad C_{ijmp}^{-1} = A_{ik}^{-1} A_{mn}^{-1} V \sum_{c} (k_{jp}^{c})^{-1} l_{k}^{c} l_{n}^{c} \quad (18) $$

In the case of random packing with sufficient number of particles, the discrete interparticle contact orientations can be approximated by a continuous function $E(\beta, \gamma)$, where $\beta$ and $\gamma$=two angles in the spherical coordinate system as shown in Fig. 1. Integration of $E(\beta, \gamma)$ over the surface of unit sphere should equal the total number of interparticle contacts, $N$ (i.e., self-consistency equation).

$$ N = \int_{0}^{\pi/2} \int_{0}^{2\pi} E(\beta, \gamma) \sin \gamma d\beta d\gamma \quad (19) $$

<figure><img src="image_4.png" /><figcaption>Fig. 3. Distribution of contact orientations with different values of $a_0$</figcaption></figure>

For any arbitrary parameter $h^c$, which is contact orientation dependent, the summation of $h^c$ over all contacts can be approximated by an expression of integral, given by

$$ \sum_{c=1}^{N} h^c = \int_{0}^{\pi/2} \int_{0}^{2\pi} E(\beta, \gamma) h(\beta, \gamma) \sin \gamma d\beta d\gamma \quad (20) $$

where $N$=total number of contact orientations; $E(\beta, \gamma)$=distribution function of interparticle contact orientations; and $h(\beta, \gamma)$=continuous expression of $h^c$.

For an orthotropic material, the distribution of contact orientations can be represented by a continuous spherical harmonic expansion in three dimensions, the truncated form of the expansion consisting of second-order terms, given by

$$ E(\beta, \gamma) = \frac{N}{2\pi} \left( 1 + \frac{a_0}{4} (3 \cos 2\gamma + 1) + 3a_{22} \sin^2 \gamma \cos 2\beta \right) \quad (21) $$

For a cross-anisotropic material, $a_{22}$ becomes equal to zero and Eq. (21) is reduced to

$$ E(\beta, \gamma) = \frac{N}{2\pi} \left( 1 + \frac{a_0}{4} (3 \cos 2\gamma + 1) \right) \quad (22) $$

and for an isotropic material,

$$ E(\beta, \gamma) = \frac{N}{2\pi} \quad (23) $$

In three dimensions, the inherent anisotropy can be represented by a distribution whose major axis often coincides with the vertical direction. An example of distributions with different values of $a_0$ is shown in Fig. 3. The anisotropy of the packing structure creates a mechanical behavior anisotropy whose directions are identical to those of the geometric anisotropy.

It was found [see Chang and Misra (1990)] that the anisotropy can be characterized by a fabric tensor that matches the coefficients of the spherical harmonic expansion, for example, in the case of an orthotropic material

$$ [F_{ij}] = \begin{bmatrix} 1 + a_0 & 0 & 0 \\ 0 & 1 - \frac{a_0}{2} + 3a_{22} & 0 \\ 0 & 0 & 1 - \frac{a_0}{2} - 3a_{22} \end{bmatrix} \quad (24) $$

In addition to the contact distribution, other parameters, which are orientation dependent, can also be expressed by the same manner. For example, the function $h(\beta, \gamma)$ for a cross-anisotropic material can be expressed as
<figure><img src="image_5.png" /><figcaption>Fig. 4. Relationship between void ratio and coordination number</figcaption></figure>

$$h(\beta, \gamma) = h_m \left( 1 + \frac{a_{0h}}{4} (3 \cos 2\gamma + 1) \right) \quad (25)$$

Note that $h_m$=mean value of $h(\beta, \gamma)$ and the coefficient $a_{0h}$ indicates the distribution of $h(\beta, \gamma)$, which needs not to have the same value as $a_0$ in the contact distribution. Using the definition of Eq. (20), the following summation can be replaced by an integral form.

$$A_{ik} = \sum_{c=1}^{N} \frac{r_{ic} r_{kc}}{2\pi} = \frac{NI^2}{2\pi} \bar{A}_{ik}$$ 
$$\text{where } \bar{A}_{ik} = \int_{0}^{\pi/2} \int_{0}^{2\pi} E(\beta, \gamma) n_i n_k \sin \gamma d\beta d\gamma \quad (26)$$

For convenience, we define a packing density index (normalized number of contacts per unit volume), $\rho=NI^3/V$. Note that the total number of contacts $N$ per volume $V$ is normalized by the cube of the mean branch length $l$. The packing density index $\rho$ is unitless. The branch length is defined as the length between centroids of two contact particles. For round particles, the mean branch length is approximately equal to the mean particle size $d$. The packing density index for an assembly of equal-sized spheres can be related to the mean coordination number $\bar{n}$ and the void ratio $e$ of the assembly (Chang et al. 1989)

$$\rho = \frac{NI^3}{V} = \frac{3\bar{n}}{\pi(1 + e)} \quad (27)$$

A relationship between $\bar{n}$ and $e$ obtained for a regular packing of spheres is plotted in Fig. 4. Using the packing density index, we can express the stiffness tensor for the kinematic method

$$C_{ijkl} = \frac{\rho}{d} \int_{0}^{\pi/2} \int_{0}^{2\pi} E(\beta, \gamma) k_{qk} n_i n_k \sin \gamma d\beta d\gamma \quad (28)$$

and flexibility tensor for the static method

$$C_{ijkl}^{-1} = \frac{d}{\rho} \bar{A}_{ik}^{-1} \bar{A}_{mn}^{-1} \int_{0}^{\pi/2} \int_{0}^{2\pi} E(\beta, \gamma) k_{jp}^{-1} n_k n_n \sin \gamma d\beta d\gamma \quad (29)$$

A numerical integration method can be used for the integrals in Eqs. (28) and (29). The integrals are computed based on a set of integration points on the surface of a sphere. Each integration point represents a contact orientation and has a weighting factor. We found that the results are more accurate if the locations of integration points on the surface of a sphere are symmetric about the x-z plane and y-z plane (i.e., fully symmetric). We studied the

performance of this model using different numbers of integration points and found that 74 integration points is adequate.

## Isotropic Loading

### Isotropic Fabric

We will at first assume an isotropic fabric for our particulate material. Under this condition, the parameters of the proposed model are the following:

* Packing density index (normalized number of contacts per unit volume): $\rho=NI^3/V$;

* Mean particle size, $d$; and

* Interparticle elastic constants: $G_g$, $\alpha$, and $n$, where the particle modulus $G_g$ and exponent $n$ have been defined in Eq. (4), and $\alpha=k_{t0}/k_{n0}$. For a representative volume consisting of a large number of round particles, in which the contact orientations are isotropically distributed [Eq. (23)], an analytical form of the shear modulus $G$ and the Poisson's ratio $\nu$ of the granular assembly can be computed based on static [Eq. (29)] or kinematic hypotheses [Eq. (28)]. With the static hypothesis, the following expressions were obtained (see Liao et al. 2000):

$$G = \frac{5\rho k_n \alpha}{6d(3 + 2\alpha)} \quad (30)$$

$$\nu = \frac{1 - \alpha}{2 + 3\alpha} \quad (31)$$

where $\rho$=packing density index. Using Eqs. (4) and (5), under an isotropic stress $\sigma_c$, the shear modulus $G$ can be expressed as

$$G = H_g \rho^{1-n} \left( \frac{\sigma_c}{p_0} \right)^n \quad (32)$$

where $p_0$=reference pressure ($p_0=1$ MPa in the computation in this paper). $H_g$=function depending on elastic property of the grain.

$$H_g = \frac{5\alpha}{2\sqrt{6}(3 + 2\alpha)} G_g \left( \frac{3p_0}{G_g} \right)^n \quad (33)$$

where $G_g$=particle modulus.

Based on the kinematic hypothesis, the shear modulus $G$ and the Poisson's ratio $\nu$ take the following expressions:

$$G = \frac{\rho k_n}{30d} (2 + 3\alpha) \quad (34)$$

$$\nu = \frac{1 - \alpha}{4 + \alpha} \quad (35)$$

An equation similar to Eq. (32) can be obtained for $G$

$$G = H'_g \rho^{1-n} \left( \frac{\sigma_c}{p_0} \right)^n \quad (36)$$

with the following expression for $H'_g$:

$$H'_g = \frac{2 + 3\alpha}{10\sqrt{6}} G_g \left( \frac{3p_0}{G_g} \right)^n \quad (37)$$

Eqs. (31) and (35) show that Poisson's ratios in both methods depend only on the interparticle stiffness ratio $\alpha$. The expressions of the shear modulus are similar in character and differ only by
<figure><img src="image_6.png" /><figcaption>Fig. 5. Shear modulus and Poisson's ratio obtained from static hypothesis compared to those obtained with the kinematic hypothesis</figcaption></figure>

the form of their dependency on α. α is theoretically related to the Poisson's ratio of the particles v_g, assuming a perfectly elastic contact between the two spherical particles (Hertz's law):

$$ \alpha = \frac{2(1 - v_g)}{2 - v_g} \qquad (38) $$

For a given value of the stiffness ratio α=k₀/k₀₀ at the particle level, the static hypothesis gives a lower value of the shear modulus (and a larger value of Poisson's ratio) of the assembly, compared to the one obtained from the kinematic hypothesis (Fig. 5). It can be noticed in particular that the kinematic hypothesis will

compute Poisson's ratio values, which are always lower than 0.25. This can be a restriction upon the use of the model derived from the kinematic hypothesis.

In fact, Eq. (38) leads to values of α in the range 0.8–0.95 for values of 0.1 < v_g < 0.3. In these conditions, the computed values of the Poisson's ratio v of the assembly using Eqs. (31) and (35) are in the range of 0.01–0.05, smaller than typical experimental results obtained from diverse granular assemblies. In real granular materials, such as natural sands, the particles have varied shapes and roughness and do not behave like ideal elastic spheres. This means that Eq. (38) from Hertz theory is not applicable in its present form. A different procedure needs to be applied, based on experimental results on granular materials. Knowing experimentally the Poisson's ratio v of the assembly, Eqs. (31) and (35) can be used to determine, for each homogenization method, the value of the ratio α=k₀/k₀₀. This means that α will have a different value in the two models, according to the following equations: α=(1-2v)/(1+3v) for the static hypothesis, and α=(1-4v)/(1+v) for the kinematic hypothesis.

If we do so, it is interesting to notice that Eqs. (32) and (36) for the shear modulus of the assembly, based on the two homogenization techniques, yield the same expression which is dependent on the Poisson's ratio v of the assembly

$$ H_g = H'_g = \frac{1 - 2v}{2\sqrt{6}(1 + v)} \qquad (39) $$

Typical values for v are around 0.2 (Hicher 1996), which lead to values of α around 0.15 for the kinematic hypothesis and around 0.45 for the static hypothesis.

The packing density index ρ can also be empirically determined from experimental results. The influences of void ratio and confining stress on the Young modulus have been experimentally studied for different granular materials (Hicher 2001). The main results are summarized in Fig. 6(a). The following two conclusions can be drawn: (1) for a given grain shape, no significant influence of the grain size on elastic modulus can be detected, and (2) the modulus is dependent on both void ratio and grain size distribution; for a given void ratio, the Young modulus is higher for soils with uniform gradation (i.e., smaller coefficient of uniformity U_c=d₆₀/d₁₀). Thus the packing density index is no longer a unique relationship with the void ratio e as it was shown in Fig.

<figure><img src="image_7.png" /><figcaption>Fig. 6. Influence of grain size and distribution on Young's modulus (E and p are in Mpa)</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Empirical relationship for packing density as a function of void ratio and of grain size distribution</figcaption></figure>

4 for an assembly of equal-sized spheres. The grain size distribution should also be a factor influencing the packing density index.

Fig. 6(b) shows results obtained on different quartzic sands (Hardin and Black 1966; Hardin and Drnevich 1972; Skoglund et al. 1976; Iwasaki and Tatsuoka 1977; Belloti et al. 1996; Hoque and Tatsuoka 1998), including some studied in this paper. One can see that the general trend is very similar to Fig. 6(a) and confirms the conclusions stated above. If we assume that the grain properties $G_g$ and the Poisson ratio $\nu$ are roughly the same for all these quartzic sands, we can back-calculate, from the experimental results in Fig. 6, the packing densities and establish an empirical relationship between the packing density index $\rho$ and the void ratio $e$ for various coefficients of uniformity $U_c$ (Fig. 7). The mean particle size, $d$, is not an explicit parameter.

Thus the parameters of the proposed model are as follows:

* Packing density index $\rho$ (obtained from void ratio $e$ and coefficients of uniformity $U_c$), and

* Particle modulus $G_g$, stiffness ratio $\alpha$, and exponent $n$.

Therefore one can theoretically determine the shear modulus of an assembly of particles, knowing the elastic properties of the particles, the grain size distribution $U_c$, and the void ratio $e$. In the subsequent study, the following phenomenological procedure was adopted to determine the interparticle parameters of the models. From test results, we can determine directly the stiffness ratio $\alpha$ from the Poisson's ratio of the assembly [Eqs. (31) and (35)]. This will lead to different values of $\alpha$ using either the static or the kinematic hypothesis. The packing density index $\rho$ can be determined from Fig. 7. The value of $G_g$ and the exponent $n$ can then be determined by matching predicted and experimental values

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Comparison between measured and predicted values of Young's modulus for glass ballotoni and Hostun Sand</figcaption></figure>

obtained for the material modulus ($E$ or $G$). The parameters $(\rho, G_g, n)$ are identical for the two homogenization models. This procedure will lead to identical numerical results for the elastic constants of the granular assembly along isotropic loading under either the kinematic or the static hypothesis.

Fig. 8 shows the comparison of experimental and numerical results of the Young modulus as a function of void ratio and isotropic stress for glass ballotini ($d_{50}=1$ mm) and Hostun Sand ($U_c=1.8$, $d_{50}=0.29$ mm). The measured values of Poisson's ratio $\nu=0.21$ for both glass ballotini and Hostun Sand (Biarez and Hicher 1994; Hicher 1998) give a value of $\alpha=0.36$ for the static hypothesis and $\alpha=0.135$ for the kinematic hypothesis. The value of $G_g=2$ MN/mm$^2$ and $n=0.5$ was determined for both glass ballotini and Hostun Sand in order to match the experimental results. By applying the chart in Fig. 7 for the packing density index $\rho$, the computed elastic moduli appear to agree very well with experimental results. Similar agreements were also found for the predicted and measured elastic moduli for Toyoura Sand, Reid Bedford Sand, Ottawa Sand, and Ticino Sand (the experimental results shown in Fig. 6). This demonstrated that elastic moduli can be accurately predicted for granular materials under isotropic loading using either static or kinematic homogenization techniques. The two homogenization models will give identical results, provided that the ratio $\alpha=k_{t0}/k_{n0}$ is determined for each model from the experimental value of the Poisson's ratio $\nu$ of the granular assembly.

## Anisotropic Fabric

Recent experimental studies (Bellotti et al. 1996; Jiang et al. 1997; Hoque and Tatsuoka 1998) have demonstrated the influence of inherent anisotropy on the mechanical response of granular materials. When prepared in the gravitational field, granular materials may have an anisotropic packing structure even for spherical grains. Stress-strain behavior reflects the inherent anisotropic fabric and therefore the elastic matrix will show cross-anisotropy with the vertical and horizontal as principal axes

$$ \left\{ \begin{array}{l} \Delta\varepsilon_{xx} \\ \Delta\varepsilon_{yy} \\ \Delta\varepsilon_{zz} \\ \Delta\gamma_{yz} \\ \Delta\gamma_{zx} \\ \Delta\gamma_{xy} \end{array} \right\} = \left[ \begin{array}{cccccc} 1/E_h & -v_{hh}/E_h & -v_{hv}/E_v & 0 & 0 & 0 \\ -v_{hh}/E_h & 1/E_h & -v_{hv}/E_v & 0 & 0 & 0 \\ -v_{hv}/E_h & -v_{hv}/E_h & 1/E_v & 0 & 0 & 0 \\ 0 & 0 & 0 & 1/G_{hv} & 0 & 0 \\ 0 & 0 & 0 & 0 & 1/G_{vh} & 0 \\ 0 & 0 & 0 & 0 & 0 & 2(1 + v_{hh})/E_h \end{array} \right] \left\{ \begin{array}{l} \Delta\sigma_{xx} \\ \Delta\sigma_{yy} \\ \Delta\sigma_{zz} \\ \Delta\tau_{yz} \\ \Delta\tau_{zx} \\ \Delta\tau_{xy} \end{array} \right\} \quad (40) $$
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Isotropic loading on anisotropic Toyoura Sand and Hime Gravel</figcaption></figure>

The decomposition of the constitutive equations along a set of planes allows us to formulate the dependency of the parameters with respect to the orientation. We can thus obtain a response, which depends on the loading orientation in relation to the direction of material anisotropy.

In order to demonstrate how the models are able to take into account an inherent anisotropy, two examples were selected from experimental results obtained by Hoque and Tatsuoka (1998) on different granular materials. Fig. 9 presents experimental and simulation results on Toyoura Sand ($U_c=1.46$, $d_{50}=0.16$ mm) and Hime Gravel ($U_c=1.33$, $d_{50}=1.73$ mm). The specimens were prepared by air pluviation and tested under air-dried conditions. Strains were measured by means of LDTs placed on the specimen inside the triaxial cell. Measurements of $E_v$ and $v_{vh}$ were obtained by applying a small vertical stress increment $\Delta\sigma_v$ ($=\Delta\sigma_{zz}$) while keeping the horizontal stress $\sigma_h$ constant (i.e., $\Delta\sigma_{xx}=\Delta\sigma_{yy}=0$)

$$E_v = \frac{\Delta\sigma_v}{\Delta\varepsilon_v} \quad \text{and} \quad v_{vh} = -\frac{\Delta\varepsilon_h}{\Delta\varepsilon_v} \qquad (41)$$

Note that $\Delta\varepsilon_v=\Delta\varepsilon_{zz}$ and $\Delta\varepsilon_h=\Delta\varepsilon_{xx}=\Delta\varepsilon_{yy}$. When applying a small horizontal stress increment $\Delta\sigma_h$ ($=\Delta\sigma_{xx}=\Delta\sigma_{yy}$) while keeping $\sigma_v$ constant (i.e., $\Delta\sigma_{zz}=0$), the following relations can be used:

$$E_h = (1 - v_{hh}) \frac{\Delta\sigma_h}{\Delta\varepsilon_h} \quad \text{and} \quad \frac{2v_{hv}}{1 - v_{hh}} = - \frac{\Delta\varepsilon_v}{\Delta\varepsilon_h} \qquad (42)$$

One can observe that the two above equations for triaxial loading are not sufficient to determine the three unknowns: $E_h$, $v_{hv}$, and $v_{hh}$. The assumption made by Hoque and Tatsuoka (1998) was $v_{hh}=v_{vh}$ under isotropic stress states. Furthermore, due to the limited accuracy of lateral strain measurements, the writers believed that the measured values of $v_{hv}$ were not reliable and thus did not report them. This shows how difficult it is to measure with sufficient accuracy all the elastic constants for an anisotropic specimen.

<table><caption>**Table 1.** Elastic Parameters for Toyoura Sand ($e=0.65$): Static Hypothesis</caption><thead><tr><th>ρ</th><th>a<sub>0</sub></th><th>G<sub>g</sub></th><th>α</th><th>a<sub>0α</sub></th><th>n</th></tr></thead><tbody><tr><td>2.5</td><td>0.24</td><td>2,100,000</td><td>0.45</td><td>0</td><td>0.5</td></tr></tbody></table>

<table><caption>**Table 2.** Elastic Parameters for Toyoura Sand ($e=0.65$): Kinematic Hypothesis</caption><thead><tr><th>ρ</th><th>a<sub>0</sub></th><th>G<sub>g</sub></th><th>α</th><th>a<sub>0α</sub></th><th>n</th></tr></thead><tbody><tr><td>2.5</td><td>0.17</td><td>2,100,000</td><td>0.3</td><td>0.27</td><td>0.5</td></tr></tbody></table>

In this study, the same assumption is made with regard to $v_{hh}$. Although it may create some uncertainty concerning the exact value of $E_h$, this assumption will not affect the discussion of the model capabilities in reproducing the overall behavior of anisotropic granular materials.

Toyoura Sand and Hime Gravel are both uniformly graded materials with $U_c<2$. The values of the vertical moduli $E_v$ (Fig. 9) are in very good accord with the moduli for Toyoura Sand and Hime Gravel presented in Fig. 6, because the moduli presented in Fig. 6 correspond to measurements of vertical stiffness. Fig. 9 shows that the values of the horizontal moduli $E_h$ are invariably smaller for Toyoura Sand and Hime Gravel, indicating the packing anisotropy. The anisotropy comes from many sources, such as the particle shape, the packing arrangement, the curvature of interparticle contact, etc. We have no detailed information on these factors. In order for the models to capture this inherent anisotropy, we lump all effects into the anisotropy in contact number $a_0=0.24$ in order to match the ratio of $E_h/E_v$. In the case of the kinematic hypothesis, the anisotropy factor was found to be slightly different: $a_0=0.17$. In both cases, the packing density index $\rho=2.5$.

A Poisson's ratio $v_{vh}=0.17$ was measured during isotropic loading and the same value was assumed for $v_{hh}$. Estimated from Eqs. (31) and (35), the corresponding value of $\alpha=0.45$ and an overall isotropy ($a_{0\alpha}=0$) was determined for the static hypothesis model, whereas a mean value of $\alpha=0.3$ and an anisotropy factor $a_{0\alpha}=0.27$ was determined for the kinematic hypothesis model. Tables 1 and 2 summarize the parameters for Toyoura Sand.

According to the void ratio value for Hime Gravel, the corresponding value for $\rho$ is 4. The measured ratio between the horizontal and vertical moduli $E_h/E_v=0.6$. To match the $E_h/E_v$ ratio, the anisotropy for the contact distribution is $a_0=1.03$ for the static hypothesis and $a_0=0.72$ for the kinematic hypothesis. The anisotropy for $\alpha$ was also introduced in order to obtain Poisson's ratio values in accordance with the measured ones. (For Hime Gravel, $v_{hh}=v_{vh}=0.15$.) Tables 3 and 4 summarize the model parameters for the two models.

Comparisons between numerical and experimental results show that the two models can accurately represent the influence of the inherent anisotropy under isotropic loading and provide identical results for two selected sets of parameters (Fig. 9). The set of parameters is different for each model and has to be determined according to the experimental results for each granular material. Therefore there is no direct physical meaning for the selected values of each parameter, in particular for those that

<table><caption>**Table 3.** Model Parameters for Hime Gravel: Static Hypothesis</caption><thead><tr><th>Material</th><th>ρ</th><th>a<sub>0</sub></th><th>G<sub>g</sub></th><th>α</th><th>a<sub>0α</sub></th><th>n</th></tr></thead><tbody><tr><td>Hime Gravel (e=0.5)</td><td>4</td><td>1.03</td><td>2,100,000</td><td>0.56</td><td>0.05</td><td>0.5</td></tr></tbody></table>
control the level of anisotropy. One can remark, however, that a constant value of $G_g$ could be assumed, identical for the two models, and that the same empirical relationship between $\rho$ and $e$ can be used, assuming either initial structural isotropy or anisotropy. The results also demonstrate that the inherent anisotropy does not change under the influence of isotropic loading. In particular, the ratio between the horizontal and the vertical moduli remains constant. We will, in a second stage, study the evolution of the elastic constants during anisotropic loading.

# Anisotropic Loading

## Isotropic Fabric

When the applied state of stress is no longer isotropic, the model will predict a stress-induced anisotropy due to the dependency of contact forces on the stiffness of interparticle contacts. For an orthotropic condition, one can then express the flexibility matrix in the axes of anisotropy as follows:

$$ \left\{ \begin{array}{c} \Delta\varepsilon_{11} \\ \Delta\varepsilon_{22} \\ \Delta\varepsilon_{33} \\ \Delta\gamma_{12} \\ \Delta\gamma_{23} \\ \Delta\gamma_{31} \end{array} \right\} = \left[ \begin{array}{ccccc} 1/E_{11} & -\nu_{21}/E_{22} & -\nu_{31}/E_{33} & 0 & 0 \\ -\nu_{12}/E_{11} & 1/E_{22} & -\nu_{32}/E_{33} & 0 & 0 \\ -\nu_{13}/E_{11} & -\nu_{23}/E_{22} & 1/E_{33} & 0 & 0 \\ 0 & 0 & 0 & 1/G_{12} & 0 \\ 0 & 0 & 0 & 0 & 1/G_{23} \\ 0 & 0 & 0 & 0 & 1/G_{31} \end{array} \right] \left\{ \begin{array}{c} \Delta\sigma_{11} \\ \Delta\sigma_{22} \\ \Delta\sigma_{33} \\ \Delta\tau_{12} \\ \Delta\tau_{23} \\ \Delta\tau_{31} \end{array} \right\} \quad (43) $$

Here we will consider only the case of a simple loading path; starting from an isotropic state of stress $\sigma_0$, the stress in direction 1 is increased, while the stresses in directions 2 and 3 are kept constant. After numerical integration of the equations at the particle level [Eqs. (28) and (29)], the elastic constitutive matrix is found to be cross-anisotropic.

For this example, we selected for each model a set of parameters, which adequately reproduces the behavior of an assembly of glass ballotini under isotropic loading. The evolution of the elastic constants is plotted in Fig. 10 as a function of the stress ratio $\sigma_{11}/\sigma_{22}$. One can see that all the components of the elastic matrix evolve with the stress ratio, albeit at different rates. Therefore the material anisotropy, represented, for example, by the ratio $E_{11}/E_{22}$, increases with stress anisotropy. From Fig. 10, one can see that the evolution of the material anisotropy is more pronounced for the kinematic homogenization than it is for the static one. $E_{11}$ ($E_{22}$) computed with the kinematic hypothesis increases at a higher (smaller) rate when $\sigma_v$ increases, compared to the values obtained with the static hypothesis. Therefore the ratio $E_{11}/E_{22}$ increases at a faster rate. The same conclusions can be drawn for the evolution of the shear moduli $G_{12}$ and $G_{23}$. The evolution of the Poisson's ratios shows a higher increase in the value of $\nu_{12}$ in the case of the kinematic hypothesis, while $\nu_{21}$ decreases at a similar rate. $\nu_{23}$ shows a slight decrease in the case of the kinematic hypothesis and a slight increase in the case of the static hypothesis.

This example shows that the two models, even if they give identical results along an isotropic stress path, produce significantly different results along more general stress paths, in particular those with a change in the principal stress ratio. In the following sections, we will compare the results of the numerical

simulations obtained by the two homogenization techniques with selected experimental results on sands with inherent anisotropic fabric subjected to various stress paths.

## Anisotropic Fabric

Toyoura Sand and Hime Gravel were also submitted to different stress paths leading to various values of the ratio $\sigma_v/\sigma_h$ (Hoque and Tatsuoka 1998). Numerical simulations were undertaken with the two models, using the same sets of parameters as the ones determined in the previous section (Tables 1–3). The experimental results showing the evolution of the elastic properties with the state of stress are given in Fig. 11 along with the results of the numerical simulations. The experimental results suggest that the Young moduli $E_v$ and $E_h$ are mainly dependent on the vertical and horizontal stress $\sigma_v$ and $\sigma_h$, respectively. Other experimental studies [see, for example, Bellotti et al. (1996) and Kuwano and Jardine (2002)] reached the same conclusions. However, one can see in Fig. 11 that the evolution of $E_v$ and $E_h$ with $\sigma_v$ and $\sigma_h$ are also influenced by the value of the stress in the perpendicular direction (i.e., $\sigma_h$ for $E_v$ and $\sigma_v$ for $E_h$). The numerical simulations show the same pattern and demonstrate the capability of the models to reproduce the observed behavior along anisotropic stress paths.

<table><caption>**Table 4.** Model Parameters for Hime Gravel: Kinematic Hypothesis</caption><thead><tr><th>Material</th><th>ρ</th><th>a<sub>0</sub></th><th>G<sub>g</sub></th><th>α</th><th>a<sub>0α</sub></th><th>n</th></tr></thead><tbody><tr><td>Hime Gravel (e=0.5)</td><td>4</td><td>0.72</td><td>2,100,000</td><td>0.42</td><td>0.85</td><td>0.5</td></tr></tbody></table>

<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Stress induced anisotropy</figcaption></figure>
<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 11. (a) Vertical modulus versus vertical stress for Toyoura Sand and Hime Gravel. For a given vertical stress, the measured and predicted range are obtained for $\sigma_h/P_a=1.0-2.0$. (b) Horizontal modulus versus horizontal stress for Toyoura Sand and Hime Gravel. For a given horizontal stress, the measured and predicted range are obtained for $\sigma_v/P_a=1.0-2.0$.</figcaption></figure>

One notices, however, that the model derived from the static hypothesis gives a more pronounced influence of the transverse stress on vertical and horizontal Young moduli, compared to the model derived from the kinematic hypothesis. As a consequence, stress induced anisotropy, which can be expressed by the influence of the $\sigma_v/\sigma_h$ ratio on the $E_v/E_h$ ratio, is less marked in the results obtained by the static model [Fig. 12(a)]. The kinematic model gives an evolution of the $E_v/E_h$ ratio, which appears to be in better agreement with the experimental evidence. The evolution of the Poisson's ratio $v_{ph}$ with the stress ratio also shows an influence of the stress state, which is well captured by the models, although with less amplitude than the measured one in the case of the static model [Fig. 12(b)].

The kinematic hypothesis appears therefore to be more suitable for modeling the evolution of the elastic constants of granular materials along anisotropic stress paths. In the following section, we will apply this approach to stress paths with rotation of the principal directions.

## Rotation of Principal Stress Axes

A recent study by Geoffroy et al. (2003) sheds new light on the evolution of the elastic constants along rotational paths by means

<figure><img src="image_14.png" /><figcaption>Fig. 12. Stress-induced anisotropy in Toyoura Sand and Hime Gravel</figcaption></figure>

of a new apparatus, a torsion hollow cylinder device. Measurements of local strains by LVDT permit us to determine the following strain tensor expressed in the axes of initial anisotropy (vertical z and horizontal r, $\theta$).

$$[\varepsilon] = \begin{bmatrix} \varepsilon_{rr} & 0 & 0 \\ 0 & \varepsilon_{\theta\theta} & \gamma_{\theta z}/2 \\ 0 & \gamma_{z\theta}/2 & \varepsilon_{zz} \end{bmatrix} = \begin{bmatrix} \varepsilon_r & 0 & 0 \\ 0 & \varepsilon_\theta & \gamma/2 \\ 0 & \gamma/2 & \varepsilon_z \end{bmatrix} \quad (44)$$

Different loading conditions can be applied, with and without stress rotation, by controlling axial, horizontal, and torsion shear stresses

$$[\sigma] = \begin{bmatrix} \sigma_{rr} & 0 & 0 \\ 0 & \sigma_{\theta\theta} & \tau_{\theta z} \\ 0 & \tau_{z\theta} & \sigma_{zz} \end{bmatrix} = \begin{bmatrix} \sigma_r & 0 & 0 \\ 0 & \sigma_\theta & \tau \\ 0 & \tau & \sigma_z \end{bmatrix} \quad (45)$$

Therefore 16 terms of the flexibility matrix can be determined as follows:

$$\begin{bmatrix} \Delta\varepsilon_r \\ \Delta\varepsilon_\theta \\ \Delta\varepsilon_z \\ \Delta\gamma/\sqrt{2} \end{bmatrix} = \begin{bmatrix} M_{rr} & M_{r\theta} & M_{rz} & M_{r\gamma} \\ M_{\theta r} & M_{\theta\theta} & M_{\theta z} & M_{\theta\gamma} \\ M_{z r} & M_{z\theta} & M_{zz} & M_{z\gamma} \\ M_{\gamma r} & M_{\gamma\theta} & M_{\gamma z} & M_{\gamma\gamma} \end{bmatrix} \begin{bmatrix} \Delta\sigma_r \\ \Delta\sigma_\theta \\ \Delta\sigma_z \\ \sqrt{2}\Delta\tau \end{bmatrix} \quad (46)$$

In their paper, Geoffroy et al. (2003) present experimental values only for the two last columns of the elastic matrix. Measurements of wave propagation velocities by bender and extender elements were used for this purpose in addition to LVDT measurements. The tests were performed on Hostun Sand prepared by air-pluviation. Two types of stress paths were applied on the specimens: a triaxial compression loading (type C) and a torsion after axial loading corresponding to a ratio $\sigma_v/\sigma_h=0.5$ (type K). The writers did not provide any data concerning the inherent anisotropy due to the mode of deposition and assumed that the specimens were close to initial isotropy. We will make the same assumption for the purpose of this study, even if the results presented in the previous section have demonstrated that the material is mainly cross-anisotropic in the case of air-pluviation preparation. Under these conditions, the parameters for Hostun Sand previously determined can be adopted.

Comparisons between experimental results and simulations using the two homogenization techniques are presented in Fig. 13. During triaxial loading (case T) the terms $M_{r\gamma}$, $M_{\theta\gamma}$, $M_{\gamma z}$, and $M_{z\gamma}$ should start null and stay null all along the loading. Despite some discrepancy in the measurements, the experimental data remain close to a mean value equal to zero in accordance with the models prediction. During a torsion shear test, these four terms will
<figure><img src="image_15.png" /><figcaption>Fig. 13. Stress-induced anisotropy during loading tests with and without stress rotation (evolution of $M_{rz}$, $M_{\theta r}$, $M_{z\gamma}$, $M_{\gamma z}$, $M_{rz}$, $M_{\theta z}$, $M_{zz}$, and $M_{rr}$)</figcaption></figure>

evolve with changes in the principal stress directions. Comparisons with experimental data show that the two models are capable of capturing reasonably well the evolution of these four parameters during torsion, considering the dispersion of experimental data.

The evolutions of the terms $M_{rz}$ and $M_{\theta z}$ during triaxial and torsion tests are also presented in Fig. 13. These two terms are always equal in the models and one can see, even though there is also here dispersion in the experimental data, that the two models give identical results and can predict correctly their evolution during triaxial tests, whereas their values remain practically constant during torsion shear tests.

Fig. 13 shows also the evolution of the terms $M_{zz}$ and $M_{\gamma\gamma}$. It can be seen that the two models can capture the decrease of these two terms during triaxial loading, which corresponds to an increase in Young and shear moduli with stress ratio. The amplitude of the evolution is higher in the case of the kinematic hypothesis in accordance with the conclusion of “Anisotropic Loading.” During shearing, these two terms slightly increase with the shear amplitude, corresponding to a decrease in Young and shear moduli. Here also, this tendency can be predicted by the two models.

In conclusion, the models are able to take into account with reasonable accuracy the evolution of the elastic matrix due to stress induced anisotropy as a function of both stress amplitude and the rotation of principal stress axes. The influence of the stress ratio is more marked with the kinematic hypothesis compared to the results obtained with the static hypothesis.

## Summary and Conclusion

Elastic models for granular materials were developed by using homogenization techniques applicable to assemblies of randomly packed particles. The stress-strain relationship for the assembly can be determined from integrating the behavior of interparticle contacts in all orientations by using either a static hypothesis

which relates the average stress of the granular assembly to a mean field of particle contact forces, or a kinematic hypothesis which estimates the interparticle displacements by assuming a linear displacement field.

It has been demonstrated that these two theories, assuming an interparticle law of the Herzian type, can provide an approximate range of the true behavior. The elastic coefficients at the level of the assembly depend on the parameters used to define the contact law at the level of the particles. These parameters comprise the contact law itself, i.e., the normal stiffness, $k_{n0}$, and the shear stiffness as a function of the normal stiffness, $k_{\theta0} = \alpha k_{n0}$, as well as the packing density index $\rho$ (normalized number of contact per unit volume) along a given plane orientation. The dependency of the normal stiffness on the normal force gives a pressure-dependent behavior for the assembly, while the interparticle stiffness ratio $\alpha$ determines the value of the Poisson’s ratio. The influence of the void ratio of the assembly on the elastic properties is taken into account through the packing density index $\rho$. An empirical relationship was proposed in order to link values of $\rho$ to the void ratio. This relationship was based on results obtained for different granular materials with different grain size distributions.

Inherent anisotropy due to the mode of deposition can be modeled by taking into account the orientation dependency of parameters. In the analysis, we assumed a fabric anisotropy by considering that the number of interparticle contacts is a function of the orientation of the contact plane. In some cases, we assume the parameter $\alpha$ is orientation dependent in order to predict the behavior of Poisson’s ratio.

Comparisons between numerical and experimental results showed that the two homogenization models can represent accurately the elastic behavior of different granular materials under isotropic loading and give identical results for two selected sets of parameters. The set of parameters is different for each model and has to be phenomenologically determined according to the experimental results for each granular material. Following this procedure, the models are capable of a very systematic account of the influence of the inherent anisotropy.
The stress-induced anisotropy along the different stress paths is also considered and the numerical simulations demonstrate the capability of the models to reproduce the observed behavior along anisotropic stress paths. The kinematic model gives an evolution of the stiffness ratio with the stress ratio, which appears to be in better agreement with the experimental evidence. The two models allow the transverse stresses to have a more pronounced influence on Young and shear moduli, a consequence for the way the local equations are integrated. This influence becomes less marked if we consider the kinematic hypothesis. It must be pointed out, however, that this discrepancy becomes noticeable only for el- evated stress ratios ($\sigma_1/\sigma_2 > 2$).

## References

Bellotti, R., Jamiolkowski, M., Lo Presti, D. C. F., and O'Neil, D. A. (1996). "Anisotropy of small strain stiffness in Ticino Sand." *Geo-technique*, 46(1), 115-131.

Biarez, J., and Hicher, P.-Y. (1994). *Elementary mechanics of soil behaviour*, Balkema, Rotterdam, 200.

Cambou, B., Dubujet, P., Emeriault, F., and Sidoroff, F. (1995). "Homo- genization for granular materials." *Eur. J. Mech. A/Solids*, 14(2), 255-276.

Chang, C. S. (1988). "Micromechanical modeling of constructive rela- tions for granular material." *Micromechanics of granular materials*, M. Satake and J. T. Jenkins, eds., Elsevier, Amsterdam, 271-279.

Chang, C. S., and Gao, J. (1996). "Kinematic and static hypotheses for constitutive modeling of granulates considering particle rotation." *Acta Mech.*, 115(1-4), 213-229.

Chang, C. S., and Liao, C. L. (1994). "Estimates of elastic modulus for media of randomly packed granules, Part 2." *Appl. Mech. Rev.*, 47(1), 197-206.

Chang, C. S., and Misra, A. (1990). "Application of uniform strain theory to heterogeneous granular solids," *J. Eng. Mech.*, ASCE. 116(10), 2310-2328.

Chang, C. S., Sundaram, S. S., and Misra, A. (1989). "Initial moduli of particulated mass with frictional contacts." *Int. J. Numer. Analyt. Meth. Geomech.*, 13(6), 629-644.

Geoffroy, H., Di Benedetto, H., Duttine, A., and Sauzeat, C. (2003). "Dynamic and cyclic loading on sands: Results and modelling for general stress-strain conditions." *3rd Int. Symposium on Deformation Characteristics of Geomaterials, ISLyon'03*, Balkema, Rotterdam.

Hardin, B. O., and Black, W. L. (1966). "Sand stiffness under various triaxial stresses." *J. Soil Mech. Found. Div.*, 92(1), 27-42.

Hardin, B. O., and Drnevich, V. P. (1972). "Shear modulus and damping in soils: Measurements and parameter effects (Terzaghi Lecture)." *J. Soil Mech. Found. Div.*, 98(6), 603-624.

Hicher, P.-Y. (1996). "Elastic properties of soils." *J. Geotech. Eng.*, 122(8), 641-648.

Hicher, P.-Y. (1998). "Experimental behaviour of granular materials." *Be- haviour of granular materials*, B. Cambou, ed., SpringerWien, New York, 1-97.

Hicher, P.-Y. (2001). "Microstructure influence on soil behaviour at small strains." *Pre-failure deformation characteristics of geomateri- als*, M. Jamiolkowski, R. Lancellotta, and D. Lo Presti, eds., Swets & Zellinger, Lisse, 1291-1297.

Hoque, E., and Tatsuoka, F. (1998). "Anisotropy in elastic deformation of granular materials." *Soils Found.*, 38(1), 163-179.

Iwasaki, T., and Tatsuoka, F. (1977). "Effects of grain size and grading on dynamic shear moduli of sands." *Soils Found.*, 17(3), 19-35.

Jenkins, J. T. (1988). "Volume change in small strain axisymmetric de- formations of a granular material. In: Micromechanics of granular materials," M. Satake and J. T. Jenkins, eds., Elsevier, Amsterdam, 143-152.

Jiang, G.-L., Tatsuoka, F., Fora, A., and Koseki, J. (1997). "Inherent and stress-state-induced anisotropy in very small strain stiffness of a sandy gravel." *Geotechnique*, 47(3), 509-521.

Kuwano, R., and Jardine, R. (2002). "On the applicability of cross- anisotropic elasticity to granular materials at very small strains." *Geo- technique*, 52(10), 727-749.

Liao, C. L., Chan, T. C., Suiker, A. S. J., and Chang, C. S. (2000). "Pressure-dependent elastic moduli of granular assemblies." *Int. J. Numer. Analyt. Meth. Geomech.*, 24, 265-279.

Liao, C. L., Chang, T. P., Young, D. H., and Chang, C. S. (1997). "Stress- strain relationships for granular materials based on the hypothesis of best fit." *Int. J. Solids Struct.*, 34, 4087-4100.

Mindlin, R. D. and Deresiewicz, H. (1953). "Elastic spheres in contact under varying oblique forces." *ASME Trans. J. Appl. Mech.*, 20, 327-344.

Rothenburg, L., and Bathurst, R. J. (1989). "Analytical study of induced anisotropy in idealized granular materials." *Geotechnique*, 39(4), 601-614.

Skoglund, G. R., Cunny, R. W., and Marcuson, W. F. (1976). "Evaluation of resonant column test devices." *J. Geotech. Eng. Div., Am. Soc. Civ. Eng.*, 102(11), 1147-1158.

Walton, K. (1987). "The effective elastic moduli of a random packing of spheres." *J. Mech. Phys. Solids*, 35, 213-226.