P.-Y Hicher

C S Chang

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

Keywords: Granular materials, Stress strain relations, Anisotropy, Constitutive models, Elasticity, Micromechanics, Homogeneity

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; [START_REF] Chang | Micromechanical modeling of constructive relations for granular material[END_REF][START_REF] Chang | Initial moduli of particulated mass with frictional contacts[END_REF][START_REF] Rothenburg | Analytical study of induced anisotropy in idealized granular materials[END_REF]. A second model was developed using a static hypothesis ͑Cambou et [START_REF] Cambou | Homogenization for granular materials[END_REF][START_REF] Chang | Kinematic and static hypotheses for constitutive modeling of granulates considering particle rotation[END_REF][START_REF] Liao | Stressstrain relationships for granular materials based on the hypothesis of best fit[END_REF], 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; [START_REF] Chang | Micromechanical modeling of constructive relations for granular material[END_REF][START_REF] Rothenburg | Analytical study of induced anisotropy in idealized granular materials[END_REF][START_REF] Jenkins | Volume change in small strain axisymmetric deformations of a granular material[END_REF][START_REF] Chang | Estimates of elastic modulus for media of randomly packed granules, Part 2[END_REF][START_REF] Cambou | Homogenization for granular materials[END_REF]. 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

⌬f i c = k ij c ⌬␦ j c ͑1͒
where k ij c = interparticle contact stiffness tensor. For two particles in contact, a local coordinate system can be constructed for each contact with three orthogonal base unit vectors: n is normal to the contact plane; s and t are tangential to the contact plane as shown in Fig. 1.

Let k n c be the compressive contact stiffness in normal direction and k 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

k ij c = k n c n i c n j c + k i c ͑s i c s j c + t i c t j c ͒ ͑ 2͒
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 ␤͒ t = ͑0,-sin ␤,cos ␤͒ ͑ 3͒
The vector s is located on the plane consisting of x and n. The vector t is perpendicular to this plane and can be obtained by the cross product of n ϫ s. The rolling resistances between two particles are not discussed in this paper.

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

k n = k n0 ͩ f n G g l 2 ͪ n ; k t = k t0 ͩ f n G g l 2 ͪ n

͑4͒

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

k n0 = G g rͩ ͱ 12 1 -g ͪ 2/3
͑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.,

u i ͑x n ͒ = u i n ͑6͒
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

u i ͑x͒ = u i + u i,j x j ͑7͒
where u i , u i,j = constants for the representative volume. The second-rank tensor u i,j has nine components. In the further derivation, for the sake of convenience, we define the branch vector l i c as the vector from the centroid of particle "a" to that of particle "b" as shown in Fig. 2.

l i c = x i b -x i a ͑8͒
Given the position vector of the contact point as x i c , the following may be written:

x i a = x i c -r i ac ; x i b = x i c -r i bc ͑9͒
Thus l i c = r i ac -r i bc . Between the two individual convex-shaped particles, as shown in Fig. 2, the interparticle compression ␦ i c at the contact point c can be obtained as follows:

␦ i c = ␦ i ab = u i ͑x b ͒ -u i ͑x a ͒ ͑ 10͒
Using the polynomial expression in Eq. ͑7͒, we can describe the interparticle compression using the continuum field of displacement.

␦ i c = u i,j l j c or ⌬␦ i c = ⌬u i,j l j c ͑11͒
The representative volume V is subjected to forces on its boundary surface S. Neglecting the body force of particles, the work done per unit volume of the discrete system can be expressed as a summation of the work done over all interparticle contacts in the unit volume.

⌬W = 1 V ͚ V f i c ⌬␦ i c = iq ⌬u q,i ͑12͒
Then, by substituting the interparticle compression ⌬␦ i c with Eq. ͑11͒, it yields

iq = 1 V ͚ V f q c l i c or ⌬ iq = 1 V ͚ V ⌬f q c l i c ͑13͒

Micro-Macro Relationship: Static Hypothesis

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

⌬u j,i = A ik -1 ͚ c=1 N ⌬␦ j c l k c ͑14͒
where the branch vector l k c is defined as the vector joining the centers of two particles, and the fabric tensor is defined as

A ik = ͚ c=1 N l i c l k c ͑15͒
Utilizing Eq. ͑12͒ and substituting Eq. ͑14͒, the mean force on the contact plane of each orientation can be obtained as

f j c = ij A ik -1 l k c V or ⌬f j c = ⌬ ij A ik -1 l k c V ͑16͒

Stress-Strain Relationship

For the case of kinematic hypothesis, the stress-strain relationship can be derived using the interparticle contact law ͓Eq. ͑1͔͒, the stress definition ͓Eq. ͑13͔͒, and the continuum expression of ⌬␦ i c ͓Eq. ͑11͔͒.

⌬ iq = C iqkl ⌬u k,l where C iqkl = 1 V ͚ c=1 N l i c l l c k qk c ͑17͒
For the case of static hypothesis, the stress-strain relationship can be derived using Eqs. ͑1͒, ͑14͒, and ͑16͒, it yields

⌬u i,j = C ijmp -1 ⌬ mp where C ijmp -1 = A ik -1 A mn -1 V ͚ c ͑k jp c ͒ -1 l k c l n c ͑18͒
In the case of random packing with sufficient number of particles, the discrete interparticle contact orientations can be approximated by a continuous function E͑␤ , ␥͒, 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͒.

N = ͵ 0 /2 ͵ 0 2

E͑␤,␥͒sin ␥d␤d␥ ͑19͒

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

͚ c=1 N h c = ͵ 0 /2 ͵ 0 2 E͑␤,␥͒h͑␤,␥͒sin ␥d␤d␥ ͑20͒
where N = total number of contact orientations; E͑␤ , ␥͒ = distribution function of interparticle contact orientations; and h͑␤ , ␥͒ = continuous expression of h c .

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

E͑␤,␥͒ = N 2 ͩ1+ a 0 4 ͑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

E͑␤,␥͒ = N 2 ͩ1+ a 0 4 ͑3 cos 2␥ + 1͒ͪ ͑22͒
and for an isotropic material,

E͑␤,␥͒ = N 2 ͑23͒
In three dimensions, the inherent anisotropy can be represented by a distribution whose major axis often coincides with the vertical direction. An example of distributions with different values of a 0 is shown in Fig. 3. The anisotropy of the packing structure creates a mechanical behavior anisotropy whose directions are identical to those of the geometric anisotropy.

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

͓F ij ͔ = ΄ 1 + a 0 0 0 0 1 - a 0 2 + 3a 22 0 0 0 1 - a 0 2 -3a 22 ΅ ͑24͒
In addition to the contact distribution, other parameters, which are orientation dependent, can also be expressed by the same manner. For example, the function h͑␤ , ␥͒ for a cross-anisotropic material can be expressed as 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.

A ik = ͚ c=1 N l i c l k c = Nl 2 2 A ¯ik
where

A ¯ik = ͵ 0 /2 ͵ 0 2 E͑␤,␥͒n i n k sin ␥d␤d␥ ͑26͒
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͒

= Nl 3 V = 3n ͑1 + e͒ ͑27͒
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

C iqkl = d ͵ 0 /2 ͵ 0 2 E͑␤,␥͒k qk n i n l sin ␥d␤d␥ ͑28͒
and flexibility tensor for the static method

C ijmp -1 = d A ¯ik -1 A ¯mn -1 ͵ 0 /2 ͵ 0 2 E͑␤,␥͒k jp -1 n k n n sin ␥d␤d␥ ͑29͒
A numerical integration method can be used for the integrals in Eqs. ͑28͒ and ͑29͒. The integrals are computed based on a set of integration points on the surface of a sphere. Each integration point represents a contact orientation and has a weighting factor. We found that the results are more accurate if the locations of integration points on the surface of a sphere are symmetric about the x-z plane and y-z plane ͑i.e., fully symmetric͒. We studied the performance of this model using different numbers of integration points and found that 74 integration points is adequate.

Isotropic Loading

Isotropic Fabric

We will at first assume an isotropic fabric for our particulate material. Under this condition, the parameters of the proposed model are the following: • Packing density index ͑normalized number of contacts per unit volume͒: = 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͒:

G = 5k n ␣ 6d͑3 + 2␣͒ ͑30͒ = 1 -␣ 2 + 3␣ ͑31͒
where = packing density index. Using Eqs. ͑4͒ and ͑5͒, under an isotropic stress c , the shear modulus G can be expressed as

G = H g 1-n ͩ c p 0 ͪ n

͑32͒

where p 0 = reference pressure ͑p 0 = 1 MPa in the computation in this paper͒. H g = function depending on elastic property of the grain.

H g = 5␣ 2 ͱ 6͑3 + 2␣͒ G g ͩ 3p 0 G g ͪ n

͑33͒

where G g = particle modulus.

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

G = k n 30d ͑2 + 3␣͒ ͑ 34͒ = 1 -␣ 4 + ␣ ͑35͒
An equation similar to Eq. ͑32͒ can be obtained for

G G = H g Ј 1-n ͩ c p 0 ͪ n ͑36͒
with the following expression for H g Ј:

H g Ј = 2 + 3␣ 10 ͱ 6 G g ͩ 3p 0 G g ͪ n ͑37͒
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͒:

␣ = 2͑1 -g ͒ 2 -g ͑38͒
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

H g = H g Ј = 1 -2 2 ͱ 6͑1 + ͒ ͑39͒
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. 6͑b͒ shows results obtained on different quartzic sands ͑Hardin and Black 1966; [START_REF] Hardin | Shear modulus and damping in soils: Measurements and parameter effects ͑Terzaghi Lecture͒[END_REF][START_REF] Skoglund | Evaluation of resonant column test devices[END_REF][START_REF] Iwasaki | Effects of grain size and grading on dynamic shear moduli of sands[END_REF]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. [START_REF] Biarez | Elementary mechanics of soil behaviour[END_REF][START_REF] Hicher | Experimental behaviour of granular materials[END_REF] 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; [START_REF] Jiang | Inherent and stress-state-induced anisotropy in very small strain stiffness of a sandy gravel[END_REF][START_REF] Hoque | Anisotropy in elastic deformation of granular materials[END_REF] 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.

Ά ⌬ xx ⌬ yy ⌬ zz ⌬␥ yz ⌬␥ zx ⌬␥ xy · = ΄ 1/E h -hh /E h -vh /E v 0 0 0 -hh /E h 1/E h -vh /E v 0 0 0 -hv /E h -hv /E h 1/E v 0 0 0 0 0 0 1/G hv 0 0 0 0 0 0 1/G vh 0 0 0 0 0 0 2͑1 + hh ͒/E h ΅ Ά ⌬ xx ⌬ yy ⌬ zz ⌬ yz ⌬ zx ⌬ xy · ͑40͒
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͒

E v = ⌬ v ⌬ v and vh = - ⌬ h ⌬ v ͑41͒ 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:

E h = ͑1 -hh ͒ ⌬ h ⌬ h and 2 hv 1 -hh = - ⌬ v ⌬ h ͑42͒
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 and2 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 and4 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 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.

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:

Ά ⌬ 11 ⌬ 22 ⌬ 33 ⌬␥ 12 ⌬␥ 23 ⌬␥ 31 · = ΄ 1/E 11 -21 /E 22 -31 /E 33 0 0 0 -12 /E 11 1/E 22 -32 /E 33 0 0 0 -13 /E 11 -23 /E 22 1/E 33 0 0 0 0 0 0 1/G 12 0 0 0 0 0 0 1/G 23 0 0 0 0 0 0 1/G 31 ΅ Ά ⌬ 11 ⌬ 22 ⌬ 33 ⌬ 12 ⌬ 23 ⌬ 31 · ͑43͒
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

͓͔ = ΄ rr 0 0 0 ␥ z /2 0 ␥ z /2 zz ΅ = ΄ r 0 0 0 ␥/2 0 ␥/2 z ΅ ͑44͒
Different loading conditions can be applied, with and without stress rotation, by controlling axial, horizontal, and torsion shear stresses

͓͔ = ΄ rr 0 0 0 z 0 z zz ΅ = ΄ r 0 0 0 0 z ΅ ͑45͒
Therefore 16 terms of the flexibility matrix can be determined as follows:

Ά ⌬ r ⌬ ⌬ z ⌬␥/ ͱ 2 · = ΄ M rr M r M rz M r␥ M r M M z M ␥ M zr M z M zz M z␥ M ␥r M ␥ M ␥z M ␥␥ ΅Ά ⌬ r ⌬ ⌬ z ͱ 2⌬ · ͑46͒
In their paper, Geoffroy et al. ͑2003͒ present experimental values only for the two last columns of the elastic matrix. Measurements of wave propagation velocities by bender and extender elements were used for this purpose in addition to LVDT measurements. The tests were performed on Hostun Sand prepared by air-pluviation. Two types of stress paths were applied on the specimens: a triaxial compression loading ͑type C͒ and a torsion after axial loading corresponding to a ratio 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 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͒.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Local coordinate at interparticle contact

Fig. 3 .

 3 Fig. 3. Distribution of contact orientations with different values of a 0

Fig. 4 .

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

Fig. 5 .Fig. 6 .

 56 Fig. 5. Shear modulus and Poisson's ratio obtained from static hypothesis compared to those obtained with the kinematic hypothesis

Fig. 8

 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

Fig. 7 .Fig. 8 .

 78 Fig. 7. Empirical relationship for packing density as a function of void ratio and of grain size distribution

Fig. 9 .

 9 Fig. 9. Isotropic loading on anisotropic Toyoura Sand and Hime Gravel

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

Fig. 11 .

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

 12 Fig. 12. Stress-induced anisotropy in Toyoura Sand and Hime Gravel

Fig. 13 .

 13 Fig. 13. Stress-induced anisotropy during loading tests with and without stress rotation ͑evolution of M zr , M r , M z␥ , M ␥z , M rz , M z , M zz , and M rr ͒

Table 1 .

 1 Elastic Parameters for Toyoura Sand ͑e = 0.65͒: Static

	Hypothesis
		a 0	G g	␣	a 0␣	n
	2.5	0.24	2,100,000	0.45	0	0.5

Table 2 .

 2 Elastic Parameters for Toyoura Sand ͑e = 0.65͒: Kinematic

	Hypothesis
		a 0	G g	␣	a 0␣	n
	2.5	0.17	2,100,000	0.3	0.27	0.5

Table 3 .

 3 Model Parameters for Hime Gravel: Static Hypothesis

	Material		a 0	G g	␣	a 0␣	n
	Hime Gravel ͑e = 0.5͒	4	1.03	2,100,000	0.56	0.05	0.5

Table 4 .

 4 Model Parameters for Hime Gravel: Kinematic Hypothesis

	Material		a 0	G g	␣	a 0␣	n
	Hime Gravel ͑e = 0.5͒	4	0.72	2,100,000	0.42	0.85	0.5

Fig. 10. Stress induced anisotropy

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

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


A DATABASE FOR CASE HISTORIES AND NUMERICAL MODELLING

This paper deals with a bibliographic database dedicated to the comparison between numerical results and in situ measurements for geotechnical structures. This database, called MOMIS, has been developed by LCPC and ECN for several years. To date, it comprises a total of 416 case histories. The generation of MOMIS has relied on a technologic watch in the field of numerical modelling. The database can be used to highlight modelling principles (in order to provide a guide for good modelling practice to users) and deviations between results given by numerical models and values measured on actual geotechnical structures. It also reflects the evolution of the accuracy of numerical computations.

INTRODUCTION

Numerical methods are now well established as an analysis tool in many engineering problems, and practical applications have been growing at prodigious rate in recent years. With the increase in modelling applications and computing codes, it becomes essential for users that the code be well validated for different point of view and the range of possible applications assessed in a reliable manner.

In the field of geotechnical engineering, validation of numerical codes is of uppermost importance for three reasons:

-results are used for the design (displacements of soil, forces and moments in structures, etc.); -complex constitutive laws are involved for modelling the behaviour of soils and rocks ; -complex computations are carried out by users. The validation of a numerical model is achieved in terms of agreement between numerical results and experimental data. Accordingly, a literature survey is indispensable. Over the past thirty years, only two qualitative compilations, by [START_REF] Duncan | The role of advanced constitutive relations in practical applications[END_REF] and [START_REF] Gens | General report: Prediction, performance and design[END_REF], and a study devoted to undergrounded facilities [START_REF] Negro | Prediction and performance: A review of numerical analyses for tunnels[END_REF] have been produced. No actual quantitative "assessment" has ever been drawn of the deviations observed between numerical modeling efforts and measurement campaigns. In light of this lack of references necessary for evaluating the capacity of both soil constitutive laws and software to reproduce complex situations, the Laboratoire Central des Ponts et Chaussées (LCPC) and the Ecole Centrale de Nantes (ECN) have instituted a technology watch program to inventory these comparisons and estimate the resultant "model error" (Mestat, 2001a).

This technology watch has covered more than four hundred case histories where comparisons between numerical predictions and field measurements have been reported (horizontal and vertical displacements, excess pore pressures, stresses, etc.). Only in situ measurements were considered. The physical models (laboratory tests 1g or centrifuged models) were not taken into account in this survey. As a means of collecting and processing these data, the MOMIS database has been developed.

MOMIS DATABASE

The bibliographic database, called MOMIS (acronym for "Modélisation des Ouvrages et Mesures In Situ"), comprises case histories originating from articles, conference papers, reports and doctoral theses identified as part of the technology watch program. The data set extends back to 1972, with entries being evenly distributed over time (figure 1).

The database currently contains 416 case histories corresponding to the modelling of embankments (84), tunnels (135), sheet-piled retaining structures (66), diaphragm walls (102) and embankments on improved soils (29). Figure 2 reveals that almost half of the cases studied originated in Europe, followed by Asia, North America, Africa and South America. 1970 1975 1980 1985 1990 1995

1

Fig. 2. Experiment locations

Following the classification by [START_REF] Lambe | Prediction in soil engineering[END_REF], the vast majority of modeling efforts recorded (382 case studies) pertain to class C predictions (ie a posteriori predictions) and only 34 cases concern class A predictions (before construction without knowledge of measures). Class A predictions are the ideal type of prediction to evaluate the performance of numerical modelling. Unfortunately they tend to be rather rare. Their limited number is due to the cost of experiments and more predominantly to their time requirements, which seem to be incompatible with the constraints imposed by construction economics.

Even though class C predictions are less demonstrative, it is valuable and useful to draw lessons from these comparisons from both an engineering and research standpoint as well as to derive recommendations for subsequent modelling set-ups (choice of model and guidelines for controlling results) and to quantify the model error.

For each type of geotechnical structure, the number of case histories would seem high enough to be able to produce a statistical overview covering a thirty-year period of numerical computation-measurement comparisons. The complete list of references, relative to embankments and tunnels, have been published (Mestat, 2001b[START_REF] Mestat | Applications de la base de données MOMIS à la validation du calcul des ouvrages souterrains[END_REF].

An electronic version using the ACCESS database management application is being implemented in coordination with the Civil Engineering Laboratory of the Ecole Centrale de Nantes [START_REF] Mestat | Database for class A and C predictions -Comparison between FEM results and measurements in excavation problems[END_REF]. This database is mainly devoted to Class A predictions. For now, the information relative to 13 blind prediction competitions is stored in 63 tables linked together in order to facilitate the data analysis.

For each case study eight families serve to organize and collate the following informations:

-type of analysis (drained vs. undrained condition, consolidation, dynamic, cyclical); -composition of the ground; -construction technique employed, actual dimensions of the structure; -constitutive laws for natural soils and construction materials, parameters; -computational model (dimensions, type of finite elements, mesh density, boundary conditions, loadings, time step, construction project phasing, interface laws, etc.); -set of numerical data for the comparison between computation results and measured values; -conclusions drawn from the comparison (maximum deviation, relative errors, etc.); -bibliographical references.

For the same experimental case, several numerical studies may be generated; such is the case for example when holding a blind prediction competition or a comparison between several constitutive laws for soils.

REVIEW OF CASE HISTORIES IN MOMIS DATABASE

The numerical studies were mainly two-dimensional finite element analyses (Table 1). Some more recent 3D analyses are included in the database (9% of the references). The quantitative comparisons between predicted results and field measurements depend on the geotechnical works.

Embankment (end of construction) -maximum settlement on center line 39 -maximum lateral displacement in depth (toe of slope) 24 -maximum excess pore pressure on center line 16

Embankment (long term) -maximum settlement on center line 38 -maximum lateral displacement in depth (toe of slope) 26 Tunnels (end of construction) -maximum surface settlement (transversal model) 120 -point of inflexion of settlement trough 87 -crown settlement 30 -maximum horizontal displacement 32

Sheet-pile walls (end of construction) -maximum horizontal displacement of the wall 69 -maximum settlement behind the wall 37 -maximum bending moment in the wall 24

Diaphragm walls (end of construction) -maximum horizontal displacement of the wall 77 -maximum settlement behind the wall 19 -maximum bending moment in the wall 18

The data contained within MOMIS database can serve to analyse, in statistical terms, the characteristics of meshes and to highlight modelling principles.

Geometrical Finite Element Model

Given that embankments, retaining structures and tunnels quite often exhibit a much greater length than width, their performance is often studied using a transversal section in plane strain. From the MOMIS data, we can justify the minimum recommendations for meshes used in both vertical and horizontal directions (Table 2). These recommendations are valuable for the numerical models which need not incorporate the presence of substratum near the surface, obstacles or other neighbouring interactions. For the mesh density, there is no real information in the references stored in the database because authors tend to offer little comment on the choices inherent in carrying out geometrical discretization. The majority of 3D models are made for tunneling problems: intersection of galleries, influence of tunnel face, construction sequences which is really 3D, etc. At present time, it is difficult to deal with 3D modeling because there are not many case histories. What we can observe is the following point : in 3D mesh, the section perpendicular to the tunnel axis has generally extensions similar to those fixed for a 2D modeling. For example, the ratio total height of mesh / tunnel axis depth is about 2 for 3D meshes reported in MOMIS database. Of course the 3D mesh density is less refined than for 2D mesh and often not sufficient. A future survey will be made with complementary publications even if there are no comparisons with measurements. It is essential to provide good recommendations for 3D numerical models.

Type of Analysis for Geotechnical Structures

Three theoretical approaches are considered : undrained conditions (29%); drained conditions (46%) and consolidation (25%). From a historical perspective, the first two analyses preceded the third. Consolidation computations with nonlinear behaviour were not possible before the existence of high-speed computers. Nevertheless, a comparison between modeling of embankments, tunnels and retaining structures shows a great difference in the types of analysis (percentages are listed in Table 3). The consolidation approach is rarer for tunnels and retaining structures. Constitutive Laws for Soils

The percentages of constitutive laws used for natural soils in the cases histories are listed in Table 4. The complexity of excavation sequences and computation time lead often to using simple constitutive laws for soil in the problems of tunnels and retaining structures: in more than 60% of case histories the constitutive law is elasticity (EL) or perfect plasticity (EP). On the other hand, the simple and monotonic loading in embankment problems allows to consider elastoplasticity with strain hardening or elasto-visco-plasticity for soft soils. Nevertheless the general trend has favored use of strain hardening models.

2D COMPUTATION MEASUREMENT COMPARISONS

For 2D FEM modelling, comparisons between measured and computed values were made at the end of construction and at long term (particularly for embankments on soft soils). For example, the analysis of surface settlement in tunnel modelling can be achieved with figure 3 (settlement) and figure 4 (characteristics of trough).

Most class C predictions have led to relatively satisfactory results for the maximum vertical settlement on the surface at the end of construction (Figure 3). The points corresponding to computation values and measured values are primarily located within the boundaries defined by the results from class A predictions. This observation was expected because class C predictions can be improved by the knowledge of the measurements.

The quality of predictions focused on the settlement trough is also acceptable (distance of the point of inflexion, figure 4), but its width is generally overestimated. The measured values of maximum lateral displacements in depth are rarer. Figure 5 shows the comparisons between computed results and measurements provided by inclinometers near the tunnel axis.

Fig. 5. Comparison in maximum lateral displacements in depth at the end of construction

Same diagrams can be drawn for embankments or retaining structures [START_REF] Mestat | Database for class A and C predictions -Comparison between FEM results and measurements in excavation problems[END_REF][START_REF] Mestat | Database for finite element prediction and performance -Application to tunnels and retaining structures[END_REF]. This analysis has enabled quantifying the model error committed during the modeling. During this year, a new synthesis has been made for all the references in the MOMIS database and new estimations of relative model errors provided on main measured variables.

ESTIMATION OF RELATIVE MODELLING ERRORS

The modelling error is estimated after numerical computations, indeed it is the "sum" of the errors related to the software, its use, the hypotheses of modelling, the approach employed to obtain the mechanical and hydraulical parameters. In deriving this error, class A predictions are obviously the best fitted. Unfortunately they are rarer. Then all the references included in MOMIS database were taken into consideration.

For each geotechnical structure and each measured variable, five intervals are distinguished for the relative error e : [0, 25%], [25, 50%], [50, 75%], [75, 100%] and [100%, + ∞[. The relative error is the absolute value of the difference between the computed value and the measured value, divided by the measured value.

Figures 6 show the distribution of relative error for the maximum surface settlement (figure 6a) and the maximum lateral displacement at the toe (figure 6b). On the same figures, the results at the end of construction and at long term are represented. The quality of predictions is satisfactory for settlements specially at long term. In contrast the modeling of lateral soil movements is often unsatisfactory ; at short term only 58 % of the predictions reveal a relative error below 50 %, and 21% have a relative error of above 100%.

For retaining structures, best results are obtained with horizontal displacements of the wall (figure 7). On the other hand, vertical displacements on soil surface behind the wall (figure 7) and heave of the bottom excavation present important relative errors : 25% of case histories reveal a relative error of above 100%.

For tunnels, the modelling of the surface settlement is generally satisfactory (figure 8) and the distance to the point inflexion of settlement trough is well described. Even the width of trough is often overestimated by FE results, it is possible to simulate in satisfactory way the settlement trough when the mechanical parameters can be estimated with a sufficient precision.

Figures 6 and9 show that FE models lead at present time to a poor simulation of lateral soil movements in the case of embankments and tunnels. In contrast figure 7 shows that horizontal movements of retaining structures are relatively well described and not the vertical soil movements. This difference confirms that the main movement of each geotechnical works is relatively well simulated by the constitutive laws used in numerical models. Unfortunately it is not sufficient because a "good numerical model" must simulate simultaneously all the aspects of soil-structure behaviour.

ELECTRONIC DATABASE FOR CLASS A PREDICTIONS

As regards the validation and demonstration of a model's capacity to predict structural behaviour, class A predictions (benchmarks) are the most incisive, yet remain the least frequent. In order to collect all the information about each benchmark and to provide better data access, an electronic database is developed [START_REF] Mestat | Database for class A and C predictions -Comparison between FEM results and measurements in excavation problems[END_REF]. At the present time its structure exclusively involves the class A predictions and the academic benchmarks (comparison of numerical calculations without measurement). The aim of this tool is to directly link the in-situ measurements with the predicted values (or predicted values with each other), the type and the parameters of the modelling. By this way, some useful indications regarding the current practice in modelling (discretization, boundary conditions, type of analysis, constitutive relationships, rheological parameters, etc.), their bias, the problems that the authors and participants encountered in these benchmarks, are expected.

For now (august 2003), 13 benchmarks, 66 soil layers characterised by benchmark authors, 184 numerical modelling reports, 369 soil layers characterised by benchmark participants, 13 constitutive models, 6400 numerical results (settlements, horizontal displacements, water pore pressure, forces, moments) have been stored in the database (Table 5).

In the database, benchmarks are distinguished by the reference data :

industrials benchmarks with in situ measurements; semi industrial benchmark with laboratory measurements; academic benchmarks without reference (measurement) to a work. The aim of these benchmarks is to compare codes and settings up of modelling. A first analysis of data contained in the database has been performed to know the use of constitutive models, the distribution of the values of some mechanical parameters, and the distribution of modelling error.

The analysis of constitutive models used in blind prediction competitions show the preponderance of Mohr-Coulomb models (54 %) and its distribution on the whole of geotechnical works seems in accordance with the idea of common practices in modelling, especially for projects (figure 10).

It is worth noting that elastoplastic models with no hardening are more currently used in class A predictions than in class C predictions. This fact is due to many reasons :

for the purpose of comparisons, some benchmarks required the use of only one model. Mohr-Coulomb model is generally imposed as the more common model; benchmark participants prefer a simple and robust model (like for a project); some laboratory or in situ tests required for complex constitutive model are not available in the geotechnical report of the benchmark; class C predictions are often performed to test new sophisticated models. The second study indicates the spread of Young's modulus values in the Mohr-Coulomb model independently of the geotechnical work (see figure 11). The Young's modulus used in these class A predictions is an "interpreted" value of soil stiffness. This value takes into account the rate of deformation in the geotechnical work. So, this figure exhibits a large scatter not representative of intrinsic elastic modulus. This distribution should be compared to the equivalent class C distribution in order to know the possible bias between "interpreted" value (class A prediction) and possible "adapted value" (class C prediction).

A similar analysis has been performed on friction angle of sand provided by participant using the Mohr-Coulomb model (see figure 12). In spite of its "more intrinsic" character, the scatter is very large: from 31 to 42 degrees. The distribution is uniform excepted for 35 and 40 degrees. if the tail of distribution in conservative area (negative values of error) is neglected, calculation results in a maximum over-estimation ratio of 4; the greatest occurrence of settlement error is located between 0 and 0.5. So, the calculation is slightly unconservative;

this greatest occurrence in the case of lateral displacement is located in conservative area. However, this displacement can be overestimated or over estimated. This error on lateral displacement is four times greater than the one on settlement

The Haarajoki competition presents a larger global modelling error. This is mainly due to the fact that the calculation concerns two areas, one of which is composed of drains, making the modelling more complicated.

- As for class C prediction analysis, the long term displacements are better estimated, but the error order is not the same. So, if class C predictions are needed for the validation and calibration of models, class A prediction are essential for the evaluation of modelling error in project conditions.

CONCLUSIONS

The "Laboratoire Central des Ponts et Chaussées" and the "Ecole Centrale de Nantes" have been conducting a technology watch mission over the past several years with respect to comparing finite element model results with in situ measurements. The primary objective is: to preserve the record of these models and their comparisons with in situ measurements, to draw lessons in the practice of geotechnical modeling, to provide orders of magnitude for computation results, and to quantify the modelling error. The information extracted from the bibliographical analysis has been combined into a database called MOMIS.

The MOMIS database is just in the nascent stages of its development. Collation of the references has served to identify a few general conclusions, such as:

-the modelling error is less than 50% for the settlement predictions at the end of tunnel or embankment construction and for the horizontal displacement of a retaining structure at the end of construction; -the modelling error remains high for the lateral displacements in tunneling and for the vertical settlement behind a retaining structure, with a significant number of relative errors surpassing 100%.

On the other hand, an electronic database is currently being developed at EC Nantes (with the ACCESS database management software) specially for the blind prediction competitions (class A predictions). The aim of this specific database is to link the in-situ measurements with the predicted values, or with the parameters of the modelling. By this way, some useful indications regarding the current practice in modelling, their bias, the problems that the authors and participants encountered in these benchmarks, are expected.

Finally it is essential to continue conducting these computation-measurement comparisons along with additional full-scale experiments. We need to ascertain the degree of realism in the models and to define modelling methodologies specific to each type of geotechnical structure.

Fig. 3 .Fig. 4 .

 34 Fig. 3. Comparison in the maximum settlement on the tunnel axis at the end of construction (squares represent class A predictions)

Fig. 7 .Fig. 8 .Fig. 9 .

 789 Fig. 6. Modelling errors for embankments -Comparison between results obtained at the end of construction (in white) and at long term (in black)

  Fig. 10. Constitutive relationships and geotechnical works modelling -Electronic database for class A predictions

Fig. 12 .

 12 Fig. 11. Young's modulus distribution, Mohr Coulomb model (electronic database, class A predictions)

  Fig. 13. Distribution of modelling error. Settlement under the embankment

Fig. 14 .

 14 Fig. 14. Distribution of modelling error. Lateral displacement under the embankment

Table 1 .

 1 2000 Number of comparisons with 2D FE analysisFig. 1. Cumulated number of references over thirty years

		tunnels : 135
	cumulated
	references	embankments : 84
		diaphragm walls : 102
	100
		sheet pile walls : 66
	50
	0
	Year
		Africa
		N. America
		S. America
		Europe
		Asia

Table 2 .

 2 Minimum dimensions for 2D plane strain meshes

	Structures	Total length	Total height
	Embankments	3 L	4 h
	Retaining structures	4 D	2 D
	Tunnels	5 d	2,5 H
	L=1/2 width of embankment base; h=embankment height;

D=height of the wall; d=tunnel diameter; H=tunnel axis depth.

Table 3 .

 3 Type of analysis (%)

	Structures	Drained	Undrained Consolidation
	Embankments	13	20	67
	Diaphragm walls	34	50	16
	Sheet-piles	48	42	10
	Tunnels	73	14	13

Table 4 .

 4 Type of constitutive laws used for soils (%)

	Structures	EL	EP	EPH	EVP HP
	Embankments	11 12,5	49	27,5	0
	Diaphragm walls	40	28	25	4	3
	Sheet-piles	27,5 51	21,5	0	0
	Tunnels	23,5 36,5	34	8	0
	EL=linear or non linear elasticity; EP=elastoplasticity without
	hardening; EPH=elastoplasticity with hardening; EVP=elasto-
	visco-plasticity; HP=hypoplasticity.

Table 5 .

 5 Blind prediction competitions stored in the electronic database (references are listed in[START_REF] Riou | A database for the performance of numerical modelling of retaining wall[END_REF]

	Geotechnical	Organisation	Place	Nbr. of
	work			Partic.
	Embankment	FNRA	Finland	5
	Embankment	DBR	Gloucester	5
			Canada
	Embankment	Instit. of Engineers	Muar Flats	4
		Road Administration	Malaysia
	Retaining wall	Univ. of Technology	Graz	15
	Anchored wall		Austria
	Retaining wall	Univ. of Technology	Graz	12
	Strutted wall		Austria
	Retaining wall	Univ. of Karlsruhe	Hochstetten	43
	Strutted wall		Germany
	Spread footing	Univ. of Technology	Graz	6
			Austria
	Spread footing	Texas A&M Univ.	USA	31
		FHWA.
	Spread footing	Texas A&M Univ.	USA	31
		FHWA.
	Spread footing	Texas A&M Univ.	USA	31
		FHWA
	Spread footing	Texas A&M Univ.	USA	31
		FHWA
	Spread footing	Texas A&M Univ.	USA	31
		FHWA
	Tunnel	Univ. of Technology	Graz	12
			Austria