Correlation and shear bands in a plastically deformed granular medium

Abstract

Recent experiments (Le Bouil et al., Phys. Rev. Lett., 2014, 112, 246001) have analyzed the statistics of local deformation in a granular solid undergoing plastic deformation. Experiments report strongly anisotropic correlation between events, with a characteristic angle that was interpreted using elasticity theory and the concept of Eshelby transformations with dilation; interestingly, the shear bands that characterize macroscopic failure occur at an angle that is different from the one observed in microscopic correlations. Here, we interpret this behavior using a mesoscale elastoplastic model of solid flow that incorporates a local Mohr-Coulomb failure criterion. This differs from the interpretation of Le Bouil et al., which is based on purely elastic considerations ignoring the potential role of local friction on deformation patterns. We show that the angle observed in the microscopic correlations can be understood by combining the elastic interactions associated with Eshelby transformation with the local failure criterion. At large strains, we also induce permanent shear bands at an angle that is different from the one observed in the correlation pattern. We interpret this angle as the one that leads to the maximal instability of slip lines.

Introduction

Plasticity is an important mechanical property in a wide variety of amorphous systems such as dense colloidal glasses, foams, emulsions, and fine-grained granular packings. It is formally defined as intense unrecoverable (shear) deformations that the material undergoes beyond its elastic limit without any crushing or crumbling. This phenomenon has been linked to the so-called yielding transition^1,2,3,4 with some universal features associated with it. Universality emerges in spite of the diversity in disordered solids – in terms of their scales, microscopic constituents, or interactions, suggesting common underlying mechanisms.

A commonly accepted picture that supports this universal character is that the bulk plastic response emerges from a collective dynamics that is not specific to the particle, but rather results from interactions mediated by the universal laws of linear elasticity. This emergent dynamics is characterized by plastic events or shear transformations that are localized in space and time, but have long-range (compared to the size of rearranging zones) elastic-type consequences⁵. In systems in which thermal fluctuations are irrelevant (which will be the case of the granular systems considered in this work), rearrangements are initially activated by external deformation, but further instability may be triggered and propagated due to non-local interactions.

In this framework, propagation of plasticity is a dynamical process which, once the characteristics of the shear transformations and the elastic properties of the medium are known, depends only on the dissipation mechanism^6,7. Near the yielding transition, the so-called avalanche dynamics may emerge in which the activation process takes place by sequentially forming clusters of all scales. In this regard, plastic yielding may be thought as a true second-order phase transition with unique characteristics such as diverging length and/or timescales and power-law distributions of avalanche sizes^8,9,10.

Due to structural heterogeneities and stochastic aspects of the shear transformation dynamics, the cascades of events described as avalanches are usually organized in a highly intermittent manner in space and time. Intermittency makes them quite distinguishable from the much longer lived localization patterns that emerge upon ultimate failure and are often known as shear bands, i.e. narrow linear (in 2 d) or planar (in 3 d) structures along which plastic activity accumulates while the bulk of the material of the system remains undeformed. Still, it is expected that pre-failure collective dynamics must have a strong connection with the formation of permanent shear bands. In fact, it has been shown in several works that elasticity-based polar features characterize the structure of correlations between plastic bursts^11,12,13, which occur preferentially at 45° for a volume conserving plastic event. This preferential direction (with respect to the principal axis of the local shear event) is also the one along which a shear band should form in a system deformed at constant volume (in particular incompressible), as it corresponds to the direction of maximum macroscopic shear stress.

However, the symmetry that Le Bouil et al.¹⁴ (see also ref.¹⁵ for the experimental setup) observed in correlation patterns of the plastic activity seems to be at odds with the morphology of the fully formed band they also observed in strongly deformed granular media. While the latter was formally described as reflecting the Mohr-Coulomb macroscopic failure angle, the former revealed a distinct orientation. Moreover, both angles differ from the “canonical” 45°. In order to rationalize the preferential direction observed in correlations, the authors proposed to take into account the possibility of a local dilation in the plastic event, a feature that was already noted to affect the preferential directions in collective bursts plastic activity¹⁶. The authors proposed a phase co-existence scenario in which recurring mini bursts persist up until failure, which may be interpreted as a signature of discontinuous first-order or spinodal transition. However, this picture does not lead to a specific prediction for the angle of global failure and emergence of macroscopic friction in the system.

In this work, our aim is to provide, based on minimal ingredients, a scenario that explains simultaneously the deviation from the “canonical” 45° direction for the different observables, and the fact that different directions are observed for the correlations in intermittent activity and in permanent shear localization. Our analysis will be based on the hypothesis that, while the interactions between plastic events are mediated by the universal laws of linear elasticity, the triggering of events depends on a local failure criterion that is completely independent of these laws, and can be quite arbitrary. In granular media, for instance, with frictional resistance being an important mechanical property, the yield criterion may be stated as a critical ratio between the resolved shear and normal stress at any arbitrary material point. As a result, the correlations in plastic activity and its propagation involve a compromise between directions favored by the elastic interactions and those favored by the local failure criteria. The results that emerge from this compromise are nontrivial, in agreement with the experimental observations.

From Le Bouil et al. perspective, friction is a topological property that emerges at global scales leading to a Mohr-Coulomb type failure; it seems to play no role in the dynamics of transient slip bands. In our approach, friction has a mesoscopic notion as well and, along with elastic couplings, governs triggering dynamics at intermediate scales. This results in a frictional flow which is collective at macroscopic scales and appears as a shear-banding phenomenon with extended length- and/or time-scales associated with it. We also find that the “bulk” friction coefficient is different than the local ingredient which suggests a separation of scales between the meso and macro view of the plastic flow.

The minimal ingredients used in our analysis are the description of local plastic events as transformations of Eshelby-type¹⁷, combined with the use of a local Mohr-Coulomb criterion that introduces a pressure sensitivity in the propagation of plasticity. It can be seen as an extension to the case of Mohr-Coulomb failure of works that involve permanent damage as a cause of strain localization^18,19,20. Remarkably, however, no permanent damage is required for inducing coexistence between transient micro events and fully developed shearing bands.

The organization of the paper is the following. In Sec. 1 we combine the hypothesis of Mohr-Coulomb local failure with the stress redistribution prescribed by Eshelby’s elastic theory, and discuss hypothetically the consequences on the correlation patterns and macroscopic failure angle. In Sec. 2 we describe the numerical model that incorporates these basic ingredients. The numerical results and conclusions are given in Sec. 3 and Sec. 4, respectively.

The Mohr-Coulomb Failure Criterion

Previous studies within the framework of mesoscopic elasto-plasticity^2,21 emphasized the role of elastic kernels in the yielding transition. The local yielding rule however, coupled with long-range elasticity, must have a strong relevance on the structure of correlations and macroscopic failure. In this section, we propose a simple theory taking into account these effects at both local and macroscopic scales. At the mesoscopic scales, we made a common assumption that, on average, the far-field stress patterns follow the classical Eshelby’s solution in response to a localized shear. The actual response, however, displays significant fluctuations due to the disordered and discrete nature of amorphous media. Interestingly, statistical quantities at large scales, i.e. plastic-activity correlations, appear to be almost insensitive to the quenched disorder and can be, therefore, replicated within the context of elasto-plastic models and continuum framework.

Meso-scale Failure Criterion and Elastic Stress Redistribution

We consider an infinite elastic matrix characterized by its bulk modulus K and shear modulus μ and an embedded inclusion going through a shear transformation ${\varepsilon }_{\alpha \beta }^{{\rm{tz}}}={\varepsilon }^{\ast }{a}^{d}({\delta }_{x\alpha }{\delta }_{x\beta }-{\delta }_{y\alpha }{\delta }_{y\beta })$ with an amplitude ε^* and microscopic volume a^d. Here δ _αβ is a Kronecker delta and we focus on the two dimensional case d = 2.

The stress tensor at a material point can be expressed as σ _αβ = −pδ _αβ + σ(δ _αx δ _βx − δ _αy δ _βy ) + σ _xy (δ _αx δ _βy + δ _αy δ _βx ) in terms of the pressure p and area preserving axial shear σ and diagonal shear σ _xy . The far-field perturbation fields given by Eshelby’s solutions at a point with polar coordinates r, θ read¹⁷

$$\begin{array}{rcl}\delta p & = & \frac{2\mu {\varepsilon }^{\ast }}{1+\frac{\mu }{K}}{(\frac{a}{r})}^{d}\,\cos \,2\theta ,\\ \delta \sigma & = & -\frac{2\mu {\varepsilon }^{\ast }}{1+\frac{\mu }{K}}{(\frac{a}{r})}^{d}\,\cos \,4\theta ,\\ \delta {\sigma }_{xy} & = & -\frac{2\mu {\varepsilon }^{\ast }}{1+\frac{\mu }{K}}{(\frac{a}{r})}^{d}\,\sin \,4\theta ,\end{array}$$

(1)

for $r\gg a$.

The power-law decay with distance r marks the non-local long-range nature of the disturbance that also has distinct angular symmetries in terms of θ. Note the four-fold symmetry in σ which is contrasted by the bi-polar structure in p as illustrated in Fig. 1(a) and (b). The negative (blue) and positive (red) lobes represent regions with decreasing and increasing stresses in space, respectively. It should be noted that the stress patterns in a disordered solid become scattered due to the spatial heterogeneity and/or anisotropy in local elastic properties. It is only beyond a characteristic length-scale that fluctuations average out spatially and the amorphous system, at macroscopic scales, has a smooth elastic response.

Figure 1

Far-field redistribution of (a) axial shear σ (b) pressure p (c) yield function f _y in an elastic matrix with $\frac{\mu }{K}=\frac{1}{2}$ in view of a shear transformation zone. Blue and red denote a decrease and increase in local stresses, respectively. The direction of maximal change in f _y is marked in (c) at ϕ = 65°.

If at a point within the bulk the shear stress on any plane becomes equal to the shear strength, failure will occur at that point. In the following, we will assume that the local shear strength can be expressed as a linear function of the confining pressure²², in accordance with the Mohr-Coulomb concept. Therefore, the distance to failure at any point can be expressed by the yield function f _y :

$${f}_{y}=|{\tau }_{m}|-(p\,\sin \,\varphi \,+\,c\,\cos \,\varphi \mathrm{).}$$

(2)

here c and ϕ are the cohesion and the internal angle of friction respectively, and ${\tau }_{m}=\sqrt{{\sigma }^{2}+{\sigma }_{xy}^{2}}$.

Assuming σ _xy = 0, the far-field perturbations in f _y to leading order following a shear transformation taking place at the origin is given by

$$\delta {f}_{y}=\delta \sigma -\delta p\,\sin \,\varphi \mathrm{.}$$

(3)

Obviously, failure is likely to localize into the overlapped sectors between positive lobes of δσ in Fig. 1(a) and negative lobes of δp in Fig. 1(b). Figure 1(c) displays perturbations in f _y for a value ϕ = 65° of the local friction. Near the shear transformation site, the pattern retains the quadrupolar shape with its positive lobes slightly tilted upwards toward regions with decreasing pressure. Note that a positive δf _y means that material points are pushed toward failure threshold. The direction of the maximum in the positive lobes (denoted by θ_max hereafter) will depend on the internal friction ϕ and is given by

$$\frac{\partial }{\partial \theta }\delta {f}_{y}{|}_{\theta ={\theta }_{{\rm{\max }}}}=0:\,\cos \,2{\theta }_{{\rm{\max }}}=-\frac{1}{4}\,\sin \,\varphi ,$$

(4)

or θ_max ≈ 45° + $\frac{\varphi }{8}$ for low values of ϕ. In the absence of friction, i.e. ϕ = 0, we recover θ_max = 45° as reported in many glassy systems²³. This implies that the change in f _y takes its maximal value along the direction of maximum change in shear stress.

In addition to distortion, we now suppose that transformation zones may undergo dilation as well, i.e. ${\varepsilon }_{\alpha \beta }^{tz}=\frac{1}{2}{\varepsilon }_{v}^{\ast }{a}^{d}{\delta }_{\alpha \beta }$ with dilatancy ${\varepsilon }_{v}^{\ast }$. This will make an anisotropic contribution with the two-fold symmetry to the shear stress $\delta \sigma =-\frac{\mu {\varepsilon }_{v}^{\ast }}{1+\frac{\mu }{K}}{(\frac{a}{r})}^{d}\,\cos \,2\theta$. Combining the effects of the local dilation and shear yields

$$\cos \,2{\theta }_{{\rm{\max }}}=-\frac{1}{4}(\frac{{\varepsilon }_{v}^{\ast }}{2{\varepsilon }^{\ast }}+\,\sin \,\varphi )\mathrm{.}$$

(5)

In comparison with Eq. 4, θ_max in the above equation has an extra ingredient that includes $\frac{{\varepsilon }_{v}^{\ast }}{2{\varepsilon }^{\ast }}$ the ratio between volumetric and shear strains. In the limit $\frac{{\varepsilon }_{v}^{\ast }}{2{\varepsilon }^{\ast }}\to 0$, the term involving friction becomes dominant and we will retrieve Eq. 4.

Macroscopic Failure: Slip-line Instability Analysis

In a highly simplified setting, a fully operating band may be idealized as a quasi-linear object made up of evenly distributed Eshelby events^16,24. Using an eigenmode based strategy, we now propose a simple theoretical argument that makes predictions on the band alignment.

Point-like events with volume a^d and amplitude ε^* are evenly positioned at points ${\overrightarrow{r}}_{i}$ on a line at an angle α with the x axis:

$${\varepsilon }_{\alpha \beta }(\overrightarrow{r})={\varepsilon }_{\alpha \beta }^{\ast }\,{a}^{d}\sum _{i}\,\delta (\overrightarrow{r}-{\overrightarrow{r}}_{i}\mathrm{).}$$

(6)

here δ(...) denotes the delta function.

Taking a Fourier transform, the q-space representation reads

$${\hat{\varepsilon }}_{\alpha \beta }(\overrightarrow{q})=N{a}^{d}{\varepsilon }_{\alpha \beta }^{\ast }\,{\delta }_{\overrightarrow{\hat{q}}.\overrightarrow{n}},$$

(7)

where $\overrightarrow{n}=(\cos \,\alpha ,\,\sin \,\alpha )$indicates the band direction, $\overrightarrow{\hat{q}}$ is the unit vector along $\overrightarrow{q}$, δ is the Kronecker symbol, and N is the total number of events. In order to keep the strain finite in the continuum limit, Na^d is constant with N → ∞ and a^d → 0. Here ${\hat{\varepsilon }}_{\alpha \beta }(\overrightarrow{q})$ takes a localized form as well but at a direction perpendicular to the band orientation in real space. Assuming pure shear distortion for the transformations, i.e. ${\varepsilon }_{\alpha \beta }^{\ast }={\varepsilon }^{\ast }({\delta }_{\alpha x}{\delta }_{\beta x}-{\delta }_{\alpha y}{\delta }_{\beta y})$, it follows that

$$\begin{array}{rcl}\delta p(\overrightarrow{q}) & = & -N{a}^{d}\frac{2\mu {\varepsilon }^{\ast }}{1+\frac{\mu }{K}}\,\cos \,2\theta \,{\delta }_{\overrightarrow{\hat{q}}.\overrightarrow{n}},\\ \delta \sigma (\overrightarrow{q}) & = & +N{a}^{d}\,[2\mu {\varepsilon }^{\ast }-\frac{2\mu {\varepsilon }^{\ast }}{1+\frac{\mu }{K}}\frac{1}{2}(1+\,\cos \,4\theta )]{\delta }_{\overrightarrow{\hat{q}}.\overrightarrow{n}},\end{array}$$

(8)

with θ the wave vector angle. The above equations follow from operating the Oseen tensor on the effective source field defined by the line of shear transformations (see Supplementary Materials).

It is clear that the solutions in Fourier space will only depend on the direction θ (not on $|\overrightarrow{q}|$) and take the exact same form as the effective source in Eq. 7, up to a normalization factor. In other words, the slip line does not induce any stress redistribution outside of the line itself, and the real-space response function is localized along the slip line at angle α

$$\begin{array}{rcl}\delta p(\overrightarrow{r}) & = & \frac{2\mu {\varepsilon }^{\ast }}{1+\frac{\mu }{K}}\,\cos \,2\alpha \,{a}^{d}\,\sum _{i}\,\delta (\overrightarrow{r}-{\overrightarrow{r}}_{i})\\ \delta \sigma (\overrightarrow{r}) & = & [2\mu {\varepsilon }^{\ast }-\frac{2\mu {\varepsilon }^{\ast }}{1+\frac{\mu }{K}}\frac{1}{2}(1+\,\cos \,4\alpha )]{a}^{d}\sum _{i}\delta (\overrightarrow{r}-{\overrightarrow{r}}_{i}).\end{array}$$

(9)

As a result, the change in the yield function δf _y is zero everywhere except on the slip line. We now argue that the most unstable lines will be those that maximize δf _y , which is positive inside the band. This corresponds to a maximum amplification of the local “damage” that is caused by the operating shear band for an infinitesimal strain. The shear band angle θ_sh is thus obtained by maximizing δf _y with respect to α:

$$\frac{{\rm{\partial }}}{{\rm{\partial }}\alpha }\delta {f}_{y}{|}_{\alpha ={\theta }_{{\rm{s}}{\rm{h}}}}=0:\,\cos \,2{\theta }_{{\rm{s}}{\rm{h}}}=-\frac{1}{2}\,\sin \,\varphi .$$

(10)

This angle lies between 45° ≤ θ_sh ≤ 60°, and for small ϕ is approximately θ_sh ≈ 45° + $\frac{\varphi }{4}$. The combination of the pure shear deformation and volumetric change in the transformation zones gives

$$\cos \,2{\theta }_{{\rm{sh}}}=-\frac{1}{2}(\frac{{\varepsilon }_{v}^{\ast }}{2{\varepsilon }^{\ast }}+\,\sin \,\varphi )\mathrm{.}$$

(11)

Note that a similar result in the case without friction was obtained in ref.¹⁶ on the basis of an energy minimization argument. Inserting Eq. 10 into Eq. 9 and integrating δp over the active domain, it becomes evident that the total change in pressure is negative. This will effectively reduce the shear strength given by the Mohr-Coulomb yield surface inside the localization zone. With such a mechanism, yielded zones become likely places where next events tend to localize, hence making possible a permanent band-like formation. This is compatible with the numerical observations presented in the following sections.

From our discussions above it follows that, theoretically, θ_max, the direction resulting from expected elastic-type correlations, will differ from the slip-line angle θ_sh. This discrepancy appears to be consistent with the observations made in the granular experiment of Le Bouil et al.¹⁴. Note that the interpretation proposed by these authors was initially based on purely elastic considerations, with no allusion to the local friction. In this analysis, a large local dilatancy had to be assumed in order to make a sensible prediction of θ_max. Our approach, by introducing the local friction angle, allows one to obtain a similar order of magnitude without invoking a large dilation. In a follow-up work²⁵, it was suggested that the material anisotropy (or its combination with the volumetric strain) may explain the characteristic angle. The latter assumption was largely attributed to the presence of force chains that build up upon loading granular solids²⁶. Our interpretation is somewhat different in that it is entirely based upon the dominant role of the friction angle on correlation patterns, and in addition makes a specific prediction for the orientation of the permanent shear bands.

Simulation Details

In order to test the predictions of the theoretical discussion above, we have made use of the Finite Elements based version of elasto-plastic models that we originally established in^6,27. In our previous work, this model was used to study the deformation of amorphous media in which the local failure is governed by a standard, maximum shear stress criterion. Here, in order to follow the hypothesis of the previous section, the microscopic failure criterion for each element will be a Mohr-Coulomb condition. Each material point is therefore assigned with the shear strength parameters c and ϕ; while the former is randomly chosen from an exponential distribution with the mean value $\bar{c}$, we allocate no disorder to the latter.

Below the failure limit at the material point in question, set by f _y = 0 in the pressure-shear stress plane, the local stress trajectory is regulated by the imposed drift and non-local elastic interactions. We further presume a linear isotropic elastic response for the pre-failure dynamics. Figure 2 plots τ _m = $\frac{1}{2}$ (σ₃ − σ₁) against $p=-\frac{1}{2}$(σ₃ + σ₁) and represents any stress state by a stress point. Here σ₁ and σ₃ denote major and minor principle stresses, respectively. Upon yielding, the stress point will take on a path which is perpendicular to the p-axis and relax visco-elastically⁶ toward a point representing the final state of stress. The released shear stress Δσ will create a localized net force that perturbs the force equilibrium in the medium. The perturbation will be of the generic form given in Eq. 1 plus a scattering term caused by currently yielded elements with decaying shear moduli²⁸. The homogeneity approximation and subsequent theoretical derivations in the preceding section appear to be legitimate at least with regard to the fluctuation patterns and shear band morphology evaluated numerically in the following sections.

Figure 2

The Mohr-Coulomb failure envelope. The red curve and white dots sketch the stress trajectory and state of stress at a material point of interest.

Simulations of shear deformations were performed by applying an area preserving axial shear rate $\dot{\varepsilon }\sim {10}^{-4}$ to an L × L periodic cell. An irregular set of triangular elements with average size $h=\frac{L}{80}$ was used to discretize the domain. The rate unit (inverse timescale) is set by the shear wave velocity ${c}_{s}={(\frac{\mu }{\rho })}^{\frac{1}{2}}$ divided by L with ρ the mass density. The rate of strain chosen is slow enough to ensure quasi-static conditions, i.e. the results are insensitive to a reduction in the value of the rate. We also set $\frac{K}{\mu }=2$ in the elastic regime, corresponding to a Poisson ratio of ν ≈ 0.33.

By setting a high value of the damping rate (in comparison with the vibrational frequency), we moreover impose an overdamped behavior during the relaxation phases, thus ensuring that inertial effects are negligible⁶. Prior to shearing, samples were prepared with random stresses assigned to each block followed by an equilibration within a purely elastic framework (no plastic events allowed) that resulted in a state of mechanical equilibrium.

Numerical Results

The results of shearing tests can be used to determine the bulk shear strength characteristics together with the structure of interactions between transient slip events. A number of tests has been performed at the same (initial) pressure $p=8\bar{c}$ and ϕ = 65°, and the resulting average stress strain curve is displayed in Fig. 3(a). For the given p and ϕ, the material shear strength has the average value of ${\sigma }_{y}\approx 8\bar{c}$. We also report ε in units of the yield strain ${\varepsilon }_{y}=\frac{{\sigma }_{y}}{{2}_{\mu }}$. A peak stress is typically reached around ε ≈ ε _y in every sample followed by a reduction in strength as the loading continues. With increasing strain, the shear strength ultimately falls to a residual value at large deformations.

Figure 3

Results of bi-axial shearing tests. (a) Stress-strain curves for three samples (denoted by different symbols) tested under the same value of pressure $p=8\bar{c}$ and ϕ = 65°. The shear stress σ is measured in units of $\bar{c}$. (b–e) Spatial maps of active sites illustrated by solid dots at $\frac{\varepsilon }{{\varepsilon }_{y}}=1.25,1.6,2.5,$ and 3.7. Each map corresponds to the full system size L.

A strong plastic activity in the form of extended linear structures is present at all times after the initial yielding for ε > ε _y . Spatial maps of active sites in Fig. 3(b–e) demonstrate the highly intermittent and non-local nature of bursts during plastic flow (see Supplementary movie). Following the stress peak, multiple shearing bands start to evolve within which most of plastic activities are taking place. The orientation of these bands are quite distinct from maximum shearing directions, 45° or 135° in our loading set-up. The quantification of the slip directions will be the subject of the next section, which presents an analysis similar in spirit to²⁹.

Structural Characterization

Here we focus on the number density of active sites and associated spatial correlations so as to quantify the spatial structure of the plastic activity. Let $\rho (\overrightarrow{r})={\sum }_{i}\delta (\overrightarrow{r}-{\overrightarrow{r}}_{i})$ where ${\overrightarrow{r}}_{i}$ is the position vector of a plastically active zone with index i and δ(...) denotes a delta function. Integrating $\rho (\overrightarrow{r})$ over the entire volume V provides the total number of plastic sites N at the current time or strain. We are interested in the variations of $\rho (\overrightarrow{r})$ about its average value $\rho =\frac{N}{V}$. These variations can be characterized by a two-point density correlation function between two different positions $\overrightarrow{r}$ and $\overrightarrow{r}^{\prime}$ in space: $S(\overrightarrow{r},\overrightarrow{r}^{\prime} )\doteq \langle [\rho (\overrightarrow{r})-\rho ][\rho (\overrightarrow{r}^{\prime} )-\rho ]\rangle$. Here the angular brackets 〈...〉 correspond to an average over different realizations and a spatial average. As a result of translational invariance $S(\overrightarrow{r},\overrightarrow{r}^{\prime} )\equiv S(\overrightarrow{r}-\overrightarrow{r}^{\prime} )$. A Voronoi cell analysis was performed that enabled interpolations of the density field onto fine regular grids. We subsequently used a 2 d Fourier transform in order to compute the correlations.

We shall naturally expect that density fluctuations are highly correlated along shear band directions, inducing strong anisotropies in $S(\overrightarrow{r}-\overrightarrow{r}^{\prime} )$. Figure 4 displays the evolution of density correlations at different loading stages averaged over 16 samples at $p=8\bar{c}$ and ϕ = 65°. Angular symmetries seen in the post-failure regime are fairly stable features showing only weak fluctuations with increasing strain. The banded regions contain correlations that are explicitly longer-ranged (as opposed to other orientations) indicating the system-spanning nature of the localization.

Figure 4

Density correlation function $S(\overrightarrow{r}-\overrightarrow{r}^{\prime} )$ at multiple strains (a) $\frac{\varepsilon }{{\varepsilon }_{y}}=1.4$ (b) $\frac{\varepsilon }{{\varepsilon }_{y}}=2.5$ (c) $\frac{\varepsilon }{{\varepsilon }_{y}}=3.7$. Here $p=8\bar{c}$ and ϕ = 65°. The scale of each density map is L. The solid and dashed-dotted lines indicate the theoretical predictions for θ_sh and θ_max, respectively.

Figure 5 displays the anisotropic part S(θ), an averaged $S(\overrightarrow{r}-\overrightarrow{r}^{\prime} )$ over different distances $|\overrightarrow{r}-\overrightarrow{r}^{\prime} |$. There are two marked maxima in every data set at 45° < θ < 135° that delineate the degree of anisotropy. We quantify the positions of these peaks by fitting the data to a sum of two Gaussian peaks, as in Fig. 5(b–e). These fits enable us to obtain the locations of the peaks–denoted by θ_sh– and to follow their evolution upon shear loading. This is shown in Fig. 5(a). Apart from the initial transient part prior to failure, the peak locations are essentially strain independent at larger strains. The dashed lines in Fig. 5(a) mark the directions of failure given by Eq. 4.

Figure 5

Anisotropic part of the density correlation function S(θ) at $p=8\bar{c}$ and ϕ = 65°. (a) Locations of the peaks in S(θ) denoted by θ_sh versus ε/ε _y . (b–e) S(θ) plotted against θ at $\frac{\varepsilon }{{\varepsilon }_{y}}=1.4,1.6,2.5,$ and 3.7. The dashed lines in the main plot designate theoretical predictions discussed in Sec. 1.2. Gaussian fits in the insets are denoted by the dashed curves.

The basis for the permanent localization observed here can be understood in simple terms. An incident avalanche tends to depressurize currently-damaged blocks, with a pressure drop proportional to sin ϕ, making them vulnerable spots against further deformations. In the standard framework of elastoplastic models, it is possible to observe a similar behavior (however with the classical 45° orientation for the shear band by adding some permanent or transient weakening mechanism to the model). Examples include models based on the damage factor¹⁸, weakened stress thresholds¹⁹, or lingered restoration time²¹, just to name a few. In our simulations, the weakening is intrinsically contained in the pressure sensitivity of the failure criterion, and depends on the local friction angle. We note that a relatively large friction angle ϕ was required to observe a permanent localization.

Mohr-Coulomb Failure Envelope

Testing several samples each under different confining pressures enables the determination of a macroscopic failure envelope and the bulk shear strength parameters. Figure 6 illustrates the results in a series of tests performed at $p/\bar{c}=4,8,12,$ and 16 and ϕ = 65°. Every set of tests was carried out on 16 independent samples. The data points in Fig. 6(b) represent the states of stress on the shear-pressure plane. These stress points correspond to the ultimate strength, marked by the symbols in Fig. 6(a). Very similar results would be obtained by using the peak values in the stress strain curves. Assuming a global Mohr-Coulomb criterion, the residual strengths are expressed by

$${\sigma }_{r}=p\,\sin \,{\varphi }_{r}\,+\,{c}_{r}\,\cos \,{\varphi }_{r},$$

(12)

where ϕ _r denotes the value of the bulk friction angle and c _r is the macroscopic cohesion.

Figure 6

Measurement of the bulk shear parameter ϕ _r . (a) stress-strain curves for multiple values of the confining pressure p and ϕ = 65°. The symbols lay out stress points used to draw the Mohr-Coulomb failure envelope. (b) stress points commensurate with the residual strengths σ _r plotted on the shear-pressure plane. The lines are Mohr-Coulomb linear fits with slopes sin ϕ _r and sin ϕ.

The parameters can be determined empirically by making linear fits to the points. The values of the measured parameters are:

$${\varphi }_{r}={40}^{\circ }\mathrm{;\ }{c}_{r}\,=\,0.5\,\bar{c}\mathrm{.}$$

As stipulated by the Mohr-Coulomb phenomenology, the theoretical angle between the major principal stress direction and the plane of failure must be ${\theta }_{{\rm{MC}}}={45}^{\circ }+\frac{{\varphi }_{r}}{2}$. Inserting ϕ _r = 40°, it follows that θ_MC = 65° which is off from θ_sh by about 6°. Therefore, we find that applying the macroscopic Mohr-Coulomb theory using the observed values of the macroscopic failure envelope tends to overestimate the actual inclination of the yield surface, in agreement with experimental observations³⁰.

Pre-failure Patterns of Plastic Activity

Beyond the yield strain, correlations in plastic activity are easily identified based on the analysis of a snapshot taken at a single time or strain, as described above using the function $S(\overrightarrow{r}-\overrightarrow{r}^{\prime} )$. This is related to the existence of well correlated linear regions of plastic activity that are operating simultaneously, eventually giving rise to localized shear bands. In the pre-failure regime, however, plastic activity is much more scattered, and the study of a single configuration does not reveal any established pattern. A calculation of $S(\overrightarrow{r}-\overrightarrow{r}^{\prime} )$ does not allow one to identify any preferred direction. The structure of the plastic activity, however, can be revealed by the study of two time correlations between configurations separated by a fixed strain interval Δε. In practice, the statistics is improved by averaging over different simulation and performing an average over a strain window for the initial configuration. We define

$$C(\overrightarrow{r}-\overrightarrow{r}^{\prime} ,{\rm{\Delta }}\varepsilon )\doteq \langle \frac{1}{{\varepsilon }_{2}-{\varepsilon }_{1}}{\int }_{{\varepsilon }_{1}}^{{\varepsilon }_{2}}d\varepsilon \rho (\overrightarrow{r},\varepsilon +{\rm{\Delta }}\varepsilon )\rho (\overrightarrow{r}^{\prime} ,\varepsilon )\rangle \mathrm{.}$$

(13)

here the brackets denote the average over different realizations. The averaging interval [ε₁,ε₂] is taken before the strain peak, ε₁ = 0.25ε _y and ε₂ = 0.7ε _y , in order to avoid the contamination of the correlation functions by the formation of permanent shear bands in the post yield regime.

Statistically speaking, $C(\overrightarrow{r},{\rm{\Delta }}\varepsilon )$ will provide spatial details about the most likely position of an event that is triggered following a local slip event at the origin. In Fig. 7(a), correlations are displayed for a small strain difference Δε = 5 × 10⁻⁴ε _y . The results exhibit a four-fold structure, which persists up to a strain interval of the order 10⁻³ before fluctuations become uncorrelated at higher strain differences. The four-fold structure is of course expected from standard elasticity theory, and has been observed in experiments and simulations of a number of glassy systems. However, we observe here a marked deviation from the 45° that would be predicted by pure elasticity, reflecting the influence of the friction angle in the yield criterion.

Figure 7

(a) $C(\overrightarrow{r},{\rm{\Delta }}\varepsilon =5\times {10}^{-4}{\varepsilon }_{y})$ at $p=8\bar{c}$ and ϕ = 65°. The dashed-dotted lines indicate the theoretical prediction for θ_max. Here the lengthscale is $\frac{L}{8}$. (b) Anisotropic portion C(θ) plotted against θ. The Gaussian fit is represented by the dashed curve.

In Fig. 7(b), the angular dependence of $C(\overrightarrow{r},{\rm{\Delta }}\varepsilon )$ is shown together with a fit to a Gaussian function, similar to the one used for the analysis of S(r, θ). The fit function has a peak centered around 52° which is in close agreement with the theoretical prediction θ_max. Note that in order to reduce the noise in the data, the integration over r was limited to small values of the distance $r < \frac{L}{16}$.

The structure of correlations preceding the ultimate failure reflects the hypothesis of our model. A localized event induces stress fluctuations in the surroundings triggering plastic rearrangements nearby in the medium. events are preferentially triggered in directions in which the changes to the yield function are maximal, as a result from non-local elastic couplings combined with a local frictional yielding rule. While one may have expected that these fluctuations would strongly influence the collective behavior that emerges upon macroscopic failure, the analysis of the shear band orientations shows that they represent a distinct phenomenon, and must be analyzed within a more collective perspective.

Conclusion and Discussion

In this work, we have taken a coarse-grained, mesoscopic view of the deformation in a disordered granular medium, by integrating the notion of a deformation taking place through local shear transformations with that of a local criterion of failure described by the Mohr-Coulomb friction condition. As a result, the mechanical response is described by two ingredients: (i) the elastic moduli of the medium which, according to the Eshelby’s picture, describe the response to a localized shear transformation and (ii) the local friction angle that characterizes the Mohr-Coulomb condition. The latter feature is, in fact, material specific and supplements Le Bouil et al.’s line of reasoning which relies entirely on non-local couplings mediated by elasticity¹⁴.

Based on this picture, we proposed a scenario in which spatial correlations of the scattered plastic activity that takes place before the yield point, and the spatial structure of the permanent failure planes that can be observed beyond the yield point are explicitly calculated as a function of these simple ingredients, but occur at different directions. The discrepancy between the morphology of permanent bands and transient correlations is in full agreement with Le Bouil et al.¹⁴, who proposed a scenario based on the phase-coexistence between correlated mini bursts and persistent shearing bands. Within this picture, the flowing, post-yielding state is intrinsically different from the states explored at small strains, and the yield transition appears similar to a first order, discontinuous transition.

The failure mechanism we have proposed can be interpreted as the maximal instability of uniform lines of slip. Localization of deformation in this form is favorable, as no further stress is built up in the surrounding medium²¹. Furthermore, the intense shearing effectively lowers the yielding strength inside the localization band giving rise to a weakening process. The maximally unstable modes are dictated by both the elastic interactions³¹ and the Mohr-Coulomb plasticity criterion (more specifically local friction) and may be viewed as a major source of mechanical instability³². It is noteworthy that the incurred damage is an ingredient of the model through the failure criterion, and does not enter as an extra material parameter.

The applicability of our analysis to a real granular medium may be questioned on two accounts: firstly, one may argue that a real granular medium is not described by an linear elastic continuum, due to the existence of force chains and micro plasticity in the form of contact breaking and formation. Second, the very existence of a Mohr-Coulomb criterion at the local scale is a rather arbitrary assumption, although in general a pressure sensitivity of the shear transformation is expected. Therefore, we have numerically tested our proposition by building a lattice-based model that incorporates rigorously, if perhaps artificially, these ingredients. The results of these simulations are in good agreement with the theoretical expectations, and, when compared to experiments, can lead to a prediction of an effective local friction coefficient.

In addition to confirming the theoretical analysis, the simulations offer the possibility to perform an analysis of the macroscopic failure envelope. This analysis yields an orientation, labeled as θ_MC, that is incompatible with θ_sh, the shear inclination. While the former is purely based on the measured bulk stresses, the later arises from kinematic considerations only. This discrepancy is in contradiction with the Mohr-Coulomb phenomenology at the global scale. Similar experimental observations were made (see³⁰ and the references herein) leading to ad-hoc remedies that incorporate extra material parameters, i.e. the dilatancy angle³³, into the theoretical framework.