Similarities between GSH, Hypoplasticity and KCR

TL;DR

This paper compares GSH, a physics-based hydrodynamic theory, with hypoplasticity and KCR, revealing their similarities in modeling elasto-plastic behavior in granular media and establishing a solid theoretical foundation for these models.

Contribution

It demonstrates the deep connections between GSH, hypoplasticity, and KCR, providing a physics-based basis for constitutive relations used in granular media modeling.

Findings

GSH and hypoplasticity show great similarity in uniform, slow elasto-plastic motion.

Proportional paths and Goldscheider rule naturally emerge from GSH equations.

GSH provides a robust physical foundation for KCR and related models.

Abstract

Accounting for elasto-plastic motion in granular media, hypoplasticity is a state-of-the-art constitutive model derived from data accumulated over many decades. In contrast, GSH, a hydrodynamic theory, is derived from general principles of physics, with comparatively few inputs from experiments, yet sporting an applicability ranging from static stress distribution via elasto-plastic motion to fast dense flow, including non-uniform ones such as a shear band. Comparing both theories, we find great similarities for uniform, slow, elasto-plastic motion. We also find that proportional paths and the Goldscheider rule used to construct barodesy, another, more recent constitutive model, are natural results of GSH's equations. This is useful as it gives these constitutive relations a solid foundation in physics, and in reverse, GSH a robust connection to reality. The same symbiotic relation exists between GSH and KCR, or Kamrin's non-local constitutive relation, a model that was successfully employed to account for a wide shear band in split bottom cells.

Equations

$${\partial }_t{\sigma }_{ij}={\mathfrak{C}}_{ij}({\sigma }_{ij},v_{k\ell },\rho )$$

$${\partial }_t{\sigma }_{ij}=H_{ijk\ell }{v}_{k\ell }+{\Lambda }_{ij}\sqrt{v_{k\ell }{v}_{k\ell }+ϵv_{\ell \ell }^2},$$

$${\pi }_{ij}\equiv -\partial w/\partial u_{ij}{|}_{\rho }.$$

$${\partial }_t{T}_g=-R_T[T_g(1-{\xi }_T^2{\nabla }_i^2)T_g-(f{v}_s{)}^2-T_a^2],$$

$$T_g=\sqrt{f^2{v}_s^2+T_a^2},$$

$${\partial }_t{u}_{ij}-v_{ij}+{\alpha }_{ijk\ell }{v}_{k\ell }=-({\lambda }_{ijk\ell }{T}_g)u_{k\ell }.$$

$${\partial }_t{u}_{ij}\to ({\partial }_t+v_k{\nabla }_k)u_{ij}+u_{ik}{\nabla }_j{v}_k+u_{jk}{\nabla }_i{v}_k.$$

$$H_{mnk\ell }=(\partial {\pi }_{mn}/\partial u_{ij})({\delta }_{ik}{\delta }_{j\ell }-{\alpha }_{ijk\ell }),{\Lambda }_{mn}=(\partial {\pi }_{mn}/\partial u_{ij}){\lambda }_{ijk\ell }{u}_{k\ell }.$$

$${\partial }_t\Delta +(1-\alpha )v_{\ell \ell }-{\alpha }_1{u}_{ij}^{\ast }{v}_{ij}^{\ast }=-{\lambda }_1{T}_g\Delta ,$$

$${\partial }_t{u}_{ij}^{\ast }-(1-\alpha )v_{ij}^{\ast }=-\lambda T_g{u}_{ij}^{\ast },$$

$${\partial }_t{u}_{ij}^{\ast }$$

$${\partial }_t\Delta$$

$${\partial }_t(\rho v_i)+{\nabla }_j({\sigma }_{ij}+\rho v_i{v}_j)=g_i\rho ,{\partial }_t\rho =-{\nabla }_i(\rho v_i),$$

$${\sigma }_{ij}=(1-\alpha ){\pi }_{ij}-{\alpha }_1{u}_{ij}^{\ast }{P}_{\Delta }+[P_T{\delta }_{ij}-{\eta }_1{T}_g{v}_{ij}^{\ast }]\text{textdagger},P\equiv {\sigma }_{\ell \ell }/3,$$

$$\Delta \equiv -u_{\ell \ell },P_{\Delta }\equiv {\pi }_{\ell \ell }/3,u_s^2\equiv u_{ij}^{\ast }{u}_{ij}^{\ast },{\pi }_s^2\equiv {\pi }_{ij}^{\ast }{\pi }_{ij}^{\ast },$$

$$w_{\Delta }$$

$${\pi }_{ij}$$

$$P_{\Delta }$$

$$\frac{4P_{\Delta }}{{\pi }_s}$$

$$u_s/\Delta \le \sqrt{2\mathcal{B}/\mathcal{A}}\text{or}{\pi }_s/P_{\Delta }\le \sqrt{2\mathcal{A}/\mathcal{B}}.$$

$$\tan {\phi }_{st}={\pi }_s^{yield}/\sqrt{2}P=\sqrt{\mathcal{A}/\mathcal{B}}.$$

$$\mathcal{B}={\mathcal{B}}_0[(\rho -\bar{\rho })/({\rho }_{cp}-\rho ){]}^{0.15},$$

$${\sigma }_{ij}={\pi }_{ij},{\partial }_t{u}_{ij}=v_{ij}.$$

$$\alpha ={\alpha }_0{r}^{a_1},{\alpha }_1={\alpha }_{10}{r}^{a_2},\lambda /{\lambda }_0={\lambda }_1/{\lambda }_{10}=r^{a_3}.$$

$$T_g$$

$$u_{ij}^{\ast }$$

$$P_c=(1-\alpha )P_{\Delta }^c,{\sigma }_c=(1-\alpha ){\pi }_c,$$

$$P_{\Delta }^c\equiv P_{\Delta }({\Delta }_c,u_c),{\pi }_c\equiv {\pi }_s({\Delta }_c,u_c),$$

$$P_{\Delta }^c/{\pi }_c=(\mathcal{B}/2\mathcal{A}){\Delta }_c/u_c+u_c/4{\Delta }_c.$$

$$u_s(t)=u_c(1-e^{-\lambda f{\epsilon }_s}),{\epsilon }_s\equiv v_{s}t,$$

$$\Delta (t)={\Delta }_c(1+f_1{e}^{-\lambda f{\epsilon }_s}+f_2{e}^{-{\lambda }_{1}f{\epsilon }_s}),$$

$$f_1\equiv \frac{{\lambda }_1}{\lambda -{\lambda }_1},f_2\equiv \frac{{\Delta }_0}{{\Delta }_c}-\frac{\lambda }{\lambda -{\lambda }_1}.$$

$$\frac{{\bar{u}}_c^2}{u_c^2}=\frac{{\bar{\Delta }}_c}{{\Delta }_c}=\frac{{\bar{\sigma }}_c}{{\sigma }_c}=\frac{1}{1+T_a^2/(f{v}_s{)}^2}.$$

$${\partial }_t\Delta /v_s=-{\lambda }_{1}f(\Delta -\eta {\Delta }_c{u}_s/u_c),$$

$${\partial }_t{u}_s/v_s=-\lambda f(u_s-\eta u_c),$$

$${\partial }_t{\sigma }_{ij}^{\ast }$$

Full text

1. I Introduction

2. II The Equations of GSH

3. II.1 The state variables and their dynamics

4. II.2 The energy

5. II.3 Two Elastic Limits

6. II.4 More on $T_{g}$ and $u_{ij}$

7. II.4.1 The granular temperature $T_{g}$

8. II.4.2 The elastic strain $u_{ij}$

9. II.4.3 Rebuttal of “Twin-Granular-Temperature Theory”

10. II.5 Summary

11. III Experiments at Given Shear Rates

12. III.1 The Critical State

13. III.1.1 A Stationary Elastic Solution

14. III.1.2 Approach to the Critical State at Constant Density

15. III.1.3 Approach to the Critical State at Constant Pressure

16. III.1.4 Perturbations and Rate-Dependence of the Critical Stress

17. III.1.5 Load and Unload

18. III.2 Hypoplasticity, Neo-Hypoplasticity and Barodesy

19. III.2.1 A Comparison of the Hypoplastic Coefficients

20. III.2.2 Elastic Energy versus the Gibbs Potential

21. III.2.3 Comparing Four Hyperelastic Models

22. III.2.4 The Third Strain Invariant

23. III.2.5 Proportional Paths and Barodesy

24. IV Experiments at Given Shear Stress

25. IV.1 Uniform Creep

26. IV.2 Narrow Shear Bands

27. IV.3 Angle of Stability and Repose

28. IV.4 Wide Shear Bands and KCR

29. V The Quasi-Elastic Regime

30. V.1 Strain Increments versus Slow Rates

31. V.2 How to Observe Quasi-Elasticity

32. V.3 Transition to Quasi-Elasticity

33. VI Conclusion

**Notations:

**

[TABLE]

[TABLE]

I Introduction

Being a subject of practical importance, elasto-plastic deformation of dense granular media has been under the focus of engineering research for many decades if not a century schofield ; nedderman ; wood1990 ; kolymbas1 ; kolymbas2 ; gudehus2010 . The state of geotechnical theories, however, remains confusing for outsiders: Innumerable models compete, employing strikingly different expressions. In a recent book on soil mechanics, Gudehus employed metaphors such as morass of equations and jungle of data to describe the situation gudehus2010 .

Granular dynamics is frequently modeled employing the strategy of rational mechanics, by postulating a function $\mathfrak{C}_{ij}$, such that the constitutive relation,

[TABLE]

• holds, whereby $\mathfrak{C}_{ij}$ is a function of the Cauchy stress $\sigma_{ij}$, strain rate $v_{k\ell}$ (the symmetric part of the velocity gradient) and density $\rho$. (More generally, ${\partial}_{t}\equiv\frac{\partial}{\partial_{t}}$ is to be replaced by an objective derivative, such as Jauman – here and below.)

Together with mass and momentum conservation, it forms a closed set of equations for $\sigma_{ij},\rho$ and the velocity $v_{i}$. Hypoplasticity kolymbas1 ; kolymbas2 is such a theory, same as barodesy barodesy . Both use a single constitutive equation, without the recourse to yield or potential surfaces. Although elasto-plastic theories are somewhat different, Einav h^2 showed that the latter can be seen as a special case of the former. All these models yield realistic accounts of the complex elasto-plastic motion, providing us with vast amount of highly condensed and intelligently organized empirical data. We concentrate on comparing GSH to hypoplasticity here, because elasto-plastic theories, with its case separations, are not easy to deal with analytically, and because a comparison to barodesy already exists GSH&Barodesy .

Typically, the hypoplastic model starts from the rate-independent constitutive relation,

[TABLE]

with $H_{ijk\ell},\Lambda_{ij}$ and $\epsilon$ functions of stress and density. Though continually improving functional dependence is being proposed for them, some known and systematic shortcomings remain:

The applicability of hypoplasticity is confined to uniform systems. Lacking a characteristic length scale, it does not account for shear band or evanescent creep flow (ie. the exponentially decaying penetration of flows into the static region). There have been two main approaches to ameliorate this: introduction of gradient terms wu2 , or addition of variables to account for the Crosserat rotation and couple stress wu1 . The latter works well, but leads to a rather more complex theory. And it provokes questions about the underlying physics: If couple stress is indeed important in the shear band, because it is fluid, why then is it not important in the uniformly fluid and gaseous state of granular media – or, coming to think of it, in nematic liquid crystals deGennes ? 2. 2.

In hypoplasticity, rate-independence is an input. Yet it is not a robust feature of granular behavior. Fast dense flow is rate-dependent, as is the critical stress in the presence of external perturbations (such as given by tapping). Rate-independence needs to be explained, understood, and the extent to where it is valid clarified. 3. 3.

Elastic waves exist in granular media. Yet if one tries to account for them employing Eq.(2), one finds falsely that they should always be over-damped. More generally, granular media appear softer and more dissipative during elasto-plastic motion than in many other circumstances – foremost elastic waves, but also incremental stress-strain relations and ratcheting. Focusing on the former, hypoplasticity struggles to account for any of the latter. (See however works on intergranular strains Herle .) 4. 4.

Fast dense flow and collisional flow are not accounted for. 5. 5.

Finally, hypoplasticity lacks energetics. One does not know whether a given state is stable because it has an energy minimum, or instable because the energy turns concave. Nor do we know how much energy is being dissipated say in the critical state. Moreover, $H_{ijk\ell}$ has 36 elements, $\Lambda_{ij}$ six, all functions of a tensor and a scalar. These are a lot of functions to choose. So it would be helpful if the second derivative of the energy were related to $H_{ijk\ell}$ and/or $\Lambda_{ij}$ (contradicting in fact the name hypoplasticity).

Granular solid hydrodynamics (GSH) is derived from basic principles of physics, including conservation laws and the second law of thermodynamics, by employing the hydrodynamic procedure Khal ; LL6 . The result is a continuum-mechanical theory granL3 ; granR2 ; granR3 ; granRgudehus ; granR4 ; GGas ; granRexp that leaves relatively little leeway, with a handful of scalar parameters to fit experiments. The most important one is the energy $w$, chosen such, first of all, that elastic waves, incremental stress-strain relation and static stress distributions are well accounted for. Second, $w$ contains the information that no stable static stress distribution exists if the density $\rho$ is too small for enduring contacts, or if the shear stress is too large for the given pressure. Third, by reducing GSH to Eq.(2) for slow, uniform shear rates, $w$’s second derivative is found to yield $H_{ijk\ell}$ and $\Lambda_{ij}$.

The thermodynamic consistency of hypoplasticity has been demonstrated independently by Fang et al, see Fang1 ; Fang2 . This was done generally, without specifying an expression for $w$.

The energy $w$ is taken as a function of the density $\rho$, granular temperature $T_{g}$ and elastic strain $u_{ij}$. The state variable $u_{ij}$ accounts for *the slowly-varying, large-scaled elastic deformations of the grains, * with the elastic stress given by the hyper-elastic relation,

[TABLE]

The state variable $T_{g}$ quantifies the fast fluctuating elastic and kinetic energy of the grains. It is frequently proportional to the shear rate, $T_{g}\propto v_{s}\,\,(\equiv\dot{\gamma})$, and a useful indicator for the system’s behavior, from the quasi-static via the dense flow to the collisional one.

For collisional and fast dense flow, with a shear rate of $10^{3}$/s and more, there is vigorous jiggling and aggitation, so $T_{g}$ is obviously important. For elasto-plastic motion (usually referred to as “quasi-static” but see below), with shear rates around $10^{-3}$/s, $T_{g}$ may seem irrelevant: Aside from an occasional slip, grains only participate in the motion dictated by the macroscopic shear rate, with no perceptible deviations that may contribute to $T_{g}$. Nevertheless, such slips lead to vibrations in the region that imply a $T_{g}$ still many orders of magnitude larger than the true temperature of the grains. More importantly, as we shall see, these vibrations are what lead to plasticity.

If $T_{g}\propto v_{s}$ is yet smaller, below $10^{-5}$/s, GSH predicts a rate-dependent transition, to irrelevance of $T_{g}$ and diminishing plasticity. This is supported by recent experiments QuasiEla . Here, the system becomes truly quasi-static, better, quasi-elastic, until statics reigns, $T_{g}\equiv 0$.

In what follows, we shall briefly present GSH in Sec.II, first its three coupled relaxation equations and an expression for the elastic energy, followed by a discussion of how best to understand $T_{g}$ and $u_{ij}$, and why a second $T_{g}$ for the fluctuating elastic energy alone is counter-productive. In Sec.III, we obtain the critical state as the steady-state solution of the three relaxation equations. The approach to the state, for both high and low initial densities, is given by a simple analytic solution. In addition, GSH accounts for the variously observed fact that disturbing the critical state with some tapping, its stress becomes highly rate-dependent. We then compare GSH in detail to hypoplasticity and neo-hypoplasticity, and discuss how the observed proportional paths and Goldscheider rule are natural results of GSH’s equations. Finally, we compare four more recently proposed elastic energies, showing that all four appear to agree well, though their second derivatives differ close to yield.

In Sec.IV, we consider remifications of GSH if the stress is held constant, stressing the difference to rate-controlled experiments, explaining why a soft spring is useful for the latter, and a triaxial apparatus for the former. Then we calculate creep, shear band, angle of repose and stability (for the infinite, uniform geometry). Finally, we demonstrate the similarity of GSH to Kamrin’s nonlocal KCR model.

In Sec.V, we discuss the transition to the quasi-elastic regime for extremely slow rates.

II The Equations of GSH

II.1 The state variables and their dynamics

A complete set of state variables is one that uniquely determines a macroscopic state, including its energy and stress. GSH’s state variables are the density $\rho$, velocity $v_{i}$, the granular temperature $T_{g}$ and the elastic strain $u_{ij}$. The latter two quantify, respectively, the temporally fluctuating elastic and kinetic energy of the grains and their slowly-varying large-scaled deformations. $T_{g}$ is increased by viscous heating $(fv_{s})^{2}$, decreased by relaxation, and its balance equation is

[TABLE]

with $R_{T}T_{g}$ the relaxation rate, and $R_{T}\xi_{T}^{2}\nabla^{2}_{i}T_{g}$ accounting for $T_{g}$-diffusion. (The strain rate is $v_{ij}\equiv\frac{1}{2}(\nabla_{i}v_{j}+\nabla_{j}v_{i})$, its traceless part $v^{*}_{ij}$, and $v_{s}\equiv\sqrt{v^{*}_{ij}v^{*}_{ij}}>0$.) $T_{a}$ is the “ambient temperature” maintained by tapping or a sound field. For a steady, uniform $v_{s}$, we have

[TABLE]

implying $T_{g}=fv_{s}$ if there is no external perturbation, and $T_{g}=T_{a}$ if the shear rate is zero. As is easy to see from Eqs.(6,8) below, it is $T_{g}=fv_{s}$ that leads to rate-independence. The easiest way to destroy it is therefore to introduce some $T_{a}$, see Sec.III.1.4.

The importance of $T_{g}$ in the collisional state is easy to accept. Recent progresses show it exists equally in the turbulent flows of rapid dry granular avalanches Fang3 . That $T_{g}$ is relevant even in elasto-plastic motion may surprise, yet as already discussed in the introduction, this is a result of $T_{g}$, many orders of magnitude larger than the true temperature, being the physics underlying plasticity, see the next few paragraphs.

The elastic strain $u_{ij}$ alone suffices to account for elasticity and plasticity. To understand this, one needs two steps: First, granular stresses relax with a variable relaxation time. Second, this relaxation is closely related to plasticity.

The free surface of a granular system at rest is frequently tilted. When perturbed, when the grains jiggle and $T_{g}\not=0$, the tilted surface will decay and become horizontal. The stronger the grains jiggle and slide, the faster the decay is. We take this as indicative of a system that is elastic for $T_{g}=0$, transiently elastic for $T_{g}\not=0$, with a stress relaxation rate $\propto T_{g}$. A relaxing stress is typical of any viscous-elastic or transiently elastic systems, such as polymers polymer-1 ; polymer-2 ; polymer-3 ; polymer-4 . In granular media, the relaxation rate is not a material constant, but a function of the state variable $T_{g}$ – a behavior that we name variable transient elasticity.

For given $\rho$, the stress $\pi_{ij}\equiv-\partial w/\partial u_{ij}$ is a monotonic function of $u_{ij}$ (as long as the energy $w$ is stable and convex, $\partial^{2}w/\partial u_{ij}\partial u_{k\ell}>0$). Therefore, there is no principle difference in considering either relaxation. Yet since strain is a geometric quantity, stress a physical one that includes material parameters such as the stiffness (stress dependent and highly anisotropic in granular media), considering strain relaxation is simpler.

We take the equation of motion for $u_{ij}$ such to reflect variable transient elasticity, finding that it suffices to capture elasto-plasticity. Starting from pure elasticity, $\partial_{t}u_{ij}-v_{ij}=0$, adding a relaxation term, $\partial_{t}u_{ij}-v_{ij}=-u_{ij}/\tau$, we allow for two tensorial coefficients, $\lambda_{ijk\ell}$ and $\alpha_{ijk\ell}$, to parameterize the efficacy of relaxation and deformation rate,

[TABLE]

If stress anisotropy changes appreciably, convective terms (or the so-called objective time derivative) become important, then one needs the substitution (see the derivation in polymer-1 )

[TABLE]

Since $\partial_{t}\pi_{mn}=(\partial\pi_{mn}/\partial u_{ij})\partial_{t}u_{ij}=-(\partial^{2}w/\partial u_{mn}\partial u_{ij})\partial_{t}u_{ij}$, and $T_{g}=fv_{s}$, see Eq.(5), this simple equation possesses the full hypoplastic structure, with

[TABLE]

Given the equality of structure, GSH is as capable of accounting for elastoplasticity as hypoplasticity, if an appropriate energy expression can be found. (Allowing the density to vary yields slightly more complicated expressions.) Next, dividing $u_{ij}$ into its trace $\Delta\equiv-u_{\ell\ell}$ and traceless part $u_{ij}^{*}$, denoting $u_{s}\equiv\sqrt{u^{*}_{ij}u^{*}_{ij}}>0,\quad\hat{u}_{ij}\equiv u_{ij}^{*}/u_{s},\quad\hat{v}_{ij}\equiv v_{ij}^{*}/v_{s}$, and specifying the two tensors with two elements each, $\alpha,\alpha_{1},\lambda,\lambda_{1}$, we arrive at

[TABLE]

which may equivalently be rewritten as two relaxation equations,

[TABLE]

Eqs.(4,11,12) are the three relaxation equations of GSH. They contain three inverse time scales: $R_{T}T_{g},\lambda T_{g},\lambda_{1}T_{g}$, a length scale $\xi_{T}$, and the parameters $f,\alpha,\alpha_{1}$. In dense media, $R_{T}T_{g}$ is large, of order $10^{3}$/s or more. With $T_{g}=fv_{s}$, the scales $\lambda T_{g},\lambda_{1}T_{g}$ are, for the shear rates typical of soil-mechanical experiments, orders of magnitude smaller, around 1/s for $v_{s}=10^{-2}$/s. The length scale $\xi_{T}$ is a few to a few tens granular diameters in dense media.

The steady state solution: $T_{g}=fv_{s}$, $u_{s}=u_{c}$ and $\Delta=\Delta_{c}$, implying constant elastic strain and stress at given shear rate, is clearly the critical state, see Sec. III.1 below.

Conservation of momentum and mass,

[TABLE]

close the set of equations, where the Cauchy stress $\sigma_{ij}$ is

[TABLE]

see granR2 for derivation. The (off-diagonal Onsager) coefficients $\alpha,\alpha_{1}$, a necessary consequence of the same coefficients in Eqs.(9,10), soften and mix the stress components. The term preceded by $\alpha_{1}$ is smaller by one order in the small elastic strain $u_{ij}$ and may usually be neglected (though not in Sec. III.2.1). $\alpha,\alpha_{1}$ are functions of the density for elasto-plastic motion, though they vanish identically for static stresses, see Eq.(77). The pressure $P_{T}\propto T_{g}^{2}$ is exerted by jiggling grains, and the last term is the viscous stress. Both are $\propto v_{s}^{2}$ for $T_{g}\propto v_{s}$ and relevant only for fast collisional and dense flows, giving rise to Bagnold scaling there.

II.2 The energy

An explicit expression for the stress via Eqs.(3,14), is delivered by the energy $w$ that we take to be $w=w_{\Delta}(\Delta,u_{s},\rho)+w_{T}(T_{g},\rho)$. The part relevant for elasto-plastic flow is $w_{\Delta}$. ($w_{T}$, important for collisional and dense flow, is neglected here.) Denoting

[TABLE]

with $u^{*}_{ij},\pi^{*}_{ij}$ the respective traceless tensors, we take $w_{\Delta}$ as

[TABLE]

Given by the energy of linear elasticity multiplied by $\sqrt{\Delta}$, this form is clearly inspired by the Hertzian contact, though its connection to granular elasticity goes far beyond that. It includes both stress-induced anisotropy and the convexity transition (see below). In addition, the expression for the elastic stress $\pi_{ij}$ has been validated for:

• •

Static stress distribution in three classic geometries: silo, sand pile, point load on a granular sheet, calculated employing $\nabla_{i}\pi_{ij}=\rho g_{i}$, see ge1 ; ge2 ; granR1 .

• •

Incremental stress-strain relation at varying static stresses KJ ; AH ; incre .

• •

Propagation of anisotropic elastic waves at varying static stresses jia2009 ; ge4 .

The GSH framework is general enough that one may employ any energy expression one deems appropriate. In fact, the above expression has been chosen for simplicity and manipulative ease, with a more realistic one given in 3inv . In Sec. III.2, we shall show and compare a number of different elastic energy expressions for realistic modeling.

Stress-induced anisotropy: In linear elasticity, $w\propto u_{s}^{2}$, the velocity of elastic waves $\propto\sqrt{\partial^{2}w/\partial u_{s}^{2}}\,$ does not depend on $u_{s}$, or equivalently, the stress. For any exponent other than two, the velocity depends on the stress, and is anisotropic if the stress is.

Convexity Transition: There are values of stress and density in granular media for which elastic solutions do not exist: First, for densities smaller than the random loose packing density, $\rho<\rho_{\ell p}$, at which the grains start to loose contact with one another, no static elastic states can be maintained. Second, static elastic states collapse when the shear stress is too large for the given pressure. On a tilted plane, the angle at which this happens is frequently referred to as the angle of stability, $\varphi_{st}$. The actually measured angle varies with the spacial dependence of the stress components and does not have a unique value. But in an infinite geometry, for uniform stresses, elastic solutions fail at a well-defined value.

Note that $\varphi_{st}$, accounting for collapses of static stresses at $T_{g}=0$, is unrelated to the critical state, which exists only in the presence of a finite granular agitation, $T_{g}=fv_{s}\not=0$. In Sec. IV.3, we shall identify the angle of repose $\varphi_{re}$ with the critical angle. Again, the measured angle may vary because the stress is not uniform.

In the space spanned by the density and stress components, the surface where the second derivative of the elastic energy changes its sign (turning from convex to concave) is a divide: On one side elastic solutions are stable, with grains at rest, on the other elastic solutions are instable, infinitesimal perturbations suffice to destroy any solution, and the grains are always agitated. The elastic energy of Eq (16) is convex only for

[TABLE]

It turns concave when this condition is violated. (The second condition may be derived by considering Eq (19), showing $P_{\Delta}/\pi_{s}=\sqrt{{\cal B}/2{\cal A}}$ is minimal for $u_{s}/\Delta=\sqrt{2{\cal B}/{\cal A}}$.)

On a plane inclined by the angle $\varphi$, with $y$ the depth of the granular layer on the plane and $x$ directed along the slope, we take the stress to be $\pi_{xx}=\pi_{yy}=\pi_{zz}=P_{\Delta}$, $\pi_{xy}=\pi_{s}/\sqrt{2}$, $\pi_{yz},\pi_{xz}=0$. Integrating $\nabla_{j}\pi_{ij}=g_{i}\rho$ assuming a variation only along $y$, we find $\pi_{xy}=g\sin\varphi\int\rho(y)dy$ and $\pi_{yy}=\pi_{xy}/\tan\varphi$. The angle of stability $\varphi_{st}$ (for infinite geometry) is reached when the energetic instability of Eqs.(20) is breached. With $\pi_{s}^{yield}\equiv P\sqrt{2{\cal A}/{\cal B}}$ denoting the yield shear stress, it is

[TABLE]

Modifications from the proximity to walls or floor are not included here.

For the parameterization of the coefficients ${\cal A},{\cal B}$, we need a second characteristic density, the random-close one, $\rho_{cp}$. It is the highest density at which grains may remain uncompressed. Assuming ${\cal B}/{\cal A}$ is density-independent (typically 5/3), and denoting $\bar{\rho}\equiv(20\rho_{\ell p}-11\rho_{cp})/9$, we take

[TABLE]

with ${\cal B}_{0}>0$ a constant. This expression accounts for three granular characteristics:

• •

The elastic energy is concave for $\rho<\rho_{\ell p}$, rendering all elastic solutions instable.

• •

The energy is convex between $\rho_{\ell p}$ and $\rho_{cp}$, ensuring the stability of elastic solutions in this region. In addition, the density dependence of sound velocities (as measured by Harding and Richart hardin ) is well rendered by an expression $\propto\sqrt{{\cal B}(\rho)}$. (The usual formula used to report the resuslts of hardin implies an unphysically concave energy.)

• •

The elastic energy diverges, slowly, for $\rho\to\rho_{cp}$, approximating the observation that the system becomes quite a bit stiffer there.

II.3 Two Elastic Limits

Because there is no space for granular rearrangements at the closed-packed density $\rho_{cp}$, plastic motions cannot take place and the system should be fully elastic. Taking $\alpha,\alpha_{1},\lambda,\lambda_{1}\to 0$ for $\rho\to\rho_{cp}$ in Eqs.(9,10,14), reduces GSH to simple elasticity,

[TABLE]

Close to $\rho_{cp}$, the dependence on $\rho_{cp}-\rho$ is more sensitive than any on $\rho$ directly. Hence we take, with $a_{1},a_{2},a_{3},\lambda_{0},\lambda_{10}>0$ and $r\equiv 1-\rho/\rho_{cp}$:

[TABLE]

Granular media should also be fully elastic for $T_{g}\to 0$, see the discussion at the end of the introduction, and in Sec.V. In this limit, the stress increases steeply with the shear rate, until the yield stress is reached, at which point the system collapses. There is no stress relaxation, hence neither a critical state.

II.4 More on $T_{g}$ and $u_{ij}$

As $T_{g}$ and $u_{ij}$ are novel variables that entail conceptual subtleties, a discussion is useful for those who want to look beyond the present equations. It is also useful because Sun Qichen and coworkers Sun have variously proposed a “Twin-Granular-Temperature Theory” that we find necessary to comment on.

II.4.1 The granular temperature $T_{g}$

In a uniform medium (such as water or copper) there are two length scales, macro- and microscopic. All degrees of freedom may be divided into one of the two groups. A hydrodynamic theory takes the degrees of freedom from the first group as explicit variables, including mass, momentum and energy density, also the strain field in the case of solids. All degrees of freedom from the second group are taken summarily, with their contribution to the energy lumped together as heat, and characterized by the temperature $T$.

In granular media, there is an intermediate group, consisting of momentum and deformation of each grain. In the kinetic theory or DEM, these are taken as independent variables. But a summary inclusion again suffices here, with the associated energy lumped into granular heat, characterized by $T_{g}$. As the kinetic energy changes rapidly into an elastic one, and back, during a collision, both types of energy interact, equilibrate, and must be lumped into one single macroscopic variable $T_{g}$. On the other hand, distinguishing $T_{g}$ from $T$ is necessary: $T_{g}$ is many orders of magnitude higher than $T$ in elasto-plastic motion, and when this is the case, granular stress relaxes, giving rise, as discussed, to plasticity.

There are also grain-sized deformations, such as in force chains, that are static or slowly varying in time. Their contribution is not included in $T_{g}$, but in that of $u_{ij}$, see below.

II.4.2 The elastic strain $u_{ij}$

Accounting for the macroscopic elastic stress that varies slowly in space, the elastic strain is also slowly varying. As energy is positive and conserved, the macroscopic energy density $w$ is the sum of all (slowly varying) mesoscopic contributions, $w^{mes}$. In other words, $w$ is the coarse-grained mesoscopic one, $w=\langle w^{mes}\rangle$, and since the latter includes the contribution from the force chains, the former also does.

It is important to realize that the elastic strain cannot be obtained by coarse-graining, $u_{ij}\not=\langle u_{ij}^{mes}\rangle$. This is because both the energy and the elastic stress are coarse-grained quantities: $w=\langle w^{mes}\rangle$ and $\pi_{ij}=\langle\pi_{ij}^{mes}\rangle$. With ${\rm d}w=\langle{\rm d}w^{mes}\rangle=-\langle\pi_{ij}^{mes}{\rm d}u_{ij}^{mes}\rangle=-\langle\pi_{ij}^{mes}\rangle{\rm d}u_{ij}$ and $\langle\pi_{ij}^{mes}{\rm d}u_{ij}^{mes}\rangle\not=\langle\pi_{ij}^{mes}\rangle{\rm d}\langle u_{ij}^{mes}\rangle$, we conclude: $u_{ij}\not=\langle u_{ij}^{mes}\rangle$.

II.4.3 Rebuttal of “Twin-Granular-Temperature Theory”

Sun and coworkers have written a number of papers Sun proposing a generalization of GSH by introducing two granular temperatures, $T^{k}$ and $T^{c}$, with $T^{k}$ taking account of the kinetic, and $T^{c}$ the elastic energy (including contributions from force chains). It should be clear from the last two sections that neither is a well-behaved macroscopic variable. The temporally fast varying elastic energy from collisional deformation is necessarily part of $T_{g}$, while the slowly varying contribution from force chains are part of $u_{ij}$.

II.5 Summary

With the differential equations given and the energy density specified, GSH is a fairly well-defined theory. It contains clear ramifications and provides little leeway for retrospective adaptation to observations. As we shall see in the rest of the paper, a wide range of granular phenomena is well encoded in these equations.

The three relaxation equations, (4,11,12), possess a steady state solution: $T_{g}=fv_{s}$, $u_{s}=u_{c}$, $\Delta=\Delta_{c}$, which is simply the critical state. With it given, we may consider the relaxation dynamics – either at given shear rate, or holding the shear stress constant. In the first case, because of the small relaxation time $1/R_{T}T_{g}$, we have $T_{g}=fv_{s}$ practically instantly, rendering $T_{g}$ dependent. Then the strain relaxation equation is manifestly rate-independent, displaying the mathematical structure of hypoplasticity. The rich phenomenology observed in triaxial apparatus is then well accounted for by hypoplasticity and GSH alike. The relaxation dynamics at constant shear stresses gives rise to phenomena such as creeping, shear band, angle of repose – all beyond hypoplasticity. They are considered in Sec.IV.

III Experiments at Given Shear Rates

III.1 The Critical State

The physics of the critical state is: Deformed grains give rise to an elastic stress. When grains are sheared past one another, they rattle and jiggle. And when they do, the deformation and stress is slowly lost, because of the loosen contact to their neighbors. As a consequence, a shear rate not only increases the deformation, as in any elastic medium, but also decreases it. A steady state exists in which both processes balance, such that the elastic deformation and stress remain constant over time. This is the critical state. Moreover, the critical stress is rate-independent, because (1) the increase in deformation is $\propto v_{s}$, (2) the relaxation is $\propto T_{g}$, and (3) $T_{g}\propto v_{s}$ in $T_{g}$’s steady state.

III.1.1 A Stationary Elastic Solution

Taking $\partial_{t}T_{g},\,\partial_{t}u_{ij}^{*},\,\partial_{t}\Delta=0$ in Eqs.(4,11,12) for $v_{ij}^{*}=$ const, $v_{\ell\ell}=0$, we have

[TABLE]

With $\Delta=\Delta_{c}$, $u_{ij}^{*}=u_{c}\hat{v}_{ij}$ and $\rho$ being constant, so is the associated elastic stress $\pi_{ij}(u_{k\ell},\rho)$. The stress value is obtained by inserting Eqs.(25,26) into Eqs.(17,18,19),

[TABLE]

The loci of the critical states thus calculated GSH&Barodesy (though employing the slightly more general energy of 3inv ) greatly resembles those calculated using either hypoplasticity or barodesy barodesy

Same as the Coulomb yield of Eq (21), the critical ratio $\sigma_{c}/P_{c}=\pi_{c}/P_{\Delta}$ is also frequently associated with a friction angle. Since the former is relevant for vanishing $T_{g}$, while the latter requires an moderately elevated $T_{g}\propto v_{s}$, it is appropriate to identify one as the static friction angle, or yield angle, and the other as the dynamic one, or the critical angle. In Sec. IV.3, we shall present reasons why the critical angle should be identified with the angle of repose $\varphi_{re}$, the angle at which granular flows will come to a halt for uniform stress.

Textbooks on soil mechanics frequently state that the friction angle is independent of the density – although they do not usually distinguish between the dynamic and the static one. We shall assume, for lack of better information, that both are.

III.1.2 Approach to the Critical State at Constant Density

Solving Eqs.(4,11,12) for $u_{s}(t),\Delta(t)$, at constant $\rho,v_{s}$, with the initial conditions: $\Delta=\Delta_{0},u_{s}=0$ (same as $P=P_{0},\sigma_{s}=0$), the relaxation into the critical state is given as

[TABLE]

This is an exponential decay for $u_{s}$, and a sum of two decays for $\Delta$. It is useful, and quite demystifying, that a simple, analytical solution exists.

Typically, we take $\lambda/\lambda_{1}\approx 3$, as the decay of $u_{s}$ and $f_{1}$ are faster than that of $f_{2}$. The prefactor $f_{2}$ may be negative, and $\Delta(t)$ is then non-monotonic. The associated pressure and shear stress are those of Eqs (27,28,29). For a negative $f_{2}$, neither the pressure nor the shear stress is monotonic. For the system to complete the approach to the critical state, the yield surface of Eq.(20) must not be breached during the non-monotonic path.

III.1.3 Approach to the Critical State at Constant Pressure

Frequently, the critical state is approached at constant pressure $P$, or a stress eigenvalue $\sigma_{i}$. The circumstances are then more complicated. As $\Delta,u_{s}$ approach $\Delta_{c},u_{c}$, the density compensates to keep $P(\rho,\Delta,u_{s})=$ const. Along with $\rho$, the coefficients $\alpha,\alpha_{1},\lambda,\lambda_{1},f,\Delta_{c},u_{s}$ (all functions of $\rho$), also change with time. In addition, with $\rho$ changing, the compressional flow $v_{\ell\ell}=-\partial_{t}\rho/\rho$ no longer vanishes (though it is still small). Analytic solutions do not seem feasible, but numerical ones are, see Fig 1. The parameters, labeled as I and II, are:

• •

${\cal B}_{0}=2,0.22$ GPa, ${\cal B/A}=5/3,8$,  $\bar{\rho}/\rho_{cp}=0.615,0.650$,

• •

$\alpha_{0}=1.04,0.85$, $\alpha_{10}=400,30$, $\lambda_{0}f=272,250$, $\lambda/\lambda_{1}=3.8,3.8$,

• •

$a_{1}=0.15,0.15$,   $a_{2}=1,0.15$,   $a_{3}=0.6,0.53$.

III.1.4 Perturbations and Rate-Dependence of the Critical Stress

If one perturbs the system, say by exposing it to weak vibrations, or by tapping it periodically, such as in a recent experiment vHecke2011 , the critical state is modified, and a strong rate-dependence of the critical shear stress is observed. The stress decreases with the tapping amplitude, and increases with the shear rate, such that the decrease is compensated at high rates. Clearly, any theory with built-in rate-independence cannot account for this observation. In GSH, on the other hand, rate-independence is a result of $fv_{s}/T_{g}=1$, and expected to be destroyed by any ambient temperature $T_{a}\neq 0$, see Eq.(5). The steady state values are then reduced to $\bar{u}_{c}\equiv(fv_{s}/T_{g})u_{c}$, $\bar{\Delta}_{c}\equiv(fv_{s}/T_{g})^{2}\Delta_{c}$, $\bar{\sigma}_{c}/\sigma_{s}={\bar{\Delta}_{c}^{0.5}}\bar{u}_{c}/{\Delta_{c}^{0.5}}u_{c}=(fv_{s}/T_{g})^{2}$, see Eqs.(11,12), implying a rate-dependence via Eq.(5),

[TABLE]

If there is no tapping, $T_{a}=0$, we retrieve the unperturbed values, $\bar{u}_{c}=u_{c}$, $\bar{\Delta}_{c}=\Delta_{c}$, $\bar{\sigma}_{c}=\sigma_{c}$. With tapping, $\bar{u}_{c},\bar{\Delta}_{c},\bar{\sigma}_{c}$ decrease for increasing $T_{a}$, and increase with increasing shear rate $fv_{s}\equiv f|v_{s}|$. see Fig 2.

(Note we have only considered the critical state at given shear rate, not the approach to it. So the result holds both at given density and pressure.)

The above consideration is the basic physics of the observation of vHecke2011 . It helps to put rate-independence, frequently deemed a fundamental property of granular media, into the proper context. A quantitative comparison is unfortunately made difficult by the highly nonuniform experimental geometry. Nevertheless, some comparison, even if unabashedly qualitative, may still be useful. In vHecke2011 , the torque $\tau$ on the disk on top of a split-bottom shear cell is observed as a function of its rotation velocity $\Omega$ and the shaking acceleration $\Gamma$. Now, $(\tau,\sigma_{s})$, $(\Omega,v_{s})$, $(\Gamma,T_{a})$ are clearly related pairs. Assuming the lowest order terms suffice in an expansion, we take $\sigma_{c}\propto\tau$, $\Gamma\propto T_{a}$. If $v_{s}$ were uniform, $\Omega\propto v_{s}$ would also hold. Since it is not, $\Omega\propto v_{s}^{n}$, $n>1$ seems plausible, because with additional degrees of freedom such as position and width of the shear band, the system has (for given $\Omega$) more possibilities to decrease its strain rate $v_{s}$. We take $\Omega=c_{2}v_{s}^{3}$ with $c_{2}=1{\rm rs}^{2}$ [replacing $T_{a}/v_{s}$ with $\Gamma/\sqrt[3]{\Omega}$ in Eq.(31)] to fit Fig 2, emphasizing that qualitative agreement exists irrespective of the fit.

III.1.5 Load and Unload

Inserting $T_{g}=fv_{s}$, $v_{\ell\ell}=0$ into Eqs (11,12), keeping the shear rate along one direction and allowing reversal, $\eta\equiv\hat{v}_{ij}\hat{u}_{ij}=\pm 1$, we have

[TABLE]

and conclude that the reason for the difference in apparent stiffness: $\partial_{t}u_{s}/v_{s}$ and $\partial_{t}\Delta/v_{s}$, between load and unload, is that $\eta$ turns negative when the shear rate $v_{ij}^{*}$ is reversed. The strain then relax towards the new stationary values, with a simple relaxation dynamics not at all history-dependent. That loading ($\eta=1$) and unloading ($\eta=-1$) have different slopes is frequently referred to as incremental nonlinearity in soil mechanics. It is the reason why no backtracing takes place under reversal of shear rate. The stress components $P,\sigma_{s}$ are calculated employing Eqs (18) for given $\Delta,u_{s}$ and $\rho$.

III.2 Hypoplasticity, Neo-Hypoplasticity and Barodesy

As discussed around Eqs.(6,8), GSK possesses the hypoplastic structure for given shear rates, $T_{g}=fv_{s}$. So, this is a lucky instance in which a comparison may be made directly, by comparing the coefficients – where one set is proposed from accumulated knowledge of granular behavior, while the other is given by the second derivative of the elastic energy. This is done in Sec III.2.1. Then we show that all recent hyperelastic models, either given as an elastic energy or a Legendre potential, are rather similar if properly transformed. The main difference in fact stems from whether or not there are elastic instabilities.

Barodesy, a recent model again proposed by Kolymbas barodesy , has a rather more complicated rate dependence than the one shared by hypoplasticity and GSH. But our comparison of GSH and barodesy did yield essentially quantitative agreement GSH&Barodesy ; P&G2009 . The idea of barodesy is to have a more modular and better organized structure than hypoplasticity, with different parts in the constitutive relation taking care of specific aspects of granular deformation, especially that of proportional paths as summed up by the Goldscheider rule. In Sec. III.2.5, we shall focus on the latter showing how the Goldscheider rule and proportional paths naturally arise from GSH.

III.2.1 A Comparison of the Hypoplastic Coefficients

Writing the hypoplastic model as

[TABLE]

we have the dimensionless coefficients $H_{1-7}$ that are functions of $\rho$ and the 3 stress invariants: $P\equiv\sigma_{nn}/3$, $\sigma_{s}\equiv\sqrt{\sigma_{ij}^{*}\sigma_{ij}^{*}}$, $\sigma_{t}\equiv\sqrt[3]{\sigma_{ij}^{*}\sigma_{jk}^{*}\sigma_{ki}^{*}}$. In Nie3 , the authors present the parameter values for more than forty different granular media, with

[TABLE]

characterized by three functions $F,f_{0},f_{d}$ and one constant $a$, given as

[TABLE]

where

[TABLE]

contain eight material constants, with $\varphi_{c}$ the critical frictional angle, and $e_{0d},e_{0c},e_{0i}$ denoting the densest, critical and loosest void ratio. For the Stuttgart sand, these are given as

[TABLE]

Note $\sigma_{ij}$ is here negative from that in Nie3 ; note also that the stress invariant $\sigma_{t}$ in $\cos 3\theta$ and $F$ will be approximated below with $\cos 3\theta=1$, ie. $\sigma_{t}^{3}=\sigma_{s}^{3}/\sqrt{6}$, valid for axial compressions.

Returning to GSH, see Eqs.(11,12) with $T_{g}=fv_{s}$, we have four dynamic coefficients,

[TABLE]

and two static coefficient $\cal B$ and ${\cal B/A}=5/3$.

Defining an operator that is valid for any function $f(\Delta,u_{s},\rho)$ such as the stress $\sigma_{ij}$,

[TABLE]

we rewrite

[TABLE]

in the form of Eq.(34) to obtain

[TABLE]

We first parameterize the static coefficients. Using the expression for the energy, Eqs.(16,22), we take ${{\cal A/B}=5/3}$ and ${{\cal B}_{0}=0.4h_{s}}$, with $h_{s}$, $e_{0i}=\rho_{G}/\rho_{lp}-1$, $e_{0d}=\rho_{G}/\rho_{cp}-1$ given by table 38. (The bulk density $\rho_{G}$ does not enter the formula.) Next, we take the dynamic coefficients to be ${\lambda f=600r^{0.5}}$, where $r\equiv(\rho_{cp}-\rho)/\rho_{cp}$; ${\lambda_{1}f=0.3\cdot\lambda f}$; and ${\alpha=(\rho_{lp}/\rho)^{2}(1-\rho_{lp}/\rho_{cp})^{-2}\,(1-r)\,r^{2}}$. Finally, we take ${\alpha_{1}=0.103\lambda f}$ (because with $\alpha_{1}/\lambda_{1}f=\Delta_{c}/u_{c}$, we obtain $\pi_{c}/P_{\Delta}^{c}=1,67\tan\varphi_{c}$, in agreement with $\varphi$ as given in table 38.) We note that the time scale $R_{T}$ and length scale $\xi_{T}$ of Eq.(4) are not contained in the hypoplastic model, and their values therefore not needed for the comparison.

Although qualitative agreement is obvious from Fig 4, a quantitative one is lacking. We do not take this as a serious defect of GSH, for three reasons. First, complete agreement of coefficients is probably an unnecessarily difficult task: Different coefficients frequently yield similar experimental curves. Second and more importantly, the parameterization of GSH, especially its energy, has been chosen stressing manipulative ease and simplicity of expressions. GSH’s robust and qualitative result is the fact that hypoplastic coefficients may be obtained from the second derivative of a material-dependent scalar potential.

Third, starting with the Gibbs potential and choosing $n=0.6,\alpha=0.1$,

[TABLE]

Niemunis, Grandas Tavera and Wichtmann most recently use its second derivative to obtain the coefficient $H_{ijkl}$ of Eq.(2) – a model that they call neo-hypoplasticity 99 . This is clearly a step that we highly welcome. Calibrating the coefficients using incremental stress-strain relations, they found good agreement with observation. The next three sections are about comparing four more recent hyperelastic models, to help optimize the choice of the energy.

III.2.2 Elastic Energy versus the Gibbs Potential

Elatic constitutive relations that possess an explicit elastic energy or Gibbs potential is usually termed a hyperelastic theory. Engineers tend to look for appropriate Gibbs potentials, because they prefer the stress as variable. On the other hand, only the energy is conserved and useful for the hydrodynamic procedure. Although one can obtain one from the other, simple energy expressions typically possess complicated Gibbs potentials, and vice versa. In this section, we collect some general results connecting both, starting with the basic one of the Legendre transformation:

[TABLE]

Assuming the elastic energy is a homogeneous function of degree $m$, we have

[TABLE]

with $B$ a constant. This leads immediately to

[TABLE]

(where prime ′ means derivative). The Legendre transformed Gibbs potential is then

[TABLE]

with $\widehat{u}_{s},\Delta$ taken as functions of $P_{\Delta},\widehat{\pi}_{s}$, especially $\Delta=\left(BP_{\Delta}\,{{\cal F}^{\prime}}/{\widehat{\pi}_{s}}\right)^{\frac{1}{1-m}}.$

A remarkble property of the above energy expression is the simple relation,

[TABLE]

As a result, the Gibbs potential also has the factorized form:

[TABLE]

where the dimensionless shearing factor $F=F\left(\widehat{u}_{s}(\widehat{\pi}_{s})\right)$ is

[TABLE]

For the Hertz contact $m=5/2$, we have $w=B^{-2/3}P_{\Delta}^{5/3}F$.

III.2.3 Comparing Four Hyperelastic Models

We compare four hyperelastic models,

[TABLE]

with ${\xi}=2{\cal B}/5{\cal A}$ a density independent constant. GSH’s elastic energy (with $a$ unspecified) was proposed and considered in 2003, see granL1 . In our later works, we concentrated on $a=1/2$, the Hertz contact value, to be definite, and because Hertz contact seems the appropriate granular picture. Three works on hyperelastic models appeared since, all of the form Eq.(45), by Houlsby, Amorosi and Rojas HousbyA , Einav and Puzrin EinavP , and, as already mentioned, as neo-hypoplasticity, by Niemunis et al 99 . We first note that the power of the “spheric part” $P_{\Delta}^{\frac{{}^{m}}{m-1}}$ of Eq.(51)) are 1.5 for both Einav-Puzrin and Houlsby-Amorosi-Rojas, 1.4 for neo-hypoplasticity, and $\frac{2+a}{1+a}$ for GSH, ie. 5/3 for a=0.5, and again 1.5 for a=1.

Calculating the stress-strain relations for the two energies and their inverse, we have

[TABLE]

The shearing factor $F$ for these models are therefore, see Eq.(52),

[TABLE]

with $F_{0}\equiv F\left(\widehat{\pi}_{s}=0\right)$ and

[TABLE]

The shearing factor of GSH explicitly shows the elastic instability that it contains: The dimensionless shear stress is restricted to $0\leq\widehat{\pi}_{s}\leq\sqrt{a\left(2+a\right){\xi}}$. The root becomes imaginary at the right edge and instability sets in. The other three models remain stable for all $\widehat{\pi}_{s}>0$.

To account for the instability of a strictly static slope that is too steep, since there is no dynamics, one has only the energy to work with. A convex-concave transition – very much in analogy to the van de Waals theory of the liquid-gas phase transition – is an appropriate candidate. Therefore, we prefer to use the elastic energy that contains these instabilities. One may alternatively use brute force and put in the end of elasticity by hand. But then an abrupt transition is forced onto the system, and one cannot account for any precursors. Circumstances are similar if one aims to account for the lack of elastic solutions at low densities, $\rho<\rho_{lp}$, relevant for the transition from elasto-plastic motion to fast flows. (We stress that the angle of repose is given by the critical value, see Sec.IV.3. If connected with the appropriate dynamics, all four models can accounted for it.)

Taking $a=1$, all four models have essentially the same dependence on $P_{\Delta}$, and the same ratio $w/\psi$. Then we only need to compare the shearing factor $F$. (We take $b=3$ for Einav-Puzrin and Houlsby-Amorosi-Rajas.) Employing appropriate values for $\xi$, we can fit the shearing factor of GSH to that of the other three, see (a) of Fig.(6), although the second derivative of GSH diverges, see (c), laying bare the precursor of instabilities.

III.2.4 The Third Strain Invariant

As stressed in the last section, an important advantage of the energy expression, Eq.(16), is its built-in shear instability, Eq.(2b-3), that accounts for granular media’s inability to sustain arbitrarily large static shear stresses. Humrickhouse observed that this instability occurs too early, as compared to the angle of stability (the angle up to which a static layer of sand will remain at rest). And he proposed adding a term given by the third invariant of the elastic strain, $u_{t}^{3}=u_{ij}u_{jk}u_{ki}$, to increase the predicted angle of stabilityhum ,

[TABLE]

Unfortunately, considering only negative values of $\chi$, he came to the erroneous conclusion that this term does not yield any improvement. We subsequently found that it is possible to yield a realistic angle of stability of approximately $30^{\circ}$ with a positive $\chi$. What is more, incremental stress-strain relations and granular acoustics also support the same $\chi$ 3inv .

III.2.5 Proportional Paths and Barodesy

Barodesy’s starting point are the proportional paths as summed up by the Goldscheider rule (GR). Taking p$\varepsilon$p and p$\sigma$p to denote proportional strain and stress paths, they are

• •

A p$\varepsilon$p starting from zero shear stress $\sigma_{s}=0$ is associated with a p$\sigma$p.

• •

The same p$\varepsilon$p starting from an arbitrary $\sigma_{ij}\not=0$ leads asymptotically to the same p$\sigma$p obtained when starting at $\sigma_{s}=0$.

Any constant strain rate $v_{ij}$ is a p$\varepsilon$p: In the principal strain axes $(\varepsilon_{1},\varepsilon_{2},\varepsilon_{3})$, a constant $v_{ij}$ means the system moves with a constant rate along its direction, with $\varepsilon_{1}/\varepsilon_{2}=v_{1}/v_{2},\,\varepsilon_{2}/\varepsilon_{3}=v_{2}/v_{3}$, independent of time. gr states there is an associated stress path that is also a straight line in the principal stress space, that there are pairs of strain and stress paths. And if the initial stress value is not on the right line, it will converge onto it.

To understand GR using GSH, we need the stationary solution for an arbitrary p$\varepsilon$p, ie. we need to generalize the stationary solution as given by Eq.(25,26) to include $v_{\ell\ell}\not=0$. Using the superscript ni to mean non-isochoric, we have

[TABLE]

where the approximate sign holds because $v_{\ell\ell}\ll v_{s}$. We also note that the deviatoric stress, elastic strain and strain rate share the same direction, $\hat{\sigma}_{ij}=\hat{\pi}_{ij}=-\hat{u}_{ij}=-\hat{v}_{ij}$, see Eqs.(14,19,26). If the strain path is isochoric, $v_{\ell\ell}=0$, these are simply the critical state, with both the deviatoric elastic strain and stress being constant. They are hence dots that remain stationary in their respective space, given by ${u_{ij}^{\ast}|_{c}}={u_{c}}(\rho)\,\hat{v}_{ij}$ and ${\sigma_{ij}^{\ast}}=-{\sigma_{s}}(\rho)\,{\hat{v}_{ij}}$. If however $v_{\ell\ell}\not=0$, with the density $\rho[t]$ changing accordingly, ${u_{c}}(\rho)$ and ${\sigma_{s}}(\rho)$ will also change, but not the direction $\hat{v}_{ij}$, making the dots walking down a straight line along $\hat{v}_{ij}$.

This would indeed be the first rule, except that GR states that it is the total (and not the deviatoric) stress that possesses a p$\sigma$p. This is a slightly more involved point. With

[TABLE]

we again have a density dependent prefactor, $[1-\alpha(\rho)]P_{\Delta}(\rho)$, multiplied by

[TABLE]

This factor is a constant direction, giving rise to a p$\sigma$p, if we assume that the system is already in the critical state – as discussed in Sec.III.1.1), $\pi_{s}/P_{\Delta}$ does not depending on the density. Otherwise, we need to rely on $\pi_{s}/P_{\Delta}$ being only weakly density dependent.

The second rule is much easier to understand. Given an initial elastic strain deviating from that prescribed by Eq (58), $u^{0}_{ij}\not=u_{c}\hat{v}_{ij}$, $\Delta_{0}\not=\Delta_{c}$, Eqs (11,12) clearly state that any deviating component will relax toward zero, implying the elastic strain and the associated stress will converge onto the prescribed line.

IV Experiments at Given Shear Stress

In this section, we examine the ramification of GSH if, instead of the shear rate, the stress or elastic strain is being held constant. As we shall see, taking $\partial_{t}u_{ij}^{*}=0$ in Eq.(11) results in a creep motion of the magnitude $fv_{s}/T_{g}=u_{s}/u_{c}$, such that an initial $T_{g}$ will co-relax with $v_{s}$, at the altered rate $R_{T}(1-\sigma_{s}/\sigma_{c})$. Moreover, $T_{g}\propto v_{s}$ relax back to the equilibrium value $T_{g},v_{s}=0$ only for $\sigma_{s}<\sigma_{c}$, growing without bounds for $\sigma_{s}>\sigma_{c}$, though a shear band solution is stable for $\sigma_{s}>\sigma_{c}$. Note hypoplasticity is not applicable here.

This behavior has implications for the angle of repose $\varphi_{re}$ (for uniform stresses): As long as the shear stress is held below the critical one, $\sigma_{s}<\sigma_{c}$, the $T_{g}\propto v_{s}$-relaxation will run its course, and the system is in a static, mechanically stable state afterwards. For $\sigma_{s}\geq\sigma_{c}$, however, the system does not come to a standstill. Therefore, defining $\varphi_{re}$ same as in Eq.(21), but employing the critical rather than yield value, we have

[TABLE]

The inequality holds because the critical state is an elastic solution, while $\varphi_{st}$ of Eq.(21) is the angle at which all elastic solutions become unstable.

Stress-controlled experiments cannot be performed in a conventional triaxial apparatus possessing stiff steel walls, because the correcting rates employed by the feedback loop to keep the stress constant are usually too strong. As a result, too much $T_{g}$ is excited that distorts its relaxation. The situation is then more one of consecutive constant rates, less of constant stress. Instead, one may employ a soft spring to couple the granular system with its driving device, to enable small-amplitude stress corrections without exciting much $T_{g}$.

Typical examples of stress-controlled experiments are creep (both uniform in space and exponentially decreasing), shear band, also flow on an inclined plane or in a rotating drum. The flow either comes to a halt, becomes jammed at the angle of repose $\varphi_{re}$, or starts flowing from that jammed state, becomes fluidized at the angle of stability $\varphi_{st}$.

IV.1 Uniform Creep

For $T_{g}=0$, the system is static, $\sigma_{s}=$ const, and there is no dynamics, only a stable, static elastic solution. But if $T_{g}\not=0$ initially, both $T_{g}$ and the elastic strain will relax according to Eqs.(4,11,12), and with them also the stress $\sigma_{s}$. Maintaining a constant $u_{s}$, or $\sigma_{s}$, therefore requires a compensating shear rate $v_{s}$, which is observed as creep. As we shall see, the characteristic time of $T_{g}$ is then $\propto(1-\sigma_{s}/\sigma_{c})^{-1}$, hence rather long close to the critical stress. So the accumulated shear strain $\varepsilon_{s}(t)=\int v_{s}{\rm d}t$ can be expected to be large.

In a recent experiment involving a fan submerged in sand, Nguyen et al. aging pushed the system to a certain shear stress at a given (and fairly fast) rate, producing an elevated $T_{g}$. Then, switching to maintaining the shear stress, they observed the accumulation of a large total strain $\varepsilon_{s}(t)$ that appears to diverge logarithmically. They called it creeping.

In the experiment, a very soft spring was used to couple the fan and the motor. This fact we believe is essential why the experiment turned out as observed. Usually, triaxial apparatus with stiff steel walls and a feedback loop is used to keep the stress constant. The correcting rates are strong, and much $T_{g}$ is excited that distorts its relaxation. The situation is hence more that of consecutive constant rates, less of constant stress. A soft spring coupling the granular system with its driving device has much lower correcting rates.

Because of the fan, the stress distribution in the experiment is rather nonuniform, rendering a quantitative comparison to GSH difficult. Nevertheless, the experiment’s new and structurally robust results are worth closer scrutiny. And we shall employ GSH assuming uniformity to yield a qualitative understanding. Also, we shall assume it is the elastic shear strain $u_{s}$ that is being kept constant, not the shear stress $\sigma_{s}\propto\sqrt{\Delta}\,u_{s}$, as both cases will turn out to be rather similar. The relevant equations are still Eqs (4,11,12) (with $T_{a}=0$).

At the beginning, as the strain is being ramped up to $u_{s}$ employing the constant rate $v_{1}$, granular temperature acquires the initial value $T_{0}=fv_{1}$. Starting at the time $t=0$, the elastic strain $u_{s}$ is being held constant. Inserting $\partial_{t}u_{s}=0$ into Eq (11), we have

[TABLE]

with $v_{s}$ the rate needed to compensate the stress relaxation. Inserting this into Eqs (4,12), we arrive at

[TABLE]

with the effective $T_{g}$-relaxation rate reduced from $R_{T}$ to

[TABLE]

Assuming uniformity, both equations may be solved analytically, if the coefficients are constant. Being functions of the density, they are if the density is. The stress $P(t),\sigma_{s}(t)$ will then change with time, as will $\Delta(t)$. The first equation accounts for the exponential decay of $\Delta$, from both below and above, to the steady state value

[TABLE]

The relaxation is faster the more elevated $T_{g}$ is. Employing the initial condition $T_{g}=T_{0}$ at $t=0$, the solution to the second equation is

[TABLE]

Because of Eq (62), the same solution also holds for the shear rate, $v_{s}={v_{0}}/(1+r_{v}v_{0}t)$, with $v_{0}\equiv(u_{s}/u_{c})T_{0}/f$, $r_{v}\equiv(fu_{c}/u_{s}){r_{T}}$, implying a logarithmically divergent total shear strain,

[TABLE]

However, we note that $\varepsilon_{s}$ does not actually diverge, because as $T_{g}\to 0$, it enters the quasi-elastic regime, see Eq.(78) below, where its relaxation becomes exponential. So even if creep is large close to $u_{c}$, it comes to a halt eventually, and the system is mechanically stable.

Assuming a large $T_{0}$, $\Delta$ is quickly relaxed. Fixing $u_{s}$ is then equal to fixing the shear stress, $\sigma_{s}\propto\pi_{s}\propto u_{s}\sqrt{\Delta}=(u_{s}^{2}/u_{c})\sqrt{\Delta_{c}}$. In addition, with $\pi_{s}/\pi_{c}=\sqrt{\Delta}\,u_{s}/\sqrt{\Delta_{c}}\,u_{c}=(u_{s}/u_{c})\sqrt{\Delta_{c}}\,u_{s}/\sqrt{\Delta_{c}}\,u_{c}=u_{s}^{2}/u_{c}^{2}$, one may rewrite the factor in $r_{T}$ as

[TABLE]

The $T_{g}$-relaxation is slower the closer $\sigma_{s}$ is to $\sigma_{c}$, infinitely so for $\sigma_{s}=\sigma_{c}$. Then $u_{s}=u_{c}$, $\Delta=\Delta_{c}$, with $T_{g}(t)=T_{0}=fv_{1}$ a constant that does not relax. This is indistinguishable from the rate-controlled critical state, which, clearly, may be maintained also at given stress.

For $\sigma_{s}>\sigma_{c}$, the relaxation rate is negative, and $T_{g}\propto v_{s}$ will grow, or explode, seemingly without bound, see the next two sections, each with a possibility of what happens next.

Only a $T_{g}$ sufficiently large will explode, destabilizing a static shear stress exceeding $\sigma_{c}$, an infinitesimal $T_{g}$ will not. This is because the critical stress diverges for $T_{g}\to 0$ [as can be seen calculating the critical stress employing Eq.(78) below], and the window between $\sigma_{c}$ and $\sigma_{s}^{yield}$ vanishes. Therefore, a static shear stress remains metastable for $T_{g}=0$, turning instable only at the yield stress, $\sigma_{s}^{yield}/\sqrt{2}P=\sqrt{{\cal A}/{\cal B}}$, as given in Eq.(21).

To summarize, as long as $\sigma_{s}<\sigma_{c}$, both the temperature and the shear flow will relax to zero, with the rate $r_{T}\equiv R_{T}[1-\sigma_{s}/\sigma_{c}]$. For $\sigma_{s}=\sigma_{c}$, we have the rate-independent critical state. For $\sigma_{s}^{yield}>\sigma_{s}>\sigma_{c}$, the system is meta-stable, and easily perturbed into developing a shear band as given in the next section, IV.2. For $\sigma_{s}>\sigma_{s}^{yield}>\sigma_{c}$, the stress necessarily has the full form of Eq.(14), containing both the gaseuous pressure, viscous stress, and also the elastic part that is always in the critical state, see Sec.IV.3.

IV.2 Narrow Shear Bands

For $\pi_{s}^{yield}>\pi_{s}>\pi_{c}$, a stable, localized steady-state solution exists that we may identify as the shear band: Taking $\partial_{t}T_{g}=0$ in Eq.(64), the steady-state solution for $0\leq x\leq\xi_{sb}$ is

[TABLE]

The velocity difference across the band is $\Delta v=\int v_{s}{\rm d}x=\int v_{s}^{0}\sin(\pi x/\xi_{sb}){\rm d}x$, hence

[TABLE]

This is to be combined with $v_{s}\propto T_{g}\equiv 0$ for $x\leq 0$ and $x\geq\xi_{sb}$. We note that although $T_{g}$ and $v_{s}$ are continuous at $0,\xi_{sb}$, neither is differentiable there. Also, the density $\rho$ must be lower in the band, because only than can the shear stress $\sigma_{s}$ (that has to be uniform) be smaller than the critical stress $\sigma_{c}(\rho)$ in the quiescent region, but larger in the band.

IV.3 Angle of Stability and Repose

If the system stays uniform at a shear stress exceeding the critical one, $\sigma_{s}>\sigma_{c}$, the variables $T_{g}\propto v_{s}$ will grow with $r_{T}$, until the gaseous pressure and viscous stress become relevant, see Eq.(14), while the elastic stress stays in the critical state,

[TABLE]

This is the stress expression for fast dense flow, again a stable and uniform solution, though no longer rate-independent. We shall not dwell on it here, see granRexp ; PG2013 for results.

On a tilted plane with a ongoing granular flow and a stress as given by Eq.(72), reducing the inclination angle will reduce both the shear stress and the shear rate $v_{s}$. And the flow will come to a stop, $T_{g}\propto v_{s}\to 0$, if the stress drops below the critical one, as discussed in Sec. IV.1. Hence the angle of repose (for uniform stress) is $\tan\varphi_{re}={\sigma_{c}}/\sqrt{2}\,P$, see Eq.(61.

Given a static elastic state on a plane inclined by the angle $\varphi$, the angle of stability $\varphi_{st}$ (again for uniform stress) is reached when the energetic instability of Eqs.(20) is breached, and the static elastic state collapses. Therefore, $\tan\varphi_{st}=\sigma_{s}^{yield}/\sqrt{2}P=\sqrt{{\cal A}/{\cal B}}$, see Eq.(21).

IV.4 Wide Shear Bands and KCR

Finally, we consider the coexistence between fast dense flow, $T_{g}=fv_{s}$, and one of static elastic stress, $T_{g},v_{s}=0$. Some $T_{g}$ will diffuse into the static region, giving rise to an exponentially decreasing creep, as observed by Komatsu et al komatsu , Crassous et al crassous .

We take this “liquid-solid boundary” to be at $x=0$, with $x>0$ being solid. In the fluid phase $x<0$, the shear rate is a constant, as is $T_{g}$, providing the boundary conditions for the solid part. Pressure $P$ and shear stress $\sigma_{s}$ are uniform, as is $v_{\ell\ell}\equiv 0$, but $\rho$ is discontinuous at $x=0$, though constant otherwise. (The liquid density needs to be lower than that of the solid one, such that the same stress exceeds the critical one in the liquid, but is below it in the solid.) The variable $T_{g}=fv_{s}$ varies along $\hat{x}$. These circumstances are similar to that of Sec IV.1, though variation is in space rather than time. Given stationarity of Eqs (11,12), we have

[TABLE]

see Eqs.(62,66). With $\Delta,u_{s}$ fixed, so are $P,\sigma_{s}$, where especially $P=P_{c}$ if $\sigma_{s}=\sigma_{c}$. Note that the above relations imply $\Delta/{u_{s}}=({\Delta_{c}}/{u_{c}})({fv_{s}}/{T_{g}})$. And since $1/\mu\equiv P/\sigma_{s}$ increases monotonically with $\Delta/{u_{s}}$, see Eq.(19), the friction $\mu$ decreases for increasing $fv_{s}/T_{g}$.

The balance equation for $T_{g}$ with $\partial_{t}T_{g}=0$, see Eq.(64), yields an evanescent creep,

[TABLE]

($T_{g}^{0}$ and $v_{s}^{0}\sqrt{}$are the liquid values at $x=0$.) That the decay length $\xi_{cr}$ diverges for $\sigma_{s}=\sigma_{c}$ is not surprising, because the whole solid region turns critical then.

Kamrin et al kamrin ; kamrin2 proposed a constitutive relation (KCR) to account for the steady flow in the split-bottom cell fenistein ; fenisteinB . A key ingredient of the theory is a new variable fluidity: $g\equiv v_{s}/\mu$, with $\mu\equiv\sigma_{s}/P$ and $\mu_{s}\equiv\sigma_{c}/P_{c}$. It is taken to obey the equation

[TABLE]

As $g_{loc}=0$ for $\mu<\mu_{s}$, this relation is rather similar to Eq.(74), with $g$ assuming the role of $T_{g}$, and the decay length diverging at the critical stress. For $\mu\geq\mu_{s}$, the system is fluid, and $g=g_{loc}$ essentially constant. Again, $T_{g}$ is constant in the liquid region.

KCR is well capable of accounting for the steady flow in the split-bottom cell. This is fortunate, because again, there is a symbiotic relation relation between GSH and KCR, similar to that between GSH and hypoplasticity. But there are also three differences.

First, KCR assumes $\sigma_{s}>\sigma_{c}$ in the liquid and $\sigma_{s}<\sigma_{c}$ in the solid, with $\sigma_{c}$ a constant. This violates momentum conservation, as the stress is continuous across the solid-fluid interface (see also the discussion in Sec.IV.2 and IV.4). Second, we have, for GSH, ${fv_{s}}/{T_{g}}=1$ in liquid, and ${fv_{s}}/{T_{g}}=\sqrt{\sigma_{s}/\sigma_{c}}$ in solid; yet $v_{s}/g=\mu\equiv\sigma_{s}/P$ on both sides for KCR. Third, Kamrin and Bouchbinder constructed a “two-temperature continuous mechanics” in a recent paper ttcm , defining a configuration temperature $\theta_{c}$ and one for the vibrational degrees of freedom $\theta_{v}$. (Microscopic degrees of freedom are not considered. There would be three temperatures if one did.) They hypothesize that the fluidity $g$ may be related to $\theta_{c}$. We disagree here. As discussed above, $g$ is essentially $T_{g}$ of GSH, which is closer to $\theta_{v}$.

V The Quasi-Elastic Regime

Quasi-static is what many in soil mechanics call the rate-independent elasto-plastic motion. The implication is: Being the slowest one possible, elasto-plastic motion stays as is however slow the rate. We believe there is, at very small $T_{g}$ and very low rates – possibly $v_{s}\sim 10^{-6}$/s, or $I\sim 10^{-8}-10^{-9}$ (where $I\propto v_{s}/\sqrt{P}$ is the inertial number of the $\mu(I)$ rheology), a rate-dependent transition to another rate-independent regime. If true, only the latter should be called quasi-static, or quasi-elastic. The reasons for our believe are the following six points [see also Sec.II.3 and the list after Eq.(19)]:

We have $T_{g}=0$ in static stress distributions, which is well accounted for by the fully elastic Eqs.(23). This should remain so for very slow stress changes. 2. 2.

The fact that cyclic load or ratcheting are less dissipative than hypoplasticity predicts Herle may well be interpreted as the onset of a transition to quasi-elasticity. One does not need to invoke intergranular strain here. 3. 3.

Quasi-static motion – a consecutive visit of static, equilibrium states – should occur at $T_{g}=0$, and be reversible. Yet elasto-plastic motion occurs at elevated $T_{g}$, is strongly dissipative and irreversible. Reactive and dissipative terms are comparable in size, being exactly equal in the critical state. (This is obvious only within the context of GSH, but one can identify the reactive and dissipative terms in hypoplasticity by comparing it to GSH. Then the same holds true there.) 4. 4.

Rate-independence is a hallmark of quasi-static motion. Yet the rate-independence of elasto-plastic motion is (as discussed in Sec.III.1.4) easily destroyed. 5. 5.

$T_{g}$ is essentially zero in elastic waves, since there is not enough time to excite appreciable amount of it. Waves are well accounted for by the elastic Eqs.(23), but not by hypoplasticity – as they should if there were only the elasto-plastic regime. Employing hypoplasticity, one finds that waves are always over-damped: Starting from momentum conservation, $\partial_{t}(\rho v_{i})=-\nabla_{j}\sigma_{ij}$, or $\,\,\partial^{2}_{t}(\rho v_{i})=-\nabla_{j}\partial_{t}\sigma_{ij}$, one inserts Eq.(2) to find the term $\propto H_{ijk\ell}$ yielding the velocity and the one $\propto\Lambda_{ij}$ damping. Both terms are of the same order in $q$, and given the values of the two tensors appropriate for elasto-plastic motion, they are also comparable in size, implying that elastic waves are always over-damped. (In solids, the damping term is an order in the wave vector $q$ higher than the velocity term, rendering it much smaller for long enough wave lengths. This is the reason elastic waves propagate there.). In fact, elastic waves are a fast yet quasi-static phenomenon, same as incremental stress-strain relations. 6. 6.

Happily and perhaps most convincingly, Peng Zhen et al. have most recently observed indications of a rate-dependent transition away from the elasto-plastic regime QuasiEla .

V.1 Strain Increments versus Slow Rates

Faced with the quasi-elastic response of sound waves and ratcheting, a frequent suggestion is to take a small strain increment to be elastic and free of dissipation, but a large one as elasto-plastic and dissipative. Unfortunately, this is not compatible with the notion of quasi-static motion, a consecutive visit of equilibrium, static states: Starting from an arbitrary static state of given stress, applying a small incremental strain that is elastic, the system ends up in another static state, with a slightly higher stress. This second state, as static as the first, is again a valid starting point – any static state is. And the next small increment must again be elastic. Many consecutive small increments yield a large change in strain, and if the small ones are not dissipative, neither can their sum be.

In GSH , it is the strain rate, not the strain amplitude, that decides whether the response is elastic or elasto-plastic. Small strain increments produced with a brief puls of shear rate will produce an elastic response, if $T_{g}$ does not have time to get to a sufficiently high value.

V.2 How to Observe Quasi-Elasticity

There are at least two ways to observe the elusive quasi-elastic behavior, both by fixing the stress rate and keeping a low $T_{g}$: Note that a given stress rate may be associated with either of two shear rates, a high, elasto-plastic one at elevated $T_{g}$ and a low, elastic one at vanishing $T_{g}$. The first method is to slowly incline a plane supporting a layer of grains. In such a situation, the shear rate remains very small, and the system starts flowing only when the yield stress is breached.

A second method is to insert a soft spring between the granular medium and the driving device. If the spring is softer by a large factor $a$ than the granular medium (which is itself rather soft), it will serve as a “stress reservoir,” and absorb most of the displacement, leaving the granular medium deforming at a rate smaller by the same factor $a$. (Employing a feedback loop in a triaxial apparatus with stiff steel walls to maintain a stress rate will not usually work, because the correcting motion typically has strain rates that are so high, that the system slips into an elasto-plastic motion quickly.)

V.3 Transition to Quasi-Elasticity

To account for the transition to quasi-elasticity, we need to specify how, in Eqs.(9,10,14), $\lambda T_{g},\lambda_{1}T_{g},\alpha,\alpha_{1}$ vanish. There are scant data to rely on, but one possibility would be

[TABLE]

For $T_{g}\gg T_{0}$, we have $\alpha=\alpha^{0},\,\,\alpha_{1}=\alpha^{0}_{1}$, ensuring rate-independence for the elassto-plastic regime; for $T_{g}\ll T_{0}$, we have $\alpha,\alpha_{1}\propto T_{g}$ vanishing with $T_{g}$.

Instead of taking $\lambda,\lambda_{1}\to 0$ with $T_{g}$, we choose to generalize $T_{g}^{2}=f^{2}v_{s}^{2}$ to

[TABLE]

because this was a natural result arising during the derivation of GSH, see granR4 ; granRexp . For $T_{g}\ll T_{1}$, we have $T_{g}=(f^{2}/T_{1})\,v_{s}^{2}$ quadratically small in $v_{s}$, and with it also the relaxation rates $\propto T_{g}$. The ratio, in Eqs.(9,10), between the irreversible, plastic term ($\propto v_{s}^{2}$) and the reversible, elastic term ($\propto v_{s}$) diminishes with $v_{s}\to 0$, as has been observed in DEM AH .

VI Conclusion

Taking the stress increment as a function of density, strain rate and the stress itself, the Eq.(1) this paper starts from to set up a constitutive model, is an appropriate and physically sound way to come to terms with the complexity of triaxial results. A rather broader view of granular behavior are captured by introducing a pair of internal state variables, the elastic strain field $u_{ij}$ and the granular temperature $T_{g}$, with the first accounting for the coarse-grained deformation of the grains, and the latter accounting for their quickly changing elastic and kinetic energy. The theory that does it is called GSH, and the phenomena accounted for include static stress distributions, creep, shear band, angle of repose and stability, also the critical state under tapping, all in addition to the trixial results.