A non-local rheology for granular flows across yield conditions

A non-local rheology for granular flows across yield conditions

Mehdi Bouzid    Martin Trulsson    Philippe Claudin    Eric Clément    Bruno Andreotti Physique et Mécanique des Milieux Hétérogènes, PMMH UMR 7636 ESPCI – CNRS – Univ. Paris-Diderot – Univ. P.M. Curie, 10 rue Vauquelin, 75005 Paris, France
July 5, 2019

The rheology of dense granular flows is studied numerically in a shear cell controlled at constant pressure and shear stress, confined between two granular shear flows. We show that a liquid state can be achieved even far below the yield stress, whose flow can be described with the same rheology as above the yield stress. A non-local constitutive relation is derived from dimensional analysis through a gradient expansion and calibrated using the spatial relaxation of velocity profiles observed under homogeneous stresses. Both for frictional and frictionless grains, the relaxation length is found to diverge as the inverse square-root of the distance to the yield point, on both sides of that point.


Granular materials belong to the class of amorphous athermal systems. Like foams BW90 (); Kabla2003 (), emulsions CCSB09 (), suspensions BDCL10 (); BGP11 (); TAC12 () or metallic glasses Xie08 (), they exhibit a dynamical phase transition between static and flowing states. Analogously to phase transitions of thermodynamic systems, this rigidity transition exhibits a divergence of correlation lengths WNW05 (); HBB10 (), revealing the presence of non-local cooperative processes called dynamical heterogeneities DJvH (). In order to describe the constitutive behavior of such systems, it is natural to adopt the Ginzburg-Landau phenomenological approach of phase transitions AT06 (); GCOAB08 (); BCA09 (); KK12 (); HK13 (). The main issue is then to identify the relevant control and order parameters. Following the now classical Liu-Nagel diagram for jamming transition LL98 () – or a revised version BZCB11 () – it is usually assumed that the solid-liquid mechanical transition is controlled by the shear stress GCOAB08 (); BCA09 (); KK12 (), which, once rescaled by its critical value, defines the yield parameter . For a granular system sheared under a fixed confining pressure , one defines the dimensionless Coulomb yield parameter , where is the friction coefficient in the zero shear-rate limit, at the jamming volume fraction .

Figure 1: (Color online) (a) Schematic of the numerical set-up. The walls, composed by the dark purple grains, are submitted to a confining pressure . In the buffer layers located close to the walls, the grains (orange) are submitted to gravity-like forces along the transverse direction . (b) Schematic profiles of the pressure (black line) and of the shear stress (red line) across the cell. In the bulk of the shear cell (white), the pressure is homogeneous: .

In a series of recent papers BCA09 (); GCOAB08 (); CMCB12 (); KK12 (); HK13 () the order parameter is a rheological quantity called the fluidity, proportional to the inverse viscosity i.e. to the ratio of the shear rate to the shear stress . Here, we consider that the relevant order parameter must be a dimensionless quantity based exclusively on state variables (which excludes ) like the shear rate rescaled by a microscopic timescale. In granular materials, the only energy scale is set by the confining pressure , so that the order parameter must be the inertial number


based on the grain diameter and on their density . compares to the microscopic rearrangement time . Considering an incompressible homogeneous flow, it can be inferred that the yield parameter is a function of noted GDRMidi (); CEPRC05 (); JFP06 (). If this local constitutive relation was still valid in heterogeneous flows, the transition between solid () and liquid () states would systematically occur at . However, different experiments have shown that the stress at a location depends on the shear rate around this point, a property called non-locality. (i) In the inclined plane geometry, thin granular layers flow anomalously DLDA06 () and stop at a yield parameter P99 (); GDRMidi (). (ii) A creeping flow is commonly observed in regions which are expected to be jammed (i.e. solid), since KINN01 (); GDRMidi (); BAU11 (). (iii) A solid plunged in grains and submitted to a force lower than the yield threshold starts moving as soon as a shear band is created far away from the solid NZBWH10 (); RFP11 ().

In this letter, we show that the liquid state continuously extends from liquid zones () into the bulk of regions that are below the yield conditions (). We furthermore find that the rheology obeys the very same non-local constitutive relation across yield conditions.

Three different pictures have emerged so far to explain non-locality A07 (). In soft amorphous systems, like foams, emulsions or glassy Lennard-Jones phases, the dynamics in the quasi-static regime is controlled by elasto-plastic events TLB06 (); LC09 (); BCA09 (); LP09 (): when sheared, energy is slowly stored and rapidly released through scale-free avalanches, in close analogy with the depinning transition of an elastic line. By contrast, the dynamics of hard non-deformable grains is essentially related to geometry: elementary plastic events are rather identified as the rapid formation of force chains followed by a slow zig-zag instability of these structures LDW12 (). Non-locality can then be related to soft modes, by essence spread in space, prescribing the cooperative motion of the particles ABH12 (). In this geometrical picture, the relevant state parameter would rather be the mean number of contacts per particle or the volume fraction H10 (). The third picture is based on an analogy with Eyring’s transition state theory for the viscosity of liquids PF09 (), where mechanical fluctuations would play the role of temperature in thermal systems. Here, we show that the non-local constitutive relation for dense granular flows can be determined from simple phenomenological assumptions, regardless the nature of the relevant dynamical mechanisms. We calibrate and test it by means of discrete element simulations.

Figure 2: (Color online) (a) Typical profiles of the yield parameter obtained numerically below (red circles, ) and above (blue circles, ) yield conditions for frictional grains. (b) Corresponding velocity profiles. Green dashed lines: predictions of the local rheology. Black solid lines: best fits by Eq. 2.

Numerical set-up — We have performed molecular dynamics simulations of massive grains, confined in a shear cell under an imposed stress field. The system is two-dimensional and constituted of circular particles of mean diameter , with a 20% polydispersity. The shear cell is composed of two rough walls moving along the -direction with opposite velocities (see Fig. 1a for notations). These walls are made of the same grains, but glued together. The walls are separated by . Their position is controlled to ensure a constant normal stress . The cell thickness then fluctuates over a fraction of grain diameter. Periodic boundary conditions are applied along the -direction. The particle and wall dynamics are integrated using the Verlet algorithm. Contact forces between particles are modeled as linear viscoelastic forces, chosen such that the restitution coefficient is . For each measurement reported here, we have varied the normal spring constant and we report the value of the plateau at large . This rigid asymptotic regime is reached in practice for . To model actual granular samples, we have considered frictional grains, which interact along the tangential direction by a Coulomb friction Cundall79 (); daCruz05 (); Luding06 () of coefficient , with a tangential spring constant . For the sake of comparison, we have also studied the same system with frictionless grains ().

What makes the set-up original is the possibility of imposing the profile of by means of gravity-like forces applied to the grains located in two buffer zones (Fig. 1a) of thickness . A grain labelled , of mass , and located at is submitted to an external force: , where is the amplitude of the localized ”gravity” field. These forces are oriented downward at the top of the cell, and upward at the bottom (Fig. 1a). The resulting pressure (Fig. 1b) starts from at the upper wall, increases due to the sum of forces applied in the buffer zone, and reaches a constant value in a central region of width , called the bulk. decreases back to at the lower boundary. By contrast, the shear stress profile is homogeneous across the cell (Fig. 1b). As a result, the stresses are homogeneous in the bulk of the shear cell: . By tuning the amplitude , can be imposed smaller or larger than . is chosen larger than so that the buffer zone remains above yield conditions, at .

Local rheology — For given stress conditions, the simulation is ran until a steady state is reached and the averaged velocity profiles are then measured. Fig. 2 shows two such profiles, one above and the other below yield conditions. One observes that the entire system always flows, even when . In the bulk, the velocity deviates from the linear profile predicted by the local rheology (as is constant, one expects to be constant as well and to vanish for ). The velocity rather tends exponentially towards such a linear profile, which suggests a linear relaxation in space. The velocity profiles inside the bulk zone (Fig. 2) are accordingly fitted with the function:


The velocity is inherited from the buffer layer, while the asymptotic shear rate (which is not the shear rate at the center of the cell) and the relaxation length are two adjustable parameters. Varying in the buffer layer but keeping constant, we systematically measured the same values of and . These two quantities characterize the bulk state and do not depend on the buffer layer characteristics.

The asymptotic shear rate provides the proper way of defining the local rheology . is indeed the inertial number selected for a certain ratio in homogeneous shear and stress conditions. It is deduced in practice from the fit of the data to Eq. 2. The resulting constitutive relations are reported in Fig. 3 for the frictional and the frictionless case. In both cases, the data is perfectly described by a law of the form


in the accessible range of (between and ): the residuals form a statistical noise. The exponent is in the frictionless case and in the frictional case, within error bars (), see also CEPRC05 (); PR08 (). We hypothesize that this difference is related to another fundamental difference between the two systems. In the frictionless situation, the jamming point ( and ) coincides with the isostatic point, while for frictional grains, the jamming point is far in the hyper static zone WNW05 (); H10 ().

Figure 3: Local rheology deduced from the fit of the bulk velocity profile to Eq. 2. Data for frictionless (a) and frictional (b) grains. Solid lines: the best fit to Eq. (3) gives for the frictionless case and for the frictional one. Bottom: as a function of in log scales for the frictionless (c) and frictional (d) cases.
Figure 4: (Color online) (a) Relaxation length as a function of , below (circles) and above (squares) yield conditions. (a) Green and yellow symbols: data for frictionless grains. (c) Blue and red symbols: data for frictional grains. The solid lines are the best fit by Eqs. (5) and (6). The values of found for frictionless () and frictional () systems are remarkably similar. (b) and (d) Log-log plots of the same quantities, revealing the divergence with an exponent .

Non-local rheology —The relaxation length is displayed in Fig. 4 as a function of . The numerical data for both frictionless and frictional systems are qualitatively similar: diverges on both sides of the critical point with an exponent . In order to account for non-local effects in the theoretical constitutive relation, we perform a gradient expansion of the functional . Assuming that non-locality results from a statistically isotropic short-range interaction between shear zones, the lowest order operator is the Laplacian . As a direct consequence, and its gradient must be continuous. Furthermore, we assume that the correction remains finite as , so that the expansion must be expressed in terms of . At the linear order in , the constitutive relation writes:


where is a phenomenological constant. is positive when the point considered is surrounded by a more liquid region (higher ). This region flows more easily than expected from the local value of , so that the corresponding shear stress is lower. is therefore positive. Importantly, our derivation does not depend on the nature of the mechanical interaction between shear zones; the reader may think of the analogy with the van der Waals gradient expansion of the Helmholtz free energy at a liquid-vapour interface Rowlinson ().

Above yielding conditions, the linearization of Eq. 4 around the bulk inertial number gives at first order a differential equation of the form , whose solutions are exponentials with a relaxation length


Below yielding conditions, as , the non-local correction is of zeroth order and Eq. 4 leads to . This gives a similar differential equation but now with a divergence of the form


As shown in Fig. 4, the measured relaxation length effectively diverges on both sides of the critical point according to the theoretical predictions (5, 6). In particular, both frictional and frictionless systems exhibit a power-law divergence as ; the multiplicative factor is the same above and below for the frictional case but differ by in the frictionless case (Fig. 4b and 4d). Note that a similar scaling was found in BCA09 () for a fluidity correlation length in a kinetic elastoplastic model.

Discussion — It must be emphasized that the creeping regime and in the fully flowing regime ) are described as a single liquid phase. The divergence of the relaxation length is indeed predicted across the yield condition by the very same non-local correction to the liquid constitutive relation. Due to the liquid boundary condition imposed by the buffer layers, the bulk is always found flowing, regardless of the stress in this zone. The liquid state persists asymptotically into the bulk, far from any direct influence of the boundary. This is consistent with our choice of order parameter: is everywhere non zero and the same constitutive relation holds in all layers. In conclusion, the liquid-solid transition is not controlled by . For smaller than static threshold value the system can be solid or fluid. If fluid then varies in space according to the non-local rheology. If solid, the stress state must be described by another constitutive relation based on elasticity.

Our derivation, based on a gradient expansion of the yield parameter , written as a functional of the order parameter , does not prejudge of any dynamical mechanisms at work at the microscopic level. The exponential behaviors and the associated length-scales hence identified are direct consequences of the linearization around the critical state. Therefore, finer investigations at the grain level must be carried out to understand the connections between the three lines of thought currently invoked to explain such a non-local rheological coupling, namely elasto-plastic TLB06 (); LC09 (); BCA09 (); LP09 (), geometrical WNW05 (); LDW12 () and stress-mediated activation PF09 (); RFP11 ().

BA is supported by Institut Universitaire de France. This work is funded by the ANR JamVibe.


  • (1) F. Bolton and D. Weaire, Phys. Rev. Lett. 65, 3449 (1990).
  • (2) A. Kabla and G. Debrégeas, Phys. Rev. Lett. 90, 258303 (2003)
  • (3) M. Clusel, E.I. Corwin, A.O.N. Siemen and J. Brujic, Nature 460, 611 (2009) .
  • (4) C. Bonnoit, T. Darnige, E. Clément and A. Lindner, J. Rheol. 54, 65 (2010)
  • (5) F. Boyer, E. Guazzelli and O. Pouliquen, Phys. Rev. Lett. 107, 188301 (2011)
  • (6) M. Trulsson, B. Andreotti and P. Claudin, Phys. Rev. Lett. 109, 118305 (2012).
  • (7) S. Xie and E.P. George, Acta Mater. 56, 5202 (2008).
  • (8) M. Wyart, S. R. Nagel and T. A. Witten, Europhys. Lett. 72, 486 (2005).
  • (9) C. Heussinger, L. Berthier and J.-L. Barrat, Europhys. Lett. 90, 20005 (2010).
  • (10) O. Dauchot, D.J. Durian, M. van Hecke, in “Dynamical heterogeneities in glasses, colloids, and granular media”, (Oxford University Press) (2011).
  • (11) I.S. Aranson and L.S. and Tsimring, Rev. Mod. Phys. 78, 641 (2006).
  • (12) L. Bocquet, A. Colin and A. Ajdari, Phys. Rev. Lett. 103, 036001 (2009).
  • (13) J. Goyon, A. Colin, G. Ovarlez, A. Ajdari and L. Bocquet, Nature. 454, 84 (2008).
  • (14) P. Chaudhuri, V. Mansard, A. Colin and L. Bocquet, Phys. Rev. Lett. 109, 036001 (2012).
  • (15) K. Kamrin and G. Koval Phys. Rev. Lett. 108, 178301 (2012).
  • (16) D.L. Henann and K. Kamrin, Proc. Natl. Acad. Sci. USA 110, 6730 (2013).
  • (17) A.J. Liu and S.R. Nagel, Nature 396, 21 (1998).
  • (18) D. Bi, J. Zhand, B. Chakraborty and R.P. Behringer, Nature 380, 355 (2011).
  • (19) S. Deboeuf, E. Lajeunesse, O. Dauchot and B. Andreotti, Phys. Rev. Lett. 97, 158303 (2006).
  • (20) O. Pouliquen, Phys. Fluids. 11, 542 (1999).
  • (21) GDR MiDI, Eur. Phys. J. E. 14, 341 (2004).
  • (22) F. da Cruz, S. Emam, M. Prochnow, J.N. Roux and F. Chevoir, Phys. Rev. E 72, 021309 (2005).
  • (23) P. Jop, Y. Forterre, and O. Pouliquen, Nature. 441, 727 (2006).
  • (24) T. S. Komatsu, S. Inagaki, N. Nakagawa and S. Nasuno, Phys. Rev. Lett. 86, 1757 (2001).
  • (25) V.B Nguyen, T. Darnige, A. Bruand and E. Clément, Phys. Rev. Lett. 107, 138303 (2011).
  • (26) K. A. Reddy, Y. Forterre and O. Pouliquen, Phys. Rev. Lett. 106, 108301 (2011).
  • (27) K. Nichol, A. Zanin, R. Bastien, E. Wandersman and M. van Hecke, Phys. Rev. Lett. 104, 078302 (2010).
  • (28) B. Andreotti Eur. Phys. Lett., 79 (2007) 34001.
  • (29) A. Tanguy, F. Leonforte & J.-L. Barrat Eur. Phys. J. E 20, 355 (2006).
  • (30) A. Lemaître & C. Caroli Phys. Rev. Lett. 103, 065501 (2009).
  • (31) E. Lerner and I. Procaccia, Phys. Rev. E 79, 066109 (2009).
  • (32) E. Lerner, G. Düring and M. Wyart, Proc. Natl. Acad. Sci. USA 109, 4798 (2012).
  • (33) B. Andreotti, J.-L. Barrat, and C. Heussinger, Phys. Rev. Lett. 109, 105901 (2012).
  • (34) M. van Hecke, J. Phys. Cond. Matt. 22, 033101 (2010).
  • (35) O. Pouliquen and Y. Forterre, Phil. Trans. R. Soc. A. 367, 5091(2009).
  • (36) P. A. Cundall and O. D. L. Strack, Geotechnique 29, 47 (1979).
  • (37) F. da Cruz, S. Emam, M. Prochnow, J.N. Roux, and F. Chevoir, Phys. Rev. E 72, 021309 (2005).
  • (38) S. Luding, Behavior of Granular Media, 137-147, Shaker Verlag, Aachen (2006).
  • (39) P.-E. Peyneau and J.-N. Roux, Phys. Rev. E 78, 011307 (2008).
  • (40) Rowlinson JS, Widom B. 1982. Molecular Theory of Capillarity. Oxford: Clarendon.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description