Hydrostatic and Shear Behavior of Frictionless Granular Assemblies under Different Deformation Conditions †

Stressand structure-anisotropy (bulk) responses to various deformation modes are studied for dense packings of linearly elastic, frictionless, polydisperse spheres in the (periodic) triaxial box element test configuration. The major goal is to formulate a guideline for the procedure of how to calibrate a theoretical model with discrete particle simulations of selected element tests and then to predict another element test with the calibrated model (parameters). Only the simplest possible particulate model material is chosen as the basic reference example for all future studies that aim at the quantitative modeling of more realistic frictional, cohesive powders. Seemingly unrealistic materials are used to exclude effects that are due to contact non-linearity, friction, and/or non-sphericity. This allows us to unravel the peculiar interplay of stress, strain, and microstructure, i.e. fabric. Different elementary modes of deformation are isotropic, deviatoric (volume-conserving), and their superposition, e.g. a uniaxial compression test. Other ring-shear or stress-controlled (e.g. isobaric) element tests are referred to, but are not studied here. The deformation modes used in this study are especially suited for the biand triaxial box element test set-up and provide the foundations for understanding and predicting powder flow in many other experimental devices. The qualitative phenomenology presented here is expected to be valid, even clearer and magnified, in the presence of non-linear contact models, friction, non-spherical particles and, possibly, even for strong attractive/ adhesive forces. The scalar (volumetric, isotropic) bulk properties, the coordination number and the hydrostatic pressure scale qualitatively differently with isotropic strain. Otherwise, they behave in a very similar fashion irrespective of the deformation path applied. The deviatoric stress response (i.e. stressanisotropy), besides its proportionality to the deviatoric strain, is cross-coupled to the isotropic mode of deformation via the structural anisotropy; likewise, the evolution of pressure is coupled via the structural anisotropy to the deviatoric strain, leading to dilatancy/compactancy. Isotropic/uniaxial over-compression or pure shear respectively slightly increase or reduce the jamming volume fraction below which the packing loses mechanical stability. This observation suggests a necessary generalization of the concept of the jamming volume fraction from a single value to a “wide range” of values as a consequence of the deformation history of the granular material, as “stored/memorized” in the structural anisotropy. The constitutive model with incremental evolution equations for stress and structural anisotropy takes this into account. Its material parameters are extracted from discrete element method (DEM) simulations of isotropic and deviatoric (pure shear) modes as volume fraction dependent quantities. Based on this calibration, the theory is able to predict qualitatively (and to some extent also quantitatively) both the stress and fabric evolution in another test, namely the uniaxial, mixed mode during compression. This work is in the spirit of the PARDEM project funded by the European Union.


Introduction and Background
Dense granular materials are generally complex systems which show unique mechanical properties different from classic fluids or solids.Interesting phenomena such as dilatancy, shear-band formation, history dependence, jamming and yield stress -among others -have attracted significant scientific interest over the past decade.The bulk behavior of these materials depends on the behavior of their constituents (particles) interacting through contact forces.To gain an understanding of the deformation behavior, various laboratory element tests can be performed (GdR-MiDi, 2004;Samimi et al., 2005;Schwedes, 2003).Element tests are (ideally homogeneous) macroscopic tests in which the experimentalist can control the stress and/or strain path.Different element tests on packings of bulk solids have been realized in the biaxial box (see Morgeneyer and Schwedes, 2003, and references therein) while other deformation modes, namely uniaxial and volume-conserving shear, have been reported in Philippe et al., 2011 andSaadatfar et al., 2012.While such macroscopic experiments are pivotal to the development of constitutive relations, they provide little information on the microscopic origin of the bulk flow behavior of these complex packings.The complexity of a granular assembly becomes evident when it is compressed isotropically.In this case, the only macroscopic control parameters are volume fraction and pressure (Göncü et al., 2010).At the microscopic level for isotropic samples, internal variables must be added to classify the microstructure (contact network), namely the coordination number (i.e. the average number of contacts per particle) and the fraction of rattlers (i.e.fraction of particles that do not contribute to the mechanical stability of the packing), see Göncü et al., 2010;Magnanimo et al., 2008.However, when the same material sample is subjected to shear deformation, not only does shear stress build up, but also the anisotropy of the contact network develops, as it relates to the creation and destruction of contacts and force chains (see Alonso-Marroquín et al., 2005;La Ragione and Magnanimo, 2012;Luding and Perdahcioğlu, 2011;Radjai et al., 1999;Walsh and Tordesillas, 2004;Yimsiri and Soga, 2010, among others).In this sense, anisotropy repre-sents a history parameter for the granular assembly.For anisotropic samples, scalar quantities are not sufficient to fully represent the internal contact structure, but an extra tensorial quantity has to be introduced, namely the fabric tensor (Goddard, 1998;Oda, 1972).To gain more insight into the microstructure of granular materials, numerical studies and simulations on various deformation experiments can be performed, see Thornton, 2010;Thornton andZhang, 2006, 2010;Yimsiri and Soga, 2010, and references therein.In an attempt to classify dif ferent deformation modes, Luding and Perdahcioğlu, 2011, listed four different deformation modes: (0) isotropic (directionindependent), (1) uniaxial, (2) deviatoric (volumeconser ving) and (3) bi-/triaxial deformations.The first two are purely strain-controlled, while the last (3) is mixed strain-and stress-controlled either with constant side stress (Luding and Perdahcioğlu, 2011) or constant pressure (Magnanimo and Luding, 2011).The isotropic and deviatoric modes 0 and 2 are pure modes which both take especially simple forms.The uniaxial deformation test derives from the superposition of an isotropic and a deviatoric test, and represents the simplest element test experiment (oedometer, uniaxial test or λ −meter as reported in Kwade et al., 1994) that activates both isotropic and shear deformation.Even though biaxial tests are more complex to realize and involve mixed stressand strain-control instead of completely prescribed strain (Morgeneyer and Schwedes, 2003;Zetzener and Schwedes, 1998), they are assumed to better represent the deformation under realistic boundary conditions -namely the material can expand and form shear bands.In this study, various deformation paths for assemblies of polydisperse packings of linearly elastic, nonfrictional cohesionless particles are modeled using the DEM simulation approach (Cundall and Strack, 1979).One goal is to study the evolution of pressure (isotropic stress) and deviatoric stress as functions of isotropic and deviatoric strain.Internal quantities such as the coordination number, the fraction of rattlers, and the fabric tensor are reported for improved microscopic understanding of the deformations.Furthermore, the extensive set of DEM simulations is used to calibrate the anisotropic constitutive model, as proposed by Luding and Perdahcioğlu, 2011;Magnanimo and Luding, 2011.After calibration through isotropic (Göncü et al., 2010) and volume-conserving pure shear simulations, the derived relations between the bulk material parameters and volume fraction are used to predict uniaxial deformations, with the goal of improving the understanding of the macroscopic behavior of bulk par ticle systems and of guiding further developments of new theoretical models that properly and predictively describe it.The focus on the seemingly unrealistic materials allows us to exclude effects that are due to friction, contact non-linearity and/or non-sphericity, with the goal of unraveling the interplay of microstructure (fabric), stress and strain.This is the basis for the present research -beyond the scope of this paper -that aims at the quantitative modeling of these phenomena and effects for realistic frictional, cohesive powders.The deformation modes used in this study are especially suited for the biaxial box experimental element test set-up and provide the fundamental basis for the prediction of many other experimental devices.The qualitative phenomenology presented here is expected to be valid, even clearer and magnified, in the presence of friction, non-spherical particles, and for strong attractive forces.This paper is organized as follows: The simulation method and parameters used are presented in section 2 while the preparation and test procedures are introduced in section 3. Generalized averaging definitions for scalar and tensorial quantities are given in section 4 and the evolution of microscopic quantities is discussed in section 5.In section 6, the macroscopic quantities (isotropic and deviatoric) and their evolution are studied as functions of volume fraction and deviatoric (pure shear) strain for the different deformation modes.These results are then used to obtain/calibrate the macroscopic model parameters.Section 7 is devoted to theory, where we relate the evolution of the stress and structural anisotropy to strain, as proposed by Luding and Perdahcioğlu, 2011;Magnanimo and Luding, 2011, and confirm the predictive quality of the calibrated model.

Simulation Method
The Discrete Element Method (DEM) has been used extensively to perform simulations in bi-and triaxial geometries (Durán et al., 2010;Kruyt et al., 2010;Luding, 2005;Sun and Sundaresan, 2011;Yimsiri and Soga, 2010) involving advanced contact models for fine powders (Luding, 2008c) or more general deformation modes, (see Alonso-Marroquín et al., 2005;Thornton, 2010;Thornton and Zhang, 2010 and references therein).However, since we restrict ourselves to the simplest deformation modes and the simplest contact model, and since DEM is other wise a standard method, only the contact model parameters and a few relevant timescales are briefly discussed -as well as the basic system parameters.

Force model
For the sake of simplicity, the linear visco-elastic contact model for the normal component of the force has been used in this work and friction is set to zero (and hence neither tangential forces nor rotations are present).The simplest normal contact force model, which takes into account excluded volume and dissipation, involves a linear repulsive and a linear dissipative force, given as (1) where is the spring stiffness, is the contact viscosity parameter and or are the overlap or the relative velocity in the normal direction .An artificial viscous background dissipation force, , proportional to the velocity of particle is added, resembling the damping due to a background medium, as e.g. a fluid.The background dissipation only leads to shortened relaxation times, reduced dynamical effects and consequently lower computational costs without a significant effect on the underlying physics of the process -as long as quasi-static situations are considered.
The results presented in this study can be seen as a "lower-bound" reference case for more realistic material models, see e.g., Luding, 2008c, and references therein.The interesting, complex behavior and nonlinearities of our most simple granular material can not be due to the contact model, but is related to the collective bulk behavior of many particles, as will be shown below.Note that the units are artificial; Luding, 2008c provides an explanation of how they can be consistently rescaled to quantitatively match the values obtained from experiments (thanks to dimensional analysis and the simplicity of the contact model used).

Simulation parameters and timescales
A typical response time is the collision duration .For for a pair of particles with masses a n d , i t i s , w h e r e is the reduced mass.The coefficient of restitution for the same pair of particles is expressed as and quantifies dissipation.The contact duration and the restitution coefficient are dependent on the particle sizes and since our distribution is polydisperse, the fastest response timescale corresponding to the interaction between the smallest particle pair in the overall ensemble is 0.228 [ ] and 0.804.For two average particles, one has 0.643 [ ] and 0.926.Thus, the dissipation timescale for contacts between two average-sized particles, = 8.37 [ ], is considerably larger than and the background damping timescale = 83.7 [ ] is much larger again, so that the particle-and contactrelated timescales are well separated.The timescale set by the maximal strain-rate (defined below) of one of our typical simulations, is = 7.2 10 -3 [s] and thus is much larger than the other timescales in the system.As usual in DEM, the integration time step was chosen to be about 50 times smaller than the shortest timescale, (Luding, 2008c).Our numerical "experiments" are performed in a three-dimensional triaxial box with periodic boundaries on all sides.One advantage of this configuration is the possibility of realizing different deformation modes with a single experimental set-up and a direct control of stress and/or strain (Durán et al., 2010;Luding and Perdahcioğlu, 2011).The systems are ideally homogeneous, which is assumed but not tested in this study.The periodic walls can be strain-controlled to follow a co-sinusoidal law such that, for example, the position of the top wall as a function of time is (2) with strain in -direction, , where is the initial box side length, is the box length at maximum strain, and is the frequency.The co-sinusoidal law allows for a smooth star t-up and ending of the motion so that shocks and iner tia ef fects are reduced.The maximum deformation is reached after half a period , and the maximum strainrate applied during the deformation at and is .Different strain-control modes are possible such as homogeneous strain-rate control for each time step (applied to all particles and the periodic walls, i.e. the system boundaries) or swelling instead of isotropic compression, as well as pressure control of the (virtual) walls.However, this is not discussed since it had no effect for the simple frictionless contact model used here and for the quasi-static deformations applied.For more realistic contact models (dating back to, e.g.Hertz, 1882;Mindlin and Deresiewicz, 1953), for friction and adhesion (Luding, 2008c and references therein) and for large strain rates, the modes of strain or stress control have to be revisited and carefully studied.

Preparation and Test Procedure
In this section, we describe first the sample preparation procedure and then the method for implementing the isotropic, uniaxial and deviatoric element test simulations.For convenience, the tensorial definitions of the different modes will be based on their respective strain-rate tensors.However, for presenting the numerical results, we will use the strain tensor as defined in section 4.2.1.

Initial isotropic preparation
Since careful, well-defined sample preparation is essential in any physical experiment to obtain reproducible results (Ezaoui and Di Benedetto, 2009), the preparation consists of three parts: (i) randomization, (ii) isotropic compression, and (iii) relaxation.All are equally important to achieve the initial configurations for the following analysis.(i) The initial configuration is such that spherical particles are randomly generated in a 3D box with a low volume fraction and rather large random velocities such that they have sufficient space and time to exchange places and to randomize themselves.(ii) This granular gas is isotropically compressed in order to approach a direction-independent configuration to a target volume fraction 0.640, slightly below the jamming volume fraction 0.665, i.e. the transition point from fluid-like behavior to solid-like behavior (van Hecke, 2009;Majmudar et al., 2007;Makse et al., 2000;OʻHern et al., 2002).(iii) This is followed by a relaxation period at a constant volume fraction to allow the particles to fully dissipate their energy and to achieve a static configuration.Isotropic compression/decompression (negative/ positive strain rate in our convention) can now be used to further prepare the system, with subsequent relaxation, so that we have a series of different initial isotropic configurations at volume fractions , achieved during loading and unloading, as displayed in Fig. 1.Furthermore, the isotropic compression can be considered as an element test itself (Göncü et al., 2010).It is realized by a simultaneous inward movement of all the periodic boundaries of the system, with strain-rate tensor where is the rate amplitude applied to the walls until the target volume fraction is achieved.A general schematic representation of the procedure for implementing the isotropic and other deformation tests is shown in Fig. 2. B represents relaxation of the system and C represents the subsequent isotropic compression up to 0.820 and then decompression.Cyan dots represent some of the initial configurations, at different , during the loading cycle, and blue stars during the unloading cycle, both of which can be chosen for further study.
The compressed and relaxed configurations can now be used for other non-volume-conserving and/ or stress-controlled modes (e.g.biaxial, triaxial and isobaric).One only has to use them as initial configurations and then decide which deformation mode to use, as shown in the figure under "other deformations".The corresponding schematic plots of deviatoric strain as a function of volumetric strain are shown below the respective modes.

Uniaxial
Uniaxial compression is one of the element tests that can be initiated at the end of the "preparation", after sufficient relaxation indicated by the drop in potential energy to almost zero.The uniaxial compression mode in the triaxial box is achieved by a prescribed strain path in the -direction, see Eq. ( 2), while the other boundaries and are non-mobile.During loading (compression), the volume fraction is increased as for isotropic compression from 0.640 to a maximum volume fraction of 0.820 (as shown in region C of Fig. 1), and reverses back to the original volume fraction during unloading.Uniaxial compression is defined by the strain-rate tensor where is the strain-rate (compression > 0 and decompression/tension < 0) amplitude applied in the uniaxial mode.The negative sign (convention) of the component corresponds to a reduction of length, so that tensile deformation is positive.Even though the strain is imposed only on the mobile "wall" in the -direction, which leads to an increase of compressive stress on this wall during compression, the nonmobile walls also experience some stress increase due to the "push-back" stress transfer and rearrangement of the particles during loading, as discussed in more detail in the following sections.This is in agreement with theoretical expectations for materials with non-zero Poissonʼs ratio (Spencer, 1980).However, the stress on the passive walls is typically smaller than that of the mobile, active wall, as consistent with findings from laboratory element tests using the biaxial tester (Morgeneyer and Schwedes, 2003;Zetzener and Schwedes, 1998) or the so-called meter.

Deviatoric
The preparation procedure as described in section 3.1 provides different initial configurations with volume fractions, .For the deviatoric deformation element test, unless stated otherwise, the configurations are from the unloading part (represented by blue stars in Fig. 1), to test the dependence of quantities of interest on the volume fraction during volume-conser ving deviatoric (pure shear) deformations.The unloading branch is more reliable since it is much less sensitive to the protocol and rate of deformation during preparation (Göncü et al., 2010, Kumar et al., 2012b).Two different ways of deforming the system deviatorically with conserved volume are used here.
The deviatoric mode D2 has the strain-rate tensor where is the strain-rate (compression in -direc-tion > 0) amplitude applied to the walls.We use the nomenclature D2 since two walls are moving while the third wall is stationary.The deviatoric mode D3 has the strain-rate tensor where is the strain-rate (> 0 for compression in -direction) amplitude applied.In this case, D3 signifies that all three walls are moving, with one wall twice as much (in the opposite direction) as the other two such that volume is conserved during deformation.Note that the D3 mode is uniquely similar in "shape" to the uniaxial mode1 , see Table 1, since in both cases two walls are controlled similarly.Mode D2 is different in this respect and thus resembles more an independent mode (pure shear), so that unless differently stated, we plot by default the D2 results rather than the D3 ones (see section 2).The mode D2 with shape factor (as defined in Table 1) = 0 is on the one hand a plane strain deformation, and on the other hand allows for simulation of the biaxial experiment (with two walls static while four walls are moving, see Morgeneyer and Schwedes, 2003;Zetzener and Schwedes, 1998).

Averaged Quantities
In this section, we present the general definitions of averaged microscopic and macroscopic quantities.The latter are quantities that are readily accessible from laboratory experiments, whereas the former are often impossible to measure in experiments but are easily available from discrete element simulations.

Averaged microscopic quantities
Here, we define microscopic parameters including the coordination number, the fraction of rattlers, and the ratio of the kinetic and potential energy.

Coordination number and rattlers
In order to link the macroscopic load carried by the sample with the microscopic contact network, all particles that do not contribute to the force network -particles with exactly zero contacts -are excluded.In addition to these "rattlers" with zero contacts, there may be a few particles with a finite number of contacts for a short time which also do not contribute to the mechanical stability of the packing.These particles are called dynamic rattlers (Göncü et al., 2010), since their contacts are transient: The repulsive contact forces will push them away from the mechanically stable backbone.Frictionless par ticles with less than 4 contacts are thus rattlers, since they are mechanically unstable and hence do not contribute to the contact network.In this work, since tangential forces are neglected, rattlers can be identified by just counting their number of contacts.This leads to the following abbreviations and definitions for the coordination number (i.e. the average number of contacts per particle) and fraction of rattlers, which must be reconsidered for systems with tangential forces or Table 1 Summary of the deformation modes and the deviatoric strain rates , as well as shape factors, , and Lode angles, , (Thornton and Zhang, 2010) for the different modes in the respective tensor eigensystem, with sorted eigenvalues , defined in section 4. : Corrected coordination number V p : Volume fraction of the particles, with V p as particle volume.
Some simulation results for the coordination numbers and the fraction of rattlers will be presented below in subsection 5.1.

Energy ratio and quasi-static criterion
Above the jamming volume fraction , in mechanically stable static situations, there exist permanent contacts between particles; hence the potential energy (which is also an indicator of the overlap between particles) is considerably larger than the kinetic energy (which has to be seen as a perturbation).
The ratio of kinetic energy and potential energy is shown in Fig. 3 for isotropic compression from = 0.640 to = 0.820 and back.The first simulation, represented by the solid red line, was run for a simulation time 5000 [ ] and the second (much slower) simulation, represented by the green dashed line was run for 50000 [ ].For these simulations, the maximum strain rates are 138 [ ] and 13.8 [ ], respectively.During compression with increasing volume fraction, the energy ratio generally decreases and slower deformation by a factor of 10 leads to more than 100 times smaller energy ratios with stronger fluctuations.Most sharp increases of the energy ratio resemble reorganization events of several particles and are followed by an exponentially fast decrease (data not shown).The decrease is controlled by the interaction and dissipation timescales and not by the shear rate; in other words the timescale of energy dissipation 10 -5 [s] is considerably smaller than the simulation timescale 7.2 10 -3 [s] or 7.2 10 -2 [s], so that kinetic energy can be well relaxed before a considerable rearrangement of particles takes place; due only to the scaling of , the decrease appears to be faster for the slower deformation.Explicitly, the rate of decay depends on material parameters only and is of the order of , however, it becomes larger the closer to jamming one gets.The large initial ratio of kinetic to potential energy indicates that the system is in the unjammed regime, whereas after some compression it enters the quasi-static regime with much smaller energy ratios (Thornton and Anthony, 1998).In this way, dynamic ef fects are minimized and the system is as close as feasible to the quasi-static state.For many situations, it was tested that a slower deformation did not lead to large, considerably different results.For the majority of the data presented, we have .Lower energy ratios can be obtained by performing simulations at even slower rates, but the settings used are a compromise between computing time and reasonably slow deformations.

Averaged macroscopic quantities
Now the focus is on defining averaged macroscopic tensorial quantities -including strain, stress and fabric (structure) tensors -that reveal interesting bulk features and provide information about the state of the packing due to its deformation.

Strain
For any deformation, the isotropic part of the in-Fig.3 Comparison of the ratio of kinetic and potential energy in scaled time for two simulations, with different period of one compressiondecompression cycle , as given in the inset.finitesimal strain tensor is defined as: (3) where dt with , or are the diagonal elements of the strain tensor in the Cartesian , , reference system.The trace integral of denoted by is the true or logarithmic strain, i.e. the volume change of the system relative to the initial reference volume (Göncü et al., 2010).Several definitions are available in literature (Imole et al., 2011;Thornton andZhang, 2006, 2010;Zhao and Evans, 2011) to define the deviatoric magnitude of the strain.For the sake of simplicity, we use the following definition of the deviatoric strain to account for all active and inactive directions in a triaxial experiment, regardless of the deformation mode, (4) Since, for our triaxial box, for all modes, the Cartesian coordinates resemble the fixed eigensystem, sor ting the eigenvalues according to magnitude leaves the eigenvalue as the maximal tensile eigenvalue, with corresponding eigendirection, and as the magnitude of the deviatoric strain 2 .The quantitative description of the tensor is completed by either its third invariant, the Lode angle (Thornton and Zhang, 2010) or, equivalently, by the shape factor , as given in Table 1.Note that the values for are during uniaxial loading where compression is performed in the -direction.The sorting will lead to different values, , when the strain is reversed for both UNI and D3 modes.

Stress
From the simulations, one can determine the stress tensor (compressive stress is positive as convention) components: (5) with particle , mass , velocity , contact , force and branch vector , while Greek letters represent components , , and (Luding, 2008a,b).The first sum is the kinetic energy tensor and the second involves the contact-force dyadic product with the branch vector.Averaging, smoothing or coarse graining (Weinhart et al., 2012) in the vicinity of the averaging volume, , weighted according to the vicinity, is not applied in this study, since averages are taken over the total volume.Since the data in this study are quasi-static, the first sum can mostly be neglected.The average isotropic stress (i.e. the hydrostatic pressure) is defined as: (6) where , and are the diagonal elements of the stress tensor in the , , and box reference system and is its trace.The non-dimensional pressure (Göncü et al., 2010) is defined as: (7) where is the mean radius of the spheres and is the contact stiffness defined in section 2. We define the deviatoric magnitude of stress (similar to Eq. ( 4) for the deviatoric strain) as: (8) which is always positive by definition.The direction of the deviatoric stress is carried by its eigen-directions (in the present case coincident with the Cartesian reference system), where stress eigenvalues are sorted like strain eigenvalues, possibly with a different sign convention, according to their magnitude.Eqs. ( 4) and ( 8) can easily be generalized to account for shear reversal using a sign taken from the orientation of the corresponding eigenvectors and eigenvalues, or from the stress shape factor, however, this will not be detailed here for the sake of brevity.
It is noteworthy to add that the definitions of the deviatoric stress and strain tensors are proportional to the second invariants of these tensors, e.g. for stress: , which makes our definition identical to the von Mises criterion (Fredlund, 1979;Hayhurst, 1972;Thornton andZhang, 2006, 2010) 3 .
2 The objective definition of the deviatoric strain defines it in terms of the eigenvalues , and , of the (deviatoric) strain tensor.However, since the global strain is given by the wall motion and the eigensystem stays practically unchanged during the deformation, the two definitions are equivalent for triaxial element tests.

Fabric (structure) tensor
Besides the stress of a static packing of powders and grains, the next most important quantity of interest is the fabric/structure tensor.For disordered media, the concept of the fabric tensor naturally occurs when the system consists of an elastic network or a packing of discrete particles.The expression for the components of the fabric tensor is: (9) where is the particle volume which lies inside the averaging volume and is the normal unity vector pointing from the center of particle to contact .
are thus the components of a symmetric ranktwo 3×3 tensor, such as the stress tensor.The isotropic fabric quantifies the contact number density as studied by Göncü et al., 2010.We assume that the structural anisotropy in the system is quantified (completely) by the anisotropy of fabric, i.e. the deviatoric fabric.To quantify it, we define a scalar similar to Eqs. ( 4) and (8) as: (10) where , and are the three diagonal components of the fabric tensor.The fabric tensor practically has only diagonal components with non-diagonal elements very close to zero, so that its eigensystem is close to the Cartesian reference system, as confirmed by eigenvalues analysis.Also for the fabric, a shape factor completes the picture.

Conclusion
Three macroscopic rank-two tensors were defined and will be related to microscopic quantities and each other in the following.The orientations of all the tensor eigenvectors show a tiny non-co-linearity of stress, strain and fabric, which we neglect in the next sections, since we attribute it to natural statistical fluctuations and consider a unique fixed eigensystem coincident with the Cartesian reference system for all our deformation modes.The shape factor, as defined for strain, can also be analyzed for stress and fabric, but this will be shown elsewhere.

Evolution of Micro-Quantities
In this section, we discuss the evolution of the microscopic quantities studied -including coordination number and fraction of rattlers -as a function of volume fraction and deviatoric strain, respectively, and compare these results for the different deformation modes.

Coordination number and rattlers
It has been observed by Göncü et al., 2010, that under isotropic deformation, the corrected coordination number follows the power law (11 where 6 is the isostatic value of in the frictionless case.For the uniaxial unloading simulations, we obtain 8.370, 0.5998 and 0.6625 as best-fit parameters.In Fig. 4, the evolution of the simple, corrected and modified coordination numbers is compared as a function of the volume fraction during uniaxial de- formation (during one loading and unloading cycle).The contribution to the coordination number originating from particles with 1, 2 or 3 is small -as compared to those with 0 -since and are very similar, but always smaller than due to the fraction of rattlers, as discussed below.The number of contacts per particle grows with increasing compression to a value of 9.5 at maximum compression.During decompression, the contacts begin to open and the coordination number decreases and approaches the theoretical value = 6 at the critical jamming volume fraction 4 after uniaxial decompression 0.662.Note that the value is smaller than 0.665 reached after purely isotropic over-compression to the same maximal volume fraction.The coordination numbers are typically slightly larger in the loading branch than in the unloading branch, due to the previous over-compression.In Fig. 5, we plot the corrected coordination number for deformation mode D2 as a function of the deviatoric strain for five different volume fractions.Two sets of data are presented for each volume fraction star ting from dif ferent initial configurations, either from the loading or the unloading branch of the isotropic preparation simulation (cyan dots and blue stars in Fig. 1).Given initial states with volume fractions above the jamming volume fraction, and due to the volume-conserving D2 mode, the value of the coordination number remains practically constant.It is only for the lowest volume fractions close to jamming, that a slight increase (decrease) in can be seen for the initial states chosen from the unloading (loading) branch of the preparation step.However, both reach similar steady-state values after large strain, as indicated by the solid lines.Hence, for further analysis and unless otherwise stated, we will only present the steady-state values of micro-and macro-quantities from deviatoric modes D2 and D3.The rearrangement of the particles during shear thus does not lead to the creation (or destruction) of many contacts -on average.There is no evidence of the change of average contacts after 10 -15 percent of strain.However, close to jamming, a clear depen-dence of on the initial state exists, which vanishes in steady state when one gets saturated values in micro-and macro-quantities after large enough strain.For the same volume fraction, we evidence a range of , where the subscripts refer to overcompressed, steady, and initially compressed states, respectively.The coordination number, or alternatively the contact number density, as related to the trace of the fabric tensor (Göncü et al., 2010), is thus a control parameter closely linked to the volume fraction that contains more information about the structure than itself (above the jamming volume fraction), see Magnanimo et al., 2008 and references therein.In Fig. 6, the corrected coordination number is shown as a function of volume fraction for the purely isotropic and for the uniaxial unloading data as well as for the large strain deviatoric deformation datasets.Different symbols show the values of for the different deformation modes for various volume fractions.Interestingly, the power law for the coordination number derived from isotropic data describes well also the uniaxial and deviatoric datasets, with coefficients given in Table 2.This suggests that (for the cases considered, when particles are frictionless) the coordination number is almost independent of the deviatoric strain in steady state, and the limit values can be approximated by Eq. ( 11) as proposed for simple isotropic deformation.The distinction between the modes at the small (isotropic) strain region ４ The value, 6, is expected since it is the isostatic limit for frictionless systems in three dimensions (Göncü et al., 2010;Maxwell, 1864), for which the number of constraints (contacts) is twice the number of degrees of freedom (dimension) -in average, per particle -so that the number of unknown forces matches exactly the number of equations.( is different from the minimal number of contacts needed for a single mechanically stable frictionless sphere in 3D).
Fig. 5 Evolution of coordination numbers with deviatoric strain for the D2 mode.Smaller symbols represent data with initial configuration from the loading branch of an isotropic simulation, while the larger symbols start from an initial configuration with the same volume fraction, but from the isotropic unloading branch.The horizontal line at the large strain of the dataset indicates an average after saturation at steady state.
is shown zoomed in the inset of Fig. 6.The mixed mode (uniaxial) is bordered on both sides by the pure modes, namely isotropic and deviatoric (D2 and D3 cannot be distinguished), indicating that the two pure modes are limit states or extrema for .Alternatively, the range in values can be seen as caused by a range in , with , which represent the maximal jamming volume fraction after previous (isotropic, strong) over-compression, the intermediate jamming volume fraction after (mixed mode) deformation, and the minimal jamming volume fraction after large deviatoric strain, respectively, with 0.6646 and 0.6602.In other words, deviatoric deformations reduce the jamming volume fraction of the packing, i.e. can disturb and dilate a dense (over-compressed) packing so that it becomes less efficiently packed.This is opposite to isotropic over-compression, where after unloading, the jamming volume fraction is higher, i.e. the system is more efficiently packed/structured.This behavior is qualitatively expected for frictional particles, however, this is to our knowledge the first time that this small but systematic range in the jamming volume fractions is reported for frictionless packings -where the most relevant and only mechanism is structural reorganization, as will be discussed further in section 6.1.1.As a related interesting microscopic quantity, we recall the analytical expression for the fraction of rattlers proposed by Göncü et al., 2010: (12) where the fit parameters for the different deformation modes are given in Table 2, and 0.6646 is obtained from extrapolation of to the isostatic coordination number 6.In Fig. 7, the evolution of the fraction of rattlers is plotted as a function of volume fraction for both isotropic and uniaxial unloading as well as for steady-state deviatoric mode simulations.We then compare these with the prediction/ fit (solid lines) from the exponential decay equation, Eq. ( 12).Interestingly, in contrast to the coordination number, the fraction of rattlers displays stronger differences at the highest volume fraction ( 0.82 in Fig. 7), and it is lower during isotropic unloading as compared to the steady-state deviatoric mode situations, and somewhat higher during uniaxial unloading.The difference between the modes is small close to jamming with largest for the UNI mode.For uniaxial simulations, at the end of unloading close to , a considerable fraction (almost 20 percent) of Table 2 Fit parameters for the analytical predictions of coordination number, fraction of rattlers, and pressure in Eq. ( 11), with 6 and Eqs. ( 12) and ( 15), respectively.
is used from the fits of to fit for the different deformation modes.The first rows of isotropic data ISOG are from Göncü et al., 2010, for   the total number of particles are rattlers that do not contribute to the stability of the network.For higher volume fractions, a strong exponential decay is evidenced5 .To better understand the peculiar behavior of the jamming volume fraction under the different modes of deformation, some macroscopic quantities are studied next.

Evolution of Macro-Quantities
In the following, we discuss the evolution of the macroscopic tensors, stress and fabric, as defined in section 4.2.For clarity, we split them into isotropic and deviatoric parts in subsections 6.1 and 6.2, respectively.

Evolution of macro-quantities: isotropic 6.1.1 Isotropic pressure
In this section, the relation between pressure and volume fraction is studied.First, we consider the contact overlap/deformation , since the force is directly related to it and stress is proportional to the force.The infinitesimal change of the normalized average overlap , can be related to the volumetric strain under the simplifying assumption of uniform, homogeneous deformation the packing.As defined in subsection 4.2.1, is the average of the diagonal elements of the infinitesimal strain tensor, and 0.425 is a proportionality constant that depends on the size distribution and can be readily obtained from the average overlap and volume fraction (data not shown), see Eq. ( 13).The integral of denoted by is the true or logarithmic volume change of the system relative to the reference volume .This is chosen without loss of generality at the critical jamming volume fraction , so that the normalized average overlap is (Göncü et al., 2010): (13) As in Eq. ( 7), see Refs.Göncü et al., 2010;Shaebani et al., 2012 for details, the non-dimensional pressure is: and the scaled pressure is: where 0.040, 0.033, and the critical volume fraction 0.6625 are fit parameters to the pressure law for uniaxial unloading.Combining the quasi-static part of Eq. ( 5) with ( 14) leads to the proportionality relation , which makes a measure for the average overlap relative to the average particle diameter.On the other hand in Eq. ( 15) scales various different pressures for different deformation modes to their respective reference jamming volume fractions and is linear for small .Note that the critical volume fraction 0.6625, as obtained from extrapolation of to the isostatic coordination number 6, is ver y close to that obtained from Eq. ( 14).When fitting all modes with pressure, one confirms again that falls in between the limits of the pure modes and (with values consistent within each mode) as summarized in Table 2.In Fig. 8, we plot the total (non-dimensional) pressure as a function of the deviatoric strain for different volume fractions in the deformation mode D2.Above the jamming volume fraction, the value of the pressure stays practically constant with increasing shear strain.A slight increase in can be seen for Fig. 7 Evolution of the fraction of rattlers as a function of volume fraction during unloading for all modes.The symbols represent the respective simulation data.The solid lines are the analytical fits of Eq. ( 12) for each mode with the values of fit parameters and for each mode shown in Table 2.The arrow indicates the unloading direction.
the lowest volume fractions when the initial states are chosen from the unloading branch of isotropic modes, whereas a slight decrease in is observed for initial states chosen from the loading branch.Independently of the initial configuration, the pressure reaches a unique steady state at large strain, similarly to what is observed for the coordination number in Fig. 6.In Fig. 9, the total pressure is plotted against volume fraction for isotropic, uniaxial and deviatoric (D2 and D3) modes, with data obtained from the unloading branch in the first two cases and after large deviatoric strain for the deviatoric modes (see Fig. 8).For these three modes, at large , the simulation data collapse on a unique curve with non-linear behavior.Due to the linear contact model, this feature can be directly related to the contact number density, i.e. the isotropic fabric, which quantifies the isotropic, direction-independent changes of structure due to rearrangements and closing/opening of contacts.When approaching the jamming transition, the pressure values diverge slightly (inset of Fig. 9) due to the difference in the critical volume fractions .In Fig. 10, we plot the scaled pressure defined in Eq. ( 15) against the volumetric strain from the same data as in Fig. 9.The three datasets almost collapse for small strain.For increasing volume fractions (larger ), the scaled pressure in the isotropic mode is considerably larger than the uniaxial and deviatoric modes, where again the uniaxial data fall in between isotropic and deviatoric values.This resembles the behavior of and is consistent with the fact that the uniaxial mode is a superposition of the purely isotropic and deviatoric deformation modes.The dependence of pressure on isotropic strain can be interpreted in relation to the sample history.The deviatoric modes (D2 or D3) lead to dilatancy and thus to higher steady-state pressure, with low ; the isotropic mode is strictly compressive, with the lowest pressure after over-compression, during unloading and the highest ; finally, the uniaxial mode is a mixed mode and thus interpolates between the two other modes.The apparent collapse of all scaled data at small strain, with similar pre-factors 0.040 is interesting since, irrespective of the applied deformation mode -purely isotropic, uniaxial, and D2 or D3 deviatoric, it boils down to a linear relation between and with a small quadratic correction -different from the non-linear power laws proposed in previous studies, e.g. in Majmudar et al., 2007.The non-linearity due to is hidden in , which is actually proportional to the isotropic fabric .

Isotropic fabric
The random isotropic orientation of the contact directions in space was studied in detail in Refs.Göncü et al., 2010 andShaebani et al., 2012, andis referred Fig. 8 Evolution of (non-dimensional) pressure, Eq. ( 7), with deviatoric strain for the D2 deformation mode, at different initial volume fractions .Small and large symbols represent simulations starting with initial isotropic configurations from the loading and unloading branch, respectively.The horizontal line at the large strain of the dataset indicates an average value of the pressure after saturation at steady state.
Fig. 9 Total (non-dimensional) pressure, Eq. ( 14), plotted as a function of volume fraction for the uniaxial and isotropic datasets during unloading, and for the D2/D3 deviatoric modes after large strain.The solid lines are the analytical fits of Eq. ( 14) for each mode, with parameters , and shown in Table 2.
to as the contact number density with , where is of order unity and depends only on the size distribution (for our case with 3, one has 1.22).Note that directly connects to the dimensionless pressure which, remarkably, hides the corrected coordination number and the fraction of rattlers in the relation , which fully determines .

Evolution of macro-quantities: deviatoric
In the following, we will show the evolution of the deviatoric stress ratio (which can be seen as a measure of stress anisotropy) and the structural anisotropy, both as a function of the deviatoric strain.
In particular, we present the raw data from deviatoric and uniaxial simulations and their phenomenology.The D2 volume-conserving simulations are used to calibrate the constitutive model, as presented in Refs.
Luding and Perdahcioğlu, 2011, Magnanimo ing, 2011 and described in section 7. We further use the fitting parameters inferred from deviatoric data to predict the evolution of stress and fabric during uniaxial deformation.

Deviatoric stress
The behavior of the deviatoric stress ratio, during deformation mode D2, is shown in Fig. 11 as a function of the deviatoric strain for various dif ferent volume fractions.The stress ratio initially grows with applied strain until an asymptote (the maximum stress anisotropy) is reached where it remains fairly constant.The asymptote is referred to as the deviatoric steady state or the macroscopic friction, where represents the mobilized friction at each step along the loading path.For lower volume fractions, higher maximum are reached and the deviatoric stress ratio increases faster, meaning a higher ratio, where is the octahedral shear modulus as defined by Barreto and OʻSullivan, 2012.This is opposite to what is expected for the shear modulus itself, being proportional to the volume fraction.Interestingly, the stress response observed for mode D3 (not shown) follows a very similar path as for mode D2, resembling independency of the results with respect to the particular deviatoric path, as will be discussed in more detail in section 7. We use the deviatoric simulations to fit the exponential relation proposed in Luding and Perdahcioğlu, 2011;Magnanimo and Luding, 2011, for the biaxial box and report the theoretical cur ves in the same Fig.11.Both the initial growth rate coefficient and the asymptotic values are inferred from the volumeconserving deviatoric data, following the procedure described in Section 7. We point out here that the softening behavior after maximal is ignored in the fitting procedure for the theoretical model, since we do not want this feature of the material to be plugged into the model as an additional element.The stress-strain behavior in the case of uniaxial compression is shown in Fig. 12 star ting from initial volume fractions 0.671, 0.695, 0.728, to a common maximum value 0.820.Unlike the volume-conserving deviatoric simulations discussed previously, the evolution of the deviatoric stress ratio during uniaxial compression leads to large fluctuations that do not allow the clear observation of a possible softening/hardening regime.This difference is because the uniaxial deformation mode has a continuously increasing density and pressure in contrast, for example, to mode D3 where is increasing and are decreasing such that the pressure remains (almost) constant.The solid lines superimposed on the data in the plot represent the predictions of the constitutive relation in Eq. ( 17), with the parameters obtained from the deviatoric modes D2 and D3, as explained in detail in section 7.Moreover, as the deviatoric strain increases during uniaxial deformation, the deviatoric stress ratio also increases with values that depend on the initial volume fraction.For lower volume fractions we observe a higher stress ratio, similar to what is obser ved in Fig. 11.Interestingly, uniaxial deformations for different initial volume fractions lead to Fig. 10 The scaled pressure plotted against the (negative) volumetric strain for the same data as presented in Fig. 9.The solid lines are the predictions from Eq. ( 15) using the fits of and for each mode.
convergence (and almost collapse) after about 7.5 percent deviatoric strain.This feature of the uniaxial simulations is also well captured by the anisotropy model in section 7.3.

Deviatoric fabric
The evolution of the deviatoric fabric as a function of the deviatoric strain is shown in Fig. 13 for mode D2 simulations and three different volume fractions.It builds up from different random, small initial values and reaches different maximum saturation values .The deviatoric fabric increases faster at lower volume fractions in a similar fashion to what was observed for the stress ratio in Fig. 11.Both the growth rate and the maximum deviatoric fabric are well defined and shown in Fig. 13 together with the fitting curves used to deduce the theoretical values and for different volume fractions (see details in section 7.2).The evolution of the deviatoric fabric for the D3 mode is not shown since it resembles the behavior of the D2 mode, implying that the fabric evolution is pretty much insensitive to the deviatoric deformation protocol employed, as was observed before also for the stress ratio.A more detailed study of the (small) differences among deformation modes with different shape factors, as predicted by Thornton and Zhang, 2010, will be reported elsewhere.Fig. 14 shows the evolution of the deviatoric fabric during uniaxial compression as presented in section 6.2.1.The deviatoric fabric builds up as the deviatoric strain (and the volume fraction) increases.It begins to saturate at 0.06 and a slight decreasing (softening) trend is seen towards the end of the loading path.The convergence of the deviatoric stress after large strain, for different volume fractions, as seen in Fig. 12, does not appear so clearly for the deviatoric fabric.The theoretical prediction of the constitutive relations from section 7 is in good qualitative agreement with the numerical data, but over-predicts the deviatoric fabric for larger strains.Their analytical form and the parameters involved will be discussed next.

Theor y: Macroscopic Evolution Equations
Constitutive models are manifold and most standard models with wide application fields such as elasticity, elasto-plasticity, or fluid-/gas-models of various kinds were applied also to granular flowssometimes with success, but typically only in a very limited range of parameters and flow conditions; for overviews see, e.g.GdR-MiDi, 2004;Luding and Alonso-Marroquín, 2011.The framework of kinetic theor y is an established tool with quantitative predictive value for rapid granular flows only -but it is hardly applicable in dense, quasi-static and static situations (Luding, 2009).Further models, such as hyper-or hypo-elasticity, are complemented by hypo-Fig.11 Deviatoric stress ratio plotted against deviatoric strain from the D2 deformation mode for initial volume fractions during unloading, from which the simulations were performed, as given in the inset.The symbols (ʻ*ʼ, ʻ×ʼ and ʻ+ʼ) are the simulation data while the solid lines through them represent a fit to the data using Eq. ( 17).

Fig. 12
Deviatoric stress ratio plotted against deviatoric strain from the uniaxial compression mode data for different initial volume fractions during unloading, from which the uniaxial deformations were initiated, as given in the inset.The symbols ( ʻ*ʼ , ʻ×ʼ and ʻ+ʼ ) are the simulation data while the solid lines represent the prediction using Eq. ( 17).plasticity (Kolymbas, 1991) and the so-called granular solid hydrodynamics (Jiang and Liu, 2009), where the latter provides incremental evolution equations for the evolution of stress with strain, and involve limit states (Mašín, 2012) instead of a plastic yield surface as in plasticity theory.A strict split between elastic and plastic behavior seems invalid in granular materials, see, e.g.Alonso-Marroquín et al, 2005.More advanced models involve so-called non-associated / non-coaxial flow rules, where some assumptions on relations between different tensors are proposed, see Thornton andZhang, 2006, 2010.While most of these theories can be or have been extended to accommodate anisotropy of the microstructure, only very few models account for an independent evolution of strain, stress and microstructure (see, for example, Thornton and Zhang, 2010;Sun and Sundaresan, 2011;Luding and Perdahcioğlu, 2011;Goddard, 2010) as found to be important in this study and many others.
In the following, we use the anisotropy constitutive model as proposed in Kumar et al., 2012a;Luding and Perdahcioğlu, 2011;Magnanimo and Luding, 2011, generalized for a dimensional system: In its simplest form, the model involves only three moduli: the classic bulk modulus (Göncü et al., 2010), the octahedral shear modulus , and the new variable "anisotropy modulus" , evolving independently of stress with deviatoric strain.Due to , the model provides a cross-coupling between the two types of stress and strain in the model, namely the hydrostatic and the shear (deviatoric) stresses react to both isotropic and deviatoric strains.
is an abbreviation for the stress isotropy with .The parameter resembles the macroscopic friction and is the growth rate of with deviatoric strain .The parameter in the evolution equation of represents the maximum anisotropy that can be reached at saturation, and determines how fast the asymptote is reached (growth rate).Both and are model parameters for the anisotropy modulus and can be extracted from fits to the macroscopic, average simulation results.Note that the evolution of is assumed to be kinematic, i.e. not explicitly dependent on pressure, but there is a possible volume fraction dependence of and , as detailed below.In the following, we test the proposed model by extracting the model parameters as functions of volume fraction from various volume-conserving deviatoric simulations.The calibrated model is then used to predict the uniaxial deformation behavior (see the previous section).The theory will be discussed elsewhere in more detail (Kumar et al., 2012a;Magnanimo and Luding, 2012).In short, it is based on the basic postulate that the independent evolution of stress and structure is possible.It comes together with some simplifying assumptions such as: the new macroscopic field is proportional to the microscopic rank-two deviatoric fabric so that they have the same non-dimensional Fig. 13 Deviatoric fabric plotted against deviatoric strain from the D2 deformation simulations of Fig. 11 The symbols (ʻ*ʼ, ʻ×ʼ and ʻ+ʼ) are the simulation data while the solid lines through them represent a fit to the data using Eq. ( 18).
Fig. 14 Deviatoric fabric plotted against deviatoric strain from the uniaxial deformation simulations in Fig. 12 The symbols (ʻ*ʼ , ʻ×ʼ and ʻ+ʼ) are the simulation data while the solid lines through them represent the prediction of the data using Eq.(18).
growth rates ; both and -to lowest order, i.e. neglecting additional (missing) terms in Eqs.(16)approach their limit states exponentially fast; only one anisotropy modulus is sufficient (valid in 2D, questionable in 3D, possibly two moduli and are needed); that lead to Eqs. ( 17) and ( 18) below.We use these two equations as empirical fit functions, since they are special cases of the complete constitutive model with anisotropy, and then use the fit result to predict another solution of the (simplified) theory for another deformation mode.

Reduced theoretical model
The reduced model consists of two evolution equations for the deviatoric stress ratio related to the mobilized macroscopic friction, and the deviatoric fabric based on DEM observations in 2D, see Luding, 2004, 2005. For pure shear, Figs. 11 and13 show that and grow non-linearly until they exponentially a constant value at steady state, with fluctuations where the material can be indefinitely sheared without fur ther change.As discussed by Luding and Perdahcioğlu, 2011, the coupled evolution equations ( 16) are (with above assumptions) consistent with approximated by: (17 where and represent the initial and maximum values of and is its growth rate.Similarly, the deviatoric fabric is approximated by:  17) and ( 18).As a final step, but not shown in this paper, in order to relate the macroscopic anisotropy (modulus) to the evolution of the deviatoric fabric , one can measure the elastic modulus directly.For this, the sample is subjected to incremental deformations (either isotropic or purely deviatoric) at various different stages along the (large strain) deviatoric paths for D2 and D3 deformations.Details of the procedure and the results will be reported elsewhere (Kumar et al., 2012a).Here, we only note that a linear relation is found such that: (19) where 0.137 is a combination of numerical constants including , .see Fig. 15.The factor decreases with increasing volume fraction and a similar trend is observed for with some larger scatter.Both and seem to saturate towards a finite limit for large volume fractions, and these values can be extrapolated by the fitting procedure described later in this sec-tion.The fit procedure applied to the deformation modes D2 and D3 leads to ver y similar results for and .This is not surprising: the same net deviatoric strain applied in the two modes leads to (almost) the same net deviatoric stress ratio response, even though the shapes of deformations are different.Fig. 16(a) shows the variation of with volume fraction for the same simulations as in Fig. 15, where the two deviatoric deformation modes D2 and D3 almost collapse on each other.decreases strongly with volume fraction for the two For higher volume fractions, the motion of spheres is more constrained by more contacts and hence the contact anisotropy developed in the system is smaller.3, and 0.6653 is chosen as the jamming volume fraction, see Table 2.For all four parameters, the values are the limit for large volume fractions, while represents the limit at and is the rate of variation (decay) with the volume fraction increasing above .We assume, as is consistent with the data, that the structural anisotropy parameters and tend towards zero as the volume fraction increases, therefore we keep = 0 in the fitting functions.Eq. ( 20) represents the solid black lines shown in Figs. 15 and 16, with coefficients given in Table 3.

Prediction of uniaxial deformation
We use the parameters determined from the deviatoric simulations presented in Table 3 to predict the behavior of uniaxial simulations in subsection 6.2,  where the volume fraction is changing with deviatoric strain and hence dependence on is needed to properly describe the deformation path.Fig. 12 above shows the normalized deviatoric stress against the deviatoric strain for uniaxial deformations star ting from three dif ferent volume fractions ( 0.671, 0.695 and 0.728), compared with the predictions of Eq. ( 17) with coefficients and taken from Eq. ( 20) and coefficients from Table 3.The proposed although in its simplified version, is able to properly capture the behavior of the material qualitatively, approaching exponentially a maximum value and then decreasing due to the volume fraction dependence of the parameters.Note that the softening present in some of the deviatoric DEM data is on purpose not plugged into the model as a constraint, which renders the softening present cal curves in Fig. 12 as a valuable prediction of the model.Fur thermore, the convergence of for uniaxial loading simulations different initial volume fraction at large strains, as discussed in section 6.2.1, is also well captured by the theoretical model with calibrated parameters from the deviatoric simulations (where this does not happen).Fig. 14 shows the evolution of the deviatoric fabric with deviatoric strain for uniaxial deformations -as above -together with the predictions of Eqs. ( 18) and ( 20).The model is still able to qualitatively describe the behavior of the deviatoric fabric, but with up to 30 percent over-prediction.For better quantitative agreement, the complete coupled model needs to be used and possibly improved as will be presented elsewhere (Kumar et al., 2012a).

Conclusions and Outlook
The discrete element method has been used to investigate the bulk response of periodic, polydisperse, frictionless sphere packings in 3D, subjected to various deformation modes, in terms of both their microand macroscopic responses.The main goal was to present a procedure to calibrate a constitutive model with the DEM data and then to use the same to predict another independent simulation (mode).The (overly) simple linear material (model) allows us to focus on the collective/bulk response of the material to different types of strain, excluding complex effects due to normal or tangential non-linearities.Therefore, the present study has to be seen as a reference "lower limit", and the procedure rather than the material is the main subject.
We focused on the strain-controlled loading and unloading of isotropic, uniaxial and two deviatoric (pure shear) type deformation modes (D2 and D3).Experimentally most difficult to realize is the isotropic deformation, while both uniaxial and deviatoric modes can be realized in various element tests where, however, often mixed strain-and stress-control procedures are applied.Both micro-mechanical and coarse-grained macroscopic properties of the assemblies are discussed and related to each other.The study covers a very wide range of isotropic, uniaxial and deviatoric deformation amplitudes and thus practically all volume fractions with mechanically stable packings -except for those very close to the jamming transition and higher than about 10 percent contact deformation, above which DEM pair contact models become questionable.

Microscopic quantities
The microscopic coordination number , defined as the ratio of the total number of contacts to the total number of particles, has been analyzed as a function of volume fraction and deviatoric strain.By disregarding particles with less than four contacts (called rattlers), the corrected coordination number is well described by Eq. ( 11) for all deformation modes (since the particles are frictionless).For the uniform size distribution used here, the fraction of rattlers shows an exponentially decaying trend towards higher volume fractions, very similar for all modes, see Eq. ( 12) and Table 2.These analytical relations provide a prediction for the coordination number that notably shows up in the macroscopic relations for both pressure and isotropic fabric, in combination with volume fraction , instead of .Note that is better accessible to theory, while is related to the wave-propagation speed, which is experimentally accessible, while both are linked by the fraction of rattlers, which was already identified as a control parameter of utmost importance (Bi et al., 2011).A small but systematic difference in and parameters appears for the dif ferent deformation modes.Most important, the jamming volume fraction is not a single, particular volume fraction, but we obser ve a range of jamming volume fractions dependent on the deformation modes, i.e. the "histor y" of the sample.Over-compression leads to an increase of , i.e. to a better, more efficient packing.Subsequent deviatoric (pure shear) deformations slightly reduce the jamming volume fraction of such a previously over-compressed packing, causing it to become less efficiently packed.Thus more/less efficient packing is reflected by a large/small jamming volume fraction and, inversely, small/large coordination numbers.The obser ved differences are more pronounced as the volume fraction becomes lower.A slight increase in the fraction of rattlers due to deviatoric deformations is also reported, as consistent with the decrease in coordination number.

Macroscopic quantities
When focusing on macroscopic quantities, an important result this study is that at small strains, the uniaxial, deviatoric and isotropic modes can be described by the same analytical pressure evolution, Eq. ( 15), with parameters given in Table 2, evidenced by the collapse of the data from these deformation modes when the scaled pressure is plotted as a (linear) function of the volumetric strain.This linearity is due to the scaling with the nonlinear terms in particular.Thanks to the linear contact model used, it allows the conclusion that the non-linear (quadratic) corrections are due to the structural rearrangements and non-affine deformations.The scaled deviatoric and uniaxial results deviate from the isotropic data.This appears at larger strains due to the build-up of anisotropy in the caused by deviatoric strain, obviously not present in the isotropic deformation mode.The good match of the data suggests an advantage of the "cheaper" uniaxial (and deviatoric) deformation modes over the experimentally difficult to realize isotropic deformation mode (three walls have to be moved simultaneously in the isotropic case, while a less complicated set-up is required for the other modes).The deviatoric stress ratio (the deviatoric stress scaled with the isotropic pressure) as a function of the deviatoric strain develops almost independently of the volume fraction when the deviatoric magnitude is defined in a similar fashion to the second deviatoric invariant (Thornton and Zhang, 2006) for all quantities studied.The deviatoric stress builds up with increasing deviatoric strain until a steady state is reached (where we do not focus on peak and softening behavior, which becomes more pronounced when the jamming volume fraction is approached).Starting from isotropic initial configurations, we also show that the slope of the normalized deviatoric stress as a function of deviatoric strain decreases with increasing volume fraction, unlike the shear modulus , which increases with volume fraction.This indicates that the pressure (and bulk modulus ) has a "stronger" dependence on the volume fraction than the shear modulus.From the macroscopic data, one observes that deviatoric and isotropic stresses and strains are crosscoupled by the structural anisotropy.The latter is quantified by the deviatoric fabric, which is proportional to the bulk anisotropy modulus/moduli , as relevant for the constitutive model.Cross-coupling means that in the presence of structural anisotropy, isotropic strain can cause deviatoric stress responses and deviatoric strain can cause isotropic stress responses (dilatancy or compactancy).The structural anisotropy response to deviatoric strain is very similar to that of the deviatoric stress ratio.The response rates of both stress and structure anisotropies with deviatoric strain are functions of volume fraction and, most important, are different from each other.

Constitutive model calibration
As a first step, the parameters of a simple constitutive model that involves anisotropy as proposed for 2D by Luding andPerdahcioğlu, 2011, Magnanimo andLuding, 2011 have been calibrated from DEM data.From the isotropic deformation mode, one can extract the bulk modulus , as was done by Göncü et al., 2010.From the volume-conserving D2 and D3 modes, by fitting the idealized evolution equations for shear stress in Eq. ( 17), the macroscopic friction and the deviatoric stress rate can be inferred as functions of the volume fraction, entering the shear modulus .Similarly, the fit of Eq. ( 18) provides a relation for the maximum fabric anisotropy at steady state and the fabric rate as functions of volume fraction.A relation between the deviatoric fabric and the anisotropy modulus/moduli in the model is finally needed to close the system and allow integration of the coupled evolution equations for stress and structure.As a second step and major result, the constitutive model calibrated on deviatoric data has been used to predict qualitatively (and to some extent also quantitatively) both stress and fabric evolution under uniaxial deformation.This is very promising, since the basic qualitative features are caught by the model, even though it was used in a very idealized and short form, with the single anisotropy modulus the only main new ingredient.Several additional terms of assumed minor impact are ignored and have to be added to complete the model, see Kumar et al., 2012a;Magnanimo and Luding, 2011;Sun and Sundaresan, 2011;Thornton and Zhang, 2010, and references therein, as postponed to future studies.(For example, an objective tensorial description of stress, strain and fabric involves also the third tensor invariants.Alternatively/equivalently, these deviatoric tensors can be completely classified by the shape factors in their respective eigen-systems, which allow us to distinguish all possible deformation and response modes in 3D.)

Outlook
In this paper, we have reviewed and presented new results for frictionless particles undergoing isotropic, uniaxial, and (pure) shear deformation.Since the particles are too idealized here, the results cannot be applied to practical systems where shape, friction, and other non-linearities are relevant.However, they form the basic reference study with details given on the calibration procedure that yields a constitutive model with satisfactory predictive quality.Therefore, the next steps in our research will involve more realistic contact models with friction, cohesion, and other physically meaningful material parameters.Only then can the validity of the analytical expressions be tested for realistic systems, to predict well the phenomenology for pressure as as the scaling arguments for the deviatoric stress and fabric.Laboratory element test experiments should also be performed with the biaxial box to validate the simulation results with realistic material properties.Macroscopic quantities that can be readily obtained experimentally -for example, the pressure-volume fraction relation and the shear stress evolution with deviatoric (pure shear) strain -can then be compared with simulation data.Moreover, the work underlines the predictive power of constitutive models with anisotropy (as in Luding and Perdahcioğlu, 2011;Magnanimo and Luding, 2011;Sun and Sundaresan, 2011) that can be further tested, validated and extended with more advanced physical and numerical experiments.
Given the detailed insights from DEM, the (missing) terms and the parameters for the constitutive models can now be further analyzed to perform the rigorous micro-macro transition.Open questions concern, among others: the validity of the 2D model in 3D, related to missing terms and parameters, the validity of global versus local coarse-graining, i.e. the scale of the micro-macro transition (Kuhn and Bagi, 2009), the microscopic (restructuring) and macroscopic (non-affine motions) origin of the peak and softening phenomenology at low volume frac-tions, possibly related to the (in)homogeneity of the packings, the validity of the model predictions for strainreversal and cyclic deformations, and the possible dependence of the moduli in the constitutive relations on other quantities (e.g.pressure) than the volume fraction, as focused on in this study.
For future application, the present calibration procedure should be checked also for other materials and applied to different element tests, among which there are (cylindrical) triaxial tests, ring-shear tests and also avalanche flow experiments like in chutes or rotating drums, all of which are more widely available than the "academic" biaxial box.In the end, the material properties and parameters should not depend on the element test chosen and the predictive value of the model(s) should be proven for more than only one validation test, be it another element test or a real-size or lab-scale process such as, e.g.granular flow in a silo or during a landslide.
kg/s].The polydispersity of the system is quantified by the width of a uniform distribution with a step function as defined in Göncü et al., 2010, where 1.5 [mm] and 0.5 [mm] are the radius of the biggest and smallest particles, respectively.

Fig. 1
Fig. 1 Evolution of volume fraction as a function of time.Region A represents the initial isotropic compression up to the initial volume fraction .B represents relaxation of the system and C represents the subsequent isotropic compression up to 0.820 and then decompression.Cyan dots represent some of the initial configurations, at different , during the loading cycle, and blue stars during the unloading cycle, both of which can be chosen for further study.

Fig. 2
Fig. 2 Generic schematic representation of the procedure for implementing isotropic, uniaxial and deviatoric deformation element tests.The isotropic preparation stage is represented by the dashed box.The corresponding plots (not to scale) for the deviatoric strain against volumetric strain are shown below the respective modes.The solid square boxes in the flowchart represent the actual tests.The blue circles indicate the start of the preparation; the red triangles represent its end, i.e. the start of the test, while the green diamonds show the end of the respective test.

Fig. 4
Fig. 4 Comparison between coordination numbers using the simple (ʻ+ʼ, blue), modified (ʻ◇ ・ʼ, green) and corrected (ʻ∇ʼ, red) definitions.Data are from a uniaxial compression-decompression simulation starting from 0.64 0.6625.The solid black line represents Eq. (11) with parameters given in the text very similar to those measured in Göncü et al., 2010, see Table 2.The compression and decompression branches are indicated by arrows pointing right and left, respectively.

Fig. 6
Fig. 6 Evolution of the corrected coordination number as a function of volume fraction during unloading for all modes.The symbols represent the respective simulation data while the solid lines represent the analytical equation according to Eq. (11) with the respective values of , , and shown in Table 2.The inset shows the corrected coordination number at low volume fractions close to jamming.
and maximum (saturation) values of the deviatoric fabric and is its rate of change.To study the variation of the parameters , , and with volume fraction , during deviatoric deformation, we perform several isochoric simulations at different volume fractions , and obtain the coefficients as shown in Figs. 15 and 16 from fits to Eqs. (

7. 2 Fig. 15
Fig. 15 Comparison of evolution parameters from Eq. (17): the maximum normalized deviatoric stress and the growth rate plotted against volume fraction for the D2 and D3 deviatoric modes.Each point represents a unique simulation; the green ʻ * ʼ s represent the D2 mode while the blue ʻ ʼ s represent the D3 mode.The solid black line is the proposed analytical form in Eq. (20), with parameters given in Table 3.(b) Fig. 16(b) shows a similar decreasing behavior of with volume fraction , where stronger scatter is seen.The analytical fits of the normalized stress parameters ( and ) are shown for comparison.A different behavior of the normalized stress and the deviatoric fabric with respect to both parameters (maximum saturation value and the evolution rate) proves that stress and fabric evolve independently of deviatoric strain (La Ragione and Magnanimo, 2012), as is the basic postulate for the anisotropic constitutive model.For the fit, we propose a generalized analytical relation for both the stress parameters , parameters with values presented in Table

Fig. 16
Fig. 16 Comparison of evolution parameters from Eq. (18): the maximum anisotropy and the growth rate plotted against volume fraction for the D2 and D3 deviatoric modes.The solid black line is the proposed theory, Eq. (20), for and , respectively, while the red lines are the corresponding parameters and in Fig. 15.
various polydispersities and also during unloading, but for different over-compression density.

Table 3
Fitting coefficients for the parameters in Eqs.