---
abstract: 'Creep mechanisms in uniaxially compressed 3D granular solids comprised of faceted frictionless grains are studied numerically using a constant pressure and constant stress simulation method. Rapid uniaxial compression followed by slow dilation is predicted on the basis of a logarithmic creep phenomenon. Micromechanical analysis indicates the existence of a correlation between granular creep and grain-scale deformations. Localized regions of large strain appear during creep and grow in magnitude and size with time. Furthermore, the accumulation of non-affine granular displacements increases linearly with local strain, thereby providing insights into the origins of plastic dissipation during stress-driven creep evolution. The prediction of slow logarithmic dynamics in the absence of friction indicates a universality in the role of plastic dissipation during the creep of granular solids.'
author:
- Ishan Srivastava
- 'Timothy S. Fisher'
bibliography:
- 'main.bib'
title: Local Plasticity as the Source of Creep and Slow Dynamics in Granular Materials
---

Long-time evolution of granular materials under external stress is of significant importance in many natural and artificial processes. Time-dependent aging and creep in soils [@Schmertmann1991] and concrete [@Vandamme2009] under uniaxial stress are critically associated with the structural stability of man-made and natural geological constructs. The arrest (*jamming*) of dry cohesionless granular materials upon densification [@OHern2003; @Majmudar2007] and/or the application of external stress [@Bi2011; @Dagois-Bohy2012] has been numerically and experimentally investigated. Concomitantly, external stress can also initiate the yielding (*unjamming*) of granular packings [@DaCruz2005]. Critical state theories in soil mechanics [@Nedderman:1992aa] have proposed rate-independent, Mohr-Coulomb-type friction laws to describe the yielding of granular solids, and more recently, rate-dependent constitutive laws [@Jop2006] have improvised a description of dense granular flows above yield stress. However, a detailed description of the internal dynamics in granular materials below the yield stress is lacking [@Balmforth2014].

Diverse experiments at various length scales—from colloidal [@Coussot2002] to bio-cellular [@Bursac2005] to geological [@Houssais2015]—on kinematically constrained, glassy materials have demonstrated aging and creep relaxation at stresses below their yield strength. Experiments on granular samples have demonstrated slow aging with time and creep when subjected to a constant external stress [@Nguyen2011; @Amon2012] or a quasistatically increasing stress in stick-slip shearing experiments [@Hartley2003; @Nasuno1998]. Logarithmically slow densification has also been observed in various experiments on compaction in confined granular systems [@Brujic2005; @Richard2005; @Knight1995].

Several explanations and models have been proposed for mechanically induced aging behavior in granular systems. Hartley and Behringer [@Hartley2003] suggested that collective irreversible granular rearrangements are responsible for slow dynamics. The emergence of such collective irreversibility in amorphous materials is often explained through long-ranged elastic coupling of localized plastic events [@Maloney2006; @Schall2007]. Macroscopically, this complex spatio-temporal phenomenon is often condensed into a dynamical internal state variable called fluidity [@Kamrin2012], based on the theoretical formulations of soft glassy rheology [@Bocquet2009]. On the other hand, other explanations attribute slowly evolving dynamics to the gradual strengthening (or aging) of the individual solid-friction contacts with a slow increase in the local coefficient of friction [@Kuhn1993].

In this letter, we simulate the temporal evolution of the micromechanics of granular packings comprised of faceted grains that are confined at a constant pressure and subjected to a constant uniaxial compressive stress below the yield limit. The simulations reveal logarithmically slow creep after initial exponentially rapid dynamics. The influence of structural transformations on creep rheology is assessed by decomposing the local grain-scale motions into elastic and plastic contributions. Initial exponential dynamics are governed by elastic deformations, and plastic transformations dominate granular motion as the system transitions into a slow creep. The accumulated local plasticity in stressed granular packings—characterized by non-affine granular motions—is observed to be strongly correlated with a slow creep at the macroscale. We hypothesize that the concert of collective granular dynamics arising from a complex interplay of local elastic and plastic deformations is the source of slow creep dynamics and it is a hallmark of the general class of amorphous solids.

#### Simulation Details.

Initially, elastic and frictionless grains were jammed under a hydrostatic pressure $P$ through a recently introduced enthalpy-based, variable-cell jamming simulation method [@Smith2014]. The system is three-dimensional and contains $4096$ monodisperse, octahedral-shaped grains, as shown in the inset of Fig. \[global\](a). A very small ratio of the hydrostatic pressure and the material elastic modulus of the grains is prescribed here to represent the limit of hard granular systems [@DaCruz2005]. Repulsive pair potentials describe the elastic force $\mathbf{F}^{\text{e}}$ and moment $\mathbf{M}^{\text{e}}$ between contacting grains (for details on force and moment calculations, see [@Smith2010]).

To investigate the creep phenomenon in jammed granular systems under a constant external uniaxial stress $\bm{\sigma}_{\text{ext}}$, we introduce dynamics into the motion of grains and the periodic cell. The magnitude of the [*[xx]{}*]{} component of the stress tensor is $\sigma_{\text{ext}}$ and the other components are zero. The geometry of the periodic cell is described by a metric tensor $\bm{g}=\bm{h}^{T}\bm{h}$, where $\bm{h}$ is a matrix whose columns are the Bravais lattice vectors of the periodic cell. Translational motion of a grain $i$ is described in lattice coordinates $\bm{s}_{i}=\bm{h}^{-1}\bm{r}_{i}$ instead of Cartesian coordinates $\bm{r}_{i}$, and the rotation $\bm{\omega}_{i}$ is described relative to its center of mass. Dissipation during granular motion is modeled using a mean-field dissipation model [@Vagberg2014] described in detail in the accompanying supplementary material [^1].

Under low external stresses, the creep progresses at low strain rates, and the inertial effects of granular and periodic cell motion are assumed to be negligible in comparison to elastic interactions and dissipation. The equations of motion of grain displacement, grain rotation and periodic cell deformation are derived from a Lagrangian formulation [@Note1]. The equations of non-affine displacement and rotation of a grain $i$ are given as: $$D_{g}\dot{\bm{s}}_{i} = \widetilde{\mathbf{F}}_{i}^{e} \qquad\text{and}\qquad D_{g}\dot{\bm{\omega}}_{i} = \mathbf{M}_{i}^{e},
\label{eq1}$$

where $\widetilde{\mathbf{F}}_{i}^{\text{e}}$ is the elastic force in lattice coordinates [@Souza1997], $\mathbf{M}_{i}^{\text{e}}$ is the net moment on grain $i$ due to all the contacting grains, and $D_{g}$ is the dissipation constant for granular motion. The equation of motion of the periodic cell in the limit of negligible inertia is defined as: $$D_{c}\dot{\bm{g}} + \left(\widetilde{\bm{\sigma}}_{\text{app}} - \widetilde{\bm{\sigma}}_{\text{int}}\right) = 0.
\label{eq2}$$ Here, $\widetilde{\bm{\sigma}}_{\text{int}}$ is the internal stress tensor [@Smith2014] in lattice coordinates and $\bm{\sigma}_{\text{app}}=\bm{\sigma}_{\text{ext}} + P\mathbf{I}$ is the total applied stress, where $\mathbf{I}$ is the identity matrix. The applied stress tensor in lattice coordinates is denoted by $\widetilde{\bm{\sigma}}_{\text{app}}$ [@Note1]. The dissipation constant for periodic cell motion is denoted by $D_{c}$.

The equations of motion, Eqs. (\[eq1\]) and (\[eq2\]), were integrated using the Heuns second order predictor-corrector method. The units of time were chosen such that all dissipation constants are unity. All displacements and volumes are scaled by grain size and volume respectively, and the stresses are scaled by the material elastic modulus of the grains.

![\[global\] (a) Evolution of deviatoric strain $\gamma(t)$ for four different uniaxial compressive stresses. The solid black curves represent fits to Eq. (\[eq3\]) \[inset: a selected realization of jammed octahedral grains. Dark lines indicate periodic-cell boundaries\]. (b) Volumetric strain evolution $\epsilon_{v}(t)$ for four different uniaxial compressive stresses $\sigma_{\text{ext}}$ (see legend in (b) for color codes; arrow represents increasing stress) \[inset: evolution of compressive strain $\epsilon_{xx}$ along the axial direction and tensile strain $\left(\epsilon_{yy}+\epsilon_{zz}\right)$ along the transverse direction of the periodic cell boundaries relative to the direction of applied external stress\].](Fig1)

#### Macroscopic Evolution.

Uniaxial compressive stress was applied to the granular packing of $4096$ octahedral grains jammed at pressure $P=10^{-4}$ and the macroscopic strain tensor $\bm{\epsilon}$ was calculated from the deformation of the periodic cell boundaries [@Smith2014]. We restrict our analysis to the first and second invariants of the tensor: volumetric strain $\epsilon_{v}=\frac{1}{3}\sum_{k}\epsilon_{kk}$ and deviatoric (shear) strain $\gamma = \sqrt{\frac{3}{2}\sum_{i,j}\left(\epsilon_{ij}-\epsilon_{v}\delta_{ij}\right)}$, where $\delta_{ij}$ is the Kronecker delta tensor. The evolution of deviatoric strain with time is well described by two-stage dynamics represented phenomenologically by: $$\gamma(t)=\gamma_{0}\left(1-e^{-t/\tau_{1}}\right) + \gamma_{1}\text{log}\left(\frac{t}{\tau_{2}}+1\right),
\label{eq3}$$ where $\tau_{1}$ is a viscoelastic time constant, and $\tau_{2}$ sets the timescale for slow dynamics. Fig. \[global\](a) displays the evolution of deviatoric strain with time (normalized by the fast dynamics time constant $\tau_{1}$) at four different compressive stresses and the same jamming pressure. Black lines indicate best fits to Eq. \[eq3\], and good agreement is observed at both long and short times. The first term of Eq. \[eq3\] denotes an initial Kelvin viscoelasticity in which $\gamma_{0}$ is the total elastic strain that is delayed by a viscous time constant $\tau_{1}$ [@Findley:1976aa]. The second term corresponds to a logarithmic creep that is caused by long-term structural rearrangements and this phenomenon has been previously observed in colloidal glasses [@Siebenburger2012] and dense granular materials [@Nguyen2011; @Amon2012]. Such distinct relaxation regimes corresponding to fast and slow dynamics has been previously reported in experiments on granular compaction [@Knight1995; @Brujic2005].

The granular packings exhibit rapid compaction at short times, whereas slow dilation occurs at longer times (see Fig. \[global\](b)). Because the variable-cell method enables full tensorial evolution of strain while maintaining a constant external stress, the temporal response of transverse and axial strains at all times can be predicted (see inset of Fig. \[global\](b)). At short times the transverse tensile strain $\left(\epsilon_{yy}+\epsilon_{zz}\right)$ is negligible, and rapid granular compaction is caused entirely by axial compression $\epsilon_{xx}$. During slow creep the tensile transverse strain competes with the compressive axial strain and results in slow dilation within the packing. This prediction correlates well with experimental observations during displacement pile setup in dense sands—comprised of angular grains—in which an initial volumetric contraction was followed by a long-term kinematically restrained dilation [@Bowman2005].

![\[local\] (a) Deviations in the compressive strain $\epsilon_{xx}$ along a [*[xy]{}*]{}-plane at three times (left to right: $t=30\tau_{1}$, $60\tau_{1}$ and $100\tau_{1}$) relative to a reference configuration at $t=20\tau_{1}$ in the slow regime. White and green circles indicate growing regions of extensional and compressive strain respectively. (b) Probability distributions of local compressive strains $\epsilon_{xx}$. The solid lines are Gaussian fits \[inset: non-gaussian parameter $\alpha$ at various time intervals\]. (c) Spatial correlation function of local deviatoric strain $C_{\gamma}(r)$ for three time intervals (see (a)) as a function of radial distance $r$ normalized by the in-sphere diameter $d_p$ of an octahedral grain \[inset: normalized correlation length $\xi$ at various time intervals\].](Fig2)

#### Heterogeneous Micromechanics.

To understand the microscopic mechanisms that govern global creep response, the complex spatio-temporal granular dynamics were decomposed in terms of local structural rearrangements. A time-dependent local strain for grain $i$ was determined by the best-fit affine matrix $\mathbf{F}_{i}$ that transforms the center-to-center vectors $\mathbf{r}_{ij}$ of $n$ contacting grains over a time step of $\Delta t$ by minimizing $D_{i}^{2} = \sum_{j=1}^{n}\left[\mathbf{r}_{ij}\left(t+\Delta t\right) - \mathbf{F}\mathbf{r}_{ij}\left(t\right)\right]^2$ [@Falk1998] . A symmetric Lagrangian strain tensor can then be calculated as $\bm{\epsilon}_{i}=\frac{1}{2}\left(\mathbf{F}_{i}^{T}\mathbf{F}_{i}-\mathbf{I}\right)$. The remaining non-affine component of motion $D_{i}^{2}$ represents a component of the granular motion in the system that is elastically irreversible.

In the regime of slow granular creep, a growing heterogeneity of the local strain field is observed. With a reference configuration at $t=20\tau_{1}$, we focus on the deviations of the uniaxial strain component $\epsilon_{xx}$ (parallel to the global uniaxial strain) along a [*[xy]{}*]{}-plane of the packing and illustrate its distribution relative to the reference configuration for $t=30\tau_{1}$, $60\tau_{1}$ and $100\tau_{1}$ in Fig. \[local\](a). The grains are color-coded based on the local uniaxial strain. White and green circles indicate regions of tensile and compressive uniaxial strain, respectively. We observe that the regions of negative and positive strain extend over many grains, and the spatial extent of these regions grows with time. Furthermore, the regions of positive (negative) strain are more positively (negatively) strained with the progress of time, as indicated within the white (green) circles. These regions indicate zones of large strain localization that are analogous to shear transformation zones in other amorphous materials [@Schall2007; @Wang2013].

The distribution of local strains broadens with time, as shown in Fig. \[local\](b). Furthermore, the strain distribution becomes less Gaussian with time, as indicated by the non-Gaussian parameter $\alpha(\Delta t)$ in the inset of Fig. \[local\](b). This phenomenon is a signature of increasingly correlated transformations wherein localized deformations are coupled across multiple grains. To quantify this coupling, we compute the spatial autocorrelation of local deviatoric strain $\gamma$, defined as: $$C_{\gamma}(\Delta r) = \frac{\langle \gamma(r+\Delta r)\gamma(r)\rangle -\langle\gamma(r)\rangle^{2}}{\langle \gamma(r)^{2}\rangle - \langle \gamma(r)\rangle^{2}},
\label{eq4}$$ where angular brackets denote ensemble averages. As shown in Fig. \[local\](c), as time progresses in the slow creep regime, the local strains become more spatially correlated with an increasing correlation length as indicated in the figure inset. The correlation length is defined as the radial distance at which the autocorrelation goes to zero (denoted with dashed lines in the figure). This correlated deformation has been previously identified as a signature of the elastic response of glassy materials to local shear transformation zones [@Jensen2014] via stress redistributions caused by Eshelby fields [@Chattoraj2013].

![\[motion\] (a) Temporal evolution of the ensemble average of local deviatoric strain rate $\langle \dot{\gamma}_{i} \rangle$ (black), non-affine velocity $\langle \dot{D}_{i} \rangle$ (green) and rotational velocity $\langle \dot{R}_{i} \rangle$ (red) as a function of time. (b) Accumulated non-affine motion $S(t)$ as a function of global deviatoric strain $\gamma(t)$ for four uniaxial compressive stresses (see legend) and a jamming pressure $P=2\times 10^{-5}$. The solid black lines are linear fits. Inset: Scatter plots of the accumulated non-affine displacements divided by the applied uniaxial compressive stress for two different jamming pressures. Red circles represent data from the main graph. Blue rectangles represent data for $P=1 \times 10^{-4}$ and applied uniaxial stresses: $\sigma_{\text{ext}} = 1\times 10^{-7}$, $2\times 10^{-7}$, $3\times 10^{-7}$ and $4\times 10^{-7}$.](Fig3)

#### Evolution of Plastic Activity.

To elucidate the correlations between macroscale creep and local transformations, we quantify the temporal evolution of local dynamics at the grain scale by calculating the ensemble average of local deviatoric strain rate $\langle\dot{\gamma}_{i}\rangle$, non-affine velocity $\langle \dot{D}_{i} \rangle$ and rotational velocity $\langle \dot{R}_{i} \rangle$ as functions of time. Granular rotation is quantified by the magnitude of the rotation of a grain about its instantaneous axis of rotation that is given by the cross product of unit vectors that connect the grain centroid to a specific vertex between consecutive time steps.

Fig. \[motion\](a) shows that at short times, when the macroscopic dynamics are exponential ($t \ll \tau_{1}$), local mechanical transformations are caused almost entirely by affine elastic deformations in the packing, as indicated by large shear strain rates and low non-affine and rotational velocities. This result indicates that the early transformations are essentially elastic in nature with little irreversibility in the granular motion. However, as the granular packing transitions into the logarithmic slow creep regime ($t \sim \tau_{1}$), it experiences increased irreversibility as evidenced by a decrease in local shear strain rates, and an increase in the non-affine and rotational velocities. During long periods of slow creep ($t \gg \tau_{1}$), both the elastic and non-affine contributions to local deformations decay slowly as the system creeps towards a metastable state. The close association between the growth of local non-affine transformations and the evolution of slow macroscopic creep indicates the presence of a causal relationship.

This relationship was further explored by correlating the temporal accumulation of squared non-affine displacements with the temporal evolution of macroscopic strain. Quantitatively, this accumulation at any time $t$ is computed as $S(t) =  S(t-\Delta t)+ \langle D_{i}^{2}\left(t-\Delta t,t\right) \rangle$, where the angular brackets denote the ensemble average of $D_{i}^{2}$ recorded between times $t-\Delta t$ and $t$. Fig. \[motion\](b) displays the relationship between cumulative non-affine displacements and macroscopic creep strain. The relationship is linear (except at short times), and the slope decreases with decreasing magnitude of applied uniaxial stress. Furthermore, all the curves collapse to a single linear curve when $S(t)$ is divided by the applied stress $\tau$ (see inset for two jamming pressures). This relationship can be expressed as $S(t) \sim \sigma_{\text{ext}} \gamma(t)$, and it reinforces the evidence that local irreversible deformations play a dominant role in the global creep process. We also note that for granular systems, jamming pressure plays a crucial role in plastic dissipation. From the inset of Fig. \[motion\](b), the scale of plastic dissipation per unit elastic relaxation is larger in systems jammed under larger hydrostatic pressure. This result is expected because the large elastic modulus of highly jammed grains increases elastic coupling between plastic events.

Lastly, we provide a mechanical definition for $S(t)$ by considering $\sigma_{\text{ext}} \gamma(t)$ to represent the elastic energy density released in the bulk during the creep relaxation process at long times. In athermal systems, elastic energy relaxation should be balanced by plastic dissipation in the system. Therefore, $S(t)$ represents the scale of plastic dissipation density in the system that is driven by external stress and manifests itself in the form of non-affine granular displacements. External stress-facilitated accumulation of non-affine displacements was also previously demonstrated for thermally driven amorphous metallic glass systems [@Wang2013].

#### Conclusion.

This letter identifies the internal structural transformations that govern the creep relaxation of uniaxially compressed granular solids. A full tensorial analysis of the system dynamics, enabled by a variable-cell simulation method, provides constitutive prediction of the temporal evolution of deviatoric and volumetric strains at all scales. The simulations predict rapid uniaxial compression at short times followed by slow dilation at long times. The accumulation of irreversible granular motion is proposed as a micromechanical explanation for slow logarithmic creep and is directly correlated with the plastic dissipation density within the system. Importantly, a simplified mean-field dissipation model captures the slow logarithmic creep dynamics observed in real granular solids [@Nguyen2011] and glassy colloids [@Siebenburger2012], thereby indicating a universality in the role of localized plasticity during slow relaxation of amorphous solids. Future work will focus on the effects of more complex grain shapes and inter-granular friction on creep relaxation processes.

#### Acknowledgements.

This work was supported by Grant No. 1344654 from the National Science Foundation Scalable Nanomanufacturing Program. We also acknowledge useful discussions with Prof M. Ashraf Alam.

[^1]: See Supplemental Material for detailed derivation and dimensional analysis of the equations of motion.
