Evaluation of Two Homogenization Techniques for Modeling the Elastic Behavior of Granular Materials

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.


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.
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.

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 microscale 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 cth contact is defined by 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 t 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 For each particle contact, the corresponding auxiliary local coordinate system is related to the global coordinate system according to ͑see Fig. 1͒ n = ͑cos ␥,sin ␥ cos ␤,sin ␥ sin ␤͒ s = ͑− sin ␥,cos ␥ cos ␤,cos ␥ sin ␤͒ 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 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ͩ ͱ 12 ͑5͒ the normal contact stiffness in Eq. ͑4͒ is equivalent to Hertz's contact formulation.
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 ͑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., where 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 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.
Given the position vector of the contact point as x i c , the following may be written: 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 ␦ i c at the contact point c can be obtained as follows: Using the polynomial expression in Eq. ͑7͒, we can describe the interparticle compression using the continuum field of displacement.
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.
Then, by substituting the interparticle compression ⌬␦ i c with Eq. ͑11͒, it yields

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͒: 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 Utilizing Eq. ͑12͒ and substituting Eq. ͑14͒, the mean force on the contact plane of each orientation can be obtained as

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 ⌬␦ i c ͓Eq. ͑11͔͒.
For the case of static hypothesis, the stress-strain relationship can be derived using Eqs. ͑1͒, ͑14͒, and ͑16͒, it yields In the case of random packing with sufficient number of particles, the discrete interparticle contact orientations can be approximated by a continuous function E͑␤ , ␥͒, where ␤ and ␥ = two angles in the spherical coordinate system as shown in Fig.  1. Integration of E͑␤ , ␥͒ over the surface of unit sphere should equal the total number of interparticle contacts, N ͑i.e., selfconsistency equation͒.
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 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 ͑3 cos 2␥ + 1͒ + 3a 22 sin 2 ␥ cos 2␤ͪ

͑21͒
For a cross-anisotropic material, a 22 becomes equal to zero and Eq. ͑21͒ is reduced to and for an isotropic material, 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 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͑␤ , ␥͒ for a cross-anisotropic material can be expressed as ͑3 cos 2␥ + 1͒ͪ ͑25͒ Note that h m = mean value of h͑␤ , ␥͒ and the coefficient a 0h indicates the distribution of h͑␤ , ␥͒, 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.
For convenience, we define a packing density index ͑normalized number of contacts per unit volume͒, = Nl 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 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 n and the void ratio e of the assembly ͑Chang et al. 1989͒ A relationship between 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 and flexibility tensor for the static method 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 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͒: = Nl 3 / V; • Mean particle size, d; and • Interparticle elastic constants: G g , ␣, and n, where the particle modulus G g and exponent n have been defined in Eq. ͑4͒, and ␣ = 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 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͒: where = packing density index. Using Eqs. ͑4͒ and ͑5͒, under an isotropic stress c , the shear modulus G can be expressed as 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.
where G g = particle modulus. Based on the kinematic hypothesis, the shear modulus G and the Poisson's ratio take the following expressions: An equation similar to Eq. ͑32͒ can be obtained for G with the following expression for H g Ј: Eqs. ͑31͒ and ͑35͒ show that Poisson's ratios in both methods depend only on the interparticle stiffness ratio ␣. The expressions of the shear modulus are similar in character and differ only by the form of their dependency on ␣. ␣ is theoretically related to the Poisson's ratio of the particles g , assuming a perfectly elastic contact between the two spherical particles ͑Hertz's law͒: For a given value of the stiffness ratio ␣ = k t0 / k n0 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Ͻ g Ͻ 0.3. In these conditions, the computed values of the Poisson's ratio 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 of the assembly, Eqs. ͑31͒ and ͑35͒ can be used to determine, for each homogenization method, the value of the ratio ␣ = k t0 / k n0 . This means that ␣ will have a different value in the two models, according to the following equations: ␣ = ͑1−2͒ / ͑1+3͒ for the static hypothesis, and ␣ = ͑1−4͒ / ͑1+͒ 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 of the assembly Typical values for 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 60 / d 10 ͒. Thus the packing density index is no longer a unique relationship with the void ratio e as it was shown in Fig.   Fig. 5. Shear modulus and Poisson's ratio obtained from static hypothesis compared to those obtained with the kinematic hypothesis Fig. 6. Influence of grain size and distribution on Young's modulus ͑E and p are in Mpa͒ 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 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 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 ͑obtained from void ratio e and coefficients of uniformity U c ͒, and • Particle modulus G g , stiffness ratio ␣, 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 ␣ from the Poisson's ratio of the assembly ͓Eqs. ͑31͒ and ͑35͔͒. This will lead to different values of ␣ using either the static or the kinematic hypothesis. The packing density index 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 obtained for the material modulus ͑E or G͒. The parameters ͑ , 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 = 0.21 for both glass ballotini and Hostun Sand ͑Biarez and Hicher 1994; Hicher 1998͒ give a value of ␣ = 0.36 for the static hypothesis and ␣ = 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 , 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 ␣ = k t0 / k n0 is determined for each model from the experimental value of the Poisson's ratio 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 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 vh were obtained by applying a small vertical stress increment ⌬ v ͑=⌬ zz ͒ while keeping the horizontal stress h constant ͑i.e., ⌬ xx = ⌬ yy =0͒ Note that ⌬ v = ⌬ zz and ⌬ h = ⌬ xx = ⌬ yy . When applying a small horizontal stress increment ⌬ h ͑=⌬ xx = ⌬ yy ͒ while keeping v constant ͑i.e., ⌬ zz =0͒, the following relations can be used: One can observe that the two above equations for triaxial loading are not sufficient to determine the three unknowns: E h , hv , and hh . The assumption made by Hoque and Tatsuoka ͑1998͒ was hh = vh under isotropic stress states. Furthermore, due to the limited accuracy of lateral strain measurements, the writers believed that the measured values of 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.
In this study, the same assumption is made with regard to 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, which is assumed orientational dependent and has the same material axis as the fabric. For Toyoura Sand, the static hypothesis provided an anisotropy factor of the distribution of 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 = 2.5.
A Poisson's ratio vh = 0.17 was measured during isotropic loading and the same value was assumed for hh . Estimated from Eqs. ͑31͒ and ͑35͒, the corresponding value of ␣ = 0.45 and an overall isotropy ͑a 0␣ =0͒ was determined for the static hypothesis model, whereas a mean value of ␣ = 0.3 and an anisotropy factor a 0␣ = 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 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 ␣ was also introduced in order to obtain Poisson's ratio values in accordance with the measured ones. ͑For Hime Gravel, hh = 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   . 9. Isotropic loading on anisotropic Toyoura Sand and Hime Gravel 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 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.

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: Here we will consider only the case of a simple loading path; starting from an isotropic state of stress 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 11 / 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 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 12 in the case of the kinematic hypothesis, while 21 decreases at a similar rate. 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 v / 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 v and 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 v and h are also influenced by the value of the stress in the perpendicular direction ͑i.e., h for E v and 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. 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 v / 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 vh 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 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, ͒.
Different loading conditions can be applied, with and without stress rotation, by controlling axial, horizontal, and torsion shear stresses Therefore 16 terms of the flexibility matrix can be determined as follows: 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 v / 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␥ , M ␥ , M ␥z , and M z␥ 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 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 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 v / P a = 1.0-2.0. Fig. 12. Stress-induced anisotropy in Toyoura Sand and Hime Gravel 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 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 ␥␥ . 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 t0 = ␣k n0 , as well as the packing density index ͑normalized number of contact per unit volume͒ along a given plane orientation. The dependency of the normal stiffness on the normal force gives a pressuredependent behavior for the assembly, while the interparticle stiffness ratio ␣ 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 . An empirical relationship was proposed in order to link values of 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 ␣ 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 elevated stress ratios ͑ 1 / 2 Ͼ 2͒.