Mechanical response of an inclined frictional granular layer approaching unjamming

Various amorphous materials, and granular media in particular, exhibit a so-called jamming transition between rigid and flowing states. The nature of this transition has been investigated during the last decade, see recent reviews [1, 2]. Most granular studies have focused on frictionless discs or spheres, typically controlled in volume fraction ϕ or in pressure P [3–5], showing that the jamming transition is critical (scaling exponents, diverging length scale) [3, 6, 7] and related to isostaticity [3, 8–11]. As the system loses its mechanical rigidity at the transition, its shear modulus G is found to vanish as a power law with respect to the distance to jamming ϕ − ϕ_c, where ϕ_c is the critical volume fraction. Vibrational mode analyses have been performed on these frictionless granular systems [3, 6, 12–14], or similarly on soft glassy materials [15–17], reporting an increasing number of soft modes, corresponding to collective low-energy motion of the particles, as the transition is approached.

Much fewer studies have dealt with frictional grains in this context and they have mostly considered homogeneous systems under isotropic pressure [11, 18–22]. In the frictional case, the Liu-Nagel jamming concept [23, 24] must be revised [25]. In particular, jamming and isostatic points do not coincide any more [1], and one thus can expect a finite shear modulus at the transition. In this letter, we consider static layers of frictional grains under gravity, by means of two-dimensional discrete element simulations (standard molecular dynamics [26]), and investigate their mechanical properties through the analysis of their stress response to a localized overload $\vec F_0$ at the layer surface. The layers are prepared at a fixed angle α with respect to the horizontal (see fig. 2 for notations), and unjamming is approached as α is close to α_c, the critical value above which static layers cannot be equilibrated at that angle and always flow. Note that this unjamming point α_c is close in spirit to the situation of a jammed solid sheared up to its yield-stress [27]. It is also close, but different, to progressively tilted granular layers, which eventually lose their mechanical stability, see, e.g., [28, 29]. In our simulations, the volume fraction in the layer is fairly uniform all through the layer depth. Its value depends very weakly on the inclination angle, varying from ϕ ≃ 0.823 at α = 0 to ϕ ≃ 0.816 at α_c (these are for the GG preparation, see below). These values are always larger —though not much— than the critical value, estimated in our system at ϕ_c ≃ 0.815 [21, 30, 31]. Then, here, the control parameter for the jamming/unjamming transition is the sole angle α. This situation is therefore qualitatively different to the homogeneous configurations submitted to isotropic pressure cited above, and is effectively closer to an experimental set-up.

The numerical model is that described in [32, 33], with N = 3600 polydisperse frictional discs coupled, when overlapping, by normal and tangential linear springs, tangential forces being limited by the Coulomb condition with a friction coefficient μ = 0.5. Two system preparations have been considered: a grain-by-grain (GG) and a rain-like (RL) deposition of the particles on a rough substrate consisting of fixed but size-distributed particles, inclined at the desired angle α. The layer is prepared when all grains have reached static equilibrium (see [32] for a description of these protocols). No external pressure applied to the topmost layer of particles, i.e., the pressure in the system is due solely to the gravitational force acting on the particles themselves. The typical thickness of the layer is ≃23 grain diameters, i.e., a system aspect ratio around 1/7. Above a certain inclination α_c, the simulations preparing the packing before the response procedure do not converge towards a static layer —the grains always flow. The "solid-liquid" transition is rather abrupt (Δα ≃ 0.5°), allowing for a well-defined value of this critical angle: α_c ≃ 20.75° for the GG and α_c ≃ 20° for the RL preparations, respectively. For the sake of concision, except otherwise stated, the displayed data are those obtained with the GG preparation.

Experimental and numerical works have shown that the linear stress response of granular systems to a point force is well described by (possibly anisotropic) elasticity [32, 34–38]. Here we show that an orthotropic elastic model is able to reproduce quantitatively the normal and tangential stress bottom profiles, and that the corresponding shear modulus G decreases with α, in parallel to a decrease of the coordination number and minor changes in the micro-structure (fabric). However, it does not vanish at the transition, but reaches a finite value as α → α_c. We finally analyze the vibration modes of the system and show that the soft modes play an increasing, though not crucial, role approaching the transition.

Once a layer is deposited, stabilized in an equilibrium state, an additional force $\vec F_0$ is applied on a grain close to the free surface, and a new equilibrium state is reached (see [32] for details). Taking the difference between the states after and before the overload, one can compute the contact forces in response to $\vec F_0$. Introducing a coarse graining length w, the corresponding stress response can be determined. Taking w of the order of few mean grain diameters (here $w = 6 \left < d \right >$) as well as an ensemble averaging of the data (here N_e ∼ 150 realisations for each tilt angle α), make the stress profiles quantitatively comparable to a continuum theory [33], such as elasticity, as discussed below. The amplitude of the overload was kept constant for all simulations: $F_0 = 1.0 \left < m \right > g$, where $\left < m \right >$ is the average mass of the grains. This value is sufficiently small to ensure a linear [39, 40] and reversible response of the system for all values of α, including close to α_c.

Three examples of stress bottom profiles σ_nn(t) and σ_tn(t) for the GG preparation are displayed in fig. 1 for two values of the slope α of the layer and two values of the angle θ that the overload force makes with the normal direction (see fig. 2). The stress profiles computed for the RL preparation have approximatively the same qualitative behavior.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In order to reproduce the stress response profile, isotropic elasticity is not enough [32], and we consider here orthotropic elasticity, characterized by a stiff axis (1) and a soft one (2) associated to two Young moduli E₁ and E₂ < E₁. Two Poisson coefficients ν₁₂ and ν₂₁ must also be introduced. However, they are not independent and verify ν₁₂/E₁ = ν₂₁/E₂. Finally, G is the shear modulus, so that the stress-strain relation, in the orthotropic axis writes

$$\begin{eqnarray} &&\left ( \begin{array}{@{}c@{}} \epsilon_{11} \\ \epsilon_{22} \\ \epsilon_{12} \end{array} \right ) = \left ( \begin{array}{@{}c@{\qquad}c@{\qquad}c@{}} \frac{1}{E_1} & -\frac{\nu_{21}}{E_2} & 0 \\ -\frac{\nu_{12}}{E_1} & \frac{1}{E_2} & 0 \\ 0 & 0 & \frac{1}{2G} \end{array} \right ) \left ( \begin{array}{@{}c@{}} \sigma_{11} \\ \sigma_{22} \\ \sigma_{12} \end{array} \right ). \end{eqnarray} \tag{ 1 }$$

Elastic energy is well defined if all moduli E₁,E₂,G are positive and 1 − ν₁₂ν₂₁ > 0. A last parameter of this modeling is the angle τ that the axes (1,2) make with (n,t) (see fig. 2).

The stress responses have been analytically computed in [41] for a semi-infinite layer (h → ∞). For a finite layer thickness, a numerical integration, in the spirit of [34], must be done. Rough bottom boundary conditions (zero displacement) are imposed. We can adjust four dimensionless parameters (two modulus ratios G/E₁ and E₂/E₁, a Poisson coefficient ν₂₁ and the orthotropic angle τ) to reproduce, for a given α, the numerical data for all θ. As shown in fig. 1 (solid lines), the fitting functions quantitatively describe the numerical data.

In fig. 3, we show the behavior of the elastic modulus ratios G/E₁ and E₂/E₁ as the slope α is larger. They both decrease, almost by a factor of two, from a fairly constant value at small α to a lower one close to α_c. In particular, G does not vanish close to the critical angle, in agreement with the observation that frictional granular systems remain hyperstatic at the unjamming transition [11, 18, 19]. Such a discontinuous behaviour has been seen in simulations by Otsuki and Hayakawa [31] investigating the rheology of sheared frictional grains close to jamming, and in experimentally created shear-jammed states reported in [25]. The ratio G/E₂ is almost constant (not shown). This behavior is also found for the RL preparation, although the variations of G/E₁ are less pronounced. This softening is finally consistent with acoustic experiments on a granular packing in the vicinity of the transition [42].

Zoom In Zoom Out Reset image size

In addition to the above global mechanical properties of the system, we have studied the evolution of various microscopic quantities with α. The first one of interest is the coordination number Z, i.e., the average number of contacts per grain, here computed in the bulk of the layer, where it is fairly uniform —it obviously drops down close to the surface. Z monotonously decreases with α (fig. 4(a)) and stays always far from the isostatic value Z_iso = 3 (for frictional grains in 2D). As evidenced by the comparison of the curves in figs. 3 and 4(a), the modulus ratio G/E₁ is not found to be a linear function of Z − Z_iso, in contrast with the finding of [18] on homogeneous frictional systems, close to isostaticity. Grains of the bulk that only carry their own weight do not contribute much to the global stability of the contact network. As for so-called rattlers in gravity-free packings (see [43], Chapt. 6), these grains can be removed from the contact counting, leading to a modified coordination number of the layer Z* (see fig. 4(a)). However, we have found that their number is roughly independent of α, so that the relation between G/E₁ and Z* is not linear either. We have also studied the friction mobilisation at the contact level. In our MD simulations, the number of contacts with a ratio of the tangential force f_t to the normal force f_n strictly equal to the microscopic friction μ is zero when static equilibrium is reached. However, some of them are effectively close to the Coulomb criterion, and we have computed the average $\langle \frac {|f_t|}{\mu f_n} \rangle$. This quantity, displayed in fig. 4(a), increases as α → α_c, but its overall variation is weak (see right y-scale).

Zoom In Zoom Out Reset image size

Finally, we have studied the contact angle distribution. Three of these distributions are represented as polar diagrams for α = 0, 10 and 20 degrees in fig. 4(b). The four strongly pronounced lobes are typical of the GG preparation [43] (Chapt. 6). The vertical and horizontal directions are always in between these lobes. When the layer is horizontal, the orthotropic stiff and soft directions are also found to be (almost) along the horizontal and vertical axis, respectively. Note that the fitting procedure effectively gives here τ = 93° in this case, while τ = 90° (or 0°) would have been expected for symmetry reasons. This indicates the typical precision we have on the measure of this orthotropic angle. Close to the critical slope, however, the orthotropic orientations are close to those of the lobes, the stiff one being in the direction of the slope. As evidenced in fig. 4(c), the transition between these two microscopic configurations occurs around α ≃ 9°, i.e., well below α_c, in correspondence with the drop by a factor of 2 of G/E₁ between 8 and 10° (see fig. 3).

For each tilt angle α, we have computed the 3N vibration modes of the layer, performing a harmonic approximation of the energy around the equilibrium point reached by deposition of the grains. We only studied layers obtained by GG deposition. They do not show contacts being at the Coulomb threshold. This avoids the issue of the relevance of the harmonic approximation for such contacts [19]. Each contact between two grains is thus being considered as two linear springs (normal and tangential stiffnesses k_n and k_t; here k_t = 0.75k_n). The reference vibration frequency is that of a single spring-mass system $\omega _0 = \sqrt {k_n/\left < m \right >}$. We denote $D(\bar \omega )=(3N)^{-1}\sum _{{\mathrm { mode}}\, i} \delta (\bar \omega - \bar \omega _i)$ and $Q(\bar \omega ) = \int _0^{\bar \omega } D(\bar \omega ')\, {\mathrm { d}}\bar \omega '$, the density and cumulative density of states, where $\bar \omega = \omega /\omega _0$ is the normalized frequency.

The function Q is displayed in fig. 5 for four values of the slope α. No excess of vanishing frequencies are observed as the transition is approached, which is expected for a system that stays hyperstatic [44]. There is no single zero-frequency mode appearing at the transition either. The criticality of the layer thus cannot be described by the emergence of an energetically costless direction in configuration space. However, some rather soft modes ($0.1 \lesssim \bar \omega \lesssim 0.5$) are enhanced when α → α_c. We have also computed the eigenfunctions corresponding to these modes. The displacement field in response to the overload with $\vec F_0$ can be decomposed on these functions, as they form a mathematical basis. We denote c_i the coefficients of this decomposition. A large part of the decomposition is owned by the  20 lowest frequency modes. Hence, we can compute the number of modes n₇₀ which are necessary to represent 70% of the displacement field, and the highest frequency $\bar \omega _{70}$ of these n₇₀ modes, by looking for the lowest frequency $\bar \omega _{70}$ such that $\sum _{i | \bar \omega _i < \bar \omega _{70}} c_i^2 > 0.70$. Both n₇₀ and $\bar \omega _{70}$ are plotted against α in the inset of fig. 5. While n₇₀ stays approximatively constant and around 20, $\bar \omega _{70}$ is found to decrease from ≃0.3 by a factor of two when α varies from zero to α_c, without however showing a clear singularity at the transition. This means that the response is progressively better represented by softer modes, even if these modes do not trigger the critical behavior of the layer close to the transition.

Zoom In Zoom Out Reset image size

To sum up, we have simulated 2D frictional and polydisperse granular layers under gravity inclined at an angle α, and investigated their mechanical and microscopic properties when the unjamming transition is approached. This work tells us what to expect in real experiments, i.e., a layer that becomes elastically softer as α → α_c, but not to the point at which the system would lose its rigidity before avalanching. Progressive enhancement of rather soft vibration modes should be observable, although they do not seem to play a particularly crucial role.

An interesting continuation of this work is to study this system with smaller and larger values of the contact friction coefficient μ. We have seen that GG layers prepared with μ = 0.1 have a significantly smaller critical angle α_c ≃ 14°, whereas those prepared with μ = 10 (α_c ≃ 21°) do not differ much from the case μ = 0.5. Similarly, the grains are more connected and more densely packed at that critical angle in layers with μ = 0.1 (Z ≃ 3.6 and ϕ ≃ 0.833), whereas the values of Z and ϕ do not change much when the system is prepared with μ = 10. However, the systematic investigation of the role of μ on the mechanical response of the granular layer, and on the elastic shear modulus in particular, is beyond the scope of the present paper. Another possible perspective could be to use granular simulations with a rolling resistance [45] in order to explore a wider range of ϕ, Z and α.

We thank B. Andreotti, E. Clément, J. Kurchan and M. Trulsson for useful discussions and a careful reading of the manuscript. This work is part of the ANR JamVibe, project No. 0430 01. APFA has been partially supported by the exchange program "Science in Paris 2010" (Mairie de Paris). APFA thanks CNPq and FAPEMIG Brazilian agencies for financial funding.