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: I. Stress-strain behavior, microstructure and irreversibility

The behavior of a model granular material, made of slightly polydisperse beads with Hertz-Mindlin elastic-frictional contacts, in oedometric compression (i.e., compression along one axis, with no lateral strain) is studied by grain-level numerical simulations. We systematically investigate the influence of the (idealized) packing process on the microstructure and stresses in the initial, weakly confined equilibrium state, and prepare both isotropic and anisotropic configurations differing in solid fraction Φ and coordination number z. Φ (ranging from maximally dense to moderately loose), z (which might vary independently of Φ in dense systems), fabric and force anisotropy parameters and the ratio, K0, of lateral stresses σ2 = σ3 to stress σ1 in the compression direction are monitored in oedometric compression in which σ1 varies by more than 3 orders of magnitude. K0 reflects the anisotropy of the assembling process and may remain nearly constant in further loading if the material is already oedometrically compressed (as a granular gas) in the preparation stage. Otherwise, it tends to decrease steadily over the investigated stress range. It is related to force and fabric anisotropy parameters by a simple formula. Elastic moduli, separately computed with an appropriate matrix method, may express the response to very small stress increments about the transversely isotropic well-equilibrated states along the loading path, although oedometric compression proves an essentially anelastic process, mainly due to friction mobilization, with large irreversible effects apparent upon unloading. While the evolution of axial strain 1 and solid fraction Φ (or of the void ratio, e = -1 + 1/Φ) with axial stress σ1 is very nearly reversible, especially in dense samples, z is observed to decrease (as previously observed in isotropic compression) after a compression cycle if its initial value was high. K0 relates to the evolution of internal variables and may exceed 1 in unloading. The considerably greater irreversibility of oedometric compression reported in sands, compared to our model systems, should signal contact plasticity or damage.

I. INTRODUCTION

The macroscopic mechanical properties of quasistatically deformed, solidlike granular assemblies are traditionally described and modeled in the realm of soil mechanics [START_REF] Wood | Soil Behaviour and Critical State Soil Mechanics[END_REF][START_REF] Biarez | Elementary Mechanics of Soil Behaviour[END_REF][START_REF] Mitchell | Fundamentals of Soil Behavior[END_REF] by phenomenological laws, often resorting to the concepts of elastoplasticity, which need nevertheless to assume complex forms if stress-strain curves are to be described with some accuracy. These laws are applied in engineering practice and have benefited, over the last decades, from sophisticated laboratory measurements [START_REF] Tatsuoka | Geotechnics for roads, rail tracks and earth structures[END_REF][START_REF]Deformation characteristics of geomaterials: Recent investigations and prospects[END_REF].

Beyond phenomenological description, investigations of the mechanics of granular materials in connection with their microscopic structural and rheophysical features are now being pursued throughout a significantly wider, multidisciplinary research community [START_REF]Physics of Dry Media[END_REF][START_REF]Powders and Grains 2009[END_REF][START_REF]IUTAM-ISIMM symposium on mathematical modeling and physical instances of granular flow[END_REF], ranging from geotechnical engineering to condensed matter physics. To this end, particle-level numerical simulations, analogous to molecular dynamics, and often referred to as "discrete element modeling" (DEM) for granular materials [START_REF] Radjaï | Discrete-element modeling of granular materials[END_REF][START_REF] O'sullivan | Particulate discrete element modeling, a geomechanics perspective[END_REF], provide very valuable information. Such approaches have successfully been employed to investigate grain-level origins of such important aspects of granular mechanics as dilatancy properties by which dense and loose configurations differ in their small strain response [START_REF] Thornton | [END_REF][START_REF] Radjaï | Granular Materials: Fundamentals and Applications[END_REF][START_REF] Suiker | [END_REF], the role of grain shape or such features as rolling resistance and angularity [14][15][16][17] in the development of internal friction. An important concept, the "critical state" -an attractor state under monotonically growing strains, as in homogeneous quasistatic shear flow, which does not depend on initial conditions -has been characterized and its properties related to micromechanical aspects [START_REF] Radjaï | Granular Materials: Fundamentals and Applications[END_REF][18][19][20].

Unlike such steady states, the configurations and the mechanical response of granular materials under small and moderate strains are sensitive to the initial material structure. In addition to their density, initial states are characterized in terms of structural anisotropy, whose importance has long been recognized in experiments [21][22][23][24][25][26], and more recently investigated by numerical means [27][28][29]. Coordination numbers (i.e., average numbers of force-carrying contacts per grain) have also been observed to vary, independently of density. As coordination numbers are mostly inaccessible to experiment, these observations were carried out in numerical studies [30][31][32], although some indirect comparisons with laboratory observations, through elastic moduli, were also proposed [33]. Compressive loads, in which stress intensities, rather than stress directions, are varied, are crucially influenced by contact deformability and elasticity. Isotropic compression [34,35] of isotropically assembled model granular materials is apparently a simple process, in which the contact network gets enriched due to a recruitement process closing gaps between neighboring grains. This simplicity is partly deceptive, though: upon unloading, while only a very small strain irreversibility is observed, the contact network may undergo profound changes; specifically, the final coordination number under low pressure may be much smaller than its initial value, before the compression cycle. Simple prediction schemes based on the assumption of homogeneous strain prove unable to capture such changes in contact number [35].

Oedometric compression, an axially symmetric process in which one principal strain component (ε 1 ) is increased, the others being maintained at zero (ε 2 = ε 3 = 0), is one of the simplest anisotropic loading process -representative of natural materials under gravity (e.g. sediments consolidating under their weight). Oedometric compression leads to transversely isotropic structures, with the symmetry of revolution about axis 1.

The present paper, the first in a set of two, reports on a numerical study of oedometric compression of a model material made of elastic-frictional spherical beads, investigating how material anisotropy, either initially present or acquired in the compression, couples to stresses and strains. We shall systematically refer to ε 1 , σ 1 as axial stress and strain, and to ε 2 = ε 3 and σ 2 = σ 3 as the lateral or transverse strains and stresses, respectively. Extending previous studies of isotropically assembled and compressed materials [30,31], a special attention is paid to the various possible initial states, which differ in density, coordination number and anisotropy, and their influence on the subsequent material response under load. As in a number of recent experimental and numerical studies [26,[36][37][38][39], the ratio K 0 of lateral to axial stresses, which results from both the initial state and the effects of the subsequent compression, is monitored, and related to internal anisotropy. The effects of unloading and compression cycles are explored, both macroscopically (strains, density, K 0 ) and microscopically (contact network and fabric). Elastic moduli are investigated in order to assess the nonelastic, irreversible nature of compression response. Elastic moduli are also measurable characteristics apt to probe material microstructure and anisotropy, but this aspect is dealt with in the companion paper [START_REF] Khalili | Submitted to Phys Rev E[END_REF].

The paper is organized in the following way. The model material and the simulation procedure are first described in Sec. II. The different numerical packing methods and the resulting initial states under low stress are presented in Sec. III. Results on the oedometric compression of the different initial states are reported in Sec. IV, both for the macroscopic behavior and the evolution of internal variables and microstructure. Sec. V discusses the mechanical response with reference to elasticity, frictional dissipation and contact network instabilities. Sec. VI then investigates the effects of unloading and compression cycles. Sec. VII finally sums up and discusses the results.

II. NUMERICAL MODEL

Based on the integration of the equations of motion for solid objects, involving linear and angular momentum, masses and moments of inertia, the DEM simulations exploited here are a standard tool in granular micromechanics, as used in many articles [START_REF] Thornton | [END_REF][START_REF] Suiker | [END_REF]27] and described in more comprehensive treatises [START_REF] Radjaï | Discrete-element modeling of granular materials[END_REF]. We therefore dispense below with a full presentation of all relevant equations, by referring adequately to previous published work in which very similar models were implemented [30,[START_REF] Peyneau | [END_REF]. Nevertheless, we need to fix notations, introduce relevant control parameters, and provide a sufficient definition of the numerical procedure.

A. Model material

We consider assemblies of spherical beads, interacting in their contacts through contact elasticity and Coulomb friction. The beads are slightly polydisperse, with diameter D distributed according to the following probability density function:

p(D) = 2D 2 1 D 2 2 (D 2 2 -D 2 1 ) 1 D 3 (D 1 ≤ D ≤ D 2 ), (1)
which ensures a uniform distribution by volume between minimum value D 1 and maximum value D 2 = 1.2 × D 1 . This distribution should avoid all crystallization phenomena (even though no nucleation tendency was detected in a previous study [30] applying a very similar treatment to monodisperse bead assemblies). For contact elasticity, the same simplified version of the Hertz-Mindlin model [START_REF] Johnson | Contact Mechanics[END_REF] is adopted as in Ref. [30], suitably adapted to a polydisperse bead collection. Specifically, considering two beads i and j, with respective centers at points r i and r j and radii R i and R j , and introducing notation h ij = ||r jr i || -R i -R j , the normal force transmitted F N ij in their contact vanishes for h ij > 0 (distant bead surfaces) and depends otherwise on deflection -h ij ≥ 0 as follows. Introducing the Young modulus E and the Poisson ratio ν of the solid material the beads are made of, and using notation Ẽ for E/(1ν 2 ), one has

F N ij = Ẽ d ij 3 |h ij | 3/2 (2) in which d ij = 4R i R j R i + R j is an effective diameter combin-
ing surface curvatures in the contact region. Due to Eq. 2 the normal stiffness expressing the response to small variations of deflection |h ij | in the contact varies as

K N ij = Ẽ d ij 2 |h ij | 1/2 = 3 1/3 2 Ẽ2/3 d 1/3 ij F N ij 1/3 . (3)
The tangential elastic force F T ij , as in [30], relates to the relative tangential displacement in the contact, δu T ij , involving a (deflection-dependent) tangential stiffness coefficient K T assumed proportional to K N :

dF T ij = K T ij d δu T ij , with K T ij = 2 -2ν 2 -ν K N ij . (4)
Tangential stiffness K T has to be suitably adapted (rescaled) whenever the normal elastic force decreases, in order to avoid spurious elastic energy creation [30,[START_REF] Elata | [END_REF].

The Coulomb condition enforces inequality ||F T ij || ≤ µF N ij , with the friction coefficient µ set to 0.3 in the present study. It is taken into account by suitably projecting F T ij onto the circle of radius µF N ij in the tangential plane, after applying incremental relation (4), whenever necessary. Normal and tangential contact force components also follow the general motion of the grain pair in order to ensure the objectivity of the model [30,44].

While we use the elastic properties of glass, E = 70 GPa and ν = 0.3 in our simulations, results, if suitably expressed in dimensionless form, exactly apply to all materials sharing the same dimensionless characteristics µ and ν. (Moreover, Poisson ratio ν only mildly affects elastic properties [33], while material properties do not vary fast with the friction coefficient in the range 0.2 ≤ µ ≤ 0.4 [45,46]).

Finally a normal viscous force is added to the elasticfrictional one, in order to ease the approach to mechanical equilibrium under static loads. The same model as in Refs. [30,47] is used, with a damping constant choice corresponding to a velocity-independent, very low coefficient of restitution in binary collisions.

B. Boundary and loading conditions

We consider cuboidal samples, periodic in all three directions. We denote as L 1 , L 2 , L 3 the dimensions of the simulation cell parallel to the three axes of coordinates, to which correspond basis unit vectors e 1 , e 2 , e 3 . In oedometric compression, L 2 and L 3 are kept fixed, while L 1 varies, either enforcing the value of strain rate ε1 = -L1 /L 1 , or requesting the system to reach an equilibrium configuration under a given level of the corresponding normal stress in direction 1, σ 1 .

In the stress-controlled case, L 1 , the cell size in direction 1, satisfies an equation of motion such that it slowly increases or decreases, according to the sign of the difference between requested and measured values of σ 1 [30,[START_REF] Peyneau | [END_REF]47]. Both strain rate controlled and stresscontrolled simulations are carried out in such conditions that inertial effects remain small, by enforcing small enough strain rates ε. Inertial effects are controlled by requesting the inertial number I [48] to remain very small. I is defined in terms of the mass m 1 of a grain of diameter D 1 and characteristic stress σ 1 as

I = ε m 1 D 1 σ 1 . (5)
We request I not to exceed I max = 10 -3 in the preparation stage, as a granular gas is gradually compressed to form the initial solid configuration. In the subsequent strain-rate controlled, quasistatic oedometric compression, we set ε to a smaller value, corresponding to I = 10 -5 ; upon unloading (as σ 1 and ε 1 decrease), even smaller strain rates are imposed, corresponding to I = 10 -6 . Such low values of I insure the absence of all rate influence on all measured quantities, and avoid, in particular, instabilities associated with spuriously large contact losses along the unloading path (as remarked in Ref. [35]).

As in [30,[START_REF] Peyneau | [END_REF], the effect of global strains is equally felt by the grains throughout the sample, upon decomposing their motion into a fluctuating, periodic part and an affine contribution; and, in the stress-controlled case, the equation for cell dimensions L α involves an acceleration term proportional to the difference between the requested value of stress σ αα and the currently measured one, using the following classical formula for stress components

σ αβ = 1 V   N i=1 m i v α i v β i + 1≤i<j≤N F α ij r β ij   . (6)
Eq. 6 expresses stress components as a kinetic term involving masses m i and velocities v i of all N grains i within sample volume V , added to a sum over pairs of interacting grains i, j transmitting force F ij (from i to j) in their contact, r ij denoting the "branch vector" pointing from the center of i to the center of j. The first (kinetic) term of the right-hand side of Eq. 6 is of course negligible at (or close to) mechanical equilibrium, but might somewhat influence the system dynamics in the initial assembling stage.

Although the present study focusses on oedometric compression, in which only L 1 varies, while L 2 and L 3 are fixed, we are also interested in the consequences of the procedure by which the initial solid configuration is assembled by compression from a loose configuration (a 'granular gas'), and we consider both oedometrically and isotropically compressed initial configurations under stress σ αβ = P δ αβ , as in Refs. [30,[START_REF] Peyneau | [END_REF], all three L α 's being simultaneously reduced in the isotropic case. Isotropic compression is, however, only applied to granular gases at the assembling stage. The resulting granular packs are then subjected to oedometric loading paths.

III. SAMPLE PREPARATION, INITIAL STATES

A. Motivation

Although widely recognized as crucially important for small strain mechanical response of granular materials, the assembling processes by which granular packs are prepared in a solid state are relatively seldom investigated, either experimentally [23][24][25]49] or numerically [28,50]. One serious difficulty in the numerical modeling of such processes is the dependence [27,28] of the final microstructure (density, coordination number) on dynamical dissipation parameters (such as restitution coefficients) which are not well known, and for which modeling choices are often guided by computational convenience as much as by physical realism.

We chose here to implement idealized assembling models, with the objective of obtaining a variety of initial structures representative of a wide range of possible material states. Although admittedly not conforming to laboratory procedures, those numerical preparation methods can be argued to exhibit some of the basic features resulting from such procedures as vibration or gravity deposition. Their main advantage is the possibility of varying, through rather wide intervals, the basic state variables: solid fraction Φ, coordination number z, and anisotropy in contact orientations. We also record the proportion of rattlers (i.e., grains that do not carry any force, in the absence of gravity), denoted as x 0 , and other data pertaining to interneighbor distances and force networks.

B. Numerical assembling process and stress control

In practice, all configurations are obtained on compressing a loose configuration (Φ 0.45, no intergranular contact, no kinetic energy) to equilibrium under low initial stress σ 0 . Specifically, one requests σ 1 = σ 0 for oedometric compression, or σ 1 = σ 2 = σ 3 = σ 0 for isotropic compression. We use σ 0 = 10 kPa, assuming the particles are glass beads. This stress level is expressed in dimensionless form using a stiffness parameter κ, defined, as in [30,[START_REF] Peyneau | [END_REF]48], by

κ = Ẽ σ 1 2/3 . ( 7
)
With glass beads ( Ẽ 77 GPa), σ 0 = 10 kPa corresponds to κ 0 39000. Definition [START_REF]Powders and Grains 2009[END_REF] is such that the typical ratio of contact deflection (-h ij in Eq. 2) to grain diameter is of order κ -1 [30].

Defining F 1 = σ 1 D 2 1 as the relevant force scale, we request, as a practical equilibrium condition, all forces to balance on each grain within tolerance 10 -4 F 1 , all torques to balance within tolerance 10 -4 F 1 D 1 , while the kinetic energy per grain should not exceed 10

C. Different initial packing states

The assembling procedure is designed such that these configurations equilibrated under σ 1 = σ 0 in the oedometric case, or σ 1 = σ 2 = σ 3 = σ 0 in the isotropic case, vary in density and in coordination number. We thus obtain 6 different initial states, as described below. Their properties are listed in Table I.

Maximally dense states with high coordination number. In order to maximize density, as in [30], we set the friction coefficient µ to zero in the assembling stage. In the isotropic case, this results in the random close packing state (RCP), as often investigated in the literature [30,[START_REF] Peyneau | [END_REF]52,53] for monodisperse spherical beads. The results obtained here are almost not affected by the slight polydispersity, as the solid fraction averages to 0.638 with 4000 beads, hardly larger than the value 0.637 reported on applying the same treatment with the same number of grains in the monodisperse case [30]. In the rigid limit, the coordination number of non-rattler grains, z * = z/(1x 0 ) approaches the isostatic value of 6 in the limit of rigid grains (κ → ∞), while x 0 remains quite small. We refer to this initial state as DHi for Dense, High coordination number, isotropically compressed in the assembling stage.

Using oedometric, rather than isotropic, compression to reach an equilibrated configuration with σ 1 = σ 0 , still without intergranular friction, we obtain the DHo state (Dense, High coordination number, oedometrically assembled). Their properties are similar to DHi, save for a slight anisotropy, evidenced in the observed value K 0 0.94 of the ratio of transverse to axial stress. As studied in Ref. [START_REF] Peyneau | [END_REF], frictionless bead assemblies may transmit anisotropic stresses, but the corresponding states (dubbed "anisotropic random close packing states" in [START_REF] Peyneau | [END_REF]) share the same density and coordination numbers as the isotropic ones.

Maximally dense packing with low coordination number. In order to mimic laboratory procedures in which dense configurations are obtained through agitation or vibration, and to obtain presumably more realistic values of coordination numbers, a very small isotropic dilation is applied to DHi configurations (multiplying bead center coordinates by 1.001), so that all contacts open; the system is then subjected (after a mixing step as in [30]) to an isotropic compression, with the final value µ = 0.3 of the friction coefficient (the one used in the subsequent study of the solid, quasistatic response), until an equilibrium state is obtained. The solid fraction is close to the initial RCP value, but the coordination number is now about 4, with a large rattler fraction. This final state is referred to as DLi (Dense, Low coordination, isotropically assembled. The same procedure may also be applied, with an oedometric final compression stage (with µ = 0.3). We note, in that case, that the mixing stage is not necessary to obtain a small final coordination number. We denote the final state as DLo (Dense, Low coordination, oedometrically assembled. Interestingly, this DLo state exhibits a rather large preparation-induced stress anisotropy, with K 0 0.51.

Looser states. States obtained on directly compressing the loose "granular gas" configurations (Φ = 0.45), with the final value µ = 0.3 of the intergranular friction coefficient (and condition I ≤ I max = 10 -3 ), are low coordination states with many rattlers, but looser than DLi or DLo, with solid fractions between 0.58 and 0.59. Depending on whether they are isotropically or oedometrically assembled, we refer to these initial static configurations as LLi or LLo, in accordance with previous notation conventions (Loose, Low coordination, isotropically or oedometrically assembled). K 0 , in the LLo case, shows a significant level of stress anisotropy, but smaller than for DLo. Although rather loose, the LLo state is likely not the loosest possible structure of a rigid bead assembly with µ = 0.3. Different procedures, possibly involving capillary cohesion in an initial stage [54], could result in lower solid fractions.

All 6 initial states of Tab. I were prepared for 3 different samples of 4000 grains, over which recorded measurements are averaged. One additional system with 13500 grains was assembled in the DLo state, in order to check for the absence of size effects. Data pertaining to this larger sample are labelled "DLo+" on some figures below.

D. Other characteristics of initial states.

The possibility to obtain a low coordination number z for a solid fraction Φ nearly equal to its maximum (random close packing) value was first pointed out in previous studies of isotropic packings [30] [START_REF]Save for the small grain size polydispersity, states DHi, DLi, LLi are identical to the states respectively denoted as A, C and[END_REF]. It is generalized here to anisotropic packings DLo. Apart from the data listed in Tab. I, initial states might also be characterized in terms of force distribution, friction mobilization and neighbor distance statistics. Some of these properties are studied in Sec. IV, in which state variable evolutions under oedometric compression are studied.

Although grains should tend to have more neighbors in denser systems, arbitrarily small displacements suffice to open contacts in the rigid limit (κ → ∞), whence the possibility of widely different z values. It is instructive to explore at which scale neighbor distance statistics are correlated to contact statistics. initial states, plots of growing function z(h), defined as the number of neighbor grains separated by a gap lower or equal to distance h, and such that z(0) = z. As in the isotropic, monodisperse sphere packings of Ref. [30], z(h) functions take larger values in denser systems, except for small values of h (say, h ≤ 0.05D 1 ) as the value of the contact coordination number is approached. Thus, direct observations of bead packs by microtomography techniques [START_REF] Aste | [END_REF] should be able to resolve distances of order D 1 /100 to provide information on coordination numbers. In Ref. [30], it is also shown that the treatment of rattlers (which would tend to rely on neighboring grains underneath in the presence of gravity) might significantly affect function z(h) at small distance. In view of the close similarity with the results obtained in isotropic systems, the structure of (orientation averaged) pair correlations is not pursued further in the present paper, which focusses more on anisotropy.

IV. OEDOMETRIC COMPRESSION

We now report on the observed material evolution in oedometric compression. Once the details of the numerical loading procedure are specified in Sec. IV A, the variation of the simplest, scalar state variables is monitored and discussed in Sec. IV B. Observations of stress anisotropy, as expressed by coefficient K 0 , are then reported and compared to literature results (Sec. IV C). The anisotropy of force networks, as investigated in Sec. IV D, is directly related to K 0 , as shown in Sec. IV E.

A. Loading process.

In order to dispel all possible confusion, let us first insist that the main objective of the present study is oedometric compression of the chosen model material with intergranular friction coefficient µ = 0.3. All samples of 4000 grains in all six initial states, even though some of them were prepared, as described above in Sec. III, without friction, and/or by isotropic, rather than oedometric, compression of a granular gas, are subjected to quasi-static oedometric compression, for which µ is set to 0.3. In such slow, quasistatic compression processes, the value of viscous damping parameters in contacts is known to be irrelevant to the material behavior [START_REF] Radjaï | Discrete-element modeling of granular materials[END_REF].

We apply a strain rate-controlled loading program [maintaining small values of I (Eq. 5), see Sec. II B]. Intermediate configurations are recorded when σ 1 reaches, for glass beads, values 31.62 kPa, 100 kPa, 316.2 kPa, 1 MPa, 3.162 MPa, 10 MPa, 31.62 MPa (increasing as a geometric progression, with factor √ 10), corresponding to κ decreasing (by constant factor 10 1/3 ) from 39000 at 10 kPa down to 181 at 31.62 MPa. Each of those intermediate configurations is subjected to equilibration under constant σ 1 -using the same numerical tolerance on equilibrium criteria as stated in Sec. III, before strain rate-controlled compression is resumed. This equilibration step is carried out in order to record accurate characterizations of contact networks. It results in a very small "creep" strain increment (typically of order 10 -5 for dense states, up to 10 -4 in looser systems), in which z increases by a small amount (from 0.5 to about 3.5 %, the highest increases corresponding to the less coordinated states) [57]. In the following, state variables are, unless specified otherwise, measured in equilibrated configurations. We checked that those measurements, and the subsequent behavior recorded on resuming constant strain rate compression, do neither depend on the chosen compression rate, provided I remains below 10 -3 , nor on changing the values of σ 1 corresponding to equilibration stages.

B. Evolution of scalar state variables

Density

Solid fractions for the different initial states are plotted versus σ 1 or κ -1 in Fig. 2. Two different sets of curves are obtained, pertaining to initially dense and initially loose systems. Quite unsurprisingly, the density increase with applied stress is of the same order as κ -1 , characterizing contact deflections, and tends to be larger in looser systems. The influence of the initial coordination number seems quite secondary, only noticeable, on the figure, for DLi and LLi systems under low stress. Note the DLo+ data points, showing the same behavior in a larger sample (N=13500 grains) -save for a very small systematic density increase, compatible with the slight size dependence (proportional to N -1/2 ) recorded for the RCP density in [30,52].

A comparison with similar isotropic compression results [35] reveals a density increase of the same order, but, as one should expect, somewhat smaller under axial stress σ 1 than under isotropic stress

P = σ 3 = σ 2 = σ 1 .
Thus dense systems reach solid fractions of about 0.658 under isotropic stress P = 31.6 MPa [35] while oedometrically compressed dense samples (DHo or DLo) do not exceed solid fraction 0.654. Axial strain ε 1 , on the other hand, is larger in the oedometric compression case, as it accounts for the whole density change ∆Φ, rather than ∆Φ/3 in isotropic compression.

Coordination number and rattler fraction

The variation of coordination number, z, under growing stress, as shown in Fig. 3, is slightly more unexpected: while it gradually increases with σ 1 in low coordinated states, whatever their density, it tends to decrease in a first stage if initially high (in DHi and DHo systems). Note the absence of sample size dependence: results for N=13500 and N=4000 grains are identical. The rattler fraction, on the other hand, steadily decreases in compression for all initial states. The nonmonotonic variation of z, if initially high, is a first clue that the oedometric compression is not always as simple and predictible as might be expected on assuming homogeneous shrinking of distances along the axial direction. Such a uniform strain assumption necessarily predicts coordination number z to increase. It would tend to explain, nevertheless, the faster increase of z in DL systems than in LL ones: the creation of new contacts in compression is achieved sooner in more densely packed structures. Quite detailed tests of the homogeneous shrinking assumption were carried out for isotropic compression in Ref. [35], revealing fairly correct predictions of coordination number increases under compression: within 20-30% of the measured values in low coordination number systems. (Yet the homogeneous strain assumption proved unable to capture the evolution of coordination numbers on unloading -see Sec. VI). As for Φ, changes in z under compression prove smaller in the present oedometric case than in the isotropic compressions of Ref. [35], with maximum values of z below 6.5 (for DHi-DHo), in range 5.8-6 (DLi-DLo states), or below 5.5 (Li-Lo states), as opposed to, respectively, 6.8, 6.6 and 5.9 under isotropic load with the same value of σ 1 .

The results of Fig. 3 also signal the enduring effects of initial anisotropy: the differences between systems DLi and DLo do not tend to vanish, even after the applied stress increased by more than 3 orders of magnitude.

Force distribution

From [START_REF]Physics of Dry Media[END_REF], in equilibrium, the average normal force F N in the contacts is readily related to the pressure

P = (σ 1 + σ 2 + σ 3 )/3, as F N = π D 3 P zΦ D . (8)
This formula, involving the first and third moments of the diameter distribution, assumes a decorrelation between normal force intensity F N ij and intercenter distance R i + R j in contacts i, j, which is satisfied in good approximation (the maximum relative error is 1.4 %, for the highest stress level). The evolution of the normal force distribution in compression is well characterized on normalizing forces by F N . The probability distribution function (p.d.f.) of f = F N / F N was observed in the isotropic case [35] to concentrate on a narrower interval about its average as compression proceeds, the faster the better coordinated the system. In the present case of oedometric compression, this effect is still present, although considerably smaller. Fig. 4 shows the p.d.f. of f for different stresses σ 1 for initial state DLo, and exhibits little change, except for large forces. The shape of force distri- butions may be characterized [35] with reduced moments

Z(α) = (F N ) α F N α (9)
As shown in Fig. 5, the reduced second moment Z(2) decreases quite slowly (except for the initial evolution of high z states, likely correlated with the nonmonotonic variation of z) as a function of σ 1 . Similarly, Z(1/3) (not shown on the figure), which from Eqs. ( 3) and ( 8) relates the average contact stiffness to P 1/3 [35] hardly changes in oedometric compression: for all systems, it stays between 0.92 and 0.95. Unlike the isotropic compression studied in [35] for similar stress levels, the oedometric compression does not cause strong changes in force distributions.

C. Stress ratio K0

Traditionally, coefficient K 0 (termed "coefficient of earth pressure at rest") is regarded as a basic characteristic of material response under oedometric load, and expected to remain constant as axial stress σ 1 increases (in a horizontal, homogeneous sand layer under its weight, both vertical and horizontal stresses thus increase proportionally to depth). However, some preparation techniques reportedly produce initial state dependent K 0 , which might also vary with σ 1 -thus raising the question of the conditions in which K 0 might indeed be regarded as a constant ratio [37,58,59]. The availability of the six widely different initial states, in the present study, provides an opportunity to investigate this issue.

Results

The variations of K 0 along the oedometric loading path are displayed in Fig. 6, in which isotropically assembled initial states are distinguished from oedometrically assembled ones, on separate plots. Systems that are first assembled by isotropic compression only gain stress anisotropy in the course of the subsequent oedometric compression, and thus exhibit K 0 values decreasing from 1, faster for a higher density, and faster for a larger coordination number. Thus K 0 , in the DHi case, reaches values slightly above 0.6 at σ 1 = 1 MPa (or κ -1 6.10 -4 ), and hardly changes under larger axial stress. K 0 , in DLi systems (dense with low initial coordination) and LLi (loose) ones, steadily decrease as functions of σ 1 , without approaching an asymptotic value, even under quite high stress levels (tens of MPa). Among the three different oedometrically assembled initial states, the dense, highly coordinated one, DHo, is in a nearly isotropic stress state (K 0 = 0.94) and close (see Tab. I) to fully isotropic state DHi. Consequently, the behavior of K 0 is quite similar in oedometric compression for DHo and DHi. However, the looser anisotropic initial state (LLo), and the dense, yet poorly coordinated one (DLo), both exhibit quite different K 0 evolutions in oedometric compression, with a remarkably constant value ( 0.5) for DLo, and a very slowly decreasing one for LLo (if a smaller interval of σ 1 is considered, K 0 0.7 might be considered constant, as a good approximation, for LLo as well). Constant K 0 values are thus observed in situations for which the anisotropy of the assembling process is similar to that of the subsequent quasistatic oedometric loading history.

Comparison to experimental and numerical literature

Okochi and Tatsuoka [58] published a detailed experimental study of factors affecting K 0 values measured in a reference, well characterized sand, subjected to many different initial treatments in the preparation stage, including a first compression (which is not oedometric as lateral strains are not set to zero), up to a relatively small stress level (about 20 kPa), in which σ 2 /σ 1 is kept constant by separately controlling axial and lateral stresses (in a triaxial cell). In the subsequent oedometric compression to higher stress values, they observed lower values of K 0 in denser systems, and, a decrease of K 0 for increasing σ 1 , from an initial isotropic state of stress. Interestingly, K 0 converges to an asymptotic value (close to 0.5) as σ 1 increases to about 200 kPa. This asymptotic value does not appear to depend on the initial stress ratio applied in the first compression. The numerical results appear to be in qualitative agreement with these observations, except that, in the numerical case, a limiting, high stress value of K 0 is not approached as soon as σ 1 is merely multiplied by 10. The experimental setup of [58] does not, however, enable the simpler procedure which consists in, first, depositing the sample under gravity, and then applying an exactly one-dimensional compression (zero lateral strain). Carried out in an oedometer, equipped with tactile pressure sensors to record lateral stresses, rather than a triaxial cell, the measurements of Gao and Wang [59] allow for such a simpler procedure to be applied. Using air pluviation (i.e., deposition under gravity, under controlled conditions) to assemble the solid sample, these authors observed a constant K 0 in the subsequent oedometric compression. Thus, regarding the oedometric assembling procedure, defining initial states LLo and DLo (Tab. I), as roughly similar to pluviation, it should be noted that simulations agree with laboratory observations in this respect: K 0 remains constant for assembling procedures resembling one-dimensional compression [60] The experiments of Lee et al. [37] also reveal a roughly constant K 0 when the tested granular materials (assembled by some unspecified process) are first oedometrically compressed from 16 to 115 kPa. Tested materials include glass beads prepared at intermediate solid fractions : Φ 0.603, for which K 0 0.55, and Φ 0.614, for which K 0 0.51 -values with which (despite possible different contact friction coefficients) our numerical results for dense system DHo approximately agree.

The experimental study of Khidas and Jia [26], carried out on glass beads in oedometric conditions, considers two different initial densities (Φ 0.605 and Φ 0.643), and aims at a characterization of anisotropic elastic properties. Values of K 0 are noted, though corresponding to a secondary compression process, after a first compression cycle -which we shall briefly discuss in Sec. VI.

On the numerical side, two recent publications are particularly relevant for comparisons with the present study, Refs. [38,39]. Both Lopera Perez et al. [38] and Gu et al. [39] consider spherical grains, with the polydispersity of the Toyoura sand particles of Ref. [58], and prepare samples by varying the friction coefficient, in the initial assembling stage by isotropic compression, from zero to its actual value used in quasistatic compression: µ = 0.5 in [39] (0.600 ≤ Φ ≤ 0.629), µ = 0.25 in [38](0.600 ≤ Φ ≤ 0.648). Those initial states thus interpolate between DHi and LLi. Investigated stress ranges in these studies are narrower than in our case, extending from 25 to 1250 kPa in [38] (where particles with elastic properties of glass beads are also simulated), and from about 130 to 1040 kPa in [39] (once stresses are rescaled in order to compare systems of equal stiffness level κ). Gu et al. also prepared samples by direct oedometric compression (similar to a series of systems interpolating between DHo and LLo).

Our results for isotropically compressed systems for DHi and LLi agree semi-quantitatively with those of [38,39] for K 0 values and trends, but the decrease of K 0 values for growing σ 1 is notably slower in our case. Similarly, our observations contradict those made by Gu et al. as regards the difference between LLi and LLo. These authors obtain roughly constant (σ 1 -independent) K 0 values for oedometrically assembled systems, and find that K 0 in isotropically assembled systems approaches this value as soon as σ 1 increases by a factor of 10 or 20, while, in the present study, K 0 values for LLi and LLo still differ after a thousandfold increase of σ 1 . One possible explanation for these discrepancies is that we could approach the quasistatic limit in compression with better accuracy: the rate of compression, as measured by the inertial number, is 250 times as small in our simulations as in [38] (while its value is left unspecified in [39]). Neither one of those two groups studied initial states of different coordination numbers for the same density, and thus our DLi and DLo results are, to our knowledge, entirely new. (Coordination numbers are not specified in [38]; values of z * specified in [39] in the looser samples approximately agree with our LLi or LLo results for the same value of κ; Gu et al. do not prepare systems as dense as our "DXy" ones).

D. Anisotropy

We now characterize anisotropy appearing in oedometrically compressed systems, both in the contact network, in the correlations between neighboring grains, and in force intensity and friction mobilization.

Contact and neighbor pair anisotropy

As a result of oedometric compression, the distribution of the orientation of unit normal vectors in contacts, on the unit sphere, ceases to be isotropic, although it remains rotationally symmetric about the compression direction (referred to as the axial direction and denoted "1" throughout the paper). Defining angle θ between direction 1 and that of normal unit vector n, with 0 ≤ θ ≤ π, the orientation distribution (or fabric) anisotropy is conveniently expressed by the probability density function (p. d. f.) of cos θ = n 1 over interval -1 ≤ n 1 ≤ 1, p(n 1 ). By construction, it is an even function (n and -n are equivalent), constant with value 1/2 in an isotropic system. p(n 1 ) might be expanded in the series of Legendre polynomials, with only terms of even order. Truncating the series after the term of order 4, one has:

p(n 1 ) = 1 + A 2 3n 2 1 -1 + A 4 35n 4 1 -30n 2 1 + 3 , (10)
in which coefficients are related to moments of the distribution: thus coefficient, A 2 , given by

A 2 = 15 4 n 2 1 - 1 3 = 15 4 1 -1 p(n 1 )n 2 1 dn 1 - 5 4 , (11)
is directly related to the difference between the second moment and its isotropic value, for which we introduce the notation

c2 = n 2 1 - 1 3 . (12)
Fig. 7 shows that expansion [START_REF] O'sullivan | Particulate discrete element modeling, a geomechanics perspective[END_REF] is already quite a good representation of the p. d. f. of |n 1 | [i.e., P (|n 1 |) = 2p(n 1 )], and that adding the term of order 4 achieves an excellent fit.

As σ 1 increases in oedometric compression, c2 evolves, for the 6 different investigated states, as displayed in Fig. 8. Isotropic packings, as well as nearly isotropic DHo ones, progressively acquire an anisotropic structure under growing oedometric load, faster in dense systems than in the LLi case: from the results of Fig. 3, contact networks undergo more changes for higher densities. Under large stress, the level of fabric anisotropy of initially isotropic systems is comparable to its value in the LLo case, which is roughly stress-independent (c 2 0.04). The larger value of c2 in the most anisotropic system, DLo, decreases slightly for the larger stress levels. The high coordination numbers reached at large σ 1 in dense systems (Fig. 3) precludes very large fabric anisotropies, as many neighbors in contact with the same central grain, by steric exclusion, tend to be more isotropically distributed at its periphery [START_REF] Radjaï | The Physics of Granular Media[END_REF].

It is worth investigating over which length scale the distribution of neighboring grains is similarly anisotropic. To this end, Fig. 9 plots function c2 (h) obtained on extending the definition of c2 (its value for h = 0) to the orientation of normal vectors joining neighbors at distance below h, both in the initial states, and under a high stress level. While the strong anisotropy of the DLo state, and the moderate ones observed under stress in DHo and initially isotropic ones, tend to vanish at distances reaching 0.05 to 0.1D 1 , the distribution of neighbor pairs in the LLo case is still notably anisotropic over a much larger range of interparticle gaps, extending to about 0.2D 1 . The anisotropic structure of loose systems, resulting from the assembling process, should be easier to detect with microtomography techniques.

Angular distribution of normal force amplitudes.

Contact force values tend to reflect stress anisotropy, resulting both from one-dimensional loading and from the initial packing process (for oedometrically assembled states). Classifying contacts by the orientation of normal vector n, some classes tend to carry larger forces than others. We denote as F(n) the average normal force am- plitude for contacts with normal direction n, normalized by the global average F N , such that its integral over the unit sphere, Σ, weighed by the orientation distribution p(n), satisfies

Σ p(n)F(n)d 2 n = 1. ( 13
)
Similarly to p(n), F, a function of |n 1 |, may be expanded in a series of Legendre polynomials. We define

f2 = 1 4π Σ F(|n 1 |)n 2 1 d 2 n - 1 3 , (14)
which vanishes in isotropic systems. Fig. 10 shows the evolution of f2 in oedometric compression. Under compression, the force anisotropy parameter f2 steadily increases in all studied systems. This increase is strikingly fast in initially isotropic states, especially those with a large coordination number. Unlike fabric, which requires changes in the contacts (presumably related to finite strains), force anisotropy might change quickly by redistributing the forces within the existing network. Such redistributions are easier in better coordinated ones, whereas force values are more strongly constrained by the network geometry in poorly coordinated states with relatively low force indeterminacy.

E. Estimation of K0 from anisotropy parameters

Relating stresses to fabric and force anisotropy parameters is quite a standard, well-known procedure in granular micromechanics [START_REF] Peyneau | [END_REF][START_REF] Azéma | [END_REF]63], which was, in particular, successfully applied to oedometric compression in the recent numerical studies discussed in Sec. IV C 2 [38,39]. We use it here, in a particularly simple form, to relate K 0 to anisotropy parameters c2 and f2 , and to discuss the roles of both kinds of anisotropies. Remarkably, the contribution of tangential forces to stresses σ 1 and σ 2 = σ 3 remains very small (below 6%) in all configurations throughout the compression cycle. Ignoring it, we obtain the desired approximative relationon truncating at the second order the expansions of p(n) and F(n), and neglect the products of anisotropic coefficients c2 and f2 (those small coefficients are dealt with to first order [64]). The normal force contribution to principal stresses, given (α = 1, 2, 3) by:

σ α = 3zΦ F N D ∂ D 3 Σ p(n)F(n))(n α ) 2 d 2 n, (15)
becomes,

σ α 3zΦ F N D π D 3 n 2 α + f (α) 2 - 1 3 , (16)
using the notation

f (α) 2 = Σ n 2 α F(n)d 2 n. Knowing that n 2 2 = n 2 3 = 1 2 (1 -n 2 1
), one may finally estimate K 0 as:

K 0 = σ 22 σ 11 2 -3(c 2 + f2 ) 2 + 3(c 2 + f2 ) . ( 17
)
This estimate of K 0 agrees quite well with measured values (see Fig. 6). Thus anisotropy parameters c2 and f2 appear combined into the sum c2 + f2 , and their relative effect is appreciated on comparing their values, as plotted in Figs. 8 and10. Force anisotropy is clearly the dominant effect in stress anisotropy (expressed by ratio K 0 , deviating from 1) in isotropically prepared systems, except for the loose initial state under high stress. In oedometrically assembled ones, with initial fabric anisotropy, both effects, of fabric and force anisotropies, are of the same order, except in the case of DHo, with its very small initial value of c2 .

V. ELASTICITY AND FRICTION

We now investigate the nature of stress-strain response in oedometric compression, discussing the role of elastic and frictional response.

A. Oedometric compression and elastic response

Elastic moduli in the six anisotropic states subjected to oedometric compression are specifically studied in the companion paper [START_REF] Khalili | Submitted to Phys Rev E[END_REF], in which their connection to density, coordination number, fabric and force anisotropies are investigated in detail. The issue we wish to address here first is whether and how the quasistatic stress-strain or stress-density curves (Fig. 2) recorded under oedometric load relate to elastic response.

Measurement of elastic moduli

An elastic response, in granular materials, is measured when small stress and strain increments about a prestressed, equilibrated configuration, are related in a reversible way, associated with an elastic potential energy. Elastic moduli may then be measured either statically, with adequate devices apt to capture very small strains, or deduced from sound wave velocities in granular materials [26,[65][66][67][68][69]. An elastic response is only observed for small strain intervals, and should in fact be viewed as an approximation, as dissipation mechanisms are always present (in particular, solid friction) and preclude the general definition of an elastic energy. The relative amount of dissipation decreases as the size of the probed strain interval approaches zero, and it is often observed, for usual conditions in which granular materials are probed, that an elastic model is satisfactory for strain increments not exceeding some upper bound of order 10 -6 or 10 -5 . For that reason, the material behavior is best characterized as "quasielastic" in that limited range. In Ref. [33], a numerical study of elastic properties of isotropic spherical bead assemblies, carried out with the same model material as the present one (except for the small polydispersity, absent in [33]), observations were made, quite similar to those of the experimental literature, as to the amplitude of the "quasielastic" domain. In simulations, an elastic model is considered for wellequilibrated configurations, in which the contact structure behaves just like a network of linear elastic springs. One may then build the stiffness matrix (also known as the "dynamical matrix") for this network, with stiffness parameters K N and K T as determined by Eqns. 3 and 4, by the procedure explained in Ref. [33], where details are provided about the necessary approximations to obtain an elastic response. The elastic moduli are then obtained by solving appropriate systems of linear equations, for the small (linear and angular) displacements of all the grains associated with global strains and stresses. We refer to [33] and to the companion paper [START_REF] Khalili | Submitted to Phys Rev E[END_REF] for details about the stiffness matrix and its treatment. In the present case, we obtain all five independent elastic moduli appropriate for a transversely isotropic material (as in [26,[START_REF] Peyneau | [END_REF]68]). Specializing to diagonal matrix components, the relation between stress increment ∆σ and strain increment ∆ε reads:

   ∆σ 1 ∆σ 2 ∆σ 3    =    C 11 C 12 C 12 C 12 C 22 C 23 C 12 C 23 C 22    •    ∆ε 1 ∆ε 2 ∆ε 3    , (18)
with a symmetric positive definite matrix of elastic moduli C αβ , 1 ≤ α, β ≤ 3, abiding by the rotational invariance about axis 1. Postponing a more complete study of the (transversely anisotropic) tensor of elastic moduli in oedometric compression to the companion paper [START_REF] Khalili | Submitted to Phys Rev E[END_REF] (in which shear moduli are also measured), we focus here on moduli C 11 and C 12 , which express the response to varying axial strain ε 1 . Those moduli increase with σ 1 in the compression, mainly due to the contact law. In view of Eqns. 3 and 4, moduli tend to scale as σ 1/3 1 . As in the isotropic case [33], they are primarily sensitive to coordination numbers, with values in poorly coordinated dense systems DLo and DLi close to the ones observed in loose systems.

Stress increments and elasticity

Note that the assumption of elastic response underlying relation (18) implies that sliding contacts are absent or have negligible effets, and that the contact network is stable. This may of course be checked by confronting the predictions of ( 18) to a complete DEM computation, in which a steadily, very slowly growing strain is applied, and the effects of friction and of network rearrangements are taken into account. Such a comparison, carried out in the isotropic case [33], showed the elastic response and the complete computation to coincide for small enough strain or stress intervals, in good agreement with laboratory results.

In the present study, a growing strain ε 1 is imposed in the axial direction, and the elastic response of an equilibrium configuration to a small increment ∆ε 1 should be:

∆σ 1 = C 11 ∆ε 1 ∆σ 2 = ∆σ 3 = C 12 ∆ε 1 . (19)
Fig. 11 compares the predictions of ( 19) for ∆σ 1 , with modulus C 11 identified from the stiffness matrix (thereby assuming an elastic behavior in all contacts), to the full DEM-computed mechanical response to small σ 1 increments, in one DLo and one DHo systems, equilibrated for different intermediate values of σ 1 along the oedometric curve. The elastic modulus correctly describes the initial slope and the first data points recorded on the curve (while the strain increment is of order 10 -6 or 10 -5 ), and then the material response turns softer. A similar comparison is made for ∆σ 2 in Fig. 12, showing similar small strain and stress intervals for which the lateral stress increments abides by the elastic prediction, using modulus C 12 . This time, the material gets stiffer as it departs from the elastic response.

The results on the 6 different systems (one sample of each type) for the amplitude of the strain interval for which the quasi-elastic model applies are gathered in Fig. 13. The convention was adopted here that the elastic response is correct as long as the relative error made on predicting stress increments with Eq. 19 remains below 5%. As noted earlier in isotropic systems [33], this quasielastic range, expressed as a strain interval, is of the same order as observed in experiments. It tends to be larger in better coordinated systems: DHo, DHi, and also DLo and DLi once z has significantly increased under compression (Fig. 3). It also increases with σ 1 , roughly proportionally to σ 2/3 1 . This exponent [33] may be regarded as a reflection of a roughly σ 1 -independent quasielastic range, if expressed in terms of relative stress increase ∆σ 1 /σ 1 . As moduli tend to scale as σ 1/3 1 , a constant ∆σ 1 /σ 1 translates into the observed scaling σ 2/3 1 for strain increments.

K0 and elasticity.

Should stress variations with axial stress ε 1 satisfy elastic behavior and relations (19), then stress ratio K 0 should be related to elastic moduli. Specifically, defining an "incremental" stress ratio K 0 = ∆σ 2 /∆σ 1 , as in Ref. [39], K 0 should be equal to C 12 /C 11 . However, as already apparent in Figs. 11 and12, stress increments differ from the predictions of the tangential quasielastic behavior and, consequently, coefficient K 0 and ratio C 12 /C 11 vary independently, as visualized in Fig. 14. These remarks raise the question of the status and validity of simulated experiments in which elastic moduli are measured, as in the case of the results shown in Figs. 11 and12. One may wonder how one can observe, e.g., a constant stress ratio K 0 along the quasistatic compression curve, on the one hand; and a different ratio of stress increments, ∆σ 2 /∆σ 1 = C 12 /C 11 for small probes applied to any intermediate equilibrium state along the curve, on the other hand. The solution to this conundrum is provided by the very small "numerical creep" phenomenon observed when well equilibrated contact networks are obtained along the primary, strain rate controlled compression curve (Sec.IV A). As the system evolves towards a well equilibrated configuration (which is necessary to build a nonsingular stiffness matrix), it is observed that the population of contacts with full friction mobilization (i.e., for which the Coulomb inequality is satisfied as an equality, ||F T || = µF N ), disappears. Instead, a number of contacts carry force values barely inside the Coulomb cone (typically, ||F T ||/µF N > 0.95). This is enough to allow for a small interval of strains within which the elastic model, assuming friction is irrelevant, applies as a good approximation. It should be recalled that experimental measurements of elastic moduli by static means, along a stress-strain curve, are often carried out in a similar way [69,70]: first equilibrated under static stresses, samples are subjected to small oscillatory probes; while the first cycles tend to cause small amounts of strain to accumulate, the subsequent ones are reproducible and quasielastic. Upon resuming the compression curve (most usually a triaxial compression test) at fixed strain rate, the initial slope of the stress-strain curve coincides with the elastic modulus.

Another way to observe an elastic response is to reverse the loading direction [33]. Anelasticity being largely due to friction mobilization, reversing the sign of strain rate ε1 tends to cause tangential relative displacements to change sign, thereby bringing back contact forces inside the Coulomb cone. Consequently, upon gradually applying a negative ε1 (with due caution, keeping accelerations very small), the obtained stress-strain curve exhibits a quasi-elastic range which is larger than in the forward direction (typically by one order of magnitude), as shown in Fig. 15.

(a) Like in experiments, on resuming the strain rate con- trolled compression after a static elastic probe, the evolution of stresses versus ε 1 in our numerical oedometric tests tends to return to the previous (nonelastic) behavior. Fig. 16 is a plot of stress increment ∆σ 2 versus ∆σ 1 following an elastic probe applied to an equilibrium configuration. The slope of this plot coincides with C 12 /C 11 for small stress increments, and gradually approaches K 0 for larger ones. (K 0 coincides with ratio ∆σ 2 /∆σ 1 along the oedometric compression curve since it is constant in good approximation for initial states DLo.)

B. Role of friction in oedometric compression

Incremental response

Results of Figs. 14 and 16 make it obvious that oedometric compression curves are not ruled by the quasielastic behavior evidenced for small stress or strain increments about a prestressed, well-equilibrated configuration. Departures from this elastic regime (as investigated previously in the isotropic case [30]) are due to frictional forces and, possibly, to contact network instabilities. As visualized in Fig. 17, the elastic response coincides, on resuming the oedometric compression from an intermediate (equilibrated) configuration, with that of a fixed contact network in which no contact creation occurs (although some of the existing contacts may open), and frictional sliding is forbidden (setting µ to an infinite value). We also investigated the response of fixed contact net-0.00 0.05 0.10 0.15 0.20 0.25 0.30 works (forbidding contact creation), and observed them not to differ from the full response on the scale of Fig. 17. The gradual departure from the elastic response is thus mainly due to frictional sliding.

Friction mobilization

Contact sliding is thus the major cause of the nonelastic nature of the mechanical response in oedometric compression. How the sliding (or full friction mobilization) status of a contact correlates to its orientation is shown in Figs. 18 and19. Fig. 18, a plot of average friction mobilization versus |n 1 | in contacts sharing normal unit vector n, shows a significantly greater proximity to the sliding limit, on average, near n 1 = 0, i.e. for nearly transversely oriented normals, close to the plane of directions 2 and 3. It also shows, as noted in isotropic systems [35], that friction mobilization tends to be larger in contacts carrying small forces. Just like in isotropic compression [35], we could also observe larger levels of friction mobilization in systems with smaller coordination numbers. Fig. 19, recording measurements carried out during controlled strain rate oedometric compression, shows moreover that the proportion of exactly sliding contacts, or even almost sliding ones (with ||F T || µF N close to 1), reaches its maximum for directions close to the transverse plane.

This angular variation of friction mobilization might seem surprising, as, from the macroscopic strain field, one does not expect any tangential displacement in the contacts oriented in the transverse plane. The assumption of uniform strain might however provide some insight. Let us write the relative displacement at the contact between grains i and j, δu ij , as

δ uij = ε • r ij = (R i + R j )ε • n ij ,
the dot denoting a derivative with respect to time.Given the uniaxial strain tensor ε = ε 1 e 1 ⊗ e 1 , the normal and tangential components of δ uij read

δ uN ij = (R i + R j ) ε1 (n (1) ij ) 2 δ uT ij = (R i + R j ) ε1 n (1) ij (e 1 -n (1) ij n ij ). ( 20
)
Thus the corresponding elastic force component derivatives are such that their normal component vanishes faster than their tangential one as n to carry tangential forces varying faster than the normal ones, and hence are likely to reach the Coulomb sliding condition more easily.

C. Discussion

To summarize our observations, we have checked that quasi-elasticity, characterizing the response of a contact network in which friction mobilization might be neglected, only applies to well-equilibrated configurations (in which full friction mobilization is lost). Anelasticity is further related to lack of reversibility in unloading in Sec. VI. On the compression curve, strains are associated to some frictional sliding, distributed among all directions of normal vector n but occurring more frequently for those close to the transverse plane. Frictional sliding reduces the apparent stiffness of the material. Unlike material deformation under deviatoric load [71], the anelastic response under oedometric compression does not appear to involve large scale internal rearrangements and failure of contact networks.

Yet, K 0 , the stress ratio, is often related to global failure conditions. First, the ratio of principal stresses, σ 1 /σ 2 , cannot exceed an upper bound 1/R a related to the internal friction angle ϕ as

1 R a = 1 + sin ϕ 1 -sin ϕ , (21)
thereby setting a lower bound R a to the possible value of K 0 (R a and R p = 1/R a are referred to as the active and passive principal stress ratios, in the context of sustaining wall engineering). ϕ (or, equivalently, R a ) is most usually, in simulations as well as in experiments, measured in triaxial compression, which consists in compressing in direction 1 while maintaining stresses σ 2 = σ 3 constant, which involves lateral expansion. In spherical grain assemblies with intergranular friction coefficient µ = 0.3, ratio 1/R a , which depends on the initial state, does not exceed 2.5 [72]. Thus K 0 should exceed 0.4 = 1/R a , which from Fig. 14 is always larger than C 12 /C 11 : the compression could not be elastic.

A second way in which K 0 is often linked to internal friction is through the Jaky relation:

K 0 = 1 -sin ϕ, (22)
which many experimental [37,58,59] and numerical [38] works attempted to check, with varying success. We did not systematically test relation 22 -which is somewhat problematic as K 0 is not constant in general (and would require a nonambiguous definition of ϕ as well). However, let us note that in state DLo for which K 0 remains, in good approximation, constant as σ 1 increases, ( 22) would yield ϕ 30 • , which is notably larger than the values recorded for the internal friction angle in spherical bead assemblies [72].

VI. UNLOADING AND COMPRESSION

CYCLES.

One major issue in oedometric compression is reversibility. It is often assumed that the compression is a plastic, irreversible process. Unloading, on the other hand -i.e., reversing the sign of ε1 or decreasing σ 1 after it has increased to some maximum value, σ p 1 -is often regarded as elastic. Strain (or density), varies less. As the initial value of σ 1 is retrieved, the system has not recovered its initial density, an irreversible density increase is observed. Then, upon reloading, the same "elastic" stress-strain path is retraced as in the previous unloading branch, until σ p 1 , the preconsolidation stress is attained. Any further stress increase beyond σ p 1 results in an additional plastic response, with a faster increase of ε 1 . Such is the classical behavior of sands and other soils in oedometric compression, as described in treatises [START_REF] Mitchell | Fundamentals of Soil Behavior[END_REF] and textbooks [START_REF] Biarez | Elementary Mechanics of Soil Behaviour[END_REF].

A. Density and coordination

To ease comparisons with the literature, in this section we describe density changes in terms of the void ratio e, defined as

e = -1 + 1 Φ . (23)
Fig. 20 displays the variations of void ratio e in the oedometric compression cycle, in which the compression described in Sec. IV is followed by a decompression, with the procedure described in Sec. II B, down to the lowest stress level σ 1 = 10 kPa. The void ratio change (Fig. 20(a)) is almost reversible, especially in dense systems -an observation which strikingly differs from the classically reported behavior of sands (see, e.g., [73]).

In fact, similar irreversible density increases under oedometric loads as in laboratory experiments on sands were to our knowledge never retrieved in DEM simulations in which grains interact merely by elasticity and Coulomb friction in their contacts. Cohesive DEM models, on the other hand, do exhibit large irreversible density increases under isotropic [54,74] or oedometric [75] loads, and behave similarly to laboratory powders, clays or sands, with the preconsolidation stress ruling the onset of further plastic compaction. Cohesionless systems, as dealt with in simulations, appear to lack some modeling ingredient to exhibit similar plasticity in compression as sands, most likely some form of plasticity or damage at the contact scale. This interpretation is confirmed by the experimental observations reported in Ref. [36]: assemblies of smooth beads, in oedometric compression, deform much less than assemblies of angular, irregular shaped particles, unless the beads break under very high stress. Contacts through small asperities tend to exhibit breakage or damage under lower stresses.

In the present numerical study, small irreversible strain changes are nevertheless observed. The lack of reversibility is also, as amply demonstrated in Sec. V, evidenced by the departure from elastic response in compression. Table II compares the total axial strain, in compression and decompression, with an elastic strain, evaluated as

ε el 1 = σ max 1 σ min 1 dσ 1 C 11 (σ 1 ) (24)
The values of C 11 have to be interpolated for all values of ) and measured in simulations along loading (ε ld 1 ) and unloading (ε uld 1 ) paths .

ε el 1 ε ld 1 ε uld 1 Lo 1.
σ 1 along the loading curve to evaluate the elastic strain according to (24). It is of the same order as the measured strain, but with a relative difference of order one. The compression cycle is most conspicuously irreversible as regards the internal microstructure of the system. Thus the coordination number (Fig. 20(b)), if initially large (as in DHo systems) decreases to values nearly as low as in poorly coordinated initial states, either dense (DLo) or loose (LLo). This behavior is quite similar to the one reported in isotropic compression [35].

This decrease of coordination number after unloading also occurs for smaller compression cycles (smaller maximum axial stress). Fig. 21 thus shows, both on the void ratio and on the coordination number, the effect of unloading from σ 1 values of 31.62 kPa, 316.2 kPa, and 3.162 MPa on the primary compression curve, in addition to the maximum stress 31.6 MPa. All compression cycles produce a decrease in z once σ 1 returns to its initial low value, the larger the wider the covered stress interval.

The first graph in Fig. 21, with a void ratio scale appropriate for DHo states, shows small density changes after a stress loop. Although in some cases the final density is smaller than the initial one (see also Tab. II), the work done in the stress loop, as evaluated by the (algebraic) area under the stress-strain curve, is of course positive, signaling energy dissipation (this is discussed in [35] in the case of isotropic compression).

B. K0 and anisotropies

In unloading, a plot of K 0 versus σ 1 , as shown in Fig. 22, first signals a gradual decrease in stress anisotropy, with K 0 increasing towards 1. In initially isotropic systems, as well as for DHo (only marginally anisotropic), the transverse directions become major principal stress directions and K 0 takes values larger than 1 (up to about 1.6 in the DHi case). DLo states, however, after a moderate decrease of K 0 , keep direction 1 as the principal stress direction (and so do LLo ones).

Loose states tend to lose their stress anisotropy for the lower stress values in decompression, which might partly be due to some instability as the load is decreased onto fragile networks.

Since formula (17) still provides a good prediction of K 0 values on the unloading branch of the cycle, the different in evolutions of K 0 upon decompressing might be ascribed to different variations of anisotropy parameters c2 and f2 : c2 changes, which request changes in the contact network, are slower than changes of f2 , which are obtained on simply redistributing forces. Fig. 23 shows that fast increases of K 0 in decompression correspond to quickly evolving force anisotropy parameters f2 in systems DHi, DLi, LLi and DHo.

C. Elastic moduli

Fig. 24 shows the variation of elastic moduli C 11 and C 12 in the compression cycle. As announced, moduli are roughly proportional to σ 1/3 1 , with, possibly a somewhat faster increase in compression associated with changes in coordination number, and some effects of fabric and force anisotropies, to be investigated in [START_REF] Khalili | Submitted to Phys Rev E[END_REF].

The evolution of elastic moduli during unloading phases reveals their dependence on the coordination number and anisotropy, rather than density. The evolution of C 11 in system DHo in loading and unloading parallels that of its coordination number z (see Fig. 20). While C 11 values on the compression branch is significantly larger for DHo systems than for DLo or LLo ones, this difference vanishes, just like the coordination number difference, upon returning to the initial stress down the decompression branch. The final value of C 11 , after the loading cycle, is even lower for DHo than for LLo, de- spite the higher density -a difference to be attributed to the different anisotropies shown in Fig. 23. Compared to coordination numbers or fabric anisotropies, elastic moduli are easier to measure in the laboratory, and some of our observations could thus be checked.

D. Further compression cycles.

Under varying axial stress σ 1 , oedometrically compressed granular materials thus undergo complex, irreversible evolutions, and one should in principle investigate the effects of arbitrary load histories, in which σ 1 may be increased or decreased, over any sequence of load intervals. We report here on a (limited) investigation of the effects of repeating the same compression cycle in systems DHo and DLo.

It is interesting to see whether the compression cycle is retraced upon compressing again: one may wonder to what extent the memory of the initial state survives re- peated cycles, and whether systems sharing the same initial density will tend to approach a common limit state. Fig 25 shows that the slight density difference between DLo and DHo does not appear to gradually vanish under repeated compression cycles. Meanwhile, the coordination number of the DHo state, after its strong decrease in the first cycle, oscillates in the following ones between values that remain somewhat larger than the ones observed in the initially poorly coordinated system DLo.

Although the difference of coordination number between DHo and DLo is greatly reduced after the first cycle, the stress anisotropy of the final states is then quite different, as noted previously (see Fig. 22), with K 0 > 1.4 for DHo, and K 0 < 0.8 for DLo. This difference in the evolution of the principal stress ratio, as shown in Fig. 26, does not tend to disappear under repeated cycles: while K 0 in DHo systems oscillates between about 0.6 at large σ 1 and nearly 1.5 under low stress, it oscillates below 1 for DL0, with a systematic increasing tendency. Thus, under repeated decompression and recompression steps, systems DLo and DHo still behave differently: while the DHo state seems to approach a limit cycle, with the same values of Φ, z and K 0 , the observation of K 0 values for DLo reveals a gradual evolution towards gradually less anisotropic, and possibly denser states.

VII. CONCLUSIONS

We now summarize and comment our observations, suggesting a few perspectives to the present study. We carried out systematic numerical simulations of quasistatic oedometric compression of model granular materials, in which contact mechanics does not involve other ingredients than (suitably simplified) Hertz-Mindlin elasticity and Coulomb friction (with friction coefficient µ = 0.3). The material is first assembled in initial states varying in density, coordination number (which might be large or small, depending on preparation, in dense systems, as in isotropic grain assemblies [30][31][32]), and anisotropy. In the compression cycle, axial stress σ 1 varies by a factor larger than 3000, corresponding, for glass beads, to the range 10 kPa ≤ σ 1 ≤ 31.6 MPa. The observed behaviors prove somewhat more complicated than superficial observations would seem to indicate, with nontrivial initial state dependence and anelastic response.

Although the strain response (or the change in solid fraction) under growing axial stress seems nearly reversible, the internal state of the material does evolve irreversibly, as apparent in the variations of coordination numbers and anisotropy parameters. Compared to the behavior of sands as described in the geomechanics literature, the stress-strain irreversibility (usually described as a plastic response) is much smaller, a difference we attribute to the absence of plasticity or damage in the implemented contact model. Just like in isotropic compression [35], coordination numbers, if ini-tially high, tend to decrease in a compression cycle, once the initial stress value is retrieved, the more the larger the maximum stress value in compression. The stress anisotropy, as expressed by ratio K 0 , is in general not constant. It varies little with σ 1 in systems assembled under similar one-dimensional compressions of gases, except in the case of a high coordination number (possibly unrealistically large), as obtained on suppressing friction in the assembling stage. We expect gravitydeposited systems, in conditions ensuring homogeneous density and microstructure, to behave similarly. Stress ratio K 0 = σ 2 /σ 1 is correctly predicted, in all configurations along the loading or unloading curves, by a formula involving the leading order anisotropic terms in Legendre polynomial expansions of normal vector-dependent contact densities and average normal force value distributions over the unit sphere. It should thus be possible to predict the mechanical response and the internal material evolution if the evolution of axial strain ε 1 could be related to anisotropy parameters. This latter task seems however arduous, given the complex non-elastic strain response of the material. Elastic moduli, although measurable upon applying small load increments onto well-equilibrated configurations along the compression path, do not correctly predict the slope of the oedometric compression curve, or the stress ratio, K 0 . This conclusion might seem paradoxical, since the compression curve is supposed to be quasistatic, i.e., consisting of a continuous sequence of equilibrium states. We attribute this apparent contradiction to the subtle role of the very small creep step preceding, in numerical simulations as well as in the laboratory, the static measurement of elastic properties. In simulations this creep stems from the small distance to equilibrium of transient configurations along a strain-rate controlled loading path. This distance decreases as the strain rate, expressed in dimensionless form by inertial number I, decreases. As well equilibrated, static states are obtained with good accuracy, such that the stiffness matrix of the contact network is well behaved, friction mobilization is lost, and for a small, non-vanishing stress increment, it is a good approximation to assume all contacts to behave elastically. In the laboratory, although strain rates are considerably smaller, and intermediate configurations are likely closer to equilibrium, some creep also takes place, possibly caused by other phenomena at contact scale, with similar results that quasielastic relations between small stress and strain increments might be subsequently measured. Creep phenomena, as observed in real laboratory materials, in general, would deserve more detailed investigations, although the elucidation of their origin is likely to involve little known micromechanical ingredients at the contact scale.

The elastic properties being easier to measure than fabric variables and coordination numbers, we provide in the companion paper [START_REF] Khalili | Submitted to Phys Rev E[END_REF] a more detailed study of the relations of all five independent moduli in the transversely anisotropic configurations obtained by oedomet-ric compression to microstructural variables. A comparison of elastic moduli in numerical and experimental works should help understanding which type of numerical preparation scenario produces initial states closer to experimental ones.

Although strains predicted by the elastic response are of the same order of magnitude as observed strains, the difference is important, and strongly affects the value of stress ratio K 0 . The anelastic response is mainly due to friction mobilization, which, although distributed over all contact orientations, surprisingly affects the most the contacts with normal direction close to the transverse plane. Contact network instabilities, avoided thanks to new contact creations, do not seem to play an important part. Detailed, strictly quasistatic analyses [51,71,76,77] of elastic-frictional response of contact networks to oedometric loads could be carried out to relate microscopic frictional sliding to macroscopic behavior.

Another remarkable result of the present numerical study is the persistent effect of the assembling process and the resulting initial state characteristics: the difference, e.g., between high coordination and low coor-dination dense systems is not lost after several compression cycles. Our numerical results, as regards the evolution of K 0 , differs somewhat from experimental observations, which sometimes report a quicker convergence to a common value of this stress ratio for different initial states. This is likely due, like irreversibility, to the absence of contact-level plasticity or damage in the numerical model. Some features of assemblies of nonspherical, rough or angular grains as probed in oedometric compression of sands might need to be modeled -with the same identification difficulties as for creep -if the effects of growing stress intensity are to be quantitatively described.

Fig. 1 FIG. 1 .

 11 FIG. 1. (Color online) h-dependent coordination number z(h)for the different initial states (σ1 = 10 kPa or κ 39000).

1 FIG. 2 .

 12 FIG. 2. (Color online) Evolution of solid fraction Φ with axial stress σ1 in oedometric compression for different initial states.

1 FIG. 3 .

 13 FIG. 3. (Color online) Coordination number z (a) and rattler fraction x0 (b) versus axial stress σ1 along oedometric loading path for the different sample preparations.

  FIG. 4. (Color online) Probability density function of f = F N / F N in oedometric compression of initial state DLo, for different axial stress levels.

1 FIG. 5 .

 15 FIG. 5. (Color online) Evolution of reduced second moment of normal contact forces, (F N ) 2 / F N 2 in oedometric compression.

1 FIG. 6 .

 16 FIG.6. (Color online) K0 versus axial stress σ1 in oedometric compression (dots joined by solid lines, as specified in legend). Dots joined by dashed lines: predictions of Eq. 17 (see Sec. IV E). (a) oedometrically assembled systems ("XYo"). (b) isotropically assembled ones ("XYi").

truncated at order 2 FIG. 7 .

 27 Fig.7shows that expansion (10) truncated at order 2

1 FIG. 8 .

 18 FIG. 8. (Color online) Anisotropy parameter c2, versus σ1 or κ -1 , for all six different initial states.

FIG. 9 .

 9 FIG. 9. (Color online) Anisotropy of orientation of pair of near neighbors as expressed by h-dependent coefficient c2(h). (a) σ1 = 10 kPa (or κ 39000, initial states). (b) σ1 = 10 MPa (or κ 390).

1 FIG. 10 .

 110 FIG. 10. (Color online) Force anisotropy parameter versus σ1 or κ -1 , for all six different initial states.

FIG. 11 .

 11 FIG. 11. (Color online) Ratio of stress increments ∆σ1 to initial stress σ0, versus strain increment ∆ε1. Dots: DEM results, after system equilibrates under application of growing ∆σ1. Dotted lines: elastic prediction (slope C11/σ0). Data recorded for different σ0 values as indicated. (a) in one DLo sample and (b) a DHo one.

FIG. 12 .

 12 FIG.12. (Color online) Ratio of stress increments ∆σ2 to initial stress σ0, versus strain increment ∆ε1. Dots: DEM results, after system equilibrates under application of growing ∆σ1. Dotted lines: elastic prediction (slope C12/σ0). Data recorded for different σ0 values as indicated, in one DHi sample

FIG. 13 .FIG. 14 .

 1314 FIG. 13. (Color online) Quasielastic range, defined as ε1 interval for which (19) holds with 5%, versus equilibrium stress in probed system. (a) axial stress response (∆σ1). (b) lateral stress response (∆σ2). Dashed lines have slopes 2/3.

FIG. 15 . 0 C 1 2 / C 1 1 FIG. 16 .

 15012116 FIG.15. (Color online) Quasielastic range (same definition as for the results of Fig.13) versus σ1 in one DLo system, for loading (positive) and unloading (negative) axial strain increments.

1

 1 FIG. 17. online) Oedometric loading curves from equilibrium state under σ1 = 31.6 kPa in a DLo sample. "Ncc, µ = ∞" labels simulations carried out without creation of any new contact, and infinite friction (no sliding). The linear elastic response of the initial contact network is shown as dashed straight line.

F

  FIG. 18. (Color online) Average level of friction mobilization, ||F T || µF N , in contacts sharing common normal orientation n, versus |n1|, in state DLo under σ1 = 1 MPa. The same average quantity is also shown separately for contacts carrying normal forces larger os smaller than the average normal force.

( 1 )FIG. 19 .

 119 FIG. 19. (Color online) Proportion of contacts for friction mobilization ||F T || µF N is equal to 1 or exceeds threshold 0.95 in contacts sharing common normal orientation n, versus |n1|, in state DLo under σ1 = 91 kPa (these data correspond to one intermediate state in the course of strain-rate-controlled compression).

FIG. 20 .

 20 FIG. 20. (Color online) Void ratio e (a) and coordination number z (b) versus σ1 in oedometric compression cycle.

1 FIG. 21 .

 121 FIG. 21. (Color online) Void ratio e (a) and coordination number z (b) versus σ1 in oedometric compression cycle on one DHo sample.

1 FIG. 22 .

 122 FIG. 22. (Color online) K0 versus σ1 in oedometric compression cycle. (a) oedometrically assembled systems ("XYo"). (b) isotropically assembled ones ("XYi").

1 FIG. 23 .

 123 FIG. 23. (Color online) c2 (a) and f2 (b) versus σ1 in oedometric compression cycle.

1 FIG. 24 .

 124 FIG. 24. (Color online) C11 (a) and C12 versus σ1, in loading cycle. Black dashed line has slope 1/3.

1 FIG. 25 .

 125 FIG. 25. (Color online) Solid fraction Φ (a) and coordination number z (b) versus σ1 in oedometric compression cycles, with three unloading and reloading steps after the first compression to maximum value of axial stress.

1 FIG. 26 .

 126 FIG. 26. (Color online) K0 versus axial stress σ1 in cycles of Fig. 25. (a) in a DHo system and (b) in a DLo one.

  -8 F 1 D 1 and controlled stress components should be measured equal to their set values, within relative error 10 -4 .

		LLo LLi DHo DHi DLo DLi
	Φ 0.584 0.589 0.639 0.638 0.634 0.637
	z	4.22 4.14 5.98 5.99 4.06 4.17
	z	4.63 4.63 6.07 6.07 4.54 4.65
	x0(%) 8.8 10.3 1.5 1.3 10.4 10.37
	K0 0.72	1	0.94	1	0.51	1
	TABLE I. Solid fraction Φ, coordination numbers of all grains,
	z, and of nonrattler ones, z * , rattler fraction x0, and stress
	ratio K0 for the initial states, according to preparation proce-
	dure. K0 = 1 in all isotropically assembled states (denoted as
	"XYi" ). Values averaged over 3 configurations of N = 4000
	grains (sample to sample differences lie below given accuracy
	level, except for x0).

  This figure makes it clear that both quantities are quite different, with much lower values of C 12 /C 11 . As remarked before, modulus C 11 overestimates the variation of axial stress, ∆σ 1 with ε 1 , while C 12 underestimates the variation of lateral stress, ∆σ 2 . Both effects entail that C 12 /C 11 is smaller than K 0 = ∆σ 2 /∆σ 1 . Furthermore, ratio C 12 /C 11 is nearly constant in systems DHi and DHo, while K 0 changes to a large extent. The opposite is true for the four other states, in which C 12 /C 11 changes much more than K 0 . Elastic moduli, in general, are thus quite uncorrelated to K 0 .

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.