How does static granular matter re-arrange for different isotropic strain rate?

The question of how soft granular matter, or dense amorphous systems, re-arrange their microstructure under isotropic compression and de-compression, at different strain rates, will be answered by particle simulations of frictionless model systems in a periodic three-dimensional cuboid. Starting compression below jamming, the systems experience the well known jamming transition, with characteristic evolutions of the state variables elastic energy, elastic stress, coordination number, and elastic moduli. For large strain rates, kinetic energy comes into play and the evolution is more dynamic. In contrast, at extremely slow deformation, the system relaxes to hyper-elastic states, with well-defined elastic moduli, in static equilibrium between irreversible (plastic) re-arrangement events, discrete in time. Small, finite strains explore those reversible (elastic) states, before larger strains push the system into new states, by irreversible, sudden re-arrangements of the micro-structure.

In the following section 2, we will first present a few typical particle simulations from loading-unloading cycles, for which we vary the compression strain rate over orders of magnitude -with the goal to zoom-in to the slow rate data to display understand what is going on during plastic, irreversible re-arrangement events.

DEM particle simulations
The discrete element method (DEM) particle simulations presented here are a simple element test in a periodic cubical cell, with only diagonal components of the strain-rate tensor active. Only isotropic loading/un-loading is presented here, from many different alternative loading paths (e.g., pure shear in plane strain or axial strain modes). Deformations are applied to all particles, each time-step, * e-mail: s.luding@utwente.nl while at the same time the periodic cell size also changes, accordingly. The reference model is using N = 4913 frictionless particles (µ p = 0), with particle diameters drawn from a random homogeneous size distribution with maximum to minimum width, w = d max p /d min p = 3, similar to the particles used in Ref. [12], and many references therein. The contact model is the simplest linear springdashpot model, set to frictionless (resulting in higher packing densities), since only this allows to focus on (irreversible) structural re-arrangements rather than contact sliding events.

Non-dimensionalization of DEM
The parameters given in the following with a prime, e.g., ρ p = 2000 or d p = 2, are used in the simulations shown in this paper. For working with units, there are two alternatives: Either one can read the numbers in chosen units 1 or the units are chosen based on physical properties to achieve non-dimensional quantities. The latter option is adopted here, i.e., the unit of length is chosen as the mean particle diameter, x u = d p = 2, so that d p = 1 is the dimensionless diameter. The second unit is the material density, ρ u = ρ p = 2000, so that one has the dimensionless density, ρ = φ, and thus the unit of mass, m u = ρ u (x u ) 3 , i.e., the particle mass, m p = (π/6). For the third unit one has several choices, where we adopt here the units of elastic stress, σ u = k n /d p , with the linear normal contact stiffness, k n = 10 5 , which yields the dimensionless stress σ = σ d p /k n , and results in the unit of time  Figure 1. Dimensionless pressure p = p d/k n plotted against volume fraction for loading (left) and unloading (right) with different strain rates indicated in the legend (every X means 10 times slower), as explained in the main text. Note that the data almost collapse, i.e., for these slow strain rates, pressure is almost independent of rate. The symbols at t 0 and t f mark beginning and ending of the simulations, respectively.
Kinetic to potential energy ratio E k /E p plotted against volume fraction, from the same simulations as in Fig. 1 with same colors/symbols. The horizontal lines are the ratio of unity (1), indicating the dynamic jamming/un-jamming transition point, and E 0 = 10 −8 , an arbitrary threshold below which the system XXS can be considered static, elastic, reversible in mechanical equilibrium.
In the chosen units, the dimensionless linear stiffness is k n = k n (t u ) 2 /m u = 1, and the linear contact viscosity, γ n = 10 3 , becomes γ n = γ n t u /m u = γ n /[k n t u ] = 4 × 10 −3 , with background viscosity, γ b = 10 2 , The consequent physically relevant properties are the restitution coefficient r = exp(−ηt c ) ≈ 0.855, with damping factor η = γ n /(2m 12 ), reduced mass, m 12 = 0.063, and contact duration, t c = π/ k n /m 12 − η 2 = 0.79, or t c = t c t u = 0.316, all considered for a contact between the largest and the smallest particle, with the larger vis- 2 The alternative dimensionless stress: σ γ = σ /[ρ p (d pγ ) 2 ], with the unit of time set by the shear rate, t u =γ , is more useful for collisional shear flows, see Ref. [23], and is thus not adopted here. cous damping time-scale, t v = 2m 12 /γ n ≈ 5, and the even larger background damping time-scale t b = 2m 12 /γ b ≈ 50. Note that this choice of units corresponds to setting t u ∝ t c .
Next, some physical quantities are plotted for loading (L) and unloading (UL). The dimensionless pressure, p = p d p /k n , in Fig. 1, shows not much differences between the different rates, even for the largest rate it is only a tiny bit larger. The energy ratio, E k /E p , in Fig. 2, is strongly rate dependent, with approximately two orders of magnitude difference between simulations that feature different rates with factors of 10 between each other.
The coordination numbers, C * , in Fig. 3, are also very similar, only the fastest simulation slightly undershoots the others, but on this scale not much differences can be seen.
The bulk (tangent) modulus, B, in Fig. 4, only shows up for the two slowest simulations, due to the small, arbitrary E 0 , and not much differences can be seen in the elastic modulus, only the peaks that indicate plastic, irreversible rearrangement events are different. Besides that, note that the peaks on the unloading branch are much less frequent due to the higher stability of previously precompressed packings, consistent with previous results, see Ref. [12] and references therein.

Rate-dependent jamming and un-jamming
In Fig. 5 we next focus on jamming and un-jamming at different rates. While faster simulations, panels (a) and (b), considerably deviate, the two slowest simulations, panels (c) and (d), are in (almost) perfect agreement, in particular concerning the unjamming density. Only the early jamming regime (at densities below 0.67) is different, with slightly lower pressure and discrete events visible for XXXS (d). Panel (a) represents the largest rate with considerable kinetic energy K = E k /E p even above jamming (here the kinetic contribution to total energy, K/(1 + K) ∈ [0 : 1] is plotted). The dynamic kick-in, left, in the yellow area, below jamming, is only a transient from the initial condition, starting at ρ 0 = 0.65778, with rather large initial kinetic energy, i.e., not to be seen as steady compression from far below jamming, rather starting from a dynamic state close to jamming.
But how to identify the jamming density? For decompression, the pressure P (blue lines) drops linearly and thus allows for a straightforward identification of φ u J , unlike during compression, when the pressure P increases rather smoothly, with some wiggles, displaying  even some drops, most clear for the smallest compression rate in panel (d). Jamming is not well defined during continuously ongoing compression, due to these plastic events overlapping; events can only be seen in P if there is enough time between them for the system to relax all fluctuation energy down close to zero, which is not possible for the faster rates. Thus, close to jamming, even the slowest simulation is not slow enough! As second indication of the jamming transition [7,24], the kinetic energy ratio K/(1 + K) (red) decays from unity down to zero in the jamming region (yellow) and reaches a low (yet rate-dependent) level with a considerable elastic stress backbone (cyan), before it drops even stronger from the right end of the cyan regime up. For ongoing compression, the kinetic energy becomes considerably smaller and smaller than the potential energy -except during plastic re-arrangement events (not visible in this scale).
The jamming transition is identified as the margin between yellow and cyan, which is placed at the transition of the coordination number C * (green) to above its isostatic value, C * ≥ C * 0 = 6 (thin black dotted line, for frictionless particles, as used here). This transition is based on a fit/smoothing of the data, and thus not considering the wiggles in the coordination number, as they are considerable during jamming.

Conclusion
To summarize, "microscopic" DEM particle simulations with different strain-rates were used to track jamming and unjamming during one loading (over-compression) and un-loading cycle. For the frictionless model system used here, a considerable difference between jamming and unjamming volume fractions is observed, which was explained previously by a series of irreversible system rearrangement events [12], -some local (involving only few particles), some global (involving many particles) -frequent and strong during loading -much less during unloading. Every event results in a change of the jamming density, increasing during loading, and only partly recovering during unloading, as described in detail in Ref. [25].
The "mesoscopic" system re-arrangement events that sometimes concern only a few particles, but sometimes also most of the system, can be tracked and identified from macroscopic information such as stress, or better kinetic energy, if the deformation rate is small enough. For larger pressure, the system is much less sensible to the rate, and -except for the fastest simulation -the system has enough time between events (discrete in time) to relax back to hyper-elastic, mechanically stable states -reversible if strain is reversed -that also have a well defined elastic bulk modulus.
However, even the smallest rate used here does not succeed to separate events very close to jamming, since their relaxation times are still too long. Thus, even slower deformations have to be performed close to jamming in order to study those events.
Further work in progress concerns the understanding of the effect of system size, and of more realistic frictional particles. Open is still how the re-arrangement events look like and how to bring those local and global events (pos-sibly a superposition of many, or an avalanche from local events) into a macroscopic continuum theory [25]. The fact that different rates lead to different packings (with different jamming density) is also relevant for packing generation procedures, and a more quantitative study of the relaxation times for different densities above jamming is ongoing and will be published elsewhere.