DEM simulation of anisotropic granular materials: elastic and inelastic behavior

In this work, Discrete Elements Method simulations are carried out to investigate the effective stiffness of an assembly of frictional, elastic spheres under anisotropic loading. Strain probes, following both forward and backward paths, are performed at several anisotropic levels and the corresponding stress is measured. For very small strain perturbations, we retrieve the linear elastic regime where the same response is measured when incremental loading and unloading are applied. Differently, for a greater magnitude of the incremental strain a different stress is measured, depending on the direction of the perturbation. In the case of unloading probes, the behavior stays elastic until non-linearity is reached.Under forward perturbations, the aggregate shows an intermediate inelastic stiffness, in which the main contribution comes from the normal contact forces. That is, when forward incremental probes are applied the behavior of anisotropic aggregates is an incremental frictionless behavior. In this regime we show that contacts roll or slide so the incremental tangential contact forces are zero.


Introduction
Granular media are complex systems widely present in civil engineering in the form of soils or granulates, in industry including chemical synthesis, food production, thermal insulation, additive manufacturing and other application consisting of granular beds. Understanding the mechanical response of granular materials is important to elucidate fundamental aspects of the behavior of these particulate systems [1,2]. To this end, numerical simulations, laboratory experiments and theoretical models have been employed. In particular, an interesting activity regards the theoretical analysis, developed in order to establish predictive models that should reproduce what seen in numerical simulations and/or laboratory experiments. There are models based upon a phenomenological approach and other based on micro-mechanics. The latter are more favorable when compared with Discrete Element Method (DEM) simulations [3] because it is possible to test not only the macroscopic response of the aggregate but also local features that characterize particle interactions. In this paper, we focus on a numerical analysis for a granular aggregate, referring to theoretical models already available in literature. In fact, it is not our goal to develop a new theory but to provide new insights of the elasto-plastic regime. We are interested to the incremental response of dense, sheared granular samples and how the response qualitatively changes as anisotropy develops in the aggregate. This is an important point, for example in the context of seismic waves that propagate in granular materials (e.g. [4,5]), geotechnical applications involving regions where deformations are small [6], the development of elasticplastic constitutive models, where elasticity needs a proper description, e.g. [7,8] or as indicator for localization [9]. Several approaches have been used to investigate elasticity numerically, e.g., dynamical unloading probes [10][11][12], response envelope [13][14][15], stiffness matrix [16,17], wave propagation [4,18,19], depending on the specific focus and final goal. In particular, procedures (and often conclusions) diverge if the interest is on "pure" elasticity at very small strains or an elasto-plastic framework.
Here we study elasticity in granular materials over a wide range of strain magnitudes, for different directions and degrees of anisotropy. We highlight that the mechanical response of anisotropic granular aggregates is different if forward or backward incremental strain are applied. Beside the well-known linear elastic regime in which the same stress is measured when both forward and backward incremental strains are applied, we recognize a second regime, associated with greater perturbations, that proceed the non-linear behavior, i.e. the stiffness depends on the strain amplitude. We identify in this regime an inelastic stiffness in which the response becomes incrementally frictionless and the major symmetry of the macroscopic stiffness is lost (e.g. [20,21]). Key parameters are the magnitude and direction of the probes applied to stressed, anisotropic states, compared with the strain under which the aggregate is initially loaded. In such framework, Froiio and Roux [17], Calvetti et al. [15], Kuhn et al. [22] use response envelopes obtained via DEM multidirectional loading probes to investigate the validity of common assumptions of elasto-plastic models for granular materials subjected to anisotropic loading paths. We, instead, operate with limited loading conditions because we have a different goal. Specifically, we focus on probes parallel and orthogonal to the initial monotonic loading in order to unravel the role of contacts elasticity, sliding and rolling in the transition of the elastic, inelastic, and plastic regime where the loss of symmetry emerges.

Numerical simulations
The DEM methods (e.g. [3,23]) are a powerful tool to study granular materials in combination with theoretical models (e.g. [24][25][26]) to predict the material response under different loading conditions (e.g. [10][11][12][27][28][29]) The code is based upon the knowledge of particles position and interaction forces. If the contact forces, acting on a particle, are known the problem is reduced to the integration of Newton's equations of motion for the translational and rotational degrees of freedom of that particle. The system considered here is a random assembly of identical, frictional, elastic spheres that interact through contacts in which gravity is neglected. A micro-mechanical analysis shows interesting features associated with the possibility that particles slide or roll.

Contact mechanics
For a given pair of particles, the interaction is represented by a non-central contact force where d i is the unit contact vector that joins the centers of contacting particles and t i is the tangential unit contact vector in the plane perpendicular to d i . The normal force F N follows the non-linear Hertz law. F T is the tangential contact force that incorporates a bilinear relationship, i.e., an elastic resistance followed by Coulomb sliding [30]. When F T ≥ F N , the tangential force in the slip direction is F T = F N [31].
The average stress ij of the aggregate, according to Cauchy [32], is given by in which the sum is extended to all N c contacts in the representative volume V. As in [33] and [28], the stress may be partitioned into its normal and tangential components and in order to capture the relative contribution of the two parts.

Preparation protocol
We employ material properties typical of glass spheres, shear modulus G = 29GPa and Poisson's ratio, = 0.2 . We use an aggregate of N = 10, 000 spheres, each with radius R = 0.1mm, randomly generated in a periodic cubic cell. Our calculations begin with a numerical protocol designed to mimic the experimental procedures used to prepare densely packed granular materials. Particles are then isotropically compressed without friction, = 0 , until a solid volume fraction slightly lower than random close-packing ≤ RCP has been reached ( RCP ≃ 0.64 for monodisperse aggregates [34]). Then, the particles are allowed to relax and reach pressure and coordination number equal to zero [11].
Furthermore, an isotropic compression is applied with friction coefficient = 0 to reach the target value of mean stress p 0 = 200kPa, followed by a new relaxation stage under which the final friction coefficient is set to be = 0.5 . In this reference isotropic configuration, the solid volume fraction reaches = 0.64 and the coordination number (the average number of contacts per particle) k 0 = 5.95 , while the volumetric strain associated with the pressure p 0 is 0 = 1 × 10 −3 . This measure of the strain is based upon a theoretical prediction proposed by Jenkins et al. [35] in which, in a succession of isotropic states, pressure and volume strain are related by: whose incremental formulation is Given the pressure p 0 and the bulk modulus of the aggregate, i.e. p 0 ∕ 0 , Eq. 6 permits to determine 0 .

Axial-symmetric compression
After the isotropic compression, an axial-symmetric deformation is applied along the direction 1 . We take the friction coefficient = 0.5 , although glass beads are characterized by a smaller value, = 0.3 , because we obtain a smoother response with a an almost identical behavior. The test is carried out at constant mean stress, p 0 = 200 kPa, by means of a servo-mechanism [28]. To ensure quasi-static conditions, the compression is performed with a sequence of small strain steps, 11 ≃ −10 −5 (compression < 0 in our convention), and relaxation steps in which particles are allowed to dissipate kinetic energy and to reach intermediate equilibrium states. At each time step, along the compression path, we measure the deviatoric stress and the normal and tangential parts q N and q T derived by means of Eqs. 3 and 4 , respectively [33]. In this work, we limit our analysis to a relative small range of deformation, ∕ 0 < 0.4 , in which = ( 22 + 33 )∕4 − 11 ∕2 deviatoric strain applied. However, in this regime, contacts already experience elastic deformation, sliding and deletion [26]. In Effective moduli A 11 and A 31 versus strain amplitude Fig. 1 we plot the normalized deviatoric stress q∕p 0 against the normalized deviatoric strain ∕ 0 . In the same figure, the partition of the deviatoric stress into the normal and tangential contributions q N ∕p 0 and q T ∕p 0 is shown as well. After an initial similar response for q N ∕p 0 and q T ∕p 0 , it is clear that the deviatoric stress is almost identifiable with its normal component, q N ∕p 0 . In this rather narrow regime of deformation the change in the coordination number (see Fig. 2), is negligible as well as in the fabric and volume. From the contact point of view, a clear picture emerges when the aggregate is axially deformed. The initial state is isotropic in terms of contact vector orientations and contact loading. When the aggregate is sheared, within the range of deformation of our interest, the contact network is still approximately isotropic, while contact loadings show an alignment with the applied deformation along 1 . Anisotropy in loading occurs and 1 becomes a preferential direction for the contact forces. Moreover, as seen in Fig. 1, the tangential component of the stress, q T ∕p 0 , reaches a plateau, e.g. [28,33], which implies that, after about ∕ 0 ≃ 0.2 , the increment in tangential forces is zero, as we will show in the present contribution.

Incremental response
At four stages along the loading path, as indicated in Fig. 1, we calculate the incremental response of the aggregate.
In order to measure the components of the stiffness matrix A ijkm , an infinitesimal strain km is applied to the aggregate and the resulting change in stress is measured after sufficient relaxation [10]:

Axial-symmetric probes
We first focus on the incremental response with perturbations that maintain the symmetry generated by the axial symmetric loading. That is, 22 = 33 = 0 , 11 ≠ 0 , so, adopting Voigt's notation, and with A 31 = A 21 for symmetry. We also distinguish between forward loading and unloading. In the former, 11 is negative as it is in the axial-symmetric loading ( 11 < 0 ), while in the latter 11 is positive ( 11 > 0). In Fig. 3 we plot the evolution with the strain of the moduli A 31 and A 11 , for the isotropic state and the anisotropic state associated with ∕ 0 = 0.24.
In the first range of perturbation, 11 ≃ 10 −6 , there is no difference between forward loading and unloading, irrespective of the state of the material, isotropic or anisotropic. That is, if the aggregate is incrementally strained with an extremely small perturbation all contacts behave elastically [16]. Instead, for slightly bigger probes, 11 ≃ 10 −5 , we see a second plateau and a difference in the response between forward and unloading probes. While unloading probes seem weakly related with anisotropy, in the case of forward loading a pronounced dependency on the stress state appears, with A 11 ( A 31 ) decreasing (increasing) with anisotropy. Details of this second regime induced by forward loading are shown in the inset.It is noteworthy to mention that, during all the applied increments, we do not see any significant change in the coordination number as we will show in details later in the section Micromechanics. An interesting feature appears in the second regime when we plot the moduli A 11 and A 31 partitioned in a normal and tangential part. A (N) 11 and A (N) 31 are the moduli inferred from the stress containing contributions from the normal component of the contact forces while A (T) 11 and A (T) 31 are related to the tangential component of the contact forces. In Fig. 4a, b we show details of the results: while in the unloading cases both moduli A 11 and A 31 have contributions from a normal and tangential part of the stress, the response to an incremental forward loading is mainly characterized by the normal contribution with a clear evidence for the anisotropic states in which q T is constant (see Fig. 1). That is, the response is incrementally frictionless. At ∕ 0 ≃ 0.2, where q T ∕p 0 has reached its plateau, both A (T) 11 and A (T) 31 are almost zero. As the slope in q N ∕p 0 changes with ∕ 0 so the moduli A 11 and A 31 vary. The transition between the two plateau, the first at very small deformation identified as the elastic regime and the second associated with the plastic regime, resembles what predicts by Rudnicki and Rice [9] in their yield-vertex constitutive model. That is, the transition represents a combination of an elastic response associated with contacts that still experience an elastic resistance and a plastic response characterized by zero incremental tangential resistance. The second plateau, also emphasized in Fig. 3, indicates instead that the response is plastic with no local, incremental tangential resistance. Something similar has been also pointed out by Tamagnini et al. [37], in a continuum model, through a bounding surface to be distinguished from the classical yield surface.

Probes with no symmetry
We now look at the incremental response of the aggregate when forward loading and unloading probes are applied along 3 , 33 ≠ 0 with 11 = 22 = 0 . We can determine: and Because of the symmetry associated with the axial loading along 1 , probes along 2 and 3 produces the same response, so and In Fig. 5 we plot the results for A 13 and A 33 . Again we note a first elastic regime, at 33 ≃ 10 −6 , in which there is no difference between incremental loading and unloading.
We recall that the axial-loading (see Fig. 1) is carried out by applying a strain 11 < 0 , with 22 = 33 > 0 to maintain a (11) constant confining pressure; so unloading probe now means 33 < 0 while forward incremental loading means 33 > 0 . These incremental conditions do not reproduce the symmetry induced by the uniaxial loading. In the Fig. 4b we note again a second plateau associated with the inelastic regime for A 13 . When we look closely at A 33 , we observe a clear elastic regime, followed by a rather narrow second plateau. In the insets, we show the dependence of the forward loading probing on anisotropy.
Furthermore, we look at the normal and tangential contributions of A 13 and A 33 versus strain. In Fig. 4a, b, we observe that the behavior of the partitioned moduli differs between A 11 and A 31 . The contribution associated with the tangential part of the contact force does not vanish because we are not deforming particles along the same path of the axial-symmetric loading. That is, some contacts that were sliding under the axial loading are now behaving elastically because the probe does not maintain the same symmetry.
Finally, Fig. 6b shows that the tangential contribution of A 33 in the inelastic regime is again negligible, suggesting an incrementally frictionless behavior.
In all cases examined, increments larger than 10 −4 produces a non-linear regime, more evident in case of anisotropy, in which the response depends on the amplitude of incremental strain ( [19,36,37]). Moreover, in the regime of deformation where the inelastic stiffness is defined, when anisotropy develops, the aggregate exhibits the loss of the major symmetry in the macroscopic stiffness, A 13 ≠ A 31 . This is a crucial condition to determine localization in a granular material [9,38,39].
Our results are in line with previous findings in [15,17] and [22]. Specifically, the authors in [15] have shown that only isotropic samples conform to the hypotheses beyond elasto-plastic constitutive models, while major deviations are observed as soon as anisotropic stress history is considered. Figs. 3, 4, 5 in the present work (a simpler context in which probes are along two directions only) provide a similar message, with the presence of an intermediate regime neither totally elastic or plastic in case of anisotropic samples. In such inelastic regime, stress and strain increments may loose alignement under uniaxial probes along y 1 and y 3 , imposed over the initial triaxial stress state. In [17] DEM simulations show the outmost importance of the rotation with respect to the axis of pre-loading on the incremental response of a granular material, in terms of the definition of a flow direction and non-associated character of the flow rule, i.e. loss of major symmetry as in our paper. However, we found the loss of major symmetry already happens in what we call an "inelastic regime", even before the friction is largely mobilized in the plastic regime. Finally, in their extensive and accurate work, Kuhn et al. [22] show that five over six principles of conventional elasto-plasticity fail when tested against preloaded anisotropic granular materials. The lost principles include direction and magnitude of the strain increments, the yield criterion as well as the separation of strain increments into elastic and plastic. Our findings in the inelastic regime conform to the analysis in [22], while it still supports the idea of a fully elastic regime for very small strain increments. All these works suggest the need of more complex constitutive models for a comprehensive description of granular materials (e.g. multi-mechanism plasticity or tangential plasticity).

Micromechanical characterization of elastic and inelastic regimes
In order to characterize the difference between the material response during probing in the inelastic regime, we investigate the micro-structure of the aggregate, with special focus on the number of contacts, mobilised friction at the contacts. In Fig. 7 we compare the contact distribution functions of the initial relaxed configuration with those after incremental inelastic forward loading and unloading, in the isotropic and anisotropic ( ∕ 0 = 0.24 ) packings. As an example we show here the response under axial-symmetric probe, 11 , as in Sect. 3.1. Interestingly, we find that the contact distribution collapses on the initial configuration, irrespective of the probing direction. That is, the contact network is not affected by either unloading or forward inelastic probes. This implies that the coordination number, i.e. the mean of the functions, as well as the fluctuation in the number of contacts coincide between the packings. Moreover, when comparing the two figures, the influence of anisotropy appears to be negligible. Furthermore, in Fig. 8, we look at the mobilized friction during incremental forward loading and unloading. We define the relative frequency R f as the number of contacts with a given ratio between tangential and normal forces, F T ∕ F N , normalized by the total number of contacts, where F T ∕ F N = 1 means sliding. The figure highlights that, despite the contact network stays unchanged under probing (see Fig. 7), slippage occurs within the contact area and contacts reach the onset of sliding. An higher percentage of contacts, about 30% , sit in proximity of sliding, in the case of incremental inelastic forward loading.

Incremental sliding and rotations
As forward step, we study the kinematics at the microscale, in terms of contact sliding and particle rotations. We focus on the axial-symmetric probes as in Sect. 3.1, and propose a micromechanical interpretation to explain why the tangential part of the moduli becomes zero under incremental forward loading, see Fig. 4a, b. We show that the tangential contribution of A 11 and A 31 disappears because particles slide and/ or roll so the incremental tangential displacement become zero. The analysis is motivated by a theoretical framework able to describe the elasticity of granular materials based upon micromechanics and the role of fluctuations, as given in [21].
Let us examine in details the incremental forward loading at ∕ 0 = 0.24 . Following [21], we can express the incremental relative contact displacement between a typical pair A and B as where c i and i are, respectively, the increment in translation and the increment in rotation of the particle. The tangential component of the contact force is in which K T is the tangential contact stiffness is the overlap between contacting particles. With Eq. (15), the incremental tangential component of the contact force becomes in which s (BA) i is the incremental tangential displacement associated with the translation of the centers. Because of equilibrium, La Ragione and Jenkins [24] show that fluctuations in spin rather than in translation play a major role. Consequently, we express s (BA) i in terms of the incremental Fig. 9 The kinematics of a typical pair in which, given the proper sign of the rotations and s , it is possible to realize a relative, tangential, contact displacement equal to zero: a projection on the plane y 1 − y 3 ; b projection on the plane y 1 − y 2 . When the incremental tangential displacement associated with the translation of the centers is equal to the corresponding contribution associated with particle rotations rolling occurs; that is: Therefore, the condition for zero incremental tangential force, F T , occurs with rolling, Eq. (20), or sliding, We want to test our findings by comparing the amount of rolling and sliding particles during incremental axial-symmetric forward loading and unloading. We recall that y 1 is the axis of anisotropy with the unit contact vector defined in terms of the polar angle, , and so d = (cos , sin cos , sin sin ). As said earlier, represents a fluctuation being the average rotation over all particles zero. We measure these fluctuations by making a partition in for all . In Fig. 9 we sketch, with the proper signs, the interaction of a typical pair.
We take M p pairs of contacting particles whose contact vectors, d i , are within a strip of width and 0 ≤ ≤ ∕2 . With incremental forward loading 11 < 0 or unloading 11 > 0 applied, we measure, for all pairs, the rotations and the corresponding average in the strip, . That is, We obtain 1 approximately zero while 2 = − 3 . The numerical results, with 11 = ±4 × 10 −5 , i.e. in the inelastic regime, for both forward loading and unloading are shown in Fig. 10. The figure shows clear differences in the microscale kinematic between the two cases. While the average fluctuation 2 is comparable for high , it becomes higher for forward loading than unloading when < 50 • .
The results can be extended, by symmetry, to the other octants, > ∕2 . We take the average of the incremental tangential displacement in each band, so Eq. (19) can be written k , we define the average over the strip centered in where we have employed the average over We could have considered an average that accounts for the fraction of contacts in each strip but the results will not differ significantly.
The difference of the averages in each strip is The corresponding numerical results, for forward probing 11 = 4 × 10 −5 , are reported in Table 1. In the strips between 0 ≤ ≤ 30 o the difference in the average, a L , is approximately zero which implies that in that range contacts experience rolling rather than sliding. For 30 o ≤ ≤ 80 o particles mostly slide. In Table 2 we report the same parameters as in Table 1 but in case of unloading. Here the parameter a U is always different from zero, implying that a contribution associated with F T is always present. This is confirmed by A T 11 and A T 31 different from zero in Fig. 4a, b. The results in Table 1 are shown in a different fashion in Fig. 11, where we plot, the relative frequency R f of sliding contacts in each strip, i.e., the number of contacts in the

Conclusions
We have analyzed the behavior of a sheared granular material via DEM numerical simulations. Anisotropy is developed in the sheared granular assembly, in a regime of deformation in which the contact network does not change. In particular, we have focused on the incremental response of the aggregate when probes, different in both direction and amplitude, are applied along the shear path.
The material behaviour depends on the smallness of the applied probes in a non-trivial way. If the amplitude of the probe is extremely small then A ijkm is essentially an elastic tensor, irrespective of forward or unloading conditions. If the probes are bigger, then a difference between unloading and forward probes occurs. Consequently the response of the aggregate, when incrementally forward strains are applied, transits from "perfectly" elastic to non-linear strain regime, through an intermediate inelastic state in which the macroscopic stiffness components, independent of strain amplitudes, are different from the elastic moduli obtained via unloading probing. Through a micro-mechanics analysis we show that the overall incremental response of the anisotropic aggregate, under particular incremental perturbations, is independent of the tangential forces between particles and mechanisms like rolling or slide occur.