Mohamed Hassan Khalili
email: mohamed-hassan.khalili@enpc.fr

Jean-Noël Roux
email: jean-noel.roux@ifsttar.fr

Jean-Michel Pereira
email: jean-michel.pereira@enpc.fr

‡ Sébastien Brisard
email: sebastien.brisard@ifsttar.fr

Michel Bornert
email: michel.bornert@enpc.fr

A numerical study of one-dimensional compression of granular materials. II. Elastic moduli, stresses and microstructure

The elastic moduli of a transversely isotropic model granular material, made of slightly polydisperse elastic-frictional spherical beads, in equilibrium along a one dimensional (oedometric) compression path, as described in the companion paper [1], are investigated by numerical simulations. The relations of the five independent moduli to stresses, density, coordination number, fabric and force anisotropies are studied for different internal material states along the oedometric loading path. It is observed that elastic moduli, like in isotropic packs, are primarily determined by the coordination number, with anomalously small shear moduli in poorly coordinated systems, whatever their density. Such states also exhibit faster increasing moduli in compression, and larger off-diagonal moduli and Poisson ratios. Anisotropy affects the longitudinal moduli, C11 in the axial direction, and C22 in the transverse directions, and the shear modulus in the transverse plane, C44, more than the shear modulus in a plane containing the axial direction, C55. The results are compared to available experiments on anisotropic bead packs, revealing, despite likely differences in internal states, a very similar range of stiffness level (linked to coordination), and semi-quantitative agreement as regards the influence of anisotropy. Effective Medium Theory (the Voigt approach) provides quite inaccurate predictions of the moduli. It also significantly underestimates ratios C11/C22 (varying between 1 and 2.2) and C55/C44 (varying from 1 to 1.6), which characterize elastic anisotropy, except in relatively weakly anisotropic states. The bulk modulus for isotropic compression and the compliance corresponding to stress increments proportional to the previous stress values are the only elastic coefficients to be correctly estimated by available predictive relations. We discuss the influences of fabric and force anisotropies onto elastic anisotropy, showing in particular that the former dominates in sample series that are directly assembled in anisotropic configurations and keep a roughly constant lateral to axial stress ratio under compression.

I. INTRODUCTION

This is the second paper in a set of two, published jointly, dealing with the macroscopic mechanical properties of a model granular material in anisotropic structural and stress states, as investigated by discrete numerical simulation. Specifically, the material is oedometrically compressed, i.e., compressed in one direction, with no strain in the orthogonal plane. Such a loading path, in which principal stress directions are fixed and strains remain uniaxial, is one of the simplest ways to obtain homogeneous anisotropic granular materials, and is often a suitable local description for layerwise assembling processes under gravity.

The first paper [START_REF] Khalili | A numerical study of one-dimensional compression of granular materials: I. Stress-strain behavior, microstructure and irreversibility[END_REF], hereafter referred to as "Paper I", studies stress-strain relations and describes how state variables evolve in oedometric compression. Paper I clearly shows that the strain-stress relation in oedometric compression or compression cycles is not elastic, but that elastic moduli express stress-strain response in very small probes superimposed on previously well-equilibrated intermediate states, provided a very small creep phase during configuration stabilization has suppressed friction mobilization. The present paper, which might be read independently, thus does not investigate the quasi-elastic response domain any further. Rather, its objective is to study how elastic moduli are related to microstructural features of the same anisotropic granular packings, and could be measured to infer useful information on such variables as coordination number and fabric. The elastic properties of anisotropic granular materials have quite often been studied experimentally, with sands [START_REF] Chen | [END_REF][3][4][5][6][7][START_REF] Geoffroy | Deformation characteristics of geomaterials[END_REF][START_REF] Duttine | [END_REF][10], less often with glass beads [7,11]. They were also recently addressed in simulations [12][13][14]. We exploit here the variety of initial material states subject to oedometric compression histories, as introduced and studied in Paper I, to compare the tensor of elastic moduli to this literature, and to test modeling schemes.

The paper is organized as follows: first we recall the basic material properties, the method by which elastic properties are computed for each equilibrated configuration (Sec. II), and the material state variation along the oedometric loading paths, explored in paper I, for the different initial packing arrangements (Sec. III). Then we present the evolution of all 5 independent elastic moduli in the transversely isotropic states along the compression curves in Sec. IV, confront our observations to experiments and previous numerical works (Sec.V), and corre-late them to state variables in Sec. VI (which involves testing the performance of simple predictive schemes), before a final discussion (Sec. VII).

II. MODEL MATERIAL AND NUMERICAL METHODS

A. Particles and contact laws

The model material, as defined in Paper I, is composed of elastic-frictional spherical beads, with diameter distribution, uniform by volume (whence D = 2D1D2 D1+D2 ), ranging from D 1 to D 2 = 1.2 × D 1 . Contact elasticity is modeled with a suitably simplified version of the Hertz-Mindlin laws [15], such that, with notation Ẽ = E/(1-ν 2 ), combining the Young modulus E and the Poisson ratio ν of the solid material the beads are made of, the normal contact force transmitted in the contact between beads i and j of respective diameters D i and D j relates to the contact deflection h ij as

F N ij = Ẽ d ij 3 h 3/2 ij (1)
with

d ij = 2D i D j D i + D j
. We use the elastic properties of glass, E = 70 GPa and ν = 0.3 in our simulations (but results, in non-dimensional form, apply to arbitrary materials). Eq. 1 defines a force-or deflection-dependent normal stiffness as

K N ij = Ẽ d ij 2 h 1/2 ij = 3 1/3 2 Ẽ2/3 d 1/3 ij F N ij 1/3 . (2)
The tangential stiffness constant K T ij , relating tangential elastic force F T ij to relative tangential displacement δu T ij , is proportional [15] to K N :

dF T ij = K T ij (h ij )d δu T ij , with K T ij = 2 -2ν 2 -ν K N ij . (3)
We refer to [15] for some necessary rescaling of K T in contacts being unloaded [16], and for implementations of contact laws abiding by the objectivity requirement [17].

The Coulomb condition, applied with friction coefficient µ = 0.3, sets the maximum magnitude of the tangential force to µF N .

B. Stiffness matrices and elastic moduli

Elastic moduli express the relations between small stress increments ∆σ and small strains ε, assuming the contact network, in equilibrium, behaves like a network of elastic springs, with stiffnesses K N and K T varying from contact to contact according to relations 2 and 3.

As we deal with systems enclosed in periodic cuboidal cells, displacements considered in elastic problems (implicitly assumed small) are conveniently parametrized, for all grains i with position r i in the simulation cell, as

u i = ũi -ε • r i , (4)
where ũi 's are fluctuating, periodic displacements, while the second, affine term represents the effect of global strain ε, a symmetric tensor (the minus sign results from our convention of positive compressive strains). Supplementing the displacements of the N particle centers, (u i ) (1≤i≤N ) with their (small) rotations θ i , this results in N f = 6N + 6 degrees of freedom, which we gather in a single displacement vector U = ( Ũ, ε), the coordinates of Ũ comprising all those of (ũ i , θ i ) 1≤i≤N . The conjugate load, denoted as F ext , contains all components of external forces F ext and torques Γ ext applied to the grains as well as the 6 independent components of external stress σ, such that

σ αβ = 1 Ω i<j F (α) ij r (β) ij . (5)
In (5), r

ij denotes the coordinate β of the vector joining the center of grain i to the center of the nearest image, by the group of translations associated to the periodic boundary conditions, of its contacting neighbor j, and Ω is the cell volume. On probing elastic properties, small stress increments ∆σ and contact force increments ∆F ij are considered, related by (5).

The small displacements associated with a load increment ∆F ext should satisfy (to first order in U)

∆F ext = K • U. (6)
(6) is the system of N f linear equations one has to solve to find the N f unknowns contained in vector U, expressing the linear elastic response about a prestressed equilibrium configuration. It involves the N f × N f symmetric, positive definite stiffness matrix K. More details on the structure of matrix K, as discussed in Ref. [18], are provided in Appendix A.

In the present case we determine the 5 independent moduli corresponding to the transversely isotropic granular systems under oedometric compression along axis 1. Those are defined by the following macroscopic relation between stress increments and small strains about an equilibrium prestressed state:

       ∆σ 11 ∆σ 22 ∆σ 33 ∆σ 23 ∆σ 31 ∆σ 12        =        C 11 C 12 C 12 0 0 0 C 12 C 22 C 23 0 0 0 C 12 C 23 C 22 0 0 0 0 0 0 2C 44 0 0 0 0 0 0 2C 55 0 0 0 0 0 0 2C 55        •        ε 11 ε 22 ε 33 ε 23 ε 31 ε 12        (7
) Eq. 7 introduces the so-called Voigt notation, in which ∆σ and ε appear as 6-dimensional vectors, and the elastic moduli are gathered in a second-rank tensor, denoted as C, writing C 11 for C 1111 , etc... Symmetries about all three planes of coordinates and isotropy within the transverse plane [START_REF] Chen | [END_REF]3) require some moduli to coincide, as already written in (7) (e.g., C 66 = C 55 ), others to vanish (e.g., C 14 = 0), and the following relation to hold:

C 22 -C 23 = 2C 44 . (8)
(8) expresses the identity of shear response along all directions within the transverse plane.

To obtain C from the numerical data, one may solve linear system (6) with appropriate values for the righthand side load vector, e.g., setting all its coordinates to zero except one component of ∆σ. This yields one line of the inverse of tensor C. One may also directly determine the elements of C on imposing appropriate values of ε, rather than imposing ∆σ. To do so, one exploits the block structure of matrix K, writing

K = K T L L k , (9)
with a 6N ×6N matrix K associated to particle displacements and rotations and a 6 × 6 matrix k associated with global degrees of freedom, while the nondiagonal block, L, of dimension 6 × 6N , couples stresses to grain displacements and rotations. Stress increments are related to imposed strains as

∆σ = L • Ũ + k • ε, with Ũ satisfying K • Ũ = -T L • ε,
so that the measured elastic moduli are given by

C = k -L • K-1 • T L. ( 10
)
Both methods yield the same results within numerical accuracy.

It should be recalled that the very existence of an elastic response involving symmetric matrix K requires several approximation steps [18]. The present paper deals with elastic moduli as identified through the stiffness matrix approach. In paper I, it was explicitly checked that this expresses the material response to small stress increments about well equilibrated configurations with very good accuracy. A dynamical simulation with all ingredients of intergranular interactions produces, for small stress increments and small strains, the same results as the elastic stiffness matrix approach.

III. EQUILIBRATED CONFIGURATIONS ALONG OEDOMETRIC LOADING PATH

We recall here the characteristics of the systems simulated in paper I, which were subjected to oedometric compression along direction 1: axial strain ε 11 , simply denoted as ε 1 , varies with no lateral strain, ε 22 = ε 33 = 0, and no shear. Equilibrated configurations are stored for different levels of axial stress σ 11 , simply denoted as σ 1 (while lateral stress σ 22 = σ 33 is denoted as σ 2 ). Equilibrated packs are labelled according to the initial state in which they were assembled, and properties are averaged over 3 statistically equivalent samples of 4000 grains each. Values of σ 1 specified in kPa correspond to glass beads with Ẽ = 76.9 GPa. Using, as a control parameter, stiffness level

κ = Ẽ σ 1 2/3
, such that typical contact deflections are of order 1/κ relative to the grain diameter [15], all results are valid for grains with arbitrary elastic properties.

A. Initial states

Initial states, prepared on gently applying a low stress σ 1 = 10 kPa (κ 39000) to a loose "granular gas" configuration until it compresses into a stable equilibrated prestressed solid pack, differ by their solid fraction Φ, coordination number z, and anisotropy. They are labelled with three letters: D or L for dense or loose, followed by H or L for high or low coordination numbers, and then by i or o for isotropic or oedometric compression in the assembling stage (e.g. systems of type DLo are prepared in a dense initial state with a low coordination number by an oedometric packing procedure). Solid fraction Φ in dense initial states (D) is close to the maximum value achieved in random close packing, from 0.634 for DLo to 0.639 for DHo (to prepare such dense packs one suppresses the effects of friction at some stage [15]). Φ is reduced to values of 0.584 (LLo) to 0.589 (LLi) in loose states. Coordination z is rather low in all loose systems, about 4.1 or 4.2 (meaning that we cannot prepare LHi or LHo states). Such loose configurations contain a significant proportion x 0 of rattlers, i.e., grains which do not carry any contact force, typically 10%. Thus the coordination number, if evaluated on the force carrying contact network, involving only N (1-x 0 ) grains, reaches somewhat larger values z * = z/(1-x 0 ) 4.6. Dense systems exhibit as many rattlers and as few contacts as loose ones, if the ultimate stages of their assembling history involves friction: values of z, z * , x 0 are similar in DLo and DLi states as in LLo, LLi ones. If, on the other hand, the effects of friction are suppressed all the way to static equilibrium, dense systems have very few rattlers (1-2 %) and very large coordination numbers (z * 6). DHo, DLo, LLo states differ from their isotropic counterparts DHi, DLi, LLi by the ratio of lateral to axial stresses, denoted as K 0 . Oedometrically assembled samples have K 0 = 0.94 for DHo, 0.51 for DLo and 0.72 for LLo, while isotropic preparations naturally enforce K 0 = 1.

B. Evolution under oedometric compression cycle

The response of all six initial states to a compression cycle in which σ 1 gradually increases to 31.2 MPa (κ = 181) and then slowly decreases back to its initial value 10 kPa was extensively studied in paper I. For the relevant variables Φ, z, and K 0 , it can be summed up as follows. Φ increases by a notable amount, especially for large σ 1 in the MPa range, reaching values larger than the initial ones by 0.01 to 0.02 (with a larger difference in loose states) under the highest axial stress. Remarkably, upon reducing σ 1 , Φ almost reverts to its initial value, meaning that the compression is to a large extent reversible. A small residual density increase is nevertheless observed, more notably in loose systems. Coordination number z increases to values above 5 (5.5 for DLo-DLi under high σ 1 ) in initially poorly coordinated systems (while x 0 decreases to 1%); it reverts to a low value as σ 1 decreases to its initial value. The behavior of DHi and DHo systems is more surprising: z is not monotonic in compression (it first decreases slightly, then increases again to reach about 6.2), and decreases to a low value along the decompression branch of the cycle (a behavior previously reported in isotropic compression [19]).

As to stress ratio K 0 , it stays roughly constant, under growing σ 1 , in states LLo and DLo, which are the most anisotropic. It steadily decreases -the system gaining anisotropy-for initially isotropic systems DLi, DHi, LLi as well as in DHo (which is nearly isotropic) down to the range 0.6-0.8 under the highest σ 1 values. Upon gradually reducing σ 1 to its lowest value 10 kPa, K 0 grows in all systems, moderately (up to 0.7-0.9) for LLo and DLo, and to a much larger extent for all three "XYi" systems and for DHo, in which case the major principal stress becomes σ 22 = σ 33 , as K 0 exceeds 1 (reaching 1.6 for DHi).

C. Fabric and force anisotropy

As in many previous numerical studies [12], the analysis of packing anisotropies in paper I distinguishes fabric from force anisotropies, which pertain to two different angular distributions. Fabric anisotropy refers to the angular distribution p(n) of the unit vectors n normal to the contacts, which solely depends on |n 1 |, or on the angle between the normal direction and the compression axis. The most relevant fabric anisotropy parameter is defined as

c2 = n 2 1 - 1 3 (11)
Force anisotropy exists because the average force carried by the contacts correlates to their orientation, it characterizes the angular dependence of F(n), defined as the average normal force amplitude for contacts with normal direction n, normalized by the global average F N . To leading order, the anisotropy of F(n) is characterized by parameter f2 :

f2 = 1 4π Σ F(|n 1 |)n 2 1 d 2 n - 1 3 . (12)
Like c2 , f2 vanishes in isotropic systems. As a very good approximation in all systems and all stages of oedometric compression cycles, stress ratio K 0 is observed in paper I to relate to fabric and force anisotropy parameters as

K 0 = σ 22 σ 11 2 -3(c 2 + f2 ) 2 + 3(c 2 + f2 ) . (13)
In both states with significant initial anisotropies, LLo and DLo, c2 and f2 first contribute about equally (with c2 f2 0.04 in LLo and c2 f2 0.08 in DLo), and then slowly change under compression. In initially isotropic (DHi, DLi, LLi) systems, both parameters increase, quite fast in the case of f2 . On reducing σ 1 from its maximum back to its initial value, c2 and f2 decrease to low values in LLo and DLo, and change sign in the four other states.

D. Implications for elastic moduli

Elastic moduli, expressing the macroscopic stiffness of a spring network, should tend to vary proportionally to the density of contacts and to the average contact stiffness. The density of contacts is conveniently written as

N c Ω = 3zΦ π D 3 . ( 14
)
The average stiffness, as in the isotropic, monodisperse case [18], can be related to the average normal contact force F N . Assuming no correlation between forces and grain diameters, F N relates to the average pressure P = trσ/3 as

F N = πP D 3 zΦ D . (15)
From ( 15) and (2), the average normal stiffness is given as (still overlooking correlations between forces and radii)

K N = (3π) 1/3 Z(1/3) 2(zΦ) 1/3 D 3 1/3 d1/3 D 1/3 Ẽ2/3 P 1/3 , ( 16
)
where we introduced the notation d for d ij as outlined in (1) and definition

Z(1/3) = F 1/3 N F N 1/3 . ( 17
)
The small level of polydispersity implies that all characteristic diameters are close to the average D = 2D1D2 D1+D2 : thus D 3 1/3 1.003× D , while d1/3 1.0006 D 1/3 .

The ratio Z(1/3), which is a characteristic of the force distribution, varies little between samples and along the loading paths, remaining between 0.92 and 0.945. We directly checked that prediction ( 16) is correct to within 1%.

A naive prediction is that moduli should be proportional to (zΦ) 2/3 [11], obtained on multiplying contact density, or zΦ, by average contact stiffness as written in (16). A power law stress dependence is also expected, as

P 1/3 = σ 1/3 1 (1 + 2K 0 ) 1/3
. Since z varies more (from nearly 4 to about 6.2 between different states and according to the stress level) than Φ (here, 0.58 ≤ Φ ≤ 0.65), one should observe, like in the isotropic systems of Refs. [18,20], that macroscopic moduli classify the different granular packings by their coordination number rather than their density.

Anisotropy should influence the tensor of elastic moduli both through contact orientations -the geometry of the elastic network is not isotropic, with more bonds being oriented nearly parallel to axial direction 1 than near the transverse plane -and through contact stiffnesses -the bonds oriented near the axial direction bear larger forces and, by Eq. 2, are consequently stiffer. To a large extent, these effects are respectively encoded in coefficients c2 and f2 . Unlike the stress ratio K 0 , which is directly related to fabric and angular force distribution, whence relation 13, elastic moduli are not expressed by such explicit formulae, and we thus have to correlate them to anisotropy parameters. We also test, in Sec. VI A, the performance of approximate expressions predicting the moduli, which could be used to directly relate them to anisotropy parameters, if sufficiently accurate.

Our study covers a rather wide variety of equilibrium states, in which the coordination number varies from 4 to 6, stress ratio K 0 from 0.5 to 1.6 and c2 reaches 0.08. This should ease comparisons with laboratory observations, as experimentally probed material states are likely to fall within the range of numerically investigated ones. This should also provide more stringent tests to theoretical attempts at predicting some of the elastic properties of spherical bead packs [13,14].

IV. ELASTIC MODULI

A. Representativity and symmetry.

For each of the six system types and each value of axial stress σ 1 , results are averaged over the 3 statistically similar samples of 4000 spherical grains obtained with the preparation and oedometric compression procedures described in paper I. Sample-to-sample differences do not exceed 5% for poorly coordinated systems, in which the fluctuations tend to be the largest, and decrease to 1% for large coordination numbers.

The full 6 × 6 matrix C as defined by ( 7) is obtained, line by line, using relation 10 and imposing independently the value of each one of the six coordinates of ε, while setting the others to zero. The 6 × 6 matrix is exactly symmetric (within numerical precision).

As to the material symmetries, i.e., transverse isotropy, imposing the matrix to depend on 5 independent coefficients as written in Eqs. 7-8 (with, e.g., C 14 = 0, or C 33 = C 22 ), it is satisfied to the extent that the sample set is statistically representative. In our case, we observe that relative differences between coefficients that should be equal (such as C 55 and C 66 ) do not exceed 2-3% in the worst cases (those, again, of low coordination systems).

To improve the accuracy and representativity of the results, averages are taken over all equivalent moduli or moduli combinations: thus C 66 , C 13 and (C 22 -C 23 )/2 are respectively regarded as independent measurements of C 55 , C 12 , and C 44 .

B. Stress dependence

Longitudinal moduli

The dependence of the longitudinal moduli C 11 and C 22 on the axial stress σ 1 in oedometric loading (leaving aside, for now, the moduli observed in the decreasing axial stress branch of the cycle) are shown in Fig. 1 for all 6 studied systems, on doubly logarithmic plots. As observed previously in isotropic systems [18], moduli are systematically larger in better coordinated systems, and increase with confining stress σ 1 , roughly as a power law, with an exponent slightly larger than the value 1/3 that results from a naive averaging of the Hertz law. A fit of the C 11 data by a power law,

C 11 ∝ σ α1 1 (18)
yield exponents α 1 , in range 0.4-0.46, depending on the system and, actually, on the stress range over which the power law is identified (there is no reason to expect an exact power law over the wide σ 1 range investigated here).

The observation of an exponent α > 1/3 is, in part, due to the increase of coordination number under compression, as new contacts are created between approaching grains. Thus, Fig. 1 shows that the moduli increase faster in systems DLo and DLi, in which the large density and the low initial coordination number entail the most significant increase in z, from nearly 4 to about 5.5 in the compression to the largest σ 1 value. In Fig. 1a the slope, on the logarithmic plot of C 11 versus σ 1 in DHi and DHo systems, appears to decrease in some interval (comprising the first 4 or 5 data points), which is to be related to the slight decrease of the coordination number in that range (see Paper I). Plots of C 22 (Fig. 1b) are farther from straight lines on the logarithmic scale, but a plot versus σ 2 = K 0 σ 1 (Fig. 2) shows that C 22 is slightly better represented as a power law of lateral stress σ 2 :

C 22 ∝ σ α2 2 (19)
This scaling of longitudinal moduli (or of longitudinal elastic wave velocities) along principal stress directions with the stress component in the same direction is reported in experimental studies on sands and spherical bead packs [6,11]. The effect of the decrease of coordination number upon unloading from the highest σ 1 value in sample series DHi is also shown in Fig. 2: initially high ( 6) coordination numbers gradually decrease to low values (hardly exceeding 4, or 4.5 if rattlers are excluded) after one compression cycle. Consequently, elastic moduli, along the unloading branch of the cycle, are close to those of poorly coordinated states.

Shear moduli

Shear moduli, as shown in Fig. 3, exhibit similar stress dependences, approximately varying as power laws of σ 1 , (a) with exponents somewhat larger than 1/3, approaching 1/2 in low coordination systems. According to the literature on sands [START_REF] Duttine | [END_REF], it could be more appropriate to express C 44 , the shear modulus in the transverse plane, as a power law function of σ 2 , the isotropic stress in this plane, while C 55 , associated to shear within planes containing the axial direction, should scale as a power of √ σ 1 σ 2 . While our data, for which stress components do not vary independently, do not enable accurate tests of those predictions, a graph of shear modulus C 44 , using a doubly logarithmic scale, does appear straighter if plotted versus σ 2 , rather than versus σ 1 .

3. Off-diagonal moduli C 12 and C 23 , coupling stresses in one direction to normal strains in orthogonal directions, vary with σ 1 somewhat similarly to diagonal elements of the elastic tensor, as shown in Fig. 4. However, it should be noted that those moduli, as opposed to diagonal terms of tensor C, tend to be smaller in systems with large coordination numbers, and that their growth with axial stress is slower than power law σ 1/3 1 .

C. Anisotropy of elastic moduli

The data of Fig. 1 also reveal the effects of anisotropy on the moduli: at low σ 1 , close to the initial states, C 22 is smaller for DLo than for DLi, while C 11 is larger. DLo and DLi have very nearly the same density and coordination number. They differ by their fabric and stress anisotropy, initially absent for DLi. Anisotropy makes DLo-type samples stiffer in direction 1, hence a larger C 11 than in isotropic systems DLi, while the depleted population of contacts oriented in the transverse directions, as well as their smaller stiffness entail smaller C 22 values. Similar differences are visible between LLo and LLi. Fig. 5 shows that the ratio of longitudinal moduli in the axial and in transversal directions reaches values larger than 2 in DLo states and about 1.4 in LLo ones, while it increases from 1 in initially isotropic systems as elastic anisotropy reflects the growing stress and microstructural anisotropies.

Anisotropy is also reflected by the shear moduli, which are larger when the shear plane contains the major principal stress direction 1. Thus ratio C 55 /C 44 is larger in more anisotropic systems, reaching nearly 1.6 in the DLo sample series, and about 1.25 for LLo systems, as shown in Fig. 6, with values in initially isotropic packings following a similar trend as for the longitudinal moduli (Fig. 5).

Remarkably, some moduli are not sensitive to stress and fabric anisotropy, such as shear modulus C 55 , shown in Fig. 3, as well as modulus C 12 , plotted in Fig. 4 (but C 44 and C 23 do exhibit anisotropy effects). The bulk modulus, expressing the response of average stress P to an isotropic strain, ε = δ1, as ∆P = 3Bδ, is given by

B = C 11 + 2C 22 + 4C 12 + 2C 23 9 . (20)
A plot of B versus average pressure P = (σ 1 + 2σ 2 )/3 does not distinguish material states according to their anisotropy. As shown in Fig. 7, it simply distinguishes large (DHi, DHo) from small (DLi, DLo, LLi, LLo) coordination numbers.

V. COMPARISONS WITH EXPERIMENTAL AND NUMERICAL LITERATURE

The six different sample series subjected to oedometric compression were chosen with the objective of exploring a rather wide variety of initial states, with large differences in density, coordination number and anisotropy. This should ease comparisons with different available experimental and numerical results, and enable tests of existing theoretical models or empirical descriptions of stress and internal state dependences of anisotropic elastic ten-sors. Theoretical attempts to relate moduli to microscopic state variables are tested in Sec.VI A. The present section first confronts our results to available numerical and experimental data.

A. Experiments on sands

Most experiments carried out on granular materials in geomechanics laboratories use devices in which stresses, rather than strains, are imposed. Thus, in the classical triaxial compression test, σ 1 and ε 1 are slowly increased, while lateral stress σ 2 = σ 3 is maintained constant (as opposed to strains ε 2 = ε 3 = 0 in oedometric compression), and lateral strain ε 2 = ε 3 is measured.

Consequently, quasistatic measurements of small strains provide direct access to the tensor of elastic compliances, the inverse of the tensor of elastic moduli, M = C -1 . Denoting as C and M the upper left blocks, of dimension 3 × 3, in C and M, the corresponding compliances are usually written in terms of Young moduli E 1 , E 2 and Poisson ratios ν 12 , ν 23 as

M =         1 E 1 -ν 12 E 1 -ν 12 E 1 -ν 12 E 1 1 E 2 -ν 23 E 2 -ν 12 E 1 -ν 23 E 2 1 E 2         (21)
Empirical formulae have been proposed to relate those compliances to material state and stresses [START_REF] Duttine | [END_REF]. They often involve a certain function f (e) of void ratio e = -1+1/Φ [4], in which the contact network properties are summed up. The correspondence between Young moduli and stress along the same direction was clearly established [6], and expressed as power laws with a single exponent, as

E 1 = f (e)σ m 1 independently of σ 2 ; E 2 = f (e)σ m 2 independently of σ 1 . ( 22
)
Experiments on sands usually record a significant elastic anisotropy under anisotropic stresses (E 1 /E 2 reaching 2) while the effect of fabric anisotropy, as created by the assembling process in isotropically compressed packs, is often smaller (causing, typically, ratios E 1 /E 2 between 0.9 and 1.1) [10]. Such pure fabric effects are ignored in simple relations (22). Such formulae, with a unique factor f (e), assume that, for a given material, the contact network density is determined by the packing density. This is contradicted by the numerical observation of systems of equal density but different coordination numbers, in the present study as well as in previous numerical investigations of isotropic bead assemblies [15,20,21].

Ref. [START_REF] Duttine | [END_REF] proposes a 3-parameter fit of elastic tensors applicable to transversely isotropic granular systems, predicting stress-independent values of Poisson ratio ν numerical study, especially in poorly coordinated systems (as in the isotropic case [18]), even in those systems (DLo, LLo) for which K 0 is approximately constant in oedometric compression -thereby contradicting the model of Ref. [START_REF] Duttine | [END_REF]. (It should be recalled, though, that most experimental studies do not explore such a wide confining stress range as the present set of numerical results.) Poisson ratios are related to moduli as

ν 12 = C 12 C 22 + C 23 and ν 23 = C 11 C 23 -C 2 12 C 11 C 22 -C 2 12 , (23)
and thus tend to decrease under compression in poorly coordinated states, just like the ratios of the off-diagonal elements of the tensor of elastic moduli, C 12 and C 23 to the diagonal ones C 11 and C 22 (see Sec. IV B 3).

B. Experiments on bead packs

Young moduli E 1 and E 2 , as obtained in our numerical simulations, expressing response to uniaxial stress variations in direction 1 or in the transverse plane, exhibit axial stress dependence along the oedometric loading path quite similar to those of longitudinal moduli, as shown in Fig. 9. They can be directly confronted to the measurements published by Kuwano and Jardine [7] on glass bead samples (Fig. 9). Those authors measured elastic moduli of bead samples initially assembled in rather loose (Φ 0.59), anisotropic states, under isotropic pressures ranging from P = 100 to P = 400 kPa. While the loading history is different, similar trends are to be expected, from relation (22). Furthermore, our results (anticipat- ing on Sec. VI) indicate that fabric anisotropy, which is present in those experimental results, is quite an important (possibly dominant) source of elastic anisotropy. Fig. 9 shows a good agreement between those experimental results and the numerical ones obtained in loose systems (or, more precisely, in poorly coordinated systems, which all share similar values of moduli), as regards the absolute values of Young moduli in the available stress range. Both data sets also show a power law increase of moduli with stress, as in (22), with some exponent m exceeding 1/3. The value of m fitted to the experimental points in [7] (0.61 for E 1 and 0.64 for E 2 ) is however larger than the one observed in simulations (about 0.42). Ratios E 1 /E 2 vary between 1.2 and 1.15 in those experiments, similar to their values in sample series LLi in the same stress range, while we observe values near 1.4 in sample series LLo (which would be more appropriate as a model for the experiment). As to shear moduli, as shown in Fig. 3, a good agreement is also to be noted between their values in poorly coordinated numerical systems and in the experiments of Kuwano and Jardine, although the experimental moduli also tend to increase a little faster with stress (with exponent 0.55, as opposed to slightly below 0.5 in simulations) and differ in terms of anisotropy (C 44 > C 55 is reported in [7], instead of the opposite inequality in the present numerical results).

(a)
Another way to obtain elastic moduli is through ultrasonic (or seismic) waves. In a transversely isotropic material, still denoting with index 1 the direction of the axis of rotational symmetry, longitudinal waves (or P waves) and transverse ones (or S waves) propagate in axial or in tranverse directions at velocities given by:

V (1) P = C 11 ρ * ; V (2) P = C 22 ρ * V (1) S = C 55 ρ * ; V (2) S = C 44 ρ * (24)
In (24), index P or S indicates the nature of the wave and index 1 or 2 its propagation direction, while ρ * is the mass density of the granular material. Khidas and Jia [11] measured all four velocities written in (24) in oedometrically compressed glass bead assemblies. They used two different modes of preparation, resulting in different densities, Φ 0.605 and Φ 0.642, and sound velocities were recorded for 100 kPa≤ σ 1 ≤ 900 kPa. In both states the stress dependence of sound velocities were fitted, for axial stresses between 300 and 900 kPa, as a power law with exponent 1/6, corresponding to moduli increasing as σ 1/3 1 . Sound velocities corresponding to the numerical results of the present study can be obtained on using relations (24), with ρ * deduced from the mass density of glass, ρ 2.5 × 10 3 kg.m -3 , as ρ * = ρΦ(1x 0 ) (the rattlers do not belong to the elastic network in which waves are propagated). Comparisons between sound velocities obtained in numerical results and reported by Khidas and Jia are shown in Figs. 10 and11. Obviously, those experimental results resemble the numerical ones obtained with highly coordinated states DHi or DHo, in two respects: the large value of wave velocities and the apparent proportionality to σ 1/3 1 . In the other sample series, Fig. 10 and Fig. 11 (which only considers the looser laboratory samples, with the better quality power law fit) clearly show that a power law fit would yield a notably larger exponent. The level of stiffness probed by the wave speed measurement in the looser state investigated by Khidas and Jia seems somewhat surprising. The discrepancy between numerical and ex-

P , normalized by σ

P , normalized by σ 1/6 2 , versus σ2. Results by Khidas and Jia [11] are shown as thick lines marked "KhJi" for the two different experimental packings.

perimental results could be due, in part, to the observed decrease of coordination number in DLi and DLo systems between the first equilibrium stress 10 kPa at assembling stage and the experimental stress range ∼ 100 kPa. On directly compressing agitated grains (the initial granular gas) under larger σ 1 , it is likely that slightly better coordinated contact networks could be observed above 100 kPa. Anyway, the comparison implies that granular samples with internal states similar to our dense, well coordinated numerical samples of types DH are observed in the laboratory. One may also note that the level of elastic anisotropy recorded in this experimental study, as expressed by ratios C 11 /C 22 (reaching up to 1.4 in denser systems, about 1.25 in looser ones) and C 55 /C 44 (with a maximum value near 1.15 for both sample types) are compatible with the numerical results on systems DHo for σ 1 in the same stress range, as shown in Figs. 56.

Khidas and Jia mainly attribute the anisotropy of elastic properties to the anisotropy of stresses. The vertical to horizontal stress ratios measured in their experiments 2 , versus σ2. Results by Khidas and Jia [11] are shown as thick lines marked "KhJi(L)" for their looser type of packing. reach 2.5 in the looser investigated state, and exceeds 3.2 in the denser one -which seems very large in an assembly of elastic-frictional beads. There are other questions arising in relating the results of Ref. [11] to micromechanics, in particular in relation to the experimental control of the stress state. Some influence of the lateral walls on stress homogeneity might be suspected, given the aspect ratio of the samples, and the apparent influence of a prior stress cycle (up to 400 kPa) on the results. Furthermore, the transversely isotropic symmetry appears not exactly satisfied, since a significant difference is observed between shear moduli C 55 and C 66 . In view of those difficulties of interpretation of the experiments, we did not strive to improve the quantitative agreement with such laboratory results, which should be deemed already satisfactory for dense systems with large coordination numbers.

C. Numerical results on anisotropic bead packs

Our results should be compared to those obtained by La Ragione and Magnanimo, who carried out simulations on anisotropic glass bead assemblies and proposed some schemes to predict elastic anisotropy. One study [13] considers different systems prepared with anisotropic fabric, with a large density and coordination numbers varying from 4.9 to 5.6 (intermediate between DLo and DHo in this respect), under isotropic compression, keeping deviator stress to zero. Fabric anisotropy parameter c2 stays constant as isotropic compression proceeds (from 50 to 500 kPa). and varies from sample to sample, between about 0.008 and 0.022, while z varies from 4.9 to 5.6. In the second study [14] a dense system with isotropic fabric (corresponding to our DHi state) is subjected to triaxial compression up to a principal stress ratio of 1.4, and its elastic moduli recorded a different stages of growing stress deviator. Those two papers thus address separately either the influence of stress anisotropy, or that of fabric anisotropy, on elastic moduli. Both advocate an approach to predict the anisotropy of the tensor of elastic moduli, based on the Voigt estimation scheme. Our results from oedometric compression inevitably mix up both sources of anisotropy. However, as our study covers a significantly wider range of state parameters (see Sec. III), the theoretical approach of Refs. [13,14] can be tested in more demanding conditions than the ones explored by their authors. The issue is discussed in Sec. VI A 1 below.

VI. ELASTIC MODULI AND INTERNAL STATE VARIABLES

We now strive to establish relations between the anisotropic tensor of elastic moduli and the internal state variables characterizing granular samples under oedometric compression, first by testing available theoretical approaches (Sec. VI A), then by resorting to simple "experimental" tests and correlations (Sec. VI B).

A. Predictive schemes

Voigt approximation

The simplest estimate for the elastic moduli is the Voigt approximation, also referred to as "effective medium theory" (EMT), in which an affine field of particle displacements is assumed, as determined by the macroscopic strain [11,18,[22][23][24]. This amounts to ignoring the non-affine contribution to displacements in Eq. 4, and discarding the second term in the right-handside of relation 10, writing, with the notation introduced in (9),

C = k. ( 25
)
To write the Voigt estimates (denoted with superscript V) for the 5 moduli of expression (7), we introduce notations f N = F N / F N for the normal contact force divided by its average,

α T = K T /K N = (2 -2ν)/(2 -ν)
(see Eq. 3) for the ratio of tangential to normal contact stiffnesses, and D for a certain averaged diameter such that (see Eq. 1)

D7/3 = 1 4 (D i + D j ) 2 d 1/3 ij . ( 26
)
With the chosen diameter distribution one has D 1.1275 D . We also introduce a certain factor C 0 , proportional to the contact density and to a typical normal contact stiffness:

C 0 = 3 4/3 2π 2/3 (zΦ) 2/3 D7/3 D 3 2/3 D 1/3 Ẽ2/3 P 1/3 , (27)
and define useful averages over all contacts, involving any coordinates α, β of unit normal vector n:

A α = f 1/3 N n 2 α ; B αβ = f 1/3 N n 2 α n 2 β . ( 28
)
Voigt estimates of elastic moduli are then written as: (no summation over repeated indices)

C αα = C 0 [(1 -α T )B αα + α T A α ] (1 ≤ α ≤ 3) (29) C αβ = C 0 (1 -α T )B αβ (1 ≤ α < β ≤ 3) ( 30
)
C 44 = C 0 (1 -α T )B 23 + 1 2 α T A 3 (31)
Beyond the simplest form of Voigt approximation, estimates of moduli should be improved [18,25,26] on imposing a suitably chosen common spin to all particles. This spin vanishes whenever the strain tensor commutes with the fabric tensor F (defined by F αβ = n α n β ). In the present case this spin effect is only present for shear modulus C 55 , and results in a modified formula

C 55 = C 0 (1 -α T )B 12 + α T A 1 A 2 A 1 + A 2 (32)
These expressions of estimated moduli rely on decoupling averages over diameters on the one hand, and over forces and fabric (which remain coupled), on the other hand. For formula (16), we could check that this does not entail any significant loss of accuracy. However, one should avoid decoupling averages written in (28), as forces and fabric are correlated: contacts oriented near the axial direction are more numerous, and also tend to carry larger forces, whence the inequality

f 1/3 N n 2 α > f 1/3 N n 2 α , (33)
for which we could check the members to differ typically by 10 to 20%. Assuming equality in (33) and decoupling averages in (28) accordingly would reduce the values of estimated moduli, and thus accidentally improve the predictions of longitudinal and shear moduli, which are too large -an improvement based on a fortuitous compensation of errors.

A more accurate treatment of averages defined by ( 28) is possible in terms of expansions of angular distributions of contacts and of f 1/3 N in Legendre polynomials, leading to some formula for Voigt-estimated moduli, analogous to expression (13) of stress ratio K 0 . We did not however, deem such a treatment justified, given the poor accuracy of the Voigt scheme to estimate moduli.

As noted in a number of previous studies [13,18,23] the Voigt approximation is quite inaccurate, and largely overestimates, in particular, shear moduli. It fails particularly badly in poorly coordinated systems, in which the fluctuating (or "non-affine") part of the displacement field contributes the most. This correlation between the inaccuracy of the Voigt estimates and the coordination number is made very conspicuous in plots of ratios of estimates to true moduli shown in Fig. 12. Only in well (a) 4.5 5.0 5.5 6.0 coordinated systems, due to a specific preparation (for DHi and DHo) or to the effect of a large stress in a dense sample (as in DLo, DLi), is the prediction of longitudinal moduli, shown in Fig. 12(a), reasonably accurate. As to shear moduli, Fig. 12(b) makes it clear that they are always severely overestimated. Another basic inadequacy of the Voigt scheme is its poor treatment of moduli C 12 and C 23 coupling normal stresses and strains along different axes. In view of expression (30), proportional to 1α T , estimates C V 12 and C V 23 are considerably too small. Poisson ratios, from Eq. ( 23), are consequently predicted to remain below 0.05, which considerably underestimates their true values shown in Fig. 8.

Despite the poor estimates of individual elements of matrix C, the Voigt estimate B V of bulk modulus B, as defined in Eq. 20, is, as in isotropic bead packs [18], in excess, but reasonably accurate (with ratios B V /B between 1 and 1.15), in all present anisotropic sample series, as shown in Fig. 13. Given its expression (20), this results

B V = C 0 3 Z (1/3) = 1 2 D7/3 Z(1/3) D 3 2/3 D 1/3 zΦ Ẽ 3π 2/3 P 1/3 , (34)
which is independent of fabric anisotropy, and coincides (save for the slight polydispersity effect) with the estimate written in the isotropic, monodisperse case [18]. While the Voigt scheme or EMT does not provide accurate predictions for elastic moduli values, it has been suggested by La Ragione and Magnanimo [13,14] that it could describe anisotropy effects. More precisely those authors concluded that EMT could satisfactorily predict the ratios of moduli between anisotropic systems and reference, isotropic ones. This would imply correct Voigt estimates of ratios C 11 /C 22 and C 55 /C 44 . Those ratios, though, as shown in Figs. 5 and6, are only accurately predicted for longitudinal moduli in moderately anisotropic systems: samples of types DHi, DLi, LLi, and DHo, under moderate axial stress σ 1 , before fabric and stress anisotropy increase due to oedometric compression. Likewise the Voigt approach only provides accurate predictions for ratio C 55 /C 44 as long as it does not exceed 1.1. The Voigt predictions for differences C 11 /C 22 -1 and C 55 /C 44 -1, which are characteristic of elastic anisotropy, remain fair (with an underestimation of ∼ 20%) in those initially isotropic systems, but typically deviate from measured values by 50% in strongly anisotropic sample series LLo and DLo. Our results therefore contradict, in part, the conclusions of Refs. [13,14], which is likely due to our investigation of larger domains of state parameters, as noted in Sec. V C.

The Voigt approach naturally incorporates the influence of stress-dependent average contact stiffness and contact density. Its failure (except for B), particularly in poorly coordinated systems, was discussed in Ref. [18], and related to the properties of granular contact networks with small force indeterminacy. The variations of shear moduli with stress, faster than the expected power law with exponent 1/3, or, correlatively, the variations, faster than the expected z 2/3 , of the ratio plotted in Fig. 12, were attributed to the tendency of moduli to vanish as the limit of vanishing force indeterminacy, occurring at z * = 4 (with a small correction due to 2-coordinated grains), is approached. This phenomenon was first observed [27] in nearly rigid (κ → ∞) frictionless systems, for which it was explained [28] in connection with the anomalous distribution of eigenmode frequencies of the stiffness matrix defined in Eq. 6 [27]. While the absence of force indeterminacy is spontaneously achieved for frictionless grain packs in the limit of large stiffness level κ [15,27,29,30], it is usually not closely approached in the presence of friction, except on setting the friction coefficient to infinity (or a very large value) [15,31]. However, anomalously low shear moduli are observed for the smallest z * values [18,32,33], which increase faster than expected from Voigt estimates as compression entails a small increase of coordination numbers. In the present study of anisotropic granular packs, similar trends are visible in Fig. 3, with a faster variation of shear moduli in poorly coordinated systems, and in Fig. 12, where the strong increase of C V 44 /C 44 for decreasing z * might signal an incipient divergence near z * = 4. The only non-anomalous modulus in poorly coordinated systems with vanishing force indeterminacy is the one expressing the response to a stress increment proportional to the preexisting stress, whence a dominant eigenvalue in the elastic moduli tensor (which is equal to the bulk modulus in isotropic conditions). This is explicitly shown and exploited in Ref. [12] for transversely isotropic assemblies of nearly rigid, frictionless beads. In such extreme situations of vanishing force indeterminacy, ratios of moduli C 11 , C 22 , C 23 , C 12 are all directly related to stress ratio K 0 , which is not the case for our data in the least coordinated systems (as observed in Paper I). The relatively larger values of off-diagonal elements C 12 , C 23 in those small z * systems, and the larger Poisson ratios, are consistent with the approach of the vanishing force indeterminacy limit.

Reuss approximation for a specific load increment

One particular compliance S P , expressing the material elastic response to a stress increment ∆σ proportional to σ, might be evaluated via the Reuss estimate, as defined in Refs. [12,18]. The Reuss approach is based on an evaluation (in excess) of the elastic energy written as a function of force increments. S P , which we refer to as the proportional load increment compliance, expressing the response to ∆σ = δ × σ, is identified through the elastic energy:

∆W = Ω S P δ 2 2 σ : σ. ( 35
)
Introducing principal stress ratio K 0 and compliance matrix elements, one has:

S P = 1 -4ν 12 K 0 E 1 + 2(1 -ν 23 )K 2 0 E 2 . ( 36
)
The Reuss estimate is based on trial force increments ∆F ij = δF ij in all contacts i-j, and reads [12,18]

S R P = 2 3π zΦ 2/3 D 3 2/3 D-1/3 D 5/3 Z(5/3) Ẽ2/3 P 1/3 (37)
Factor Z(5/3) is defined from the force distribution and the friction mobilization as in [18]: approached by its estimate S R P for poorly coordinated systems, since with little force indeterminacy the trial force increments used in the Reuss approximation become closer to the real ones. Remarkably [18], (S R P ) -1 and B V (Eq. 34) only differ by factors of order 1, related to polydispersity and force distribution. The Reuss estimate of S P is accurate with a relative error below 16%, as shown in Fig. 14. In isotropic systems (K 0 = 1), upper bound S V P to S P = 1/B provides a lower bound to bulk modulus: B ≥ 1/S R P . For anisotropic systems, estimates B V for B and S R P for S P are in general the only available quantitatively accurate predictions.

Z(5/3) = F 5/3 N (1 + r 2 T α T ) F N 5/3 , (38)
More sophisticated estimation schemes, in which fluctuations about the main strain field are dealt with selfconsistently, have been developed for isotropic packings [34] (and corresponding ideas, or generalizations thereof, tested by numerical means [35]) with some partial success for shear moduli in well-coordinated systems [18]. Such schemes are complex, and have yet to be generalized to anisotropic systems.

B. Correlation of elastic, force and fabric anisotropies

For lack of adequate predictive schemes, we resort to systematic investigations of correlations between the anisotropies of the tensor of elastic moduli, and structural and mechanical anisotropy parameters c2 and f2 (Sec. VI B 1). In order to sort out the effects of stiffness (or force) anisotropy from those of fabric anisotropy, we also compute the elastic moduli, keeping the same contact network but assuming linear contact elasticity instead of the Hertz-Mindlin behavior (Sec. VI B 2).

Correlations to c2 and f2

Fig. 15 shows how ratios C 11 /C 22 and C 55 /C 44 correlate to anisotropy parameters f2 for forces and c2 for fabric. All values of those ratios of moduli are represented, in oedometric compression as well as in decompression, the latter leading to negative values of both anisotropy coefficients. c2 varies between -0.05 and 0.08, while f2 covers a somewhat larger interval, from -0.06 to about 0.12. However, the gradient of the level of elastic anisotropy, as expressed by those ratios, clearly tends to be along the c2 axis of the graphs, while a number of equal elastic anisotropy (coded as color) data point clusters tend to be oriented along the f2 axis. This implies that elastic anisotropy is to a greater extent determined by fabric anisotropy parameter c2 than by force anisotropy parameter f2 . Note that some data points correspond to states of "reversed anisotropy", i.e. situations in which c2 and f2 are negative, and ratios C 11 /C 22 and C 55 /C 44 smaller than 1. Those data correspond to low stress states after an oedometric compression cycle has been applied to initially well coordinated states DLo and DLi (see Paper I).

Elastic moduli obtained with linear contact elasticity

One more direct way to assess the relative effects of fabric and force anisotropy on moduli is to compute the elastic response of the same contact networks with a linear elastic behavior in the contacts, i.e., on attributing to stiffness K N the same value in all contacts, instead of using the Hertz law (2). As to ratio K T /K N , it is maintained to the same fixed value α T as in the simplified Hertz-Mindlin model used so far in the stiffness matrices. Fabric anisotropy affects the elastic properties of both the Hertzian and the linear models. However, force anisotropy, by which Hertzian contacts oriented near the axial direction tend to be stiffer, does not influence the properties of the network of linear elastic contacts. A quantitative comparison of the results of the Hertzian and the linear models is presented in Figs. 16

c αβ = C αβ d K N (1 ≤ α ≤ 3, 1 ≤ β ≤ 3), (39)
thereby eliminating the influence of the average stiffness (which depends on P as in Eq. 16 in the Hertzian case). The naive expectation that moduli increase linearly with contact density appears to be nearly fulfilled (although of course a proportionality of reduced moduli to zΦ, as would be predicted by the Voigt approach, is not satisfied). As to the values of Poisson ratios and their evolution in oedometric compression, Fig. 8 shows them to be rather well reproduced by the linear model, which yields nevertheless a smaller difference between ν 12 and ν 23 , especially in highly coordinated systems DHi and DHo.

VII. CONCLUSIONS

Using model granular systems initially differing by density, coordination number and anisotropy of fabric and stress, we have been studying how elastic properties reflect the evolution of their internal state under oedomet- ric compression. Although such a loading history does not enable a strict separation of stress anisotropy from fabric anisotropy effects, our study includes a large variety of possible microstructures and leads to a number of quite clearcut conclusions. First, it is amply demonstrated (Figs. 1 to 4) that the elastic moduli are mostly sensitive to coordination number z. Poorly coordinated dense granular materials are not stiffer than loose ones. With low z, moduli tend to increase somewhat faster with confining stress, as a power law with apparent exponents approaching 0.4 for longitudinal moduli (Figs. 1 and2) or 0.5 for shear moduli (Fig. 3). Limited evidence suggest that longitudinal moduli corresponding to compression along principal stress directions are better fitted by a power law of the stress in the same direction (Fig. 2), and that the modulus in the transverse plane is better described as growing like some power of the lateral stress (Fig. 3b). Anomalously low shear moduli are controlled by the rattlercorrected coordination number, z * , and tend to vanish as the value signaling the disappearance of force indeterminacy is approached (Fig. 12b). Poisson ratios (Fig. 8)are larger than in well-coordinated states, and decrease under compression as the system state moves farther away from this singular limit. The Voigt estimates of moduli are not accurate, except for the bulk modulus B (correctly predicted within 15%). Longitudinal and, especially, shear moduli are largely overestimated, especially in states with low coordination numbers. Predictions for off-diagonal moduli or Poisson ratios are much smaller than measured values. A Reuss estimate can be written which correctly predicts (with the same accuracy as the Voigt estimate for B) the proportional load increment compliance. All those observations are direct generalizations of the results previously obtained on isotropically assembled and compressed granular packs [18].

The anisotropy of the tensor of elastic moduli stems from the anisotropies of both the contact network (the fabric) and the angular force distribution. In the set of equilibrium configurations along the oedometric compression curves studied here, both effects tend to vary simultaneously. As a direct consequence of the anisotropic assembling process in dense (DLo) or looser (LLo) configurations with small coordination numbers, the relative difference of longitudinal moduli C 11 /C 22 -1 may be observed in the range 0.5-1.2 in low stress states (Figs. 15 and 18); the same difference decreases to about -0.4 on reversing the anisotropy after one oedometric compression cycle, from initially well coordinated states (Fig. 15). Similar observations apply to the relative difference of shear moduli, C 55 /C 44 -1 (for which values reach up to about 0.25). Comparisons with moduli computed with linear contact elasticity and direct correlations to anisotropy parameters indicate that, among the range of investigated material states of the present study, the larger elastic anisotropies are determined by fabric more than by force anisotropies. In moderately anisotropic systems (say, for |c 2 | ≤ 0.02) both origins of elastic anisotropy appear to have similar contributions, and the Voigt estimates can be used, as a rough approximation, to predict values of C 11 /C 22 -1 and C 55 /C 44 -1 (but not the Poisson ratios). Those differences would be strongly underestimated for larger anisotropy parameters.

To reach more complete conclusions on possible anisotropic states and connections between elasticity, fabric and force anisotropies, it would be necessary to explore different states and microstructures, and to vary stresses and fabric independently, using different (e.g., triaxial) loading paths. The observed strong, possibly dominant effect of geometric (fabric) anisotropy, is encouraging as to the possibility of inferring the fabric anisotropy from appropriate ratios of moduli, in the practically relevant case of transversely anisotropic materials. Interestingly, some of the results, if the influence of the Hertz nonlinearity is secondary to the one of the network fabric, could bear some relevance for other types of contact elasticity (e.g., through conical asperities [36]).

Comparisons with experimental data reveal interesting convergences, but also raise a number of questions. The relations, used in sands, between moduli and stresses, appear more suited to the cases in which stress anisotropy is the major cause of elastic anisotropy, and ignore the possibility, evidenced in numerical simulations, of obtaining very different coordination numbers with the same density. We compared our numerical data to those of two experimental studies of elastic properties of glass bead assemblies. Interestingly, the first one [7] measures moduli in a loose state (similar to numerical state LLo), while the second one [11] deals with dense, well coordinated systems (analogous to numerical states DHi or DHo): we thus found laboratory analogs for very different numerical systems, among the wide range which we numerically explored. While numerical and experimental values of moduli are in quite good agreement, some features of the laboratory results, especially regarding anisotropy, remain unexplained, such as the differences between all three different longitudinal moduli [11], or the sign of the difference between both shear moduli [7]

FIG. 1 .

 1 FIG. 1. (Color online) Longitudinal moduli: (a) C11 and (b) C22 versus σ1 in systems DHo, DLo, LLo, DHi, DLi, and LLi. Dashed lines in plot (a) have slopes 1/3 and 0.42.

FIG. 2 .

 2 FIG. 2. (Color online) C22, in DHi, DLi and LLi systems (same data as in Fig. 1b) versus σ2. Data corresponding to unloading of DHi systems (as shown by arrows) are added in this graph. Straight line slopes are 1/3 and 0.4.

  FIG. 4. (Color online) Modulus C12 versus σ1. Dashed line has slope 1/3.

1 FIG. 5 .

 15 FIG. 5. (Color online) C11/C22 versus σ1 or κ -1 in initially anisotropic (a) or isotropic (b) sample series. Dashed lines depict the prediction of the Voigt approximation to be discussed in Sec. VI A.

1 FIG. 6 .

 16 FIG. 6. (Color online) C55/C44 versus σ1 or κ -1 in initially anisotropic (a) or isotropic (b) sample series. Dashed lines: Voigt approximation (Sec. VI A).

FIG. 7 .

 7 FIG. 7. (Color online) Bulk modulus B versus average stress P . Dashed line slopes: 1/3 (top) and 0.4 (bottom).

1 FIG. 8 .

 18 FIG. 8. (Color online) Poisson ratios (a) ν12 and (b) ν22 versus σ1 or κ -1 . Data points connected by dashed lines show results of linear contact model, discussed in Sec. VI B 2.

FIG. 9 .

 9 FIG. 9. (Color online) Young moduli: (a) E1 versus σ1, and (b) E2 versus σ2. Results of Ref. [7] are shown as thick lines marked "KuJa".

FIG. 10 .

 10 FIG. 10. (Color online) (a) Velocity of longitudinal sound waves propagating in axial direction, V

1 / 6 1

 16 , versus σ1; and (b) velocity of longitudinal sound waves propagating in transverse directions, V

FIG. 11 . 1 S , normalized by σ 1 / 6 1

 11116 FIG. 11. (Color online) (a) Velocity of transverse sound waves propagating in axial direction, V 1 S , normalized by σ 1/6 1 , versus σ1; and (b) velocity of transverse sound waves propagating in transverse directions, V 2 S , normalized by σ 1/6

FIG. 12 .

 12 FIG. 12. (Color online) Ratios of Voigt estimates to true moduli C11 (a) and C44 (b) versus rattler-corrected coordination number z * .

1 FIG. 13 .

 113 FIG. 13. (Color online) Ratio of Voigt estimate to measured value of bulk modulus in all sample series, versus σ1, in oedometric compression. form the underestimation of C 12 and C 23 compensating the overestimation of C 11 and C 22 . Note that, on combining Eqs. 20, written for Voigt estimates of moduli, and 28, a simpler expression of B V emerges,

FIG. 15 .

 15 FIG. 15. (Color online) Color or grey-level-coded values of ratios of moduli, characterizing elastic anisotropy, in plot with coordinates c2, f2, with one round dot per simulated state, in oedometric compression and unloading of all six sample series. (a) ratio of longitudinal moduli, C11/C22; (b) ratio of shear moduli, C55/C44.

FIG. 16

 16 FIG. 16. (Color online) Reduced modulus c11 versus zΦ in (a) isotropically assembled and (b) oedometrically assembled sample series. Data points joined by solid lines: Hertz model; data points joined by dashed lines: linear model.

1 FIG. 18 .

 118 FIG. 18. (Color online) Difference C11/C22 -1 versus σ1 or κ in oedometrically compressed, isotropically (a) or oedometrically (b) assembled sample series. Data points joined by continuous lines for Hertzian contacts, by dashed lines for linear elastic ones.

  and 17 for longitudinal moduli C 11 and C 22 . Those figures show the variations, versus zΦ (proportional to the contact

  Results pertaining to highly coordinated systems DLo and DLi stand apart, because coordination number does not grow monotonically with axial stress in those cases. Results of Hertzian and linear elastic contact networks are quite close, showing that a fair approximation of moduli would be obtained on attributing the average stiffness to all contacts: C 22 would be almost exact, and C 11 somewhat underestimated. Comparing ratios C 11 /C 22 for Hertzian and linear elastic contacts, as shown in Fig.18, shows that similar values are reached in both cases for systems LLo and DLo, which exhibit the largest levels of elastic anisotropy. In those strongly anisotropic systems, the difference (C 11 /C 22 ) -1 is a little smaller, typically by 10% or 20%, with linear, com-

		0.9
		0.8	LLi DHi
		0.7	DLi
	c 22	0.6
		0.5
		0.4
		2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 4.0 4.2 0.3
			zΦ
	(a)
		0.9
		0.8	LLo DHo
		0.7	DLo
	22	0.6
	c	0.5
		0.4
		0.3
		2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 4.0 4.2 0.2
			zΦ
	(b)
	FIG. 17. (Color online) Analog of Fig. 16 for reduced modulus
	c22.
	pared to Hertzian, contact elasticity. It is typically 50%
	smaller in the other sample series, for which anisotropy
	is acquired in the course of oedometric compression from
	an isotropic (or nearly isotropic) initial states. Those
	are the systems in which stress anisotropy is gradually
	induced by compression. As discussed at the end of
	Sec. VI A 1, the prediction of the Voigt approximation
	for (C 11 /C 22 ) -1 in those systems, unlike in sample se-ries LLo and DLo (which are strongly anisotropic from
	the initial state), achieves about 20% accuracy.
	Quite similar conclusions can be drawn for the
	anisotropy of shear moduli, as expressed by difference
	(C 55 /C 44 ) -1.

  . The agreement between numerical and experimental results for the level of anisotropy, as expressed by differences (C 11 /C 22 ) -1 or (C 55 /C 44 ) -1, is only semi-quantitative at best. More simulations investigating a wider range of different states are certainly necessary, but it should be pointed out that the role of the assembling procedures in laboratory tests should also be clarified. For instance, as remarked in[10], samples prepared by gravity deposition (often referred to as 'air pluviation') are observed with greater stiffness in the vertical direction (C 11 > C 22 ) in some experimental studies, while the opposite sign of the anisotropy (C 11 < C 22 ) is obtained in others, carried out by different groups. A fuller understanding of anisotropic initial states thus requires more accurate experimental as well as numerical investigations.

ACKNOWLEDGMENTS

This work was partly supported by a French government grant managed by ANR within the frame of the national program Investments for the Future, ANR-11-LABX-022-01.

Appendix A: Structure of stiffness matrix

We provide here additional information on the definition and structure of stiffness matrix K, as introduced in (6) and discussed in Ref. [18]. Contact force F ij transmitted from i to j is split into its normal and tangential components as

Considering an equilibrium configuration with N c force-carrying contacts, the static contact law relates the 3N c -dimensional contact force increment vector ∆f , formed with the values ∆F N ij , ∆F T ij of the normal and tangential parts of all contact force increments, to relative displacements, gathered in 3N c -dimensional vector δu as ∆f = K • δu.

(A1)

This defines the block diagonal local stiffness matrix K, of dimensions (3N c × 3N c ). K does not couple different contacts. The 3 × 3 block of K corresponding to contact i, j, K ij , provided friction is not fully mobilized reads:

Stiffnesses K N (h ij ) and K T (h ij ) are given by ( 2) and ( 3). This simple form of K ij is an approximation, relying on simplifications of the Hertz-Mindlin laws (see [18]). The rigidity matrix [15], which should not be confused with the stiffness matrix, relates the relative displacements to the N f degrees of freedom. Normal unit vector n ij pointing from i to j (along the line joining centers for spheres), the relative displacement δu ij , for spherical grains with radii R i , R j writes

r ij denoting the vector pointing from the center of the first sphere i to the nearest image (by the periodic translation group of the boundary conditions) of the center of the second one j. The normal part δu N ij of δu ij is the increment of normal deflection h ij in the contact.

The rigidity matrix G is 3N c × N f -dimensional, it is defined by the linear correspondence expressed by relation (A3), which transforms U into the 3N c -dimensional vector of relative displacements at contacts δu:

The equilibrium equations -the statements that contact forces f balance load F ext -can simply be written with the transposed rigidity matrix, as [15]

Stiffness matrix K results from Eqs. A4, A1 and A5:

Note that the symmetry and the positiveness of the stiffness matrix is immediately apparent in (A6), given the diagonal form of K in (A2), and the positive signs of stiffnesses K N , K T . In Ref. [18], a correction to K is identified, which stems from the preexisting forces prior to application of stress increments. This correction, which is not symmetric in general, is however negligible as long as K is positive definite, which is equivalent to the nonexistence of a non-trivial kernel of G, i.e. "mechanism motions" or "floppy modes". Such free motions are absent [15] in equilibrated packs of spherical balls, save for global translations of all grains in the periodic cell, and for (rare) two-coordinated balls. Global rigid-body translations are readily eliminated in matrix computations on assuming one grain to be elastically tied to a fixed "ground", and two-coordinated grains are also innocuous (see [15,18]).