Elastic model for partially saturated granular materials

This paper presents the development of an elastic model for partially saturated granular materials based on micromechanical factor consideration. A granular material is considered as an assembly of particles. The stress-strain relationship for an assembly can be determined by integrating the behavior at all interparticle contacts and by using a static hypothesis, which relates the average stress of the granular assembly to a mean field of particle contact forces. As for the nonsaturated state, capillary forces at grain contacts are added to the contact forces created by an external load. These are then calculated as a function of the degree of saturation, depending on the grain size distribution and on the void ratio of the granular assembly. Hypothesizing a Hertz-Mindlin law for the grain contacts leads to an elastic nonlinear behavior of the particulate material. The prediction of the stress-strain model is compared to experimental results obtained from several different granular materials in dry, partially saturated and fully saturated states. The numerical predictions demonstrate that the model is capable of taking into account the influence of key parameters, such as degree of saturation, void ratio, and mean


Introduction
A granular material can be considered as a collection of particles of different sizes and shapes.This discontinuous medium can be represented by a fictitious continuous medium whose mechanical properties should depend on the geometric arrangement and contact interactions between interacting particles.For dry granular materials, the interparticle forces are related to the applied external stresses.When a certain amount of water is added to the grain assembly, significant changes, such as an increase in stiffness and strength, can be observed.These changes can be explained by the formation of water menisci between neighboring particles, which creates capillary forces.The amplitude of these capillary forces depends upon the degree of saturation and the morphology at particle level ͑grain and pore sizes͒ ͑Bishop and Blight 1963; Hicher 1998͒.A number of studies have addressed the determination of the capillary forces for a two-particle system ͑Fredlund and Rahardjo 1993; Lian et al. 1993;Taibi 1994;Cho and Santamaria 2001;Molenkamp and Nazemi 2003a͒.However, fewer results are available for more complex particle packing ͑Molenkamp and Nazemi 2003b; Dangla and Coussy 1998͒.Macroscopic analyses of the mechanical behavior of unsaturated granular materials in the elastic domain can be classified into two categories.Most of the data are analyzed on the basis of the double-variable theory by considering that the strain tensor is governed by the total net stress tensor ij − p g ␦ ij , p g being the gas pressure, and the negative pore water pressure or suction s ͑Fredlund and Rahardjo 1993; Alonso et al. 1990͒.However, several authors have shown that the experimental results could be interpreted using the concept of a generalized effective stress Ј by introducing a capillary stress c , function of the suction s ͑Biarez et al. 1993͒ The above equation can be considered as a general formulation of the equation proposed by Bishop and Blight ͑1963͒ for unsaturated soils Ј = − p g I + X͑S l ͒͑p g − p l ͒I ͑2͒ where p l = liquid pressure; and X͑S l ͒ϭfunction of the degree of saturation S l of the liquid phase ͑X = 0 for a dry material, X = 1 for a fully saturated material͒.Dangla and Coussy ͑1998͒ have confirmed theoretically the validity of the effective stress approach in the elastic domain.
An interesting method for granular material modeling is to develop homogenization techniques in order to derive a stressstrain relationship from forces and displacements at the particle level.Two different approaches can be followed based on the so-called kinematic hypothesis, which states that every particle displaces in accordance with a uniform deformation field ͑Chang 1988͒, or on the static hypothesis, which relates the average stress of the granular assembly to a mean field of particle contact forces ͑Cambou et al. 1995;Chang and Gao 1996;Liao et al. 1997͒.Hicher and Chang ͑2005͒ demonstrated the capability of these two approaches to provide a good approximation of the nonlinear elastic behavior of dry granular materials along different loading paths.In this study, the static hypothesis was considered and the capillary forces were introduced at each grain contact in order to extend the model to unsaturated granular media.The determina-tion of the elastic parameters of the granular assembly based on the characteristics of the particle contacts will be discussed and the predictions of the model along different stress loading will be compared to experimental results obtained from different granular materials, in dry, partially, or fully saturated states.

Micromechanical Approach
In this model, we envision a granular material as a collection of particles.The deformation of a representative volume of this material is generated by the mobilization of contact particles inside the whole volume.Thus, the stress-strain relationship can be derived as an average of the mobilization behavior at all interparticle contact planes in the representative volume.At the ␣th contact plane, the local forces f i ␣ and the local movements ␦ i ␣ can be denoted as follows: , where the subscripts n, s, and t represent the components in the three directions of the local coordinate system.The direction normal to the plane is denoted as n; the other two orthogonal directions, s and t, are tangential to the plane ͑Fig.1͒.
The forces and movements at all contact planes are suitably superimposed to obtain the macroscopic stress-strain tensors.The macroscopic stiffness tensor is obtained by suitably superimposing the stiffness properties at all interparticle contacts.In such a formulation, it has usually been assumed that the microstructure is statically constrained, which means that the forces acting on each contact plane are assumed equal to the resolve components of the macroscopic stress tensor.

Contact Forces and Capillary Forces
For dry samples, the contact forces are directly determined from the external stresses applied on the granular assembly.In the case of wetted samples, different stages of saturation can be identified.The fully saturated regime corresponds to a two-phase material with water filling completely the voids between grains.The water pressure p w can be positive or negative ͑suction͒, but in both cases the effective stress concept ͑Terzaghi 1925͒ can be applied and the contact forces determined by considering the effective stresses Ј as the external stresses ͑De Buhan and Dormieux 1996; Hicher 1998͒ In the case of partially saturated samples, the liquid phase is distributed in menisci located between close grains and as a consequence, capillary forces are applied on the grains and are added to the contact forces defined above.When the water content decreases inside a saturated granular sample, the air breaks through at a given point.The negative water pressure or suction corresponding to that point is called the air-entry pressure.Its value depends on the pore sizes.Afterwards, the sample becomes unsaturated with the water phase being connected inside the pores.This state is called the funicular regime, in which a constant decrease of the degree of saturation corresponds to a gentle decrease in pore water pressure, which is homogenized along the continuous water phase.The menisci start to form between two grains, not necessarily in contact, and are connected to each other.The pendular regime starts when the water menisci become disconnected.The water is no longer continuous and equilibrium in the water pressure is obtained by the vapor pressure.At this stage, the capillary forces start to increase drastically according to the pore size.If the drying process continues, the water bridges begin to fail, at first with the noncontacting grains, until a completely dry state is achieved and no capillary forces are left inside the granular assembly.
Several studies have shown that the attractive capillary forces between two grains connected by a water bridge are a decreasing function of the distance between the grains until the bridge fails ͓see, for example Biarez et al. ͑1993͒, Taibi ͑1994͒, and Chateau et al. ͑2002͔͒.This function depends on the volume of liquid found between the grains.Different mathematical expressions have been proposed for these capillary forces, which are the sum of the pressure forces exerted by the liquid on the wetted area and the surface tension forces acting along the wetted perimeter.In this study, we retained the following expression: where f cap = capillary force between two neighboring grains, not necessarily in contact; f max = value of f cap for two grains in contact; and R = mean grain radius; d = distance between two grains and is equal to l −2R, l being the branch length given as a distribution function of the grain size and the void ratio; c = material parameter, dependent on the grain morphology and on the water content; f max depends on the capillary pressure defined as the pressure jump across the liquid-air interface, on the liquid-air interface surface tension, as well as on the geometry of the menisci governed by the solid-liquid contact angle and the filling angle.In this study, a simplified approach was to consider an empirical relation between f max and the degree of saturation S r , without taking into account the hysteresis along drying and wetting paths ͑Hicher 1998͒ where f 0 and S 0 = material parameters; f 0 depends on the grain size distribution; S 0 = degree of saturation at which any further drying of the specimen will cause substantial breaking of the menisci in the pendular domain; S 0 depends on the nature of the granular material.Since the menisci are not necessarily all formed in the funicular regime, Eq. ͑5͒ may not be applicable for high degrees of saturation.However, in this first approach, we decided to extend it to the whole range of saturation, considering that the amplitudes of capillary forces were small for degrees of saturation higher than 80% and could, therefore, be approached with sufficient accuracy by using the same equation.

Interparticle Behavior
The contact stiffness of an orientation includes normal stiffness k n ␣ and shear stiffness k t ␣ of the contact plane.The elastic stiffness tensor is defined in where k ij ␣ = interparticle contact stiffness tensor; ␦ j ␣ = relative displacement at the interparticle contact.
Let k n be the compressive contact stiffness in normal direction and k s be the shear contact stiffness.Assuming the shear contact stiffness to be the 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 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 ␤͒ 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 crossproduct 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 can be adopted ͑Chang et al. 1989͒, given by where f n = contact force in normal direction; f ref = reference force to make the units consistent.In this paper, f ref is selected to be 1 N. k n0 , k t0 , and n = material constants.For two spherical particles, if the parameters are chosen to be where G g and g = respectively, the shear modulus and the Poisson's ratio for the grains; and r = mean radius of the grains.Eq. ͑9͒ is 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 interparticles at all contacts.In the integration process, a micro-macro relationship is required.
For simplicity, the effects of particle spin and the change of interparticle orientations are not considered during a small strain increment ͑i.e., finite-strain condition is not considered͒.Using the static hypotheses proposed by Liao et al. ͑1997͒, we obtain the relation between the strain and interparticle displacement where the branch vector l k ␣ is defined as the vector joining the centers of two particles; ␦ j ␣ = relative displacement at the interparticle contact ͓see Eq. ͑6͔͒; and the fabric tensor is defined as The length of branch vectors for the unsaturated case may be greater than the summation of the two particle radius due to the film of water between particles.Thus, the branch vectors do not have the same value for all particle pairs.It is noted that in Eqs.͑11͒ and ͑12͒, the summation is not limited to direct contacts.The summation should also include indirect contacts of neighboring particles within a Voronoi polygon, as discussed by Cambou et al. ͑2000͒.
Using Eq. ͑11͒ and invoking the condition that the rate of energy dissipation expressed in terms of the assembly stress and strain must be equivalent to that expressed in terms of interparticle forces and movements, we are seeking to establish the expression of interparticle forces in terms of assembly stress.It is obvious that there exist an infinite number of solutions that can satisfy the conditions stated above.However, without invoking further constraints, it can be proved that the following relationship is always one of the possible solutions ͑Liao et al. 1997͒: Eq. ͑13͒ also satisfies the stress definition in terms of interparticle contact forces ͑Cauchy theory of stress mean͒ Because of its approximation nature, Eq. ͑13͒ can be viewed as an averaged solution.The interparticle force can be regarded as the mean value for forces on all contact planes of the same orientation.For convenience, we also regard the branch vector as the mean value for all contact planes of the same orientation.
Using the definition of Eq. ͑14͒, the stress induced by capillary forces can be computed and is termed as capillary stress, given by It is noted that this term is not analogous to the usual concept of capillary pressure or suction, which represents the negative pore water pressure inside the unsaturated material ͑Bishop and Blight 1963͒.In this model, the capillary stress depends on the geometry of the pores and is a tensor rather than a scalar.Only for an isotropic distribution of the branch lengths l ␣ , the capillary stress can be reduced to an isotropic tensor.This can be the case for an initially isotropic structure during isotropic loading, but during deviatoric loading, an induced anisotropy is created and the capillary tensor is no longer isotropic.
By adding the contact forces due to the external stresses ͓Eq.͑14͔͒ and the capillary forces ͓Eq.͑15͔͒, we can define a generalized effective stress tensor * as

͑16͒
This equation represents a generalization of Eq. ͑2͒ proposed by Bishop, in which the capillary stress is reduced to an isotropic tensor.Dangla and Coussy ͑1998͒ demonstrated the validity of the effective stress approach in elasticity by means of an energy approach.They obtained an expression similar to Eq. ͑2͒, but with an additional term corresponding to the work of the interfaces.As pointed out before, the capillary forces in our model, and as a consequence, the capillary stresses depend on the negative pore water pressure, or suction, and on the water-air interface surface tension.A similar approach can be found in the work of Fleureau and et al. ͑2003͒ who obtained an explicit expression of the capillary stress as a function of the suction for regular arrangements of spherical grains by neglecting the surface tension.The definition of the capillary stress tensor in Eq. ͑15͒ can, therefore, be considered as an extension of the results obtained from these previous studies to cases of nonisotropic granular assemblies.

Stress-Strain Relationship
Using Eqs.͑6͒, ͑11͒, and ͑13͒, the following relationship between stress and strain can be obtained: For elastic behavior, Eq. ͑17͒ can lead to a closed-form solution in the case of a random packing of equal-size particles ͑Liao et al. 1997͒.It is noted that for contacts in the same orientation, the branch vector and contact stiffness are their mean values.
The orientation of each interparticle contact can be expressed by two angles ␤ and ␥ in a spherical coordinate system as shown in Fig. 1.The orientational distribution of interparticle contacts can be described by a probability density function E͑␤ , ␥͒.In the case of random packing with sufficient number of particles, the discrete probability density function can be approximated by a continuous probability density function.Integration of E͑␤ , ␥͒ over the surface of unit sphere should equal to one ͑i.e., selfconsistency equation͒ For isotropic packing structure, E͑␤ , ␥͒ =1/ 2. For any arbitrary parameter h c , which represents interparticle contact property, the summation of h c over all interparticle contacts can be approximated by an expression of an integral, given by 1 where N = total number of interparticle contacts; and h͑␤ , ␥͒ = continuous expression of h c .Using the definition of Eq. ͑19͒ and assuming an isotropic packing structure, 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.
A relationship between n ¯and e obtained for a regular packing of spheres is plotted in Fig. 2. For an isotropic packing structure, using the packing density index, we can express the flexibility tensor Although the packing structure is isotropic, the contact stiffness can be orientational dependent because the contact forces are orientational dependent ͓see Eq. ͑9͔͒.Thus, the assembly can exhibit anisotropic behavior due to the anisotropy of applied stresses.
A numerical integration method can be used for the integrals in Eqs.͑17͒ and ͑22͒.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 weighing factor.We found that the results were more accurate by using a set of fully symmetrical integration points.A study of the performance of using different numbers of orientations showed 37 points to be adequate ͑Chang and Hicher 2005͒.

Elastic Properties of Dry Granular Materials
As a first step, we will describe the procedure used to model the elastic behavior of dry granular materials.If one assumes an isotropic fabric for the material, the parameters of the proposed model are the following: • packing density index: Nl 3 / V. • mean particle size, D.
• interparticle elastic constants: k n0 , k t0 , and n.The unit volume V and the mean particle size D are in fact related, so that only one of these two values has to be considered.
The number of contacts per unit volume governs the density of the packing.If we assume that the same relationship as the one presented in Fig. 2 for a regular packing of spheres can be roughly applied to irregular packing ͑Liao et al. 2000͒, a relationship can be obtained between N / V and void ratio e, for R =1.
This relationship, applicable to spherical grains, has to be adapted for grains with irregular shapes, which can only be done by an empirical method.Based on the study of different granular materials behavior ͑Hicher 2001͒, a chart, giving an empirical relationship between N / V and void ratio e for various coefficients of uniformity Uc = d 60 / d 10 , was proposed ͑Hicher and Chang 2005͒ and is presented in Fig. 3.
One can, therefore, determine theoretically the shear modulus of an assembly of particles, knowing the elastic properties ͑G g , g ͒ of the particles, the grain size distribution, and the void ratio ͓Eq.͑22͔͒.In fact, Chang and Gao ͑1996͒ have demonstrated that static approaches underestimate the true stiffness of the packing.As a consequence, the approach hereby adopted was to determine the parameters of the model through the following procedure.From the test results, we could directly determine the coefficient of nonlinearity n as well as the ratio ␣ = k t0 / k n0 , which is directly related to the Poisson's ratio of the assembly.The ratio N / V could be determined from Fig. 3.The value of k n0 could, therefore, be determined by fitting predicted and experimental values obtained for the material stiffness ͑E or G͒.A typical value of k n0 = 2,420 was derived by fitting predicted and experimental values in the case of a glass ballotini assembly as well as for quartzic sands.Using this procedure, Hicher and Chang ͑2006͒ showed the capacity of the present model to reproduce with very good accuracy the elastic behavior of different granular materials.
The same procedure was used in this work to determine the parameters of the model for the studied granular materials in dry state.The grain size distribution curves of the studied materials are presented in Fig. 4.

Elastic Properties of Unsaturated Granular Materials
Resonant column tests made on fine-grained cohesionless soils by Wu et al. ͑1984͒ showed the influence of the capillary stresses on the values of the shear modulus G.The unsaturated samples were compacted in a mold that was placed on the pedestal of the resonant column device, at a desired density and water content.Dry samples were prepared with a vibrator to obtain the required density.Fully saturated samples were prepared by introducing distilled water into prepared dry samples.For all the tested soils, the authors found a significant increase in the shear modulus with the decrease of the degree of saturation.A maximum influence was obtained for degrees of saturation between 5 and 20%, according to the nature of the soil.For lower degrees of saturation, a decrease of the shear modulus was observed, and completely dry samples had the same modulus as the fully saturated samples.The authors proposed a correlation of the optimum degree of saturation, giving the maximum shear modulus, with the grain size, represented by the d 10 value in mm ͑S r ͒ opt = ͑− 6.5 log͑d 10 ͒ + 1.5͒/100 ͑23͒ in which ͑S r ͒ opt is equivalent to S 0 in Eq. ͑5͒, d 10 being the effective grain size in mm.Fig. 5 presents the results obtained for the glacier way silt for two different confining pressures ͑see grain size distribution curve in Fig. 4͒.The value of G was the same for S r = 0% and S r = 100%.A well-marked peak for the two confining pressures could be observed, which occurred for a value of S r around 17%.The maximum value of G is 1.8 times its value for dry and saturated states at a confining pressure of 50 kPa and 1.6 times at a confining pressure of 100 kPa.The capillary effect appears more pronounced at lower confining pressures.The general trend showing the evolution of the shear modulus with the degree of saturation is in agreement with the model concepts.That implies an identical value of G at dry and fully saturated states in the absence of capillary forces and a maximum influence of these forces at an optimum degree of saturation, as given by Eq. ͑5͒ in which S 0 corresponds to ͑S r ͒ opt in Eq. ͑23͒.To determine the capillary forces, Eq. ͑4͒ includes two material parameters c and d.Very little information could be gathered on these two parameters by simply using the test results.Therefore, we decided to take a constant value c = 8 and a value d function of the mean grain size d 50 : d = 1.1d50 .These values give a standard evolution for the capillary forces, function of the distance between particles, as well as an initial isotropic distribution of these forces.For the numerical simulations, the following parameters had to be determined: • The parameters of the elastic nonlinear law: k n0 , k t0 , and n.
• The parameters of the capillary forces equation: f 0 and S 0 .
k n0 and n were determined from the values of G obtained at dry ͑or saturated states͒.A ratio k t0 / k n0 = 0.35 was assumed, leading to a Poisson's ratio n = 0.2, which corresponds to a standard value for this type of material ͑Hicher 1996͒.
The value of N / V was determined by using the chart in Fig. 3.For a material with U c Ͼ 5, a void ratio e = 0.58 corresponds to N / V = 0.31.
S 0 was considered equal to ͑S r ͒ opt in Eq. ͑23͒ for d 10 = 0.0024 mm.f 0 was computed in order to fit the maximum value of G for the confining pressure of 100 kPa.The material parameters for Glacier Way Silt are summarized in Table 1.
The numerical results are presented in Fig. 5 together with the experimental ones.A good overall agreement was obtained for the two confining pressures.Starting from the fully saturated state, the shear modulus values increase first slowly for degrees of saturation varying between 100 and 70% and then faster up to the maximum values obtained for S r = S 0 .The computed peak value of G was close to the experimental one for the confining pressure of 50 kPa, showing a more significant effect of the capillary forces at low confining pressure.For degrees of saturation lower than S 0 , the shear modulus values decrease toward the values obtained for dry states.
Cho and Santamaria ͑2001͒ also measured the elastic properties of granular materials, using, however, a different preparation technique.The samples were prepared in a oedometric cell.They were mixed with distilled water at saturated conditions and gradually poured into the cell.They were then subjected to gradual drying in an incubator at 50°C.Shear wave velocity was measured continuously during drying by means of bender elements placed on top and bottom platens.In this study, we concentrated on two of the tested granular materials: monosize glass beads and sandboil sand ͑see grain size distribution curves in Fig. 4͒.The results obtained on glass beads are presented in Fig. 6.They show a continuous increase of the shear modulus with the decrease of the degree of saturation up to a peak value obtained for S r = 0.007 and then a sharp decrease as drying continues.The G value at the dry state remains slightly higher than the one obtained at a fully saturated state.The authors explain the difference by an increase in internal coordination during drying and residual compressive stress.The value of S r at peak shear modu-lus was found to be equal to 0.007, which differs from the value computed from Eq. ͑23͒: ͑S r ͒ opt = 0.053.This could be due to the mode of preparation, which is different in the two cases and could lead to different specimen states at small degrees of saturation.For a specimen prepared at a given low water content, menisci might not form between two wet particles because the water film is not thick enough, whereas during drying, the menisci may subsist longer before starting to fail when the water continues to evaporate.
For the numerical simulations, k n0 and n were determined from the values of G obtained at saturated state, using previous studies on elastic properties on glass beads ͑Hicher and Chang 2006͒.S 0 was chosen equal to 0.007 and f 0 was determined by fitting the G peak value ͑Table 2͒.
The experimental testing was performed under a very small applied vertical stress: 1.5 kPa.Numerical simulations were performed at a higher stress ͑10 kPa͒ to prevent numerical errors.As a result, we decided to present the evolution of the ratio G / G 0 , rather than the absolute values of G, as function of S r .A comparison of experimental and numerical results in Fig. 6 shows an overall agreement in the pendular domain ͑S r Ͻ 5%͒ except for the very dry state where the computed ratio G / G 0 converges toward 1, while it remains around 1.4 in the test results.However, the numerical values are smaller than the experimental ones in the funicular domain.We can see in Fig. 6 that the difference is mainly due to the fact that the experimental values increase significantly between S r = 100% and S r = 80%: for S r = 80%, G / G 0 = 1.31, while this increase remains more limited between 80 and 10%: for S r = 10%, G / G 0 = 1.51.A similar trend was obtained in the numerical simulations in this intermediate range of degrees of saturation, while the initial strong increase failed to appear.Its physical meaning lies in the drying process, which follows a mechanism often termed "ink bottle."When water is progressively removed from the sample, the pore space remains saturated as long as the negative water pressure is not big enough to allow the formation of a meniscus.The air entry point corresponds to the suction value for which an air bubble will form inside the pore.When the material is highly heterogeneous, this phenom-  enon is progressive, as water is removed first from the biggest pores and then from smaller and smaller pores during the drying process.On the contrary, for monosized grain assemblies, the air entry point is well marked and corresponds to a given value of the negative water pressure and function of the grain size.In the studied case, the mean grain size was equal to 0.32 mm with a maximum value of 0.5 mm and a minimum value of 0.2 mm.A theoretical study by Taibi ͑1994͒ gave the following relation between the pressure of desaturation and the grain size for monosized assemblies: where u wd = pressure of desaturation; d = diameter of the grains; and b = parameter that depends on the density of the assembly.For a void ratio e = 0.601 and 200 m Ͻ d Ͻ 500 m, we obtain: 1.2 kPaϽ −u wd Ͻ 3.1 kPa.This small value of the suction is large enough to affect the material stiffness, due to the very small value of the confining stress applied to the specimen.Since we performed our numerical simulations at larger confining stress, the effect cannot be as marked, even if we replace the net applied stress by the effective stress, adding to the effect of the suction.Therefore, we decided to perform another series of numerical simulations, starting from the value of the shear modulus corresponding to S r = 80%, as a reference point.The f 0 value was readjusted ͑Table 3͒ and the numerical results are presented in Fig. 6.They are, this time, in much better agreement with the experimental values in the range 0 Ͻ S r Ͻ 80%.
The same kind of test was conducted on a natural sand: The sandboil sand.For this material, the shear modulus increased continuously and showed no drop even at perfectly dry conditions ͑Fig.7͒.The authors explained the absence of drop at dry states as a consequence of chemical reactions, such as salt precipitation and bonding, at particles contacts.Our model is not elaborate enough at this stage to be able to deal with this kind of phenomenon.Therefore, we decided to choose a small value for S 0 , typically S 0 =10 −3 and to investigate only the influence of degrees of saturation higher than S 0 .The same procedure, as explained before, was used to determine the material parameters, which are presented in Table 4. Fig. 7 shows the comparison between experimental and numerical results for sandboil sand.A good agreement was obtained, in particular a gentle increase of G / G 0 for degrees of saturation higher than 0.2 in the funicular regime and then a strong increase for lower values of S r in the pendular regime.The value of the negative pore water pressure at the air entry point was probably small enough in this case, due to the existence of bigger pores created by a more heterogeneous size distribution, to prevent any noticeable influence on the stiffness values at high degrees of saturation.
Hadiwardoyo ͑2002͒ reported results on an Indonesian silty sand.The material is classified as a silty sand with spread distribution ͑see grain size distribution curve in Fig. 4͒.The samples were all prepared according to the same procedure: The soil was mixed with a given quantity of water and then compacted to the chosen water content and density.The elastic properties were measured in a precision triaxial cell, equipped with local sensors for measuring axial and radial strains.The enhanced precision of the transducers allows the measurement of strains as small as 10 −6 .Tests were performed at different confining pressures on samples having different void ratios and water contents.Fig. 8 presents the evolution of the Young's modulus E with the confining pressure p 0 for dry specimen with void ratio equal to 0.74 and unsaturated specimen with void ratio varying from 0.48 to 0.7 and water content from 8.3 to 16%, corresponding to various degrees of saturation from 33 to 90%.
The elastic parameters were determined by using the experimental relation between E and p for dry samples ͑Table 5͒.Comparison between experimental and computed values is shown in Fig. 8.The test results did not allow us to determine S 0 .Using the results by Wu et al. ͑1984͒ and Eq.͑23͒, a typical value of S 0 = 0.26 was selected for d 10 = 0.2 m. f 0 was then determined from the test results on unsaturated specimen at a given confining pressure ͑p 0 = 52 kPa͒.The set of parameters for the Indonesian silty sand is presented in Table 5. Numerical results, as well as experimental ones, are presented in Fig. 8.A good overall agreement could be obtained.For unsaturated samples at a given void  ratio and water content, the relation between Young's modulus and confining pressure in a bilogarithmic coordinate plane is linear.The slope value decreases with the degree of saturation down to 0.3 for S r = 33% and remains smaller than the one obtained on dry samples, which was equal to 0.67.This decrease is due to the increased influence of the capillary forces compared to the applied forces when S r decreases.Using Eq. ͑15͒, we can calculate the capillary stress tensor for each initial state.The material is assumed to have an isotropic structure and, therefore, the capillary stress tensor is also isotropic.The capillary stress is a decreasing function of the degree of saturation, which depends on the grain size distribution curve.
Using the definition of the generalized effective stress tensor in Eq. ͑16͒, we can now plot the relation between the Young's modulus E and the generalized mean effective stress for all the tests.The results are presented in Fig. 9.One can see that, in this plane, the relationship between log͑E͒ and log͑p*͒ is linear and depends only on the void ratio e.This relationship can be written n being equal here to 0.67.The Young's modulus appears, therefore, to be a function of the generalized mean effective stress and of the void ratio.
In conclusion, the model is able to take into account the evolution of the elastic properties of a granular material with the degree of saturation varying from a perfectly dry state to a fully saturated state.The determination of the material parameters requires to know the elastic coefficients at dry or saturated state, assuming that they are equal in both cases.If we examine the model parameters obtained for each material in Tables 1-5, it appears that the parameters of the contact law k n0 , k t0 , and n have very similar values, except for n in the case of the Indonesian silty sand, which has a significantly higher value than the standard value n = 0.5 obtained for the other materials.This difference can be explained by the higher content of clay material ͑20% in weight Ͻ2 m͒.
The model also requires the determination of the value of the parameters controlling the capillary force amplitude S 0 and f 0 .These can be obtained by knowing the values of the elastic coef-ficients of the granular material at a given degree of saturation different from 0 and 1, and the degree of saturation corresponding to the peak stiffness value.If this last parameter is not determined experimentally, an approximate value can be derived from the grain size distribution using Eq.͑23͒.We have seen, however, that this parameter can be sensitive to the mode of preparation.For the two materials prepared by compaction at a given water content, i.e., Glacier Way silt and Indonesian silty sand, Eq. ͑23͒ can be applied and leads to values of S 0 function of the grain size of the smallest particles d 10 .In these conditions, f 0 is also a function of the grain size distribution and increases when d 10 decreases: f 0 = 0.02 N for d 10 = 0.0024 ͑glacier way silt͒ and f 0 = 0.1 N for d 10 = 0.2 mm ͑Indonesian silty sand͒.This shows the influence of size of the smallest pores on the amplitude of the capillary forces.The same conclusions can be drawn by comparing the S 0 and f 0 values determined for the two materials prepared by progressive drying.Eq. ͑23͒ cannot be applied to determine S 0 , since the degree of saturation corresponding to the maximum stiffness appears to be smaller in this case than the computed one.However, one can see that S 0 is here also an increasing function of d 10 : S 0 = 0.007 for d 10 = 0.26 mm ͑glass beads͒ and S 0 = 0.001 for d 10 = 0.17 mm ͑sandboil sand͒.In the same way, f 0 is a decreasing function of d 10 : f 0 = 0.025 N to 0.045 N for d 10 = 0.26 mm and f 0 = 0.051 N for d 10 = 0.17 mm.
In the case of uniform size grain material, such as the assembly of glass beads, the negative pore pressure corresponding to the air entry point needs to be considered and the calculation of the elastic coefficients has to be performed using the effective stress as long as the sample remains saturated.

Summary and Conclusion
An elastic model for granular materials has been developed based on micromechanics considerations.The material is considered 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, given a static hypothesis, which relates the average stress of the granular assembly to a mean field of particle contact forces.For each contact plane, an elastic behavior based on the Hertz-Mindlin contact formulation was assumed.Under these conditions, the stress-strain relationship leads to an elastic nonlinear behavior of the 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 t0 = ␣k n0 , as well as the number of contacts per unit volume N / V 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 ratio ␣ determines the value of Poisson's ratio.The influence of the density or void ratio of the assembly on the elastic properties is taken into account during the integration procedure by the ratio N / V.An empirical relationship has been proposed in order to link N / V values to the void ratio.This relationship was based on results obtained for different granular materials with different grain size distributions.
For dry samples, the contact forces are directly determined from the external stresses applied on the granular assembly.In the case of partially saturated samples, the liquid phase is distributed in menisci located between close grains and, as a consequence, capillary forces are applied on the grains and added to the contact forces.The capillary forces were computed as a function of the pore morphology and of the water content, using a simple empirical relationship.Further improvement of the model will require a better understanding of the capillary force amplitude and evolution, based more closely on the physics of the phenomenon.
Comparison between numerical and experimental results on several granular materials demonstrated the ability of the model to reproduce the main characteristics of the elastic behavior of an unsaturated granular material.In particular, the model was able to reproduce the evolution of the material stiffness with the degree of saturation.The integration of the capillary forces over all spatial orientations permitted a definition of a capillary stress tensor.A generalized effective stress tensor could, therefore, be determined by adding the total stress and the capillary stress tensors.It was demonstrated that the evolution of the Young's modulus obtained on a given granular material at different water contents was controlled by the generalized effective stress.

Fig. 2 .
Fig. 2. Relationship between void ratio and coordination number

Fig. 3 .Fig. 4 .Fig. 5 .
Fig. 3. Empirical relationship between contact density N / V and void ratio as a function of grain size distribution

Fig. 6 .
Fig. 6.Evolution of shear modulus with degree of saturation for glass beads

Table 1 .
Model Parameters for Glacier Way Silt

Table 2 .
Model Parameters for Glass Beads ͑Simulation 1͒

Table 3 .
Model Parameters for Glass Beads ͑Simulation 2͒

Table 4 .
Model Parameters for Sandboil Sand Fig. 7. Shear modulus versus degree of saturation for sandboil sand

Table 5 .
Model Parameters for Indonesian Silty Sand Fig. 9. Young's modulus versus generalized effective stress for Indonesian silty sand