# An anisotropic non linear elastic model for particulate materials

Pierre-Yves Hicher, Ching S. Chang

► To cite this version:

Pierre-Yves Hicher, Ching S. Chang. An anisotropic non linear elastic model for particulate materials. Journal of Geotechnical and Geoenvironmental Engineering, 2006, 132 (8), pp.1052-1061. 10.1061/(ASCE)1090-0241(2006)132:8(1052) . hal-01004982

HAL Id: hal-01004982

https://hal.science/hal-01004982

Submitted on 12 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>
# Anisotropic Nonlinear Elastic Model for Particulate Materials

P.-Y. Hicher<fnref n="1" /> and C. S. Chang<fnref n="2" />

This paper presents the development of an elastic model for particulate materials based on micromechanics considerations. A particulate material is considered as an assembly of particles. The stress-strain relationship for an assembly can be determined by integrating the behavior of the interparticle contacts in all orientations and using a static hypothesis which relates the average stress of the granular assembly to a mean field of particle contact forces. Hypothesizing a Hertz-Mindlin law for the particle contacts leads to an elastic nonlinear behavior of the particulate material, we were able to determine the elastic constants of the granular assembly based on the properties of the particle contacts. The numerical predictions, compared to the results obtained during experimental studies on different granular materials, show that the model is capable of taking into account both the influence of the inherent anisotropy and the influence of the stress-induced anisotropy for different stress conditions.

## Introduction

Particulate materials, such as soils or powders, can be considered as a collection of particles of different sizes and shapes. These discontinuous media can be represented by fictitious continuous media whose mechanical properties should depend on geometric arrangement and contact interactions between interacting particles. It has been experimentally shown that the elastic properties of such materials depend on many parameters, such as the material type, the stress condition, the void ratio, the coordination number or the packing structure (Hardin and Black 1966; Hardin and Drnevich 1972; Iwasaki and Tatsuoka 1977; Chung et al. 1984; Lade and Nelson 1987; Hicher 1996). However, it is not easy to predict elastic moduli for particulate material. The main difficulty lies in the complex behavior of stress dependency leading to nonlinear behavior and stress-induced anisotropy. Furthermore, inherent anisotropic behavior occurs due to the mode of deposition. Recent developments in measuring techniques have allowed us to investigate more accurately the behavior of particulate media at very small strains, giving a more comprehensive insight into their elastic properties along different stress paths.

In previous works on granular mechanics, several empirical equations were proposed in the prediction of elastic moduli (Hicks and Monismith 1971; Hardin and Drnevich 1972; Boyce 1980; Brown and Pappin 1981; Biarez and Hicher 1994; Santha 1994). These equations were derived from curve fitting to

measured data. Therefore, the equations usually work well only in a limited range of stress conditions.

Another approach to the prediction of elastic moduli is based on micromechanics in which the granular material is considered as an assembly of particles. The moduli can be predicted directly from the properties of the constituent particles, for example, as in the work by Chang et al. (1989) and Ng and Petrakis (1996). This approach provides a better understanding of the pressure sensitive behavior of granular media. However, under deviatoric stress conditions, the prediction is computationally intensive. It does not provide a straightforward equation as a prediction tool.

In this study, a constitutive model is developed from micro-scale behavior of interacting particles 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. Comparisons with discrete element models for particle assemblies have shown that the static hypothesis leads to representative results. Using a static hypothesis, the particle displacements are considered to fluctuate around a mean displacement field that represents the least square fit for the actual displacements (Liao et al. 1997). This leads to a more relaxed kinematic condition in comparison to the kinematic hypothesis.

The determination of the elastic constants of a granular assembly based on the properties of interparticle contacts will be discussed and the predictions of soil modulus under different loading conditions will be compared to experimental results obtained from different granular materials.

## 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; Jenkins 1988; Rothenburg and Bathurst 1989; Chang and Liao 1994; Cambou et al. 1995). All these studies assume fundamentally the following two relationships:

1. Interparticle contact law—This provides the relationship between force and displacement at an interparticle contact; and

<footnote n="1">¹Research Institute in Civil and Mechanical Engineering, UMR CNRS 6183, Ecole Centrale Nantes, Univ. of Nantes, Nantes, France.</footnote>

<footnote n="2">²Professor, Dept. of Civil and Environmental Engineering, Univ. of Massachusetts, Amherst, MA 01002 (corresponding author). E-mail: chang@ecs.umass.edu</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1. Local coordinate at interparticle contact</figcaption></figure>

2. Link between micro- and macrovariables. There are two micro-macro relationships:

a. Relationship between interparticle contact forces and the assembly stress; and

b. Relationship between interparticle contact displacements and the assembly strain.

## Interparticle Behavior

The contact stiffness of an orientation includes normal stiffness, $k_n^\alpha$, and shear stiffness, $k_t^\alpha$, of the contact plane. The elastic stiffness tensor is defined in

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

where $k_{ij}^\alpha$=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$ be the compressive contact stiffness in normal direction and $k_s$ be the shear contact stiffness. Assuming the shear contact stiffness is same in **s** and **t** directions and that there is no coupling effect between normal and shear directions, the contact stiffness tensor $k_{qk}^c$ can then be expressed in terms of the unit vectors **n**, **s**, and **t**, as

$$k_{ij}^\alpha = k_n^\alpha n_i^\alpha n_j^\alpha + k_t^\alpha (s_i^\alpha s_j^\alpha + t_i^\alpha t_j^\alpha) \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 on the plane defined by **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 Hertz-Mindlin's formulation (Mindlin 1969). For sand grains, a revised form was 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 quartz mineral, the Poisson ratio is about 0.1 and $k_{n0}$ can be approximated by

$$k_{n0} = \sqrt{3/2} G_g l \quad (5)$$

For two spherical particles, the branch length is the particle diameter $l=d$. Using $n=1/3$, Eqs. (5) and (4) are equivalent to Hertz-Mindlin's contact formulation.

## Macro-Micro Relationship

The stress-strain relationship for an assembly can be determined from integrating the behavior of interparticle contacts in all orientations. During the integration process, a relationship is required to link the macro- and microvariables. Using the static hypotheses proposed by Liao et al. (1997), we obtain the relation between the macrostrain and interparticle displacement (here, we do not consider the finite strain condition)

$$\dot{u}_{j,i} = A_{ik}^{-1} \sum_{\alpha=1}^{N} \dot{\delta}_{j}^{\alpha} l_{k}^{\alpha} \quad (6)$$

where $\dot{\delta}_{j}^{\alpha}$=relative displacement between two contact particles and the branch vector $l_k$=vector joining the centers of two contact particles. It is noted that particle contacts include both direct contacts and indirect contacts of neighboring particles associated with a Voronoi polyhedron as discussed by Cambou et al. (2000). For convenience, we let $N$ be the total number of contact orientations. The variables $\dot{\delta}_{j}^{\alpha}$ and $l_k^{\alpha}$ are defined, respectively, as the averaged values of $\dot{\delta}_{j}^{\alpha}$ and $l_k^{\alpha}$ for all contacts belonging to the $\alpha$th orientation. The fabric tensor in Eq. (6) is defined as

$$A_{ik} = \sum_{\alpha=1}^{N} l_i^{\alpha} l_k^{\alpha} \quad (7)$$

Using the principle of energy balance and using Eq. (6), the mean force on the contact plane of each orientation is

$$\dot{f}_{j}^{\alpha} = \dot{\sigma}_{ij} A_{ik}^{-1} l_{k}^{\alpha} V \quad (8)$$

In Eq. (8), the stress increment $\dot{\sigma}_{ij}$ can be obtained by the contact forces and branch vectors for contacts in all orientations (Christofferson et al. 1981; Rothenburg and Selvadurai 1981)

$$\dot{\sigma}_{ij} = \frac{1}{V} \sum_{\alpha=1}^{N} \dot{f}_{j}^{\alpha} l_{i}^{\alpha} \quad (9)$$

Applying the defined contact force in Eq. (8), Eq. (9) is unconditionally satisfied.

## Stress-Strain Relationship

The problem is defined as follows. Initially, we know the global variables ($\sigma_{ij}$ and $\varepsilon_{ij}$) for the assembly and the local variables ($f_j^\alpha$ and $\delta_j^\alpha$) for each contact orientation. For a given loading increment, which can be stress control, strain control or mixed mode, out of the 12 global variables ($\Delta\sigma_{ij}$ and $\Delta\varepsilon_{ij}$), six of them are unknown. The objective is to determine all global vari-
ables ($\sigma_{ij}$ and $\epsilon_{ij}$) and local variables ($f_j^\alpha$ and $\delta_j^\alpha$) at the end of load increment. For a system with $N$ interparticle orientations, the number of unknowns is $3N$ for $f_j^\alpha$ and $3N$ for $\delta_j^\alpha$. The total number of unknowns is $3N+3N+6$.

The following constraints must be satisfied:

1. The local constitutive equation, i.e., Eq. (1). Since there are three equations for each contact plane orientation, the total number of equations is $3N$, $N$ being the total number of interparticle orientations.

2. Static hypothesis between global stress and local forces, i.e., Eq. (8): the number of equations is $3N$.

3. Strain definition between global strain and local displacement, i.e., Eq. (6). The number of equations is 6 (strain is symmetric).

The total number of unknowns is the same as the total number of equations. Therefore, the solution can be determined.

Using Eqs. (1), (6), and (8), the following relationship between stress and strain can be obtained:

$$ \dot{u}_{i,j} = C_{ijmp} \dot{\sigma}_{mp} $$

where

$$ C_{ijmp} = A_{ik}^{-1} A_{mn}^{-1} V \sum_{\alpha=1}^{N} (k_{jp}^{ep})^{-1} l_k^{\alpha} l_n^{\alpha} \quad (10) $$

The summation in Eq. (10) can be replaced by an integral over orientations. The integral can lead to a closed-form solution for the elastic modulus of randomly packed equal-size particles (Chang et al. 1995). However, due to the nonlinearity nature of the local constitutive equation, a numerical calculation with iterative process is necessary to carry out the summation in Eq. (10) (Chang and Hicher 2005). In order to facilitate the numerical calculation, the orientations are selected to coincide with the locations of Gauss integration points in a spherical coordinate. Summation over these orientations with the Gauss weighting factor for each orientation is equivalent to determining the integral over orientations. The results were more accurate by using a set of fully symmetric integration points. From a study of the performance of using different numbers of orientations, we found $N \ge 74$ to be adequate (Chang and Hicher 2005).

## Particulate Materials with Isotropic Fabric

### Isotropic Loading

Extensive studies of the elastic properties of soils have been done in the last 20 years, thanks to the development of new experimental techniques, which measure with a higher degree of accuracy very small deformations of particulate materials. The results have shown that elastic behavior exists only for very small strains, typically smaller than $10^{-5}$; the elastic behavior is nonlinear, with a dependency upon both the actual state of stress and the actual void ratio.

We will at first assume an isotropic fabric for our particulate material. In these conditions, the parameters of the proposed model are as follows:

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

* Mean particle size $d$; and

* Interparticle elastic constants: $G_g$, $\alpha$, and $n$.

Note that the total number of contacts $N$ per volume $V$ is normalized by the cube of mean branch length $l$. The packing

<figure><img src="image_3.png" /><figcaption>Fig. 2. Relationship between void ratio and coordination number</figcaption></figure>

density $\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 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{Nl^3}{V} = \frac{3\bar{n}}{\pi(1 + e)} \quad (11) $$

A relationship between $\rho$ and $e$ obtained for a regular packing of spheres is plotted in Fig. 2.

For a representative volume consisting of a large number of round particles, in which the contact orientations are isotropically distributed, an analytical form of the shear modulus $G$ and the Poisson's ratio $\nu$ of the granular assembly can be obtained (Liao et al. 2000)

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

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

where $\alpha = k_{t0}/k_{n0}$ and $\rho$=packing density. Poisson's ratio depends therefore only on the interparticle stiffness ratio $\alpha$. 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 (14) $$

where $p_0$=reference pressure ($p_0$=1 MPa in the computations in this paper) and $H_g$=function depending on the elastic properties of the grains

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

where $G_g$ is defined in Eq. (4).

Eq. (14) was derived from spherical grains. In a subsequent discussion, the relationship is extended to sands by empirically determining the packing density from experimental results. The influences of void ratio and confining stress on Young modulus have been experimentally studied for different granular materials (Hicher 2001). The main results are summarized in Fig. 3(a). We can draw the following two conclusions: (1) for a given grain shape there is no influence of the grain size on elastic modulus. This implies that the modulus is dependent on the packing density
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Influence of grain size and distribution on Young modulus ($E$ and $p$ are in megapascal)</figcaption></figure>

(normalized by particle size), not on the particle size. (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_{60}/d_{10}$). Thus the packing density is no longer a unique relationship with the void ratio $e$ as it was shown in Fig. 2 for an assembly of equal-sized spheres. The grain size distribution should also be a factor influencing the packing density.

Fig. 3(b) shows results obtained on different quartzic sands, including some studied in this paper. One can see that the general trend is very similar to Fig. 3(a) and confirms the conclusions stated above. If we assume that the grain properties are roughly the same for all these quartzic sands, we can back-calculate, from the experimental results in Fig. 3, the packing densities and establish an empirical relationship between the packing density $\rho$ and the void ratio $e$ for various coefficients of uniformity $U_c$ (Fig. 4). Mean particle size, $d$, is not an explicit parameter.

Thus, the parameters of the proposed model are as follows:

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

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

One can therefore determine theoretically 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$.

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Empirical relationship for packing density as function of void ratio and of grain size distribution</figcaption></figure>

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

The approach that we adopted in the subsequent study was to determine the parameters of the model through the following procedure. From the test results, we can determine the ratio $\alpha=k_{t0}/k_{n0}$ directly from the Poisson's ratio of the assembly [Eq. (12)]. The packing density $\rho$ can be determined from Fig. 4. The value of $G_g$ and the coefficient of nonlinearity $n$ can therefore be determined by fitting predicted and experimental values obtained for the material modulus ($E$ or $G$).

Fig. 5 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 (Hicher 1998), giving a value of $\alpha=0.36$. The value of $G_g=2,000,000$ N/mm$^2$ and $n=0.5$ was used for both glass ballotini and Hostun sand. By using the chart in Fig. 4 for packing density $\rho$, the computed elastic modulus has a very good agreement with experimental results. Similar agreement was 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. 3). Results on Toyoura sand and Ticino sand will be discussed in the following sections. This demonstrated that elastic moduli can be predicted accurately for granular materials under an isotropic loading, provided that one can assume the initial isotropy of the packing. This hypothesis will be discussed.

## Anisotropic Loading

When the applied state of stress is no longer isotropic, the model will predict a stress-induced anisotropy due to the dependency

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Stress induced anisotropy</figcaption></figure>
of the interparticle contact stiffness with the contact forces. The flexibility matrix can be expressed in the axes of anisotropy as follows:

$$ \begin{Bmatrix} \varepsilon_{11} \\ \varepsilon_{22} \\ \varepsilon_{33} \\ \gamma_{12} \\ \gamma_{23} \\ \gamma_{31} \end{Bmatrix} = \begin{bmatrix} \frac{1}{E_{11}} & -\frac{\nu_{21}}{E_{22}} & -\frac{\nu_{31}}{E_{33}} & 0 & 0 & 0 \\ -\frac{\nu_{12}}{E_{11}} & \frac{1}{E_{22}} & -\frac{\nu_{32}}{E_{33}} & 0 & 0 & 0 \\ -\frac{\nu_{13}}{E_{11}} & -\frac{\nu_{23}}{E_{22}} & \frac{1}{E_{33}} & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{1}{G_{12}} & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{1}{G_{23}} & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{1}{G_{31}} \end{bmatrix} \begin{Bmatrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{33} \\ \tau_{12} \\ \tau_{23} \\ \tau_{31} \end{Bmatrix} \quad (16) $$

We will consider here 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 remain constant. After numerical integration of all interparticle stiffness, the elastic constitutive matrix is found to be cross anisotropic.

The evolution of the elastic constants is plotted in Fig. 6 as a function of the stress ratio $\sigma_1/\sigma_2$. All the components of the elastic matrix increase with the stress ratio, albeit at different rates. The material anisotropy, represented for example by the ratio $E_1/E_2$, increases with the stress anisotropy.

## Particulate Materials with Anisotropic Fabric

One of the main advantages of the model is its capability of taking into account the structural anisotropy. In fact, 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 the case of an orthotropic material, the anisotropy can be characterized by a fabric tensor, defined by

$$ [F_{ij}] = F_{\text{ave}} \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 (17) $$

where $F_{\text{ave}}=(F_{11}+F_{22}+F_{33})/3$. $a_0$ and $a_{22}$=fabric constants describing the packing anisotropy (Chang et al. 1995). By representing the distribution of contact orientations with a continuous spherical harmonic expansion in three dimensions, it was found [see Chang and Misra (1990)] that the fabric tensor matches the distribution described by a truncated form of the expansion consisting of second-order terms, given by

$$ \bar{E}(\gamma, \beta) = \frac{F_{\text{ave}}}{4\pi} \left[ 1 + \frac{a_0}{4}(3 \cos 2\gamma + 1) + 3a_{22} \sin^2 \gamma \cos 2\beta \right] \quad (18) $$

The meaning of angles $\beta$ and $\gamma$ is given in Fig. 1.

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

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

$$ \bar{E}(\gamma, \beta) = \frac{F_{\text{ave}}}{4\pi} \left[ 1 + \frac{a_0}{4}(3 \cos 2\gamma + 1) \right] \quad (19) $$

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. 7. The anisotropy of the packing structure is responsible to the anisotropy of the mechanical behavior.

When deposited in the gravitational field, granular materials comprise an anisotropic packing structure (or anisotropic fabric), often termed inherent anisotropy. Recent experimental studies (Belloti 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. The elastic matrix due to this inherent anisotropy can be considered as cross anisotropy. Assuming the vertical and horizontal axes are two principal material axes, Eq. (16) becomes

$$ \begin{Bmatrix} \Delta\varepsilon_{xx} \\ \Delta\varepsilon_{yy} \\ \Delta\varepsilon_{zz} \\ \Delta\gamma_{yz} \\ \Delta\gamma_{zx} \\ \Delta\gamma_{xy} \end{Bmatrix} = \begin{bmatrix} \frac{1}{E_h} & -\frac{\nu_{hh}}{E_h} & -\frac{\nu_{vh}}{E_v} & 0 & 0 & 0 \\ -\frac{\nu_{hh}}{E_h} & \frac{1}{E_h} & -\frac{\nu_{vh}}{E_v} & 0 & 0 & 0 \\ -\frac{\nu_{hv}}{E_h} & -\frac{\nu_{hv}}{E_h} & \frac{1}{E_v} & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{1}{G_{vh}} & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{1}{G_{vh}} & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{2(1 + \nu_{hh})}{E_h} \end{bmatrix} \times \begin{Bmatrix} \Delta\sigma_{xx} \\ \Delta\sigma_{yy} \\ \Delta\sigma_{zz} \\ \Delta\tau_{yz} \\ \Delta\tau_{zx} \\ \Delta\tau_{xy} \end{Bmatrix} \quad (20) $$

### Isotropic Loading

In order to demonstrate the capability of the model to take into account an inherent anisotropy, three examples were selected
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Isotropic loading on anisotropic Toyoura sand, Hime gravel, and Chiba gravel</figcaption></figure>

from experimental results obtained by Hoque and Tatsuoka (1998) and Jiang et al. (1997) on different granular materials. Fig. 8 presents experimental and simulation results on Toyoura sand ($U_c=1.46$, $d_{50}=0.16$ mm), Hime gravel ($U_c=1.33$, $d_{50}=1.73$ mm), and Chiba gravel ($U_c=10$, $d_{50}=8$ mm). The specimens were prepared by air pluviation and tested under air-dried conditions. Strains were measured by means of LVDTs 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$ while keeping the horizontal stress $\sigma_h$ constant

$$E_v = \frac{\Delta\sigma_v}{\Delta\varepsilon_v}$$

and

$$v_{vh} = -\frac{\Delta\varepsilon_h}{\Delta\varepsilon_v} \qquad (21)$$

When applying a small horizontal stress increment $\Delta\sigma_h$ [i.e., $\Delta\sigma_{xx}$ in Eq. (20)] while keeping $\sigma_v$ constant [i.e., $\Delta\sigma_{zz}=0$ in Eq. (20)], Eq. (20) leads to the following relations:

$$E_h = (1 - v_{hh}) \frac{\Delta\sigma_h}{\Delta\varepsilon_h}$$

and

$$\frac{2v_{hv}}{1 - v_{hh}} = - \frac{\Delta\varepsilon_v}{\Delta\varepsilon_h} \qquad (22)$$

One can observe from Eq. (22) that the two equations 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, Hoque and Tatsuoka believe that the measured values of $v_{hv}$ were not reliable, and thus were not reported. This shows the difficulties of measuring with sufficient accuracy all the elastic constants for an anisotropic specimen.

In this study, the same assumption is made concerning $v_{hh}$. Although it leads to some uncertainty concerning the exact value of $E_h$, it will not affect the discussion of 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 vertical moduli $E_v$ (Fig. 8) are in very good accord with the results presented in Fig. 3. The values of horizontal moduli $E_h$ are invariably smaller for Toyoura sand, Hime gravel, and Chiba gravel. One can use the chart in

<table><caption>**Table 1.** Elastic Parameters for Toyoura Sand ($e=0.65$)</caption><thead><tr><th>ρ</th><th>a<sub>ρ</sub></th><th>G<sub>g</sub></th><th>α</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.5</td></tr></tbody></table>

Fig. 4 to determine the packing density corresponding to the void ratio of the tested specimens ($e=0.65$). If we do so, a packing density $\rho=2.5$ is obtained for Toyoura sand.

A Poisson's ratio $v_{vh}=0.17$ was measured during isotropic loading and the same value was assumed for $v_{hh}$. Estimated from Eq. (12), the corresponding value of $\alpha=0.45$. The ratio between horizontal and vertical moduli $E_h/E_v$ was found constant during isotropic loading and equal to 0.9, indicating the packing anisotropy. The anisotropy comes from many sources: the particle shape, the packing arrangement, the curvature of interparticle contact, etc. We have no information on these factors. In order for the model to capture this inherent anisotropy, an orientational dependent $\rho$ was assumed, which has the same material axis as the fabric. For a cross-anisotropic material, Eq. (19) can be used to define the orientational dependent $\rho(\gamma, \beta)$, in which $F_{ave}=\rho$ and $a_0=a_\rho$. The anisotropy of the packing density $a_\rho=0.24$ was determined in order to match the ratio $E_h/E_v$. Table 1 summarizes the parameters for Toyoura sand.

According to void ratios for Hime gravel and Chiba gravel ($U_c<2$), the values of $\rho$ are 4 and 4.2, respectively. The measured ratio between the horizontal and vertical moduli were $E_h/E_v=0.6$ for Hime gravel and 0.45 for Chiba gravel. To match the $E_h/E_v$ ratios, the anisotropy for the packing density was $a_\rho=1.03$ for Hime gravel and 1.5 for Chiba gravel. The anisotropy for $\alpha$ is also introduced in order to obtain anisotropic Poisson's ratio values in accordance with the measured ones (for Hime gravel, $v_{hh}=v_{vh}=0.15$, and for Chiba gravel $v_{hh}=v_{vh}=0.25$). For a cross-anisotropic material, Eq. (19) can be used to define the orientational dependent $\alpha(\gamma, \beta)$, in which $F_{ave}=\alpha$ and $a_0=a_\alpha$. Table 2 summarizes the model parameters for these two materials.

Comparisons between numerical and experimental results show that the model can represent accurately the influence of inherent anisotropy under isotropic loading. The results also demonstrate that the inherent anisotropy does not change during isotropic loading. In particular, the ratio between horizontal and vertical moduli remains constant. We will, in a second stage, study the evolution of the elastic constants during anisotropic loading.

## Anisotropic Loading

The three previous materials were also subjected to different anisotropic loading. For Toyoura sand and Hime gravel, test programs corresponding to different stress paths lead to various values of the ratio $\sigma_h/\sigma_v$ (Hoque and Tatsuoka 1998). The experimental results showing the evolution of the elastic properties with the state of stress are given in Fig. 9 along with the results of the numerical simulations. The parameters determined from isotropic loading were used for these calculations (Tables 1 and 2). The experimental results suggest that the Young moduli $E_v$ and $E_h$ are

<table><caption>**Table 2.** Model Parameters for Hime Gravel and Chiba Gravel</caption><thead><tr><th>Material</th><th>ρ</th><th>a<sub>ρ</sub></th><th>G<sub>g</sub></th><th>α</th><th>a<sub>α</sub></th><th>n</th></tr></thead><tbody><tr><td>Hime gravel (e=0.5)</td><td>0.53</td><td>1.03</td><td>2,100,000</td><td>0.56</td><td>0.05</td><td>0.5</td></tr><tr><td>Chiba gravel (e=0.49)</td><td>0.36</td><td>1.5</td><td>2,100,000</td><td>0.33</td><td>-0.07</td><td>0.52</td></tr></tbody></table>
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Stress dependency on the horizontal and vertical Young's modulus for Toyoura sand and Hime gravel</figcaption></figure>

mainly dependent on the vertical and horizontal stress $\sigma_v$ and $\sigma_h$ respectively. Other experimental works [see for example Bellotti et al. (1996) and Kuwano and Jardine (2002)] reached the same conclusions. However, one can see in Fig. 9 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$). Numerical simulations show the same pattern and demonstrate the capability of the model to reproduce the observed behavior along anisotropic stress paths. One can notice, however, that the model gives a more pronounced influence of the transverse stress on vertical and horizontal Young moduli. This influence is due to the mode of integration along a set of planes of different orientations. Each plane contributes to the assembly behavior and, therefore, the stresses in all directions affect this behavior. Similar results were obtained on Chiba gravel along different stress paths. As a consequence, the 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 model compared to experimental evidence [Fig. 10(a)]. The evolution of Poisson's ratio $v_{vh}$ also shows an influence of the stress state which is captured by the model, although with less rate of change than the measured one [Fig. 10(b)].

An extensive study of the elastic behavior of Ticino sand was performed by Bellotti et al. (1996). The tests were performed in a large calibration chamber and elastic properties were measured by means of propagating seismic body waves inside the specimens. With this technique, it is possible to measure wave propagation velocities in more than two directions (vertical and horizontal), which opens the way to determine all the elastic

<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Stress induced anisotropy in Toyoura sand and Hime gravel</figcaption></figure>

<table><caption>**Table 3.** Model Parameters for Ticino Sand along Isotropic Loading Path</caption><thead><tr><th>ρ</th><th>a<sub>ρ</sub></th><th>G<sub>g</sub></th><th>α</th><th>n</th></tr></thead><tbody><tr><td>0.37</td><td>-0.4</td><td>2,100,000</td><td>0.38</td><td>0.45</td></tr></tbody></table>

constants. In comparison with the previous measurement technique, the advantage is that all the elastic constants can be directly measured. The specimens were prepared by air pluviation, leading to an initial fabric which can be considered, as for the examples examined before, as cross anisotropic in vertical and horizontal axes. Several constant stress ratio tests were performed over a range of values of stress ratio $K=\sigma_h/\sigma_v$ from 0.5 to 2. Results under isotropic loading conditions confirmed that the inherent anisotropy can be considered as cross anisotropic with isotropy in the horizontal plane since the wave velocities propagating in this plane were independent of the direction of propagation. The wave velocities either propagated along the z axis or polarized in the vertical plane were always lower than those confined in the horizontal plane. These last results lead to lower vertical than horizontal stiffness, contradicting the results discussed previously. Explanation for this phenomenon must await further study.

Considering the physical characteristics of Ticino sand: $U_c=1.6$ and $d_{50}=0.55$ mm, values of $E_v$ obtained along isotropic stress paths are in agreement with the results presented in Fig. 3. We can therefore apply the chart in Fig. 4 to determine $\rho=1.2$ for a given mean unit weight $1.505$ MN/mm$^3$ (relative density $D_r=41\%$, $e=0.79$). Values of $E_h/E_v=1.21$ and $n=0.45$ could be determined from the test results. Measured values of $v_{hh}$ and $v_{hh}$ were used to determine the values of $\alpha$ in vertical and horizontal directions. Table 3 summarizes the set of model parameters.

Comparison between experimental results and numerical simulations is presented in Fig. 11(a). An overall agreement was obtained for the different elastic constants. Constant values of Poisson's ratios are predicted by the model along the isotropic loading path, contrary to the experimental results which showed a slight increase of $v_{hh}$ and a slight decrease of $v_{vh}$ and $v_{hv}$. This

<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Comparison between measured and predicted elastic constants for Ticino sand</figcaption></figure>
<table><caption>**Table 4.** Sets of Parameters for Ticino Sand for Different Stress Ratios</caption><thead><tr><th>K</th><th>ρ</th><th>a<sub>ρ</sub></th><th>α</th><th>a<sub>α</sub></th><th>n</th></tr></thead><tbody><tr><td>0.5</td><td>0.32</td><td>0.11</td><td>0.5</td><td>0.48</td><td>0.45</td></tr><tr><td>1.5</td><td>0.34</td><td>-0.65</td><td>0.48</td><td>0.67</td><td>0.45</td></tr><tr><td>2.0</td><td>0.40</td><td>-0.78</td><td>0.38</td><td>0.23</td><td>0.45</td></tr></tbody></table>

could be the result of a certain evolution of the fabric during isotropic loading due to strain induced anisotropy, even if the authors estimated that strain induced anisotropy could be neglected in this study. This issue cannot be answered by the model in its present form.

When examining experimental data obtained for different stress ratios, one cannot observe any clear trend in the evolution of elastic constants (see for example in Fig. 11, the evolution of Poisson's ratios for different stress ratios), which could also be due to a change of fabric anisotropy. Therefore, the model parameters governing the initial anisotropy were assigned for each loading condition, given in Table 4. The evolution of the packing density anisotropy $a_ρ$ is an indicator of the evolution of the level of anisotropy with the stress ratio.

The results are presented in Fig. 11. Satisfactory agreement is found between experimental results and numerical simulations for moduli. However discrepancies are found in the comparisons for Poisson's ratio. One should be aware that Poisson's ratios for an anisotropic material cannot be experimentally determined with accuracy as discussed previously. Therefore, comparisons of measured Poisson's ratio with numerical simulations are difficult to assess. In any case, the comparisons show that the model captures the general trend for the evolution of elastic constants along stress paths with constant principal stress ratios.

## Rotation of Principal Stress Directions

In all examples treated above, the principal stress directions were maintained constant and coincided with the axes of inherent anisotropy. We will now study the evolution of elastic constants along rotational stress paths. A recent study by Geoffroy et al. (2003) sheds new light on this evolution by means of a new apparatus, a torsional 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, θ):

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

Controlling axial, horizontal, and torsion shear stresses allowed different loading conditions to be applied with and without stress rotation

$$[\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 (24)$$

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

$$\left\{ \begin{matrix} \varepsilon_r \\ \varepsilon_\theta \\ \varepsilon_z \\ \frac{\gamma}{\sqrt{2}} \end{matrix} \right\} = \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} \sigma_r \\ \sigma_\theta \\ \sigma_z \\ \frac{\sqrt{2}}{\tau} \end{bmatrix} \quad (25)$$

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_h/\sigma_v=0.5$ (Type K). The authors did not give 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. In this condition, the parameters for Hostun sand were adopted.

In their paper, Geoffroy et al. present experimental values only for the two last columns (eight elastic constants) of the elastic matrix in Eq. (25). Comparisons between experimental results and simulations are presented in Fig. 12, where the elastic constants are normalized to the value under initial loading condition. For the triaxial loading case, the elastic constants are plotted against normalized stress $\sigma_z/\sigma_{z0}$. For the torsional tests, the elastic constants are plotted against stress ratio $\tau/p_{mean}$. During triaxial loading (Case C) the terms in the last column $M_{r\gamma}$, $M_{\theta\gamma}$, $M_{\gamma z}$, and $M_{z\gamma}$ should be 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 model prediction (Fig. 12). During a torsion shear test, these four terms will evolve with the change in the principal stress directions. In the model, these two relations are always verified: $M_{r\gamma}=M_\theta$, $M_{\gamma z}=M_{z\gamma}$. Comparisons with experimental data show that the model is 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. 12. These two terms are always equal in the model prediction and one can see that, despite the dispersion of the experimental data, the model can predict correctly their evolution during triaxial tests, while their values remain practically constant during torsion shear tests.

Fig. 12 shows the evolution of the terms $M_{zz}$ and $M_{\gamma\gamma}$. It can be seen that the model can capture the decrease of these two terms during triaxial loading, which corresponds to an increase in Young and shear moduli with stress ratio. During shearing, these two terms slightly increase with the shear amplitude, corresponding to a decrease in Young and shear moduli. Here this tendency can also be predicted by the model.

In conclusion, the model is able to take into account with reasonable accuracy the evolution of the elastic matrix due to stress induced anisotropy, which may be caused by either stress amplitude or rotation of principal stress axes.

## Summary and Conclusion

An elastic model for particulate materials was developed based on micromechanics considerations. A particulate material is consid-
<figure><img src="image_13.png" /><figcaption>Fig. 12. Stress-induced anisotropy during loading tests with and without stress rotation (evolution of $M_{zr}$, $M_{\theta r}$, $M_{zy}$, $M_{yz}$, $M_{rz}$, $M_{\theta z}$, $M_{zz}$, and $M_{rr}$)</figcaption></figure>

ered as an assembly of particles. The stress–strain relationship for the assembly can be determined from integrating the behavior of interparticle contacts in all orientations using a static hypothesis which relates the average stress of the granular assembly to a mean field of particle contact forces. This stress–strain relationship leads to an elastic nonlinear behavior of the particulate material.

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_{r0}=\alpha k_{n0}$, as well as the packing density $\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 Poisson’s ratio. The influence of the void ratio of the assembly on the elastic properties is taken into account through the packing density $\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 assessed by taking into account a dependency of the parameters defining the normalized number of contacts in a given plane on the orientation of this plane. The nonlinear character of the contact stiffness leads to a stress-induced anisotropy with regard to the overall behavior of the grain assembly. This stress induced anisotropy includes the effects of the principal stress amplitudes as well as the rotation of the principal stress directions.

The predictions of the model for different stress conditions were compared to the results obtained during several experimental studies on different granular materials. The comparisons showed that the model is capable of taking into account very precisely the influence of the inherent anisotropy. The stress induced anisotropy along the different stress paths is also taken into account, but the model allows the transverse stresses to have a

more pronounced influence on Young and shear moduli. This is the result of the way the local equations are integrated. However, it must be pointed out that this discrepancy becomes noticeable for elevated stress ratios. In this regard, several studies have demonstrated that the change in the fabric due to strain-induced anisotropy could no longer be neglected (Kuwano and Jardine 2002; Dano and Hicher 2003). Under these conditions, a mechanism of strain-induced anisotropy needs to be introduced in the model, which was not the subject of this study. This could be achieved in future studies, given that an elastoplastic model based on the same numerical concepts as the ones exposed here has recently been developed (Chang and Hicher 2005), which allows nonreversible strains to be computed that could be used to take into account the fabric evolution and the related evolution of the elastic properties.

## References

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

* Biarez, J., and Hicher, P.-Y. (1994). *Elementary mechanics of soil behavior*, Balkema, Rotterdam, The Netherlands.

* Boyce, H. R. (1980). "A nonlinear model for the elastic behaviour of granular materials under repeated loading." *Proc., Int. Symp. Soils under Cyclic and Transient Loading*, Vol. 1, Swansea, U.K., 285-294.

* Brown, S. F., and Pappin, J. W. (1981). "Analysis of pavement with granular bases, layered pavement system." *Transportation Research Record 810*, Transportation Research Board, Washington, D.C., 17-23.

* Cambou, B., Dedecker, F., and Chaze, M. (2000). "Relevant local variables for the change of scale in granular materials." *Constitutive modeling of granular materials*, Dimitrios Kolymbas, ed., Springer, Berlin, 275-290.

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

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

Chang, C. S., Chao, S. J., and Chang, Y. (1995). "Estimates of elastic moduli for granular materials with anisotropic random packing structure." *Int. J. Solids Struct.*, 32(14), 1989-2008.

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 Hicher, P.-Y. (2005). "An elastoplastic model for granular materials with microstructural considerations." *Int. J. Solids Struct.*, 42, 4258-4277.

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

Chang, C. S., and Misra, A. L. (1990). "Application of uniform strain theory to heterogeneous granular solid." *J. Eng. Mech.*, 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.

Christofferson, J. S., Mehrabadi, M. M., and Nemat-Nassar, S. (1981). "A micromechanical description on granular material behavior." *ASME J. Appl. Mech.*, 48, 339-344.

Chung, R. M., Yokel, F. Y., and Drnevich, V. P. (1984). "Evaluation of dynamic properties of sands by resonant column testing." *Geotech. Test. J.*, 7(2), 60-69.

Dano, C., and Hicher, P.-Y. (2003). "Behavior of uncemented soils and grouted soils before maximum shear strength." *Soils Found.*, 43(4), 13-20.

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." *Proc., 3rd Int. Symp. on Deformation Characteristics of Geomaterials, ISLyon'03*, Lyon, France.

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." *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." *Behaviour of granular materials*, B. Cambou, ed., Springer, New York, 1-97.

Hicher, P.-Y. (2001). "Microstructure influence on soil behavior at small strains." *Prefailure deformation characteristics of geomaterials*,

M. Jamiolkowski, R. Lancellotta, and P. Lo Presti, eds., Swets & Zellinger, Lisse, The Netherlands, 1291-1297.

Hicks, R. G., and Monismith, C. L. (1971). "Factors influencing the resilient response of granular materials." *Transportation Research Record. 345*, Transportation Research Board, National Research Council, Washington, D.C., 15-31.

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 deformations of a granular material." *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." *Geotechnique*, 52(10), 727-749.

Lade, P. V., and Nelson, R. B. (1987). "Modelling the elastic behavior of granular materials." *Int. J. Numer. Analyt. Meth. Geomech.*, 11(5), 521-542.

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. (1969). "Microstructure in linear elasticity." *Arch. Ration. Mech. Anal.*, 16, 51-78.

Ng, T. T., and Petrakis, E. (1996). "Small-strain response of random arrays of spheres using discrete elements method." *J. Eng. Mech.*, 122(3), 239-244.

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

Rothenburg, L., and Selvadurai, A. P. S. (1918). "Micromechanical definitions of the Cauchy stress tensor for particular media." *Mechanics of structured media*, A. P. S. Selvadurai, ed., Elsevier, Amsterdam, 469-486.

Santha, B. L. (1994). "Resilient modulus of subgrade soils: Comparison of two constitutive equations." *Transportation Research Record. 1462*, Transportation Research Board, Washington, D.C., 213-226.

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