Macroscopic model with anisotropy based on micro–macro information

Physical experiments can characterize the elastic response of granular materials in terms of macroscopic state variables, namely volume (packing) fraction and stress, while the microstructure is not accessible and thus neglected. Here, by means of numerical simulations, we analyze dense, frictionless granular assemblies with the final goal to relate the elastic moduli to the fabric state, i.e., to microstructural averaged contact network features as contact number density and anisotropy. The particle samples are first isotropically compressed and then quasi-statically sheared under constant volume (undrained conditions). From various static, relaxed configurations at different shear strains, infinitesimal strain steps are applied to “measure” the effective elastic response; we quantify the strain needed so that no contact and structure rearrangements, i.e. plasticity, happen. Because of the anisotropy induced by shear, volumetric and deviatoric stresses and strains are cross-coupled via a single anisotropy modulus, which is proportional to the product of deviatoric fabric and bulk modulus (i.e., the isotropic fabric). Interestingly, the shear modulus of the material depends also on the actual deviatoric stress state, along with the contact configuration anisotropy. Finally, a constitutive model based on incremental evolution equations for stress and fabric is introduced. By using the previously measured dependence of the stiffness tensor (elastic moduli) on the microstructure, the theory is able to predict with good agreement the evolution of pressure, shear stress and deviatoric fabric (anisotropy) for an independent undrained cyclic shear test, including the response to reversal of strain.


Introduction
Granular materials behave differently from usual solids or fluids and show peculiar mechanical properties like dilatancy, history dependence, ratcheting and anisotropy [24,25,28,30,39,40,60,70,75,76,87]. The behavior of these materials is highly nonlinear and involves plasticity also for very small strain due to rearrangements of the elementary particles [4,15,22]. The concept of an initial purely elastic regime (small strain) for granular assemblies is an issue still under debate in the mechanical and geotechnical communities. On the other hand, approaches that neglect the effect of elastic stored energy are also questionable, i.e., all the work done by the internal forces is dissipated. Features visible in experiments, like wave propagation, can hardly be described without considering an elastic regime. In a general picture, both the deformations at contact and the irrecoverable rearrangements of the grains sum up to the total strain. The former represents the elastic, reversible contribution to the behavior of the material. That is, for very small strain, the response of a finite granular system in static equilibrium can be assumed to be linearly elastic [20,42,59,70], as long as no irreversible rearrangements take place.

Numerical simulation
The Discrete Element Method (DEM) [2,30,47] helps to study the deformation behavior of particle systems. At the basis of DEM are laws that relate the interaction force to the overlap (relative deformation) of two particles. Neglecting tangential forces, if all normal forces f i acting on particle i, from all sources, are known, the problem is reduced to the integration of Newton's equations of motion for the translational degrees of freedom: with the mass m i of particle i, its position r i , velocity v i (=ṙ i ) and the resultant force f i = c f i c acting on it due to contacts with other particles or with the walls, and the acceleration due to gravity, g (which is neglected in this study). The force on particle i, from particle j, at contact c, has normal and tangential components, but the latter are disregarded in this study to focus on frictionless packings.
For the sake of simplicity, the linear visco-elastic contact model for the normal component of force is used, where k is the spring stiffness, γ is the contact viscosity parameter, δ = d i + d j /2 − r i − r j ·n is the overlap between two interacting species i and j with diameters d i and d j ,n = r i − r j / r i − r j andδ is the relative velocity in the normal direction. In order to reduce dynamical effects and shorten relaxation times, an artificial viscous background dissipation force f b = −γ b v i proportional to the moving velocity v i of particle i is added, resembling the damping due to a background medium, as e.g., a fluid. The standard simulation parameters are N = 9,261(=21 3 ) particles with average radius r = 1 (mm), density ρ = 2,000 (kg/m 3 ), elastic stiffness k = 10 8 (kg/s 2 ), particle damping coefficient γ = 1 (kg/s), background dissipation γ b = 0.1 (kg/s). Note that the polydispersity of the system is quantified by the width (w = r max /r min = 3) of a uniform size distribution [24], where r max and r min are the radii of the biggest and smallest particles, respectively.

Coordination number and fraction of rattlers
In order to link the macroscopic load carried by the sample with the active microscopic contact network, all particles that do not contribute to the force network are excluded. Frictionless particles with <4 contacts are thus "rattlers," since they cannot be mechanically stable and hence do not contribute to the contact or force network [24,30,39]. The classical definition of coordination number is C = M/N , where M is the total numbers of contacts and N = 9,261 is the total number of particles. The corrected coordination number is: C * = M 4 /N 4 , where, M 4 is the total number of contacts of the N 4 particles with at least four contacts. Moreover, we introduce here the reduced number of contacts M p 4 , where contacts related to rattlers are excluded twice, as they do not contribute to the stability of both the rattler and the particle in contact with it. Hence, M The total volume of particles is N P=1 V P = (4/3)π N r 3 , where r 3 is the third moment of the size distribution [24,39], and the volume fraction is defined as ν = (1/V ) N P=1 V P , where V is the volume of the periodic system.

Macroscopic (tensorial) quantities
Here, we focus on defining averaged tensorial macroscopic quantities-including strain-, stress-and fabric (structure) tensors-that provide information about the state of the packing and reveal interesting bulk features.
By speaking about the strain-rate tensorĖ, we refer to the external strain that we apply to the sample in a time interval dt. The isotropic part of the infinitesimal strain tensor ε v [24,30,39] is defined as: where ε αα =ε αα dt with αα = x x, yy and zz are the diagonal components of the tensor in the Cartesian x − y − z reference system. The trace integral of 3ε v is denoted as the volumetric strain ε v , the true or logarithmic strain, i.e., the volume change of the system, relative to the initial reference volume, V 0 . On the other hand, from DEM simulations, one can measure the "static" stress in the system [9] as averaged over all the contacts in the volume V of the dyadic products between the contact force f c and the branch vector l c , where the contribution of the kinetic fluctuation energy has been neglected [30,47]. The isotropic component of the stress is the pressure P = tr(σ σ σ )/3. In order to characterize the geometry/structure of the static aggregate at microscopic level, we will measure the fabric tensor, defined as where V P is the volume of particles P, which lie inside the averaging volume V , and n c is the normal unit branch vector pointing from center of particle P to contact c [40,47,86]. We want to highlight that a different convention for the fabric tensor involves only the orientation of contacts as follows [61,69,87]: where N c is the total number of contacts in the system. An approximated relationship between Eqs. (5) and (6) can be derived as: with tr(F o ) = 1. This relation is exactly equal for monodisperse assemblies but largely deviates for assemblies with high polydispersity (see further discussion in Sect. 3). The difference also becomes more significant when the jamming volume fraction [52,79] is approached. In the following, when not explicitly stated, we will refer to Eq. (5), since we combine the effects of volume fraction and number/orientation of contacts, both relevant quantities when the elastic moduli are considered [24]. In a large volume with a given distribution of particle radii, the isotropic fabric, i.e., the trace of F, is proportional to the volume fraction ν and the coordination number C, see Refs. [24,30,39], as where C, C * and φ r have been introduced in previous Sect. 2.1, and g 3 ≈ 1.22 for polydispersity w = 3, being only a weighted factor of order unity of the non-dimensional moments size distribution [24,39,71].

Isotropic and deviatoric parts
We choose here to describe each symmetric second order tensor Q, in terms of its isotropic part (first invariant) and the second, , and third, invariants of the deviator, with Q D 1 , Q D 2 and Q D 3 eigenvalues of the deviatoric tensor Q D = Q − (tr(Q)/3)I. We use the following definition (of the Euclidean or Frobenious norm) to quantify with a single scalar the magnitude of the deviatoric part [39,40] of Q: where Q x x , Q yy and Q zz are its diagonal, Q xy , Q yz and Q zx its off-diagonal components and the deviators ε dev , σ dev and F dev refer to strain E, stress σ σ σ and fabric F, respectively. Fsgn (Q) is the sign function that relates the tensorial quantity to be measured, Q, with the reference tensor that describes the (strain-or stress-controlled) path applied to the sample, H 0 : with ":" being the inner product between the two tensors. For a given, complex deformation path, the reference tensor H 0 must be chosen in a convenient way, in order to take into account both the actual loading path and/or the previous deformation history of the sample. In the special case of an undrained shear test, as introduced later in Sect. 3, we use as reference the director of the strain-rate H 0 =D(−Ė) ∝ (−1, 1, 0), where only the diagonal terms are given (and the normalization can be ignored), see Eq. (10) below, so that Fsgn simplifies to with x-wall expanding, y-wall compressing and z-wall non-mobile [39]. We want to point out here that, during a deformation, the response of stress σ σ σ and fabric F is opposite in sign to applied strain-rateĖ. Unless mentioned explicitly, we will be using a sign convention for strain (isotropic δε v = −(1/3)tr(δE) and deviatoric δε dev = −δ E dev ), such that consistently a positive strain increment leads to a positive stress and fabric response. Finally we note that in this work we will use k * = k/ (2 r ) to non-dimensionalize pressure P and deviatoric stress σ dev to give P * and σ * dev , respectively, and will be referring to deviatoric stress as shear stress. 1

Volume-conserving (undrained) biaxial shear test
In this section, we first describe the sample preparation procedure and then the details of the numerical shear test.
The initial configuration is such that spherical particles are randomly generated, with low volume fraction and rather large random velocities in a periodic 3D box, such that they have sufficient space and time to exchange places and to randomize themselves. This granular gas is then compressed isotropically, to a target volume fraction ν 0 = 0.640, sightly below the isotropic jamming volume fraction [30,39,52] ν c ≈ 0.658 and then relaxed to allow the particles to fully dissipate their potential energy [30,39]. 2 The relaxed state is then compressed (loading) isotropically from ν 0 to a higher volume fraction of ν = 0.82, and decompressed back (unloading) to ν 0 [30,39].
The preparation procedure, as described above, provides many different initial configurations with volume fractions ν i , each one in mechanical equilibrium. Starting from various ν i chosen from the unloading branch [30,39], the samples are then sheared keeping the total volume constant, that is with a strain-rate tensoṙ whereε yy = 2.84 10 −5 [s −1 ] is the strain-rate (compression > 0) amplitude applied to the moving xand y-walls, while the third z-wall is stationary. Our shear test, where the total volume is conserved during deformation, resembles an undrained test typical in geotechnical practice [87]. The chosen deviatoric path is on the one hand similar to the pure shear situation, and on the other hand allows for simulation of the biaxial element test [56,66] (with two walls static, while four walls are moving, in contrast to the more difficult isotropic compression, where all the six walls are moving). Pure shear is here used to identify constant volume deviatoric loading with principal strain axis keeping the same orientation as the geometry (cuboidal) of the system for the whole experiment. In this case, there is no rotation (vorticity) of the principal strain (rate) axis, and no distortion/rotation of the sample due to shear deformation. Different types of volume-conserving deviatoric deformations can be applied to shear the system, but very similar behavior has been observed [30], in terms of shear stress.

Evolution of stress
The evolution of non-dimensional pressure P * with deviatoric strain ε dev is presented in Fig. 1a during undrained shear tests for some exemplary volume fraction. For frictionless systems analyzed here, only a slight variation of the pressure is observed at the beginning of the test, due to the development of anisotropy in the sample, after which P * remains constant. 3 Both the (small) initial pressure change and the final saturation value vary with the vicinity of ν to the jamming volume fraction ν c . Interestingly, depending on the volume fraction, some of the samples show increase in the pressure (dilatancy) with respect to the initial value and some other decrease (compactancy), as shown in Fig. 1. This supports the idea of a certain threshold value ν p d = 0.79, as shown in Fig. 2a, where the pressure of the system changes behavior, similarly to the switch between volumetric dilation and contraction visible in triaxial tests. The evolution of the (non-dimensional) shear stress σ * dev during shear, as function of the deviatoric strain ε dev , is shown in Fig. 1b, for the same simulations as in Fig. 1a. The stress grows with applied strain until an asymptote (of maximum stress anisotropy) is reached where it remains fairly constant, with slight fluctuations around the maximum σ * dev [10]. The growth rate and the asymptote of σ * dev , both increase with ν.

Evolution of fabric
Complementary to stress, in this subsection we study the evolution of the microstructure in the sample during the volume-conserving shear test. Figure 3a shows that the isotropic fabric F v behaves in a very similar fashion as P * , with a slight increase/decrease at the beginning, followed by saturation stage, whose value increases continuously with ν. Figure 2b shows that the difference between the initial value of F v and its saturation value changes sign when a certain volume fraction, ν F d = 0.755, is reached. Note that ν F d = ν p d , that further confirms the independent evolution of F F F and σ σ σ .
From Eq. (8) F v is proportional to the product of volume fraction ν, that remains unchanged during deviatoric deformations, and coordination number C, that varies only slightly for sheared frictionless systems [30]. Note that as C = C * (1 − φ r ), knowing the (empirical) relations of C * and φ r with volume fraction, as presented in Refs. [30,39], we can fully describe the isotropic fabric state. In this study, we assume F v to stay constant during the shear test. This assumption will be used later in Sect. 5 for the prediction of a cyclic shear test. However, the small changes in F v or P * can be associated with a (small) change in the jamming volume fraction [41].  The evolution of the deviatoric fabric, F dev , as function of the deviatoric strain is shown in Fig. 3b during shear for five different volume fractions. It builds up from different random small initial values (due to the initial anisotropy in the sample that develops during preparation) and reaches different maxima. The deviatoric fabric builds up faster at lower volume fractions but the maximal values are higher for smaller volume fractions, qualitatively opposite to the evolution of σ * dev [10]. As mentioned in Sect. 2.2 the validity of Eq. (7), that relates the two different definitions of fabric, depends on polydispersity. In order to check the relation, in Fig. 4 the evolution of the three eigenvalues of the fabric tensor is plotted, for both definitions, Eqs. (5) and (6), during the volume-conserving shear test, for three different values of polydispersity w = 1, 2 and 3. For all polydispersities, the chosen volume fraction ν = 0.685 is close to the jamming points, that slightly vary with w [39]. The difference between the definitions of fabric becomes higher for higher polydispersity w = 3, as in Eq. (6) the contribution of each particle is weighted by its surface area (via its number of contacts), whereas in Eq. (5) it is weighted by the volume. Only for the monodisperse case, the relation is exact, as can be seen in Fig. 4a. The differences are considerable for w = 2 and w = 3, for both compressive and tensile direction, while the non-mobile direction is not affected. Note that the difference of the two fabrics will be smaller for denser systems.

Elastic moduli
In this section, we focus on the evolution of the elastic properties of the material and neglect the plastic contribution to the granular behavior that will be superimposed to the present analysis later in Sect. 5. We first describe the numerical procedure to measure the elastic moduli of the anisotropic aggregate, and later analyze the data and their relation with stress and fabric.

Numerical probes
In a general framework, a possible description for the incremental, elastic behavior of an anisotropic material is where the isotropic and deviatoric components of stress have been isolated and are expressed as functions of ε v and ε dev via a non-dimensional stiffness matrix [26] (by multiplying the moduli with k * , the real stiffnesses can be extracted). B is the classical bulk modulus, and G oct the octahedral shear modulus. The anisotropy moduli A 1 and A 2 provide a cross-coupling between the two parts (isotropic and deviatoric) of stress and strain increments. Equation (11) provides a partial description for the evolving stress and stiffness of a sheared material, as it applies to a triaxial-box configuration (with eigen-system coincident with the axes of the box), where no off-diagonal terms are measured and stress and moduli are assumed to be collinear. Moreover, the increase in stress and stiffness in the out-of-plane direction (z-direction here) due to the non-planar (triaxial) stress state associated with the plane deformation mode, is not independently accounted for. These are rather hidden in the expression for deviatoric stress as proposed in Eq. (9) and used in Eq. (11). However, we have chosen this representation, since advantages are obtained by investigating the elasticity of a granular material (e.g., soil), not by its resistance to direct stresses expressed by Young's modulus and Poisson's ratio, but rather in terms of (purely volumetric and deviatoric) stress response to volume and shape changes, as described by the bulk modulus B and the octahedral shell modulus G oct . This aspect will be further addressed in Sect. 5, where Eq. (11) will be included in the theoretical model.
To study the evolution of the effective moduli during shear, we choose different initial states (forty) as shown in Fig. 5, and apply sufficient relaxation, so that the granular assemblies dissipate the kinetic energy accumulated during the original shearing path. When the states along the shear path are relaxed, a much higher change is visible in σ * dev rather than in F dev , see Fig. 5. This shows that the contact network remains almost intact and F dev does not change; on the other hand, the average particle overlap is more sensitive to the relaxation stage and decreases (in steady state), leading to a finite drop in σ * dev . Then we perform a small strain perturbation to these relaxed anisotropic states, i.e., we probe the samples, and measure the incremental stress response [40,50]. Finally, the elastic moduli are calculated as the ratio between the measured increment in stress and the applied strain. We obtain all the different moduli in Eq. (11), by applying an incremental pure volumetric or pure deviatoric strain and measuring the incremental volumetric or shear stress response: The system is allowed to relax again after the incremental strain is applied, that is both the stress and the change in stress are measured after relaxation [50,52]. Since the numerical probe experiments are carried out with zero contact friction, we are measuring the resistance of the frictionless material [40], where only normal forces are involved. The first big question concerns the amplitude of the applied perturbation to get the elastic response [6,21,72].

How small is small?
In this section, we discuss the amplitude of the perturbations applied to measure the elastic stress response of the granular material. Also, we will discuss the results for larger amplitudes and the threshold between elastic and plastic regimes. Figure 6 (column 1 and 2) shows the changes in non-dimensional pressure δ P * , non-dimensional shear stress δσ * dev , isotropic fabric δ F v and deviatoric fabric δ F dev for different amplitudes of the isotropic perturbation 3δε v , applied to two relaxed states that have been sheared until ε dev = 0.0065 (nearly isotropic configuration: column 1) and ε dev = 0.31 (steady state configuration: column 2), respectively. The data correspond to the shear test with ν = 0.706. The linear elastic response is also plotted (red solid curve) in the whole strain range, as derived from the incremental behavior at very small strain, to give an idea of the deviation from elasticity when strain increases.

Effect of isotropic perturbations 3δε v
δ P * initially increases linearly and smoothly with 3δε v , in agreement with the prediction of linear elasticity. Also the difference between the two initial states (near isotropic and steady state as shown in Fig. 6a, b, respectively) is minimal, meaning that the bulk modulus B (slope of δ P * with 3δε v in the elastic regime) is almost constant. This is not surprising, as we expect B to be dependent on isotropic quantities that stay mostly unchanged during the shear deformation, as discussed in Sect. 4.3. δσ * dev behaves similar as δ P * for small strain, but shows several sharp drops for large strain. These correspond to sudden changes in the coordination number δC * (see Fig. 7a, b), due to rearrangements in the system during the probe. For the nearly isotropic state (Fig. 6e), the ratio of δσ * dev with 3δε v in the linear elasticity regime, i.e., A 2 , is small when compared with the steady state (Fig. 6f). This clearly tells that A 2 evolves during the shear deformation for a given volume fraction, and must be linked with deviatoric quantities.
δ F v increases with 3δε v , with more fluctuations compared with δ P * , for both states considered here, ε dev = 0.0065 (nearly isotropic state Fig. 6i) and ε dev = 0.31 (steady state, Fig. 6j). Moreover, the prediction using Eq. (8) for F v matches the dataset very well. F dev does not change (δ F dev = 0) with increasing 3δε v , until the first rearrangement in the structure occurs (see Fig. 7c, d). After this δ F dev starts to decrease with increasing amplitude 3δε v , faster in the steady state (Fig. 6m) than in the near isotropic state, see Fig. 6n. We note here that, when a non-incremental volumetric strain (3δε v > 10 −4 ) is applied, the system moves from a volume-conserving to a new non-volume-conserving deformation path. As this system is already anisotropic, this leads to a decrease (δ F dev < 0) in deviatoric fabric F dev , opposite to the increase (δσ dev > 0) in deviatoric stress, see Fig. 6e, f, higher in the steady state ( Fig. 6n) than in the nearly isotropic state (Fig. 6m). The last observation suggests that the distance between the volume-conserving and non-volume-conserving configurations increases with ε dev .
Hence, during isotropic compression (increasing 3δε v ) of a pre-sheared (anisotropic) state, both the pressure P * and shear stress σ * dev increase, with pressure increasing much faster leading to a decrease in deviatoric stress ratio s dev = σ * dev /P * . The deviatoric fabric F dev also decreases with isotropic compression of a pre-sheared state, and the decrease is initially faster than the exponential decay of F dev (see Sect. 5 below) with volume fraction ν, as seen in Fig. 6n. This decrease in F dev becomes slower for large strain, as also seen in Fig. 6m. These observations are consistent with the findings of Imole et al. [30], where the authors noticed a decreasing steady state deviatoric fabric and deviatoric stress ratio with increasing volume fraction, or ε v . Figure 6 (column 3 and 4) shows the changes in the same quantities as before for different amplitudes of the deviatoric perturbation δε dev , applied to a relaxed state with volume fraction ν = 0.706 that has been sheared until ε dev = 0.0065 (nearly isotropic configuration: column 3) and ε dev = 0.31 (steady state configuration: column 4).

Effect of deviatoric perturbations δε dev
δ P * increases linearly with δε dev (the slope in the elastic regime is A 1 ), with A 1 much smaller for the nearly isotropic state (Fig. 6c) than for the steady state (Fig. 6d). This shows that A 1 evolves during the shear deformations, like A 2 , for a given volume fraction, and must be linked with the deviatoric state of the system. Moreover, after large deformation, both states show drops in δ P * , which can be linked to the particle rearrangements (see Fig. 7c, d). A nonlinear, irregular behavior shows up for δε dev > 10 −4 , with δ P * positive in case of loose sample (present sample) and negative for dense samples (data not shown), in agreement with 10 -7 10 -6 10 -5 10 -4 10 -3 10 -2 10 -1 10 -7 10 -6 10 -5 10 -4 10 -3 10 -2 10 -1 10 -7 10 -6 10 -5 10 -4 10 -3 10 -2 10 -1 10 -7 10 -6 10 -5 10 -4 10 -3 10 -2 10 -1 nearly isotropic The green line in m and n represents the evolution of change in deviatoric fabric δ F dev in critical state using parameters from Table 3 of Ref. [30], with the assumption that the new state after volumetric deformation is also in critical state. The green line in o and p represents Eq. (18) from Ref. [30] when subtracted from its initial value F 0 dev = 0.03 for o and F 0 dev = 0.113 for p, with the growth rate β F = 39 and F max dev = 0.12 (color figure online) the observations in Fig. 2a. δσ * dev also increases linearly with δε dev , with G oct (slope of the line) slightly higher for the near isotropic state (Fig. 6g) than for the steady state (Fig. 6h). Again, similar to δ P * , δσ * dev shows drops after large deformations, which can be linked to the particle rearrangements (see Fig. 7c, d). In the steady state, the incremental stresses (δ P * and δσ * dev ) increase linearly for very small strain, as the relaxed configuration, starting point for the probes, has lower stress than the main deviatoric path (see Fig. 5a) and the system tends to regain the "missed" stress, when the shear restarts. After the first elastic response, δ P * and δσ * dev fluctuate around zero for larger amplitudes (Fig. 6d, h), as no change in stress is expected with increasing deviatoric strain in the steady state.
δ F v stays mostly zero when small δε dev is applied for both near isotropic and steady state configurations (Fig. 6k, l). With increasing strain amplitude, δ F v increases in the case of a loose sample close to the isotropic state (Fig. 6k), and decreases for denser samples (data not shown), in agreement with the behavior in Fig. 2b. In Fig. 6o, δ F dev for the nearly isotropic state stays zero for δε dev < 10 −4 , when no rearrangements happen and the behavior is elastic, while it reaches a positive finite value for larger amplitude (that coincides with the slope of the curve in Fig. 3b). This finite value drops to zero in the steady state, where no variation of deviatoric fabric is expected with further applied deviatoric strain (see Fig. 6p). When compared with the model predictions in Ref. [30], the simulation data for F dev well match with the theoretical line, where F dev increases due to shear for the near isotropic state, and does not change for the steady state simulation.

Discussion and comparison
Since we are interested in measuring the pure elastic response of the material, we take care that no rearrangements happen in the system during the numerical probe, that is, 3δε v and δε dev are applied only up to 10 −4 (with very slow wall movement ≈ 10 −5 [s −1 ], i.e., smaller than for the main large shear strain preparation experiment). Looking at Fig. 6, we note that much bigger drops appear in the deviatoric response when the isotropic perturbation is applied. Vice versa, the fluctuations/drops are much larger in pressure rather than in shear stress, when we deal with deviatoric perturbations. It is worthwhile to mention here that we have tested our method by applying strain perturbations in opposite directions i.e., 3δε v and −3δε v , δε dev and −δε dev . This does not lead to any difference in the elastic response, as long as we stay in the limit of elastic perturbations.
We test the rearrangements argument in Fig. 8, by plotting the calculated bulk modulus B and octahedral shear modulus G oct against the amplitude of the applied isotropic 3δε v and deviatoric δε dev strain, respectively, for states at ε dev = 0.0065 and 0.31 (nearly isotropic and steady state configurations, respectively) of the main deviatoric experiment. Both B and G oct stay practically constant for small amplitudes, and we can assume the regime to be linear elastic [10]. At 3δε v 10 −4 , the first change in the number of contacts happens (Fig. 7a, b), and B starts to increase non-linearly. Similarly, when ε dev 10 −4 , the first change in the number of contacts happens (Fig. 7c, d), and G oct starts to decay. It is interesting to notice that for both B and G oct the elastic regime shrinks when the main deviatoric strain ε dev increases (Fig. 8) and, also, when the volume fraction reduces, going toward the jamming volume fraction (data not shown). A fabric modulus may be plotted as δ F dev /δε dev that, due to the finite size of the system, would be identically zero, until the first change of structure (fabric) occurs (see Fig. 6).
We further check the elasticity of the probe by reversing the incremental strain. We plot the stress responses to volumetric/deviatoric strain in Fig. 9 and compare loading and unloading probes for different volume fractions (ν = 0.706 and 0.812) and amplitudes. Looking at Figs. 6, 7 and 9 together, three regimes seem to appear. The first one at very small strain (on the order of 10 −5 ) is characterized by no opening and closing of contacts and shows perfect reversibility of the data, i.e., elasticity, in Fig. 9a-d. This is related with the finite size of the system. The second regime in Fig. 9(e-h) shows some weakly irreversible behavior, but only for the smaller Evolution of the incremental non-dimensional pressure P * , non-dimensional shear stress σ * dev during small (a-d), medium (e-h), and large (i-l) perturbations in the loading (symbols) and then unloading (solid lines) direction. Red plus symbols represents loading and the green line represents unloading for ν = 0.812. Similarly, blue asterisk symbols represents loading and the black line represents unloading for ν = 0.706. The deformation is applied to the state corresponding to ε dev = 0.31 (steady state configuration) of the main deviatoric experiment (color figure online) volume fraction and a mixed perturbation mode, see the sample at ν = 0.706 in Fig. 9f; we associate this behavior to minor contact changes, as visible in Figs. 6 and 7, but no large scale rearrangements occur. Finally, in the third regime, for perturbations two orders of magnitude higher (on the order of 10 −3 ), a residual strain after reversal shows up for both volume fractions and all types of perturbations, see Fig. 9i-l, proving also that plasticity is much more pronounced in the deviatoric modes than in isotropic ones. We claim that small drops are related to local (weak, almost reversible) re-structuring, while in the last case, the whole system (or a big portion of it) is involved in the collapse of the structure, with a more pronounced effect for samples close to the jamming volume fraction [35,49].
For granular materials, the strain cannot be split in elastic and plastic contributions by "trivially" referring to the residual deformation like in classical solids: as soon as we are out of the elastic range, rearrangements happen during loading and (even though less probably) during unloading, and most likely no original position is recovered for any particle. Finally, we note that the results shown here are valid for finite-size systems; for much larger (real) samples of much smaller particles, we expect the first elastic regime to reduce to much smaller strains. The transition between the second and third regime is an issue for further research [68].

Evolution of the moduli
Using five packings at different ν i , we next determine which variables affect the incremental response of the aggregates at different deviatoric strains along the main path. In order to understand the role of the microstructure, i.e., the fabric tensor F, the volumetric and deviatoric components, F v and F dev , are considered. We postulate that the incremental response of the granular material can be uniquely predicted, once its fabric state (along with the stress state) is known, irrespective of the path that the system experienced to reach that state. In this sense the fabric tensor can be referred to as a state variable.

Bulk modulus B
In Fig. 10a, we plot the incremental non-dimensional pressure δ P * against the amplitude of the applied isotropic perturbation 3δε v for one volume fraction, ν = 0.706, and various initial anisotropic configurations. The slope of each line is the bulk modulus of that state. It practically remains unchanged for different states and suggests that B is constant for a given volume fraction.
In Fig. 10b, we plot the variation of the bulk modulus B with the isotropic fabric F v for packings with different volume fractions ν i . B increases systematically when the five different reference configurations are compared, and it is related to the value of F v constant at a given ν i [24,40,70]. As expected B is a purely volumetric quantity and varies with changes in the isotropic contact network. The inset in Fig. 10b shows that the bulk modulus remains almost constant with applied shear during a single deviatoric experiment [40], behaving qualitatively similar to pressure P * and isotropic fabric F v , see Figs. 1a and 3a, respectively. That is, the contact orientation anisotropy, F dev , which changes during the main deviatoric deformation path (see Fig. 3b) does not affect it. In agreement with observations on the volumetric fabric in Sect. 3.2, also B shows a slight increase/decrease in the first part of the deviatoric path, more pronounced for loose samples, as clearly seen in Fig. 2b. The trend of B slightly deviates from F v in the low strain regime, while the dependence is well captured in the steady state, after large strain. The relation between bulk modulus and fabric was given in Ref. [24] as: where p 0 , γ p and the jamming volume fraction ν c are fit parameters presented in Table 1. 4 g 3 ≈ 1.22 is dependent on the particle size distribution as presented in Refs. [24,30,39], see Sect. 2. For a given volume fraction, the above relation only requires the knowledge of the isotropic fabric where the empirical relations for C * (ν) and φ r (ν) are taken from Refs. [30,39], see Sect. 2. The numerical data show good agreement with the theoretical prediction presented in [24] and reported in Fig. 10b. The minimum F v is obtained at the jamming volume fraction, with ν c = 0.658, C * = 6, and φ r = φ c = 0.13, leading to F min v ≈ 4.2. At the jamming transition, we can extrapolate a finite value of the bulk modulus B min = 0.22, while it suddenly drops to zero below ν c [14,31,55,[62][63][64][65]85]. The discontinuity of B is related to the discontinuity in F v , that jumps from zero to a finite value at ν c due to equilibrium requirements.

Anisotropy moduli A 1 and A 2
In Fig. 11a, we plot the non-dimensional pressure increment δ P * against the strain amplitude, when the material is subjected to small deviatoric perturbations δε dev , to measure the first anisotropy modulus A 1 as defined in Eq. (12), in given anisotropic configurations, as in Fig. 10a. Since the material is in an anisotropic state, an increment in deviatoric strain leads to a change in volumetric stress, along with shear stress. The slope of the curves, A 1 , increases with the previous shear strain the system has experienced, from small values in the initial isotropic configuration, to an asymptotic limit.
We are interested in the ratio A 1 /B. In this ratio, the dependence of isotropic fabric F v cancels out, all that remains is a pure dependence on F dev . In Fig. 11b, we plot the variation of A 1 /B, with F dev for packings with different volume fractions ν i as shown in the inset. Besides the fluctuations, the data collapse on a unique curve irrespective of volume fraction and pressure, that is, once a state has been achieved, a measurement of the overall anisotropy modulus is associated with a unique F dev . An increasing trend of A 1 /B with the deviatoric fabric shows up. As the deviatoric fabric decreases with volume fraction (see Fig. 3b), this leads to lower values of the scaled anisotropy modulus for denser systems. In conclusion, we have a linear relation between the first anisotropy modulus A 1 and fabric: where B is the bulk modulus, F dev is the deviatoric part of fabric, and a I ≈ 1.0 is a fit parameter presented in Table 1.
In Fig. 12a, we plot the stress response of the material δσ * dev to isotropic perturbation 3δε v , for the same anisotropic initial configurations as in Fig. 10a, to measure the second anisotropy modulus A 2 as defined in Eq. (12). Similarly to A 1 , the slope of the elastic curves, i.e., A 2 , increases with the previous shear strain the system has felt, starting form zero until an asymptotic limit is reached. In Fig. 12b, we plot the variation of A 2 /B, with F dev for different volume fractions ν i as shown in the inset. Data show a very similar trend as in Fig. 11b, and besides the fluctuations, a collapse of data is observed. 5 Again we can relate A 2 to deviatoric fabric as The equality between the two fitting constants a I ≈ a II ≈ 1 (see Table 1), establishes the symmetry of the stiffness matrix in Eq. (11). Equations (14) and (15) provide an interesting, novel way to back-calculate the deviatoric structure in a granular sample via F dev = A/B, where A and B can be inferred from wave propagation experiments, while the direct measurement of fabric is still an open issue [29,34,36,74]. Table 1 Summary of fit parameters extracted from the small perturbation results in Eqs. (13), (14), (15) and (17) Modulus Parameters Bulk modulus B p 0 = 0.0425, γ p ≈ 0.2, ν c = 0.658 First anisotropy modulus A 1 a I = 1 ± 0.01 Second anisotropy modulus A 2 a II = 1 ± 0.02 Octahedral shear modulus G oct g I = 130 ± 3, F α v ≈ 1.9

Octahedral shear modulus G oct
In Fig. 13a, we plot the shear stress response δσ * dev of the material when the initial configurations in Fig. 10a are subjected to small deviatoric perturbations δε dev . The octahedral shear modulus G oct is then measured, as defined in Eq. (12). The slope of the curves for different initial configurations slightly decreases with the deviatoric state of the system, and saturates for high deformation ε dev , when the steady state is reached. Figure 13b shows the variation of G oct against shear strain ε dev . G oct starts from a finite value in the initial configuration, related to the isotropic contact network, and slightly decreases with increasing strain, with different rates for different volume fractions. The behavior of G oct differs from that observed for the bulk modulus in the inset of Fig. 10b: the shear resistance consistently decreases with shear strain, and no transition between initial decrease/increase is observed, meaning that a factor other than F v influences the change of G oct during the deviatoric path. Similarly to what done for A 1 and A 2 , we look at the ratio of the shear modulus with the bulk modulus G oct /B plotted against the isotropic fabric F v in Fig. 13c. The ratio increases with increasing F v , with higher values in the initial state than in the steady state (data are averaged over shear strain ε dev ≤ 0.0065 to get the initial value and in the steady state to get the final one). The isotropic ratio G oct /B ini increases with F v , following the empirical relation: where G oct /B max ∼ 0.51 represents the maximum value of the ratio G oct /B for large F v (or volume fraction), F min v ∼ 4.2 (see Sect. 4.3.1) and F α v ∼ 1.9 are the volumetric fabric and the rate of growth of (G oct /B) ini at jamming transition. This is in agreement with previous studies that find an upper limit equal 0.5 for the ratio between the shear and bulk moduli [18,38,50,73]. In the limit of high F v , the granular assembly becomes highly coordinated and practically follows the affine approximation that predicts a constant value for the ratio G oct /B [78]. Here, a qualitatively similar behavior is observed for the values in the steady state, approaching a saturation ratio lower than the isotropic one.
Next, in Fig. 13d, we subtract the initial value G oct /B ini from G oct /B and assume that F v does not change during the deviatoric deformation. Interestingly, we find that in this case the deviatoric microstructure alone is not able to capture the variation of the modulus along the shear path, but both stress σ σ σ and fabric F seem to influence the incremental shear response, in agreement with findings in [87]. We relate the decrease in G oct to the deviatoric components of stress and fabric via: where σ * dev is the non-dimensional shear stress, F dev is the deviatoric fabric, and g I ≈ 130 is a fit parameter reported in Table 1. Two contributions of the fabric to the shear stiffness can be recognized-isotropic and deviatoric. The overall contribution is multiplicative proportional to B, related to the isotropic contact network and changing very little with deviatoric strain. In the bracket, the first term gives the resistance of the material in the initial isotropic configuration, whereas the second part only depends on the deviatoric (state) variables and characterizes the evolution of the shear modulus with deviatoric strain. Given the initial isotropic configuration, the corresponding G oct is known [16,50,78]; on the other hand, the anisotropic network of such a configuration uniquely defines the reduction in the shear stiffness. The product of deviatoric stress and fabric σ * dev F dev as proposed in [77,87] is able to capture the evolution of the ratio of the elastic moduli along the whole undrained (c) (d) Fig. 13 a Evolution of the incremental change in shear stress δσ * dev versus strain amplitude during purely deviatoric perturbations δε dev for different states, with volume fraction ν = 0.706, along the main path as shown in the inset. b Evolution of octahedral shear modulus G oct along the main deviatoric path ε dev for five different volume fractions as shown in the inset. The corresponding lines passing through the data represent Eq. (17). c Evolution of ratio of octahedral shear modulus and bulk modulus, i.e., G oct /B with isotropic fabric F v , together with the averaged values at the initial (near isotropic state averaged over shear strain ε dev ≤ 0.0065) and the steady state (averaged dataset in the steady state), as given in the legend. Note that the difference between initial and steady state increases with denser systems. The solid orange line passing through the isotropic dataset represents Eq. (16). d Evolution of the ratio of octahedral shear modulus and bulk modulus when its initial value, i.e., G oct /B − G oct /B ini is subtracted, plotted using Eq. (17), for five different volume fractions as shown in the inset (color figure online) path, as seen in Fig. 13d. 6 No other relation with volumetric quantities needs to be considered, as the evolution of σ * dev F dev depends on the volume fraction of the sample ν i . Note that when G oct /B is plotted against F v in Fig. 13c, a deviation from the fitting law is observed for each volume fraction, showing that extra correction terms might be needed for a more accurate description. This is neglected in this preliminary work. It is interesting to point out that the isotropic fabric has different effects in case of the anisotropy moduli A 1 , A 2 and G oct , as in the former two cases F v , through B, is multiplied to F dev and contributes to the growth of the moduli from zero to the asymptotic values, while in the latter case F v defines mostly the initial values of G oct via the bulk modulus, but does not affect the further decrease.

Prediction of undrained cyclic shear test
In this section, the constitutive model is presented, involving the elastic moduli measured and calibrated in Sect. 4, and the plastic response of the material under large strain. The model is then used to predict the material response under cyclic shear, involving reversal.

Calibration: constitutive model with anisotropy
We introduce here a constitutive model as proposed in Refs. [30,37,39,48,51], extended to three dimensions, that takes into account the evolution of fabric, independently of stress: In its simple, reduced form, the model involves only three moduli B, A and G oct , defined in the previous section in Eqs. (13)- (17). Due to A, the model provides a cross-coupling between the two types of stress and strain in the model, namely the isotropic stress P * and shear stress σ * dev reacting to both isotropic (ε v ) and deviatoric (ε dev ) strains. F dev evolves differently from stress, as the rate of change with deviatoric strain can be (and in many cases is) different from the respective rate for the shear stress evolution. Note that additional terms (cross-coupling of fabric with strain) might be needed for the incremental evolution of δ F dev in Eq. (18), due to the observations from Fig. 6m, k, where F dev and F v change also with ε v and ε dev , respectively. However, both cross-terms appear to be mostly activated in the highly anisotropic state, with values of the out-of-plane fabric considerably smaller than out-of-plane stress-but this has to be confirmed by other deformation paths also, i.e., we claim that some features are related to the specific deformation path proposed here. If a dependence between stiffness and fabric similar to that proposed in Eqs. (14) and (15) is assumed, previous arguments also lead to the conclusion that the out-of-plane stiffness terms developing during plane strain and neglected in Eq. (11) must be small compared with B, A and G. For the sake of simplicity, both evolution of cross-coupling fabric terms and out-of-plane stiffness are neglected in the present work, and postponed to future investigations, for the description of arbitrary deformation paths. The use of non-frictional particles is another possible reason for the simplest model to work astonishingly well-so the general model is expected to show all contributions for arbitrary deformation, in the presence of friction.
S σ = S/S I σ , with S = (1 − s dev /s max dev ) is a measure of the stress isotropy with normalized shear stress ratio s dev = σ * dev /P * , and S I σ is the initial stress isotropy at the start of a new deformation direction and/or after relaxation. 1 − S σ is the measure for the probability of plastic events. Similarly, S F = (1 − F dev /F max dev )/S I F is the fabric isotropy, and S I F is the initial fabric isotropy at the start of a new deformation direction and/or after relaxation. s max dev and F max dev represent the maximum (saturation) values of normalized shear stress ratio s dev and deviatoric fabric F dev , respectively, and β F is the rate of change in F dev at smaller strains (as shown in Fig. 3b).
It is worthwhile to point out that the definitions of S σ and S F are different from those used in Refs. [48,51], as both S σ and S F are now scaled by the initial reference value and can take values between 0 and 1. Due to S σ and S F , the incremental response of the material is purely elastic, after relaxation or at strain reversal, with the elastic moduli evolving, as given by Eqs. (13)- (17), as functions of the momentary stress and structure states. At reversal, the probability for plastic deformation drops to zero, and plastic events-as related to the approach to steady state-only occur after relatively large strain, that is, the reversal stiffness is not affected. Due to S σ and S F , the incremental response of the material in the large strain steady state (S = 0) becomes elastic (S = 1), just when the strain is reverted, or after relaxation (which is allowed before the probes). Due to the dependence of the elastic moduli on the stress/fabric state, the model involves nonlinear elasticity in its present form (without contact nonlinearities), while plasticity due to rearrangements is entirely associated with S σ . On the other hand, the equation that describes the evolution of fabric is "purely plastic," as there is no change in fabric (δ F v = 0), in the elastic regime, when no contact opening/closing and no multi-particle rearrangement happen. 7 Thus, the rate β F is associated with changes of structure with deviatoric (shear) strain amplitude (not rate); changes are becoming more and more probable when approaching the steady state.  [30,39]. The red circle data points are the DEM simulation data over which the calibration of moduli was done, while the green asterisk data points represent unloading (reversal) and re-loading. The solid line is the prediction of the DEM observations using Eq. (18) (color figure online) Now, we can predict an independent experiment, by using Eq. (18), and the relations between the three moduli B, A and G oct with microscopic state variables given by Eqs. (13)- (17) with the numerical scaling factors from Table 1 (starting from B, we express the other moduli using their ratio with B). Moreover, three other parameters s max dev , F max dev and β F are needed to fully solve the coupled Eq. (18). The dependence of s max dev , F max dev and β F on volume fraction ν is well described by the exponential decay relation proposed in Refs. [30,39], where constant values, as given in Fig. 14 are used, as the volume is conserved during the cyclic shear test, as discussed next.

Prediction: (undrained) cyclic shear test
We choose an initial isotropic configuration, with volume fraction ν = 0.711 and apply deviatoric (volume conserving) shear for one cycle: loading, unloading and final re-loading, to recover the initial box configuration. Footnote 7 continued points in Fig. 5, the incremental measured elastic response (moduli) is different between states A and B as it depends on stress and fabric, that is, the stiffness matrix in Eq. (11), varies non-linearly with ε dev . On the other hand, when the incremental strain δε dev is applied to each state (e.g., A or B), the incremental response is linearly elastic (by definition of incremental) and becomes plastic for high δε dev , when rearrangements happen and the moduli in that given state go from elastic to plastic. Figure 14 shows the evolution of pressure P * , shear stress σ * dev , shear stress ratio s dev and deviatoric fabric F dev with deviatoric shear strain ε dev for one cycle, compared with the prediction using Eq. (18). Since the initial configuration is isotropic, the shear stress σ * dev and F dev start from zero and approach saturation values (with fluctuations) at large strains. During reversal, both drop from their respective saturation value and decrease with unloading strain, crossing their zero values at different strain levels, and finally reach their steady state with negative signs. This supports the need of independent descriptions for the evolution of stress and fabric. Finally, re-loading is applied to reach the initial box configuration.
The qualitative behavior of pressure P * is similar in simulations and model, going from a finite initial value to saturation with much less pronounced variations, since the deformation path is volume conserving. It is also interesting that the final state after the complete cycle, which corresponds to the initial box configuration, is highly anisotropic (non-zero stress σ * dev and deviatoric fabric F dev ). Both, the shear stress σ * dev and deviatoric fabric F dev , as well as their responses during strain reversal are well predicted by the model. P * increases during loading ε dev by ∼ 9% and saturates at large strains. After reversal, P * drops because of opening and release of contacts and then increases again with unloading strain. Although P * is not quantitatively predicted by Eq. (18), the qualitative behavior is captured by the model, which requires a correction as proposed by Krijgsman and Luding [37]. The concept of a history-dependent jamming point, introduced by Kumar and Luding [41], is capable of capturing the behavior of P * quantitatively; however, this goes beyond the scope of this study.
Equations (18) are able to describe the volumetric/deviatoric behavior of a granular assembly, in terms of stress and fabric. Once the initial state and the deformation path are defined, the evolution of isotropic fabric can be determined (using the coordination number and the fraction of rattlers) along the deformation path. The knowledge of isotropic and deviatoric fabric and the incremental relations in Eq. (18) allow for the definitions of the moduli at each incremental step. Given also the probabilities for the plastic events (1 − S σ and 1 − S F ), the coupled system can be solved. That is, the characterization of the initial state is the information needed to fully describe the behavior of the material along a general deformation path, defined in terms of strain, since the incremental evolution equations for both stress and structure are given.
In the case of granular matter, the concept of a (homogeneous) material point in a continuum model is debated, and many studies have been devoted to the introduction of a length scale in the constitutive model, starting from the Cosserat brothers, see [11,44,53] among others. Here we limit ourselves and state that a finitesize system is always needed, in order to calibrate a continuum model. Every model interpretation works only between the upper/lower bounds of infinite system and particle scale. When a finite-size system is considered an elastic range can always be detected, such that rearrangements happen (see Sect. 4.2) with negligible(tiny) probability for very small strain, and an elasto-plastic framework could then make sense. In this work we introduce a local rate-type model in Eq. (18), and identify elasticity as the unique initial, static, configuration, from which the (incrementally irreversible) evolution of stress and structure follows. Our choice is to reduce elasticity to a "punctual range," as plastic deformations (which include irreversible opening/closing of contacts by large scale rearrangements) will dominate for large deformations. Dynamics and kinetic fluctuations, leading to relaxation, are not considered here, but also need to be taken into account, see e.g., [33].

Summary and outlook
In a triaxial box, the four elastic moduli that describe the incremental, elastic constitutive behavior of an anisotropic granular material in terms of volumetric/deviatoric components, namely the bulk modulus B, the two anisotropic moduli A 1 , A 2 and the octahedral shear modulus G oct , can be measured by applying small strain perturbations to relaxed states that previously experienced a large strain along a volume-conserving (undrained) shear path. A connection between the macroscopic elastic response and the micromechanics is established, by considering both stress and fabric tensors, σ σ σ and F, respectively. While the bulk modulus B depends on the isotropic contact network F v , the deviatoric component of the fabric tensor F dev is the fundamental state variable determining the ratios between the (cross-coupling) anisotropic and bulk moduli. When the deviatoric stress and strain are appropriately scaled (normalized), we find that the moduli reduce to three relevant ones, i.e., A = A 1 = A 2 = a F dev B, with a ≈ 1. The anisotropy moduli are related to both deviatoric and isotropic fabric, as the whole contact network determines how the system will react to a perturbation. Surprisingly, when the shear resistance G oct is considered, both the contact network and the deviatoric stress determine the incremental behavior of the assembly. When the initial response is subtracted, the residual ratio G oct /B − G oct /B ini scales with the deviatoric state of the system, by the product σ * dev F dev . For strain amplitudes larger than about 10 −4 , rearrangements in the sample take place, and the behavior deviates from elastic (reversible). The effect of increasing amplitude of isotropic/deviatoric strain perturbations on isotropic/deviatoric stress and fabric is investigated, in various cases from nearly isotropic states to steady states at various different densities. For very small strain, the initial (linear) elastic regime, visible in the stress response, is associated with strictly no change in fabric. For higher strain amplitude applied to a nearly isotropic state, plasticity comes into play, and the incremental stress-strain relation deviates from linear as soon as the contact network changes. When the system is strained up to the steady state, the elastic regime is tiny (compared to the nearly isotropic states), and both stress and fabric display first elastic and then plastic response to smaller deviatoric strains. In other words, when the system was strained up to the steady state, the elastic regime becomes smaller, and both stress and fabric display no change associated to saturation. Large volumetric strain induces substantial modifications, as the sample previously subjected only to volumeconserving deformation experiences now large volume changes. In the limit of large strain, the tangential moduli of the stress-strain and fabric-strain curves (see Fig. 5) are recovered. The relation between particle rearrangements and macro-scale plasticity is a present object of investigation, as well as the transition between local/global plastic regimes. As an important result, our study provides a new way to indirectly characterize the granular structure. Once the moduli in a given isotropic/anisotropic configuration have been measured by wave propagation experiments, they can be uniquely associated with the internal fabric. However, we do not expect the proportionality factors to remain constant for different materials. As further step, a simple constitutive model is introduced that involves anisotropy, as proposed in Refs. [48,51]. The nonlinear elastic behavior is established, and the irreversible/ plastic contribution is introduced via empirical probabilities for plastic events, that require more research and theoretical support. The dependence of the model parameters on volume fraction and polydispersity has been analyzed in previous extensive work [30,39]. Here, by using the new relations for the elastic moduli, we are able to integrate the increments at each state along a generic deformation path. Hence, we can predict the evolution of pressure, shear stress and fabric for large strain, and also at and after reversal. The method is first calibrated and then applied to a volume-conserving (undrained) shear cycle. When the prediction is compared with numerical simulations, quantitative agreement is found for the deviatoric field variables. The most notable feature of soft but different response of shear stress and fabric to reversal of strain is well captured; the pressure response amplitude is underestimated by the present model.
This study concerns a seemingly unrealistic material of spheres without friction and interacting with linear contact forces to exclude effects that are due to contact nonlinearity, friction and/or non-sphericity. This allows to unravel the peculiar interplay of stress with microstructure. However, the work should be extended to more realistic cases involving particle shape, friction, and non-linear contact behavior. We expect that friction will not completely change the observed relations between stiffness and fabric state, but possibly will add new effects to be explored in the future; the deviatoric fabric and the moduli are expected to change quantitatively when tangential forces are included. The contacts nonlinearity will introduce an extra pressure-dependence for the moduli, as already shown by many authors (see e.g., [7,16,50,78] in the case of Hertzian interactions). Speculating about the effects of shape goes beyond the scope of this study. A similar analysis is already in progress to check the influence of polydispersity on the relation between elastic stiffness and microstructure, as polydispersity affects the contact network, the structure, and the orientation of contacts [24,25,39].
Future work will focus on the extension of our small perturbation approach to elasto-plasticity, by using concepts like e.g., the Gudehus response envelope [27,54]. Other theoretical approaches involve ideas proposed by Einav [17], or by Jiang and Liu [33], for which our results can provide a microscopically-based calibration of parameters, but details are not discussed here. The information obtained for the pure elastic range can then be used to decouple the plastic contribution, associated with rearrangements, and to study the flow rule. The validation of the present analysis with experimental data is another important goal. Nevertheless, the issue of measuring fabric from laboratory experiments is far from solved, even though big advances have been made in recent years using photoelasticity and microtomography CT-scans [5,29,34,74]. A partial validation is anyway possible when measuring the residual dependence of the elastic response from variables other than stress and porosity [19], by means of acoustic measurements [36]. The behavior after more than one cycle deserves further investigation, from both simulational and theoretical points of view, to detect features like creep, liquefaction and ratcheting, analyzed in preliminary works [51] with constant elastic moduli and for many cycles [41].
Finally, a general tensor formulation that allows for highly different orientations of strain-rate, stress and fabric is an open issue, but can be inspired by the works of Thornton [77] and Zhao and Guo [87].