Non-locality in Granular Flow: Phenomenology and Modeling Approaches

This paper reviews the emergence of non-local flow phenomena in granular materials and discusses a range of models that have been proposed to integrate an intrinsic length-scale into granular rheology. The frameworks discussed include micro-polar modeling, kinetic theory, three particular order-parameter-based models, and strongly non-local integral-based models. An extensive commentary is included discussing the current capabilities of these existing models as well as their implementational ease, physical motivation, and breadth of predictive ability.


INTRODUCTION
Let us define a local rheology as a constitutive model whose well-developed strain-rate response depends only on the stress and no higher-gradients of kinematic quantities or stress. Such models include relations between the stress tensor and the tensors of strain and/or strain-rate and state variables evolving locally within each material element. In granular flow modeling, local models can provide qualitative and sometimes quantitative predictions for flows and stresses [1][2][3][4][5][6][7][8][9][10][11][12][13]. However, it is now widely accepted that such models lack robustness in their ability to predict all flow phenomena and, in many cases, predictions of local models disagree with experiments by nontrivial amounts [14][15][16][17][18][19]. This is true even though such models predict uniform flows well, such as steady simple shearing. It is the fact that these models succeed in homogeneous cases but break down in the presence of spatial inhomogeneity that the origin of this difference can be attributed to a non-local effect. In such models, a microscopic length-scale emerges on dimensional grounds. In the case of granular media, it is evident that the mean grain size d provides this microscopic length unit, and the observed length-scale of non-local effects should be some multiple of d.

COMMONLY OBSERVED MANIFESTATIONS OF NON-LOCALITY IN GRANULAR FLOWS
There are a few particular manifestations of non-locality in granular flow that bear mentioning as canonical examples. All of these effects are evident in systems of the simplest grains-round, stiff, frictional particles, with a restitution coefficient.
First, there is the class of problems wherein the size of shear features is influenced by the grain size, not just the size scale of stress gradients.
• In geometries with wall-located shear bands at steady state, such as flow in an annular couette cell, the shear stress decays gradually toward zero as one moves away from the rapid zone near the inner wall. Therefore, any local concept of a yield stress would predict that the velocity field should vanish in some extended domain outside the main shear band. However, slow, creeping flow is seen throughout the entire geometry. The velocity field has an exponential-type decay moving away from the rapid flow zone, with decay length governed by the grain size [20][21][22][23].
• A similar phenomena is observed in steady flow geometries where the zone of flow is not concentrated near walls, such as silo flows [24][25][26][27], split-bottom couette flows [16,28], or sub-surface flows in heaps/rotating-drums [19,29]. In these cases, despite the fact that there are spatial zones where the stress varies below what a local model may predict to be the yield stress, non-zero flow is observed everywhere. In triaxial compression of a highly packed column of grains, an unsteady flow case, a shear-band of grain-size-dependent width emerges within the sample to mitigate the boundary displacement [30,31].
Another non-local effect is connected to the reverse problemrather than material flowing even though the stress seems too small, there are cases where material does not flow even though the stress appears to be above the yield limit. In these cases, non-locality makes thinner domains act stronger. There are two common examples of this effect.
• In silos, hoppers, and hourglasses, it is common experience that the nozzle arch-jams when the opening size is on the order of a few d [26,[32][33][34]. The critical opening size plays a key role in models that predict silo outflow rates such as the seminal Beverloo correlation [33]. However, local continuum models cannot capture this effect. Nozzle jamming requires an intrinsic length-scale in order to give a critical opening size. • In flows down rough inclined surfaces, classical understanding gives an angle of repose characteristic of the material; material flows until the tilt is reduced below this angle. However, as layers get thinner it is observed that the repose angle is not constant; it is thickness dependent with thinner layers acting stronger [14,35,36]. This effect is often described inversely through the function H stop (θ )/d, which maps tilt angles to the critical thickness H stop at which such a flow would arrest.
The third non-local effect we discuss is the "secondary rheology, " where flow somewhere causes the rest of the domain to lose its yield stress. This more recently studied non-local effect can be seen in tests where a constant-force probe is placed in approximately rigid subdomains of a flowing collection of grains. For example: • A heavy ball that is placed on the surface of a cup of grains will sit on the free surface after some mild sinkage [37,38]. However, if any part of the cup's boundary is moved-such as by replacing the bottom of the cup with a rotating turn-tablethe ball will sink even though the material near it is otherwise quiescent. This is true even if the free surface of the grains is far from the moving bottom boundary of the container. The notion that perturbations from a boundary cause rheological changes away has been observed as well in suspensions of fine particles [39]. • In annular shear of a granular media, the material far from the rotating wall appears to be still. However, if a probe is placed in this part of the domain and a constant force applied to it, it is observed to move through the material even as the force on the probe shrinks to zero [18]. When the rotation of the inner wall stops, the same probe only moves if the force on it exceeds a critical value; i.e., the yield stress is "restored" once the far-away boundary stops moving.

NON-LOCAL MODELING APPROACHES
Below is a non-exhaustive summary of non-local granular flow models with a brief description of the model assumptions and how the length-scale is introduced. The discussion is limited to models that address (at least) steady-state flow behavior.

Cosserat Continuum
One of the most straightforward and oldest ways to introduce a size-scale is to presume the continuum is a "micropolar" or "Cosserat" continuum [40][41][42][43]. In such models it is assumed that a body posseses not just a stress field σ and velocity field v but also a couple-stress field M and intrinsic spin field ω (see Figure 1). Like a standard continuum, on any internal surface with outward normal n the force per area (traction) on the surface is given by σ n. However, Cosserat material elements are presumed to also have a non-vanishing torque-per-area given by Mn. As a result, while standard continua satisfy torque balance as long as one imposes σ = σ T , a Cosserat continuum requires an additional, non-trivial differential equation for angular momentum balance invoking the divergence of the tensor M, possible asymmetry in σ , and the rate-of-change of ω. It is understood that ω represents the spin of the individual microscopic components of the media, which can differ from the vorticity, ∇ × v, which measures the spin of a collection of those components as seen from following their center-of-mass velocities. An intrinsic length-scale arises in such models when imposing a constitutive relation, in our case a yield criterion-since the couple-stress and stress differ by a unit of length, a yield criterion mixing these two fields requires a length-scale [44]. The motivation to connect Cosserat continuum assumptions to granular media stems from the clear possibility that the microscopic components are the grains themselves, and that their individual spins give rise to ω. Thus, the needed length-scale is some multiple of the grain size d. One might visualize a "rolling resistance" mechanism, where couple-stress per d governs the onset of rolling, which blends with sliding to produce a joint yield condition.

Kinetic Theory
The kinetic theory of gases has a well-known extension to the case of granular media, which ultimately brings a dependence on granular temperature (given microscopically by particle velocity fluctuations) into the granular constitutive relations [45][46][47][48][49]. While it may at first seem odd to describe a thermal flow model as non-local, note that it passes our definition from the introduction-by letting the rheology depend on the temperature field, the length-scale controlling the spread in temperature, a multiple of d, can now influence the flow field.
Unlike an ideal gas, granular media lose energy every collision. Thus, energy balance can be expressed via where, up to prefactors involving packing fraction and restitution, the viscosity η and conductivity K are ∼ d √ T, and the dissipation rate Ŵ is ∼ 1 d T 3/2 . At steady state the temperature field finds a distribution whereby the thermal energy generated by the work term (the first on the right side) is diffused by the conduction term (second on the right side), and reduced by the dissipation term (last term). Therefore, thanks to the non-vanishing conductivity K, non-locality emerges through a self-heating process; work done deforming a material element somewhere produces heat that spreads to material nearby and changes its rheology (see Figure 2). As will be seen in the upcoming sections, this notion that fluctuations spread into neighboring material and alter the flow rule is important in a variety of non-local theories.

Order Parameter Models
Some models append a scalar order parameter to the system of equations, which indicates how "fluid-like" the material is. The idea to use an order parameter is connected to the notion that the transition from solid to flowing is a type of phase change. The order-parameter is presumed to obey a diffusion-like equation and its value affects the rheology, thereby imbuing a length-scale into the flow. While these features are similar to kinetic theory, different order-parameter models differ on what precise variable is diffusing and what its diffusion equation looks like.

Partial Fluidization Theory
One of the first granular models of this type was the partial fluidization theory of Aranson and Tsimring [50,54] and Volfson et al. [55]. Their order parameter, ρ, varies from 1 (solid) to 0 (fluid). In its common usage, the material's viscosity is presumed to obey η f /(1 − ρ) so that the viscosity diverges as the material becomes solid and reverts to a (constant) fluidlike viscosity η f when ρ goes to zero. The order parameter is assumed to participate in a free-energy functional F ∼ FIGURE 2 | The basic "self-heating" picture common to many recent non-local models [50][51][52][53]. Zones whose stress ratio µ are beneath the local friction coefficient µ s can be made to flow in inhomogeneous stress environments by diffusive passage of a state field (e.g., temperature, order-parameter, force-fluctuations) sourced in higher stress media nearby, which spreads over some grain-size dependent correlation length. The influx elevates the value of the state field in the low-stress media, which reduces its flow resistance and allows flow to occur. The newly "heated" material then "heats" neighboring material triggering a cascade that causes all material to flow. l 2 |∇ρ| 2 + f (ρ, σ ) dV and its evolution is obtained assuming purely dissipative dynamics,ρ = δF/δρ. The evolution equation that emerges looks like a basic diffusion equation but with a source term given by the choice of f . Their suggested choice of f gives the final result: where δ ≡ (µ − µ d )/(µ s − µ d ). The constants µ s and µ d represent static and dynamic friction coefficients, respectively, and µ is the ratio of local shear stress to pressure. The cubic nature of the source term causes stability switches that give rise to frictional hysteresis. By setting l ∼ d, the equation forces an inherent grain-size dependent length-scale into the model's flow predictions, and causes flow stability to be size-sensitive, giving rise to size-dependent flow start/stoppage.

Non-local Granular Fluidity
The Non-local Granular Fluidity (NGF) model is a more recent approach, developed by the author and coworkers (primarily Henann and Koval) [51,56,57]. Inspired by a similar model for emulsions [58,59], the NGF model presumes a scalar-valued order-parameter-like field called the granular fluidity, g, whose value ranges from zero to infinity. Whereas, the fluidity in the emulsion flow model operates as an inverse-viscosity-like field in their flow rule, the granular fluidity g enters the flow rule througḣ γ = gµ and hence may be interpreted as a pressure-weighted inverse viscosity field. The fluidity field, in steady flow, is then posited to obey The diverging of the correlation length, ξ , where µ = µ s has been shown directly in simulations [51,56] and experiments [23], with A an order-one material constant. The function g loc is what the fluidity would be if the rheology were local. The local rheology, obtained in uniform simple shearing, is the µ(I) rheology [1], where the stress ratio µ = τ/P is given by the normalized shear rate I =γ d 2 ρ s /P. The fluidity at a point is then viewed as having a contribution from local effects (g loc ) and a non-local contribution due to "disturbances" from neighboring material (ξ 2 ∇ 2 g). Out of the many choices of "flow variable"/"load variable" one could choose for the role of fluidity, g =γ /µ was initially selected for practical reasons as it was the only one that produced quantitative flow predictions for a range of data from several different flow geometries [51,56]. This choice gained additional support when it was shown that g has a clear kinematic meaning as a function of particle velocity fluctuations and packing fraction; other possible fluidity definitions do not [52]. This suggests the fluidity field measures the degree of fluctuation-based activation in the material. Parameter studies indicate the NGF model is robust to material changes-varying particle properties such as surface friction affects the parameters, including A, but the model itself maintains its accuracy [60]. The model has been extended into a transient form that evolves g over time, which can be used to model the smaller-is-stronger effect as an instability process [57,61]. The more primitive transient form can be justified thermomechanically under the virtual power principle with g as an energetic order parameter [61].

I-gradient Model
Shortly after the NGF model, an alternative model for steadystate non-local flow was proposed in Bouzid et al. [62,63]. The proposed form is a more direct non-local expansion of the µ(I) rheology: The above does not propose an independent order parameter with its own equation, but rather treats I itself as a fluidity-type field within a gradient expansion of the flow rule. If anİ term were added to account for unsteady cases (like in [64]), the model would claim direct diffusion of the I field. For constant ν ∼ d 2 , the form above has the primary feature that effective friction µ is reduced when neighboring material is flowing faster, and increased when neighboring material flows slower. In doing so, the model can, by design, replicate the effect of creeping flow beneath the local yield criterion, µ s , as well as rate independence of the non-local response due to the scaling of ∇ 2 I by 1/I. Like NGF, the model replicates the effect of a diverging length-scale about µ s .

Integral Equation Approaches
While the previous models produce non-locality by letting a rheological variable obey an additional PDE, non-locality can also be achieved through an integral equation, which describes the flow response at a point in terms of influences coming from surrounding material finite distances away. In simplified 1D, these "strongly non-local" models generally take the forṁ where H can be the entire domain, or be limited to a horizon a finite distance away. It is well known that one can Taylor expand in x ′ under the integral and integrate through to obtain an approximate model in the form of a PDE. One model of the strongly non-local type is the self-activation model of Pouliquen and Forterre [53], in which the function F is selected to represent the product of the frequency of fluctuations emitted from position x ′ and the probability that a fluctuation from x ′ will result in a forward vs. backward shear event at x. Another strongly non-local model is the recent model of Nott [65], in which the integral is a weighted spatial convolution of what the local flow-rate formula would give in the surrounding media. A strongly non-local eddy-viscosity-type model has also been written [66], in which the shear-rate has a contribution from redistribution of vorticity occurring through the geometry. The corresponding PDE approximation of this model takes a form similar to an order-parameter model [66,67].

COMMENTARY
All models discussed could be evaluated based on three criteria: (i) ease of implementation, (ii) physical motivation, and (iii) breadth of predictive capability. Below, we summarize where each model stands, and conversely where models need further understanding.
On issue (i), clearly, PDE's are simpler to solve than integral equations. Kinetic theory, partial fluidization, and NGF propose one additional PDE involving one new state variable. These can be solved numerically using the finite-element method (see for example, the implementation in [68]) and sometimes analytically in unidirectional flow cases [44,54,57,62,69]. Based on its mathematical form, it is not clear what the difficulty of numerical integration of the I−gradient model would be in arbitrary cases. However, the fact that it requires no new state variables and is a straightforward extension of µ(I) is a plus. Another issue is that of (Hadamard) well-posedness, which determines if small wavelength flow perturbations grow unboundedly in the linear regime, an issue known to arise for basic local models [70,71]. Hadamard unstable models may lack existence of a solution, and/or have grid-size dependent or unstable numerical solutions. It is known that gradient corrections have the potential to resolve the ill-posedness [72]. It has been directly shown that NGF and certain strongly non-local models [65] are Hadamard wellposed, while the I−gradient model has been proven Hadamard ill-posed [73].
On issue (ii), kinetic theory and Cosserat continua have a strong microscopic basis, with their additional PDE's representing energy balance and moment balance, respectively. The physics of the boundary conditions are, hence, fairly clear in these models. While derived for dilute systems, the kinetic theory has gone through a number of recent enhancements to bring it closer to representing dense media [74,75]. With Cosserat theory, to this author's knowledge, stress asymmetry, couple stresses, and spin-differences are notably small in tests of simple particles excluding cases of non-negligible rolling resistance or directly wall-adjacent behavior [15,76,77]. The partial fluidization and NGF models have each identified (postfacto) kinematic observables that match the behavior of their hypothesized order parameter, i.e., static contact fraction for partial fluidization [55], and weighted velocity fluctuations for NGF [52]. The full physical meaning of these models' phenomenological PDE's is not understood, and selection of boundary conditions can be unclear. That said, NGF has been shown to follow a virtual power principle, which aids in boundary condition selection [61]. The I-gradient model does not require a new state variable, though it is also conceivable, in cases of high ∇ 2 I, that the model produces negative µ values, violating nonnegative dissipation. Strongly non-local models have the benefit of absolving the question of boundary conditions [65]. Moreover, the strongly non-local self-activation model is obtained entirely from statistical and mechanical propositions at the small scale.
On issue (iii), the extent to which the different models have been tested varies quite a bit. NGF, with its own dedicated numerical solver [68], has been validated in the most geometries and conditions [78]. The NGF and the I−gradient models have each shown quantitative agreement with non-local creep in multiple 2D experiments [23,51,62], capturing proper flow spreading and the emergence of rate-independence. NGF also produces quantitative predictions in 3D cases such as flows in split-bottom cells [56] and wall-bounded heap flows [79]. It also captures smaller-is-stronger phenomenology in multiple geometries including H stop [57,69], as well as all features of secondary rheology [80] including the effect's direction sensitivity [81]. Partial fluidization captures these different phenomena to a qualitative level [50], and additionally is equipped to model stickslip and start-stop differences due to the inclusion of hysteresis (although our understanding of the form for frictional hysteresis is evolving [82]). The strongly non-local self-activation model has matched an experimental H stop curve and produces reasonable non-local flow-spreading [53]. Cosserat and kinetic theories produce non-local flow spreading as well [44,83], and certain kinetic theories can also produce an H stop effect [84]. It is not clear if these predictions are quantitative across many geometries. Similarly, the strongly non-local plasticity model of Nott [65] possesses qualitative elements of non-local flow spreading.
In closing, to those interested in new model development, it should be emphasized that the universe of non-local models is not unconstrained. The reader is directed to Gurtin et al.
[85] and Maugin [42] for a review of the constraints of thermodynamic consistency and existing modeling formalisms.

AUTHOR CONTRIBUTIONS
The author confirms being the sole contributor of this work and has approved it for publication.