Modeling the nonlocal behavior of granular flows down inclines

Ken Kamrin and David L. Henann

Received Xth XXXXXXXXXX 20XX, Accepted Xth XXXXXXXXX 20XX

First published on the web Xth XXXXXXXXXX 20XX

DOI: 10.1039/b000000x

Flows of granular media down a rough inclined plane demonstrate a number of nonlocal phenomena. We apply the recently proposed nonlocal granular fluidity model to this geometry and find that the model captures many of these effects. Utilizing the model’s dynamical form, we obtain a formula for the critical stopping height of a layer of grains on an inclined surface. Using an existing parameter calibration for glass beads, the theoretical result compares quantitatively to existing experimental data for glass beads. This provides a stringent test of the model, whose previous validations focused on driven steady-flow problems. For layers thicker than the stopping height, the theoretical flow profiles display a thickness-dependent shape whose features are in agreement with previous discrete particle simulations. We also address the issue of the Froude number of the flows, which has been shown experimentally to collapse as a function of the ratio of layer thickness to stopping height. While the collapse is not obvious, two explanations emerge leading to a revisiting of the history of inertial rheology, which the nonlocal model references for its homogeneous flow response.

## 1 Introduction

^{†}

^{†}footnotetext: Department of Mechanical Engineering, MIT, Cambridge, MA, USA. E-mail: kkamrin@mit.edu

^{†}

^{†}footnotetext: School of Engineering, Brown University, Providence, RI, USA. E-mail: david_henann@brown.edu

Nonlocal effects in granular media manifest in myriad different ways. At the origin of the nonlocality is the finite size of the grains themselves, inducing cooperative behaviors that defy local rheological description. Examples include grain-size-dependent shear features in the steady flow profiles of granular media. A local rheology can be extracted from uniform simple shearing data of a granular media ^{1}; however, nonuniform steady flows of the same material can be seen to violate such a relation ^{2, 3, 4}, as the grain-size sets up an internal length-scale that effectively penalizes variations in flow-rate over space. A more recently observed nonlocal manifestation is the mechanically-induced creep effect, also known as “secondary rheology” ^{5, 6, 7}. A highlight of this phenomenon is that flow anywhere in a granular media removes the yield condition everywhere, and the rheology of probes in the body varies with distance from the source of motion. Some of the most compelling demonstrations of nonlocality in granular media can be observed in the behavior of grains on an inclined surface ^{8, 9, 10, 11}. Contrary to local rheological models, which predict a thickness-independent repose angle, experiments and discrete simulations verify that granular layers have a critical “stopping height” proportional to the grain size — layers thinner than this value come to a stop, whereas thicker layers admit steady flow down the incline.

Recently, a nonlocal rheological model based on the concept of “granular fluidity” has shown itself able to reconcile the issues of grain-size dependent shear features as well as the mechanically-induced creep effect in granular media ^{12, 13, 14, 15, 16}. While these demonstrations pertain primarily to the way in which grain size influences spatial fields in materials that are driven to flow, nonlocality as studied in inclined plane flows is of a relatively different nature, particularly the notion of a stopping height, i.e., whether a material flows at all. Despite this distinction, in this paper we shall show that the nonlocal fluidity model captures the phenomenology observed in the inclined plane geometry. Doing so provides a new test of the framework, in a family of problems distinct from those previously studied.

Herein, we compare theoretical results to known results (i.e., experimental or DEM) on the issues of (i) the presence of a stopping height and how it varies with inclination, (ii) the variation in shape of the flow profile as inclination increases above the stopping height, and (iii) the dependence of mean flow speed on the ratio of the thickness to the stopping height. These three phenomena have been well-studied (see aggregated data in MiDi ^{11}) and constitute the major manifestations of nonlocality in inclined plane flows.

## 2 From local to nonlocal granular rheology

We begin by describing the inertial granular rheology and the nonlocal granular fluidity model (per the review in Kamrin and Koval ^{14}). The inertial rheology for steady flow of dense granular media relates the local stress state and the strain-rate ^{1, 17, 18} and is expressed via the dimensionless relationship:

(1) |

In the above, is the ratio of shear stress and normal pressure , and is the inertial number, where is the shear rate, is the mean grain diameter, and is the grain density. The function for is empirically fit and is typically characterized by a yield coefficient such that . Under quasi-static circumstances, for which does not increase substantially above , a linear fit is often used,

(2) |

where is a dimensionless material parameter. As increases further above , curvature of the function is observed with an asymptote at an upper-limiting value of , denoted . A common fit is^{17}

(3) |

where and .

The inertial rheology works well in describing uniform flows (e.g., planar shearing) over a wide range of flow rates ^{1}. However, in non-uniform flow geometries in zones of slowly-flowing, quasi-static material the one-to-one inertial relation between and is violated ^{3}. The behavior displayed in these zones is definitively nonlocal; the bulk stress/strain-rate behavior at steady-state changes as the macroscopic geometry varies, even when the local kinematics are identical.

The “nonlocal granular fluidity” (NGF) model may provide the solution to the above issues. It has demonstrated the ability to quantitatively predict granular flows in many disparate geometries, with predictions verified in 2D and 3D as compared to both discrete-particle simulations ^{12} and experiments ^{13}. It is the first continuum model to quantitatively describe all experimentally-obtained flow data in the complex split-bottom family of geometries ^{13}. It has also been shown to correctly capture other nonlocal phenomena such as nonlocally induced material weakening, whereby the motion of a boundary removes the yield strength of the material everywhere, permitting far-away loaded objects to creep through the grains when otherwise they would remain static ^{16}.

Our exisiting work has implemented the NGF equations in a reduced, steady-state-only form

(4) |

where the linearized version of , Eq. (2), was used since flows of interest were all close to quasi-static. The field is a state parameter called the granular fluidity, and is the plastic cooperativity length, which is proportional to . Note that in planar shear, flow gradients vanish and the above reduces appropriately to the local rheology, but in the presence of gradients, the Laplacian term “spreads” fluidity based on . In our previous work ^{12, 13}, we have verified that is, in fact, not a constant length but satisfies

(5) |

roughly in-line with prediction of the kinetic-elasto-plastic (KEP) theory on which other fluidity models are based ^{19}. The parameter , the dimensionless nonlocal amplitude, is the only new parameter in the model, which quantifies the spatial extent of cooperativity in the flow.

### 2.1 Dynamical system for fluidity

In Henann and Kamrin ^{15}, we describe how the steady-state-only NGF model – i.e., Eqs. (1), (2), (4), (5) – is obtained in its entirety from the steady flow limit of a thermomechanically derivable dynamical system for . One treats the fluidity and its gradient as kinematic state variables with separate free-energy contributions. By selecting the corresponding free-energies in a fashion that preserves the linear inertial rheology in uniform flow conditions, Eq. (2), we obtain the following system:

(6) |

where is a constant time-scale. The steady-state-only model arises as the approximate solution of the stable, steady behavior of in the above dynamical partial differential equation (PDE). The approximation is valid as long as for some small function , where emerges as the stable solution when the Laplacian term above vanishes. To switch the local response to the more robust form, Eq. (3), the same analysis can be reapplied giving

(7) |

When the above is reduced to a steady-state only model, the nonlocal system obtained matches the aforementioned one in the appropriate limit of near , as it should, with maintaining the same inverse square-root divergence behavior^{*}^{*}*The precise form of in the steady-state system corresponding to Eq. (7) is
(8)
Note that the limit of Eq. (8) as goes to infinity corresponds to Eq. (5)..

The dynamics of presented above have yet to account for bistability of granular flows nor do they account for the totality of transient effects that occur in a sheared granular media. To the former point, recent experimental evidence suggests non-monotonicity of the local response (contrary to Eqs 2 and 3) may exist giving rise to bistable flow hysteresis in certain circumstances ^{20}. To the latter point, Reynolds dilation during shear initiation induces transient variations in material strength, which are commonly described using critical-state models ^{21} but which are not yet included in our fluidity dynamics. Instead, our model above intends to provide is an accurate long-term dynamical behavior of the material when passing through a sequence of developed flowing states. We note that the process of obtaining a steady-state-only relation from a fully dynamical form has a history within fluidity-based modeling for other amorphous media ^{19}. The dynamical form of the nonlocal fluidity model shows mathematical similarities to order-parameter-based rheological approaches, which also feature a diffusing state variable ^{22}.

## 3 Strengthening due to thinness

In the inclined plane geometry, a layer of thickness is inclined to an angle . Assuming a packing fraction of and acceleration of gravity , the stress distribution in the layer obeys

(9) |

See Fig. 1(a) for reference. Accordingly, if a local relation were applied, either Eq. (2) or (3), one would predict a universal angle of repose , i.e., any layer tilted above would be predicted to flow.

A signature of the cooperativity of granular motion is the fact that this universal repose angle is contradicted in experiments. As shown initially by Pouliquen ^{8} and verified by others ^{9, 10, 11}, the angle at which an initially flowing layer of grains comes to a stop, , depends sensitively on the thickness of the layer when is small. Inverting, one can extract a function for every granular media and substrate, which represents the critical thickness at which a flowing layer at a certain angle would arrest.

Unlike our past studies, which have focused on nonlocal effects in media that is necessarily flowing, the problem of size-dependent strengthening of thin layers often presents conditions that do not satisfy the approximation necessary to reduce the mathematics to steady-state-only form; i.e., material stops () although the inclination can be steep enough for to be notably greater than zero. To study inclined flows, we must revert to the primitive, dynamical form of the nonlocal fluidity model, Eq. (7), to ensure accuracy of solutions, which should be valid up to and including steep inclination angles.

We compare the results of our model to the experiments of Pouliquen ^{8}, where glass beads were used as the media. The local continuum parameters for glass beads have been calibrated in Jop et al. ^{17} and are , , , and kg/m. The nonlocal amplitude for glass beads was calibrated in Henann and Kamrin ^{13} to be . We can therefore apply our model in this geometry without additional parameter fitting, as long as we can identify accurate boundary conditions for the fluidity and basal slip velocity. We have tried two different fluidity boundary conditions in our past studies, chosen primarily based on simplicity, with the understanding that the choice is less important far from system boundaries due to the source-diffusion nature of the fluidity system. However, due to the thinness of the layers we wish to consider here, the accuracy of the boundary condition has a much larger influence on the flow behavior. To remove this issue, we opt instead to extract the fluidity boundary conditions directly from existing flow data. In discrete particle simulations of Silbert et al. ^{9}, who also used spherical particles, it is apparent that adjacent to a fully rough boundary, the shear-rate (and velocity fluctuations) approach an approximately vanishing state, and the mean velocity vanishes, indicating no slip. At or near the free surface, there is usually (though not always) a zone where the shear-rate appears to level off and a vanishing shear-rate gradient is observed. Likewise, we will assume and along the rigid surface at and take at the free surface, . Since these boundary conditions are homogeneous, they always permit both flowing and static solutions.

The phase diagram of flowing and non-flowing states according to the nonlocal model can be obtained by determining the conditions necessary for the global solution to be stable, as this decides if a system perturbed to flow returns to a static state. In Eq. (7), a perturbation about renders the term negligible. The remaining PDE is linear and, substituting , it is solved by

(10) |

where is a constant prefactor. The solution decays to as long as the term in the exponent is negative. Hence, a perturbed layer solidifies if

(11) |

A similar mathematical technique was utilized in ^{22}. The above shows that for tall layers, , the material is predicted to have a common repose angle of . However, for thinner layers, an inclination higher than this value is needed to stop flow. This applies up to an upper limit of at which all layers, independent of height, are predicted to flow. The existence of two critical angles having these properties has been verified as a common trait of data in inclined plane flows in multiple experiments and particle simulations involving different materials and substrates ^{11}. When it comes to specific functional forms for data, several fit curves have been proposed, a topic we shall return to later, but all invoke critical and values as described here. It bears noting that the derived relation, Eq. (11), is linearly proportional to the cooperativity length, Eq. (8), i.e., . This is in line with experimental observations of Pouliquen ^{23}, who reported that velocity correlation lengths in inclined plane flows vary with (and hence ) in a similar manner as . We emphasize that this is a derivable consequence – not an underlying assumption – of the dynamical system, Eq. (7). Further, one will obtain the same result when working with the simplified dynamical system, Eq. (6), albeit with slightly modified functional forms.

For a direct comparison, in Fig. 1(b), the above predicted form for , using the parameters for glass beads, is compared against Pouliquen’s experimentally determined values using glass beads with a fully rough base. It is worth pointing out that the model’s result was obtained using the same continuum parameters that were used to successfully predict steady flow fields of glass beads in split-bottom cells and other geometries ^{13}. Yet here, the question is of a different nature, one of predicting input conditions for flow stoppage rather than velocity profiles in a flowing body.

We note that a size-dependent starting height – the height at which a static layer initiates flow – is also observed experimentally. Its curve depends on the preparation of the static layer and differs somewhat from , though it shares the same qualitative features (e.g., and values where the curve asymptotes and vanishes respectively). Our analysis above is tailored to since it reflects the limiting behavior of a flowing solutions as (or ) is decreased, however we note that without a bistable term in our dynamics, a distinct curve does not arise theoretically.

## 4 Velocity profiles

For angles between and , and exceeding , the material has a steady flow state down the incline. As reported collectively in GDR MiDi ^{11} and explored directly in Silbert et al. ^{9}, the shape of the flow profile varies from concave up, to roughly linear, to concave down as is increased. The concave down shape observed for large is well-fit by the so-called ”Bagnold profile” obtainable by integrating the strain-rate predicted by the inertial rheology through the thickness of the layer. The result is of the form

(12) |

Figure 2(a) shows the discrete particle simulation data of Silbert ^{9} for inclination angle and various values of . We have computed numerical solutions to our model at the same inclination angle and several layer thicknesses. We note that the stopping curves of the simulated material and glass beads are different; at , the simulated particles have an value more than twice that of the experiment. To make an appropriate albeit qualitative comparison, we select our values so as to (roughly) match the relative heights used in the discrete simulations. The results are in Fig. 2(b). The model predictions are calculated by evolving Eq. (7) using finite-differences in MATLAB, using and the Bagnold profile as the initial guess^{†}^{†}†In the numerics, we allow a small pressure on the top surface corresponding to the weight of a layer of thickness, i.e., , to ensure that the coefficient of the term in Eq. (7) remains bounded..

The model shows that for near the profiles display a concave-up feature in agreement with the particle simulations. It could be said that the concave up appearance is due to the comparative unimportance of the term in Eq. (7) when flows are slow, which, if neglected, gives a cosinusoidal solution for and . As increases to be much greater than , the profiles approach the Bagnold profile, which agrees with the discrete simulations and is sensible since the relative importance of the particle-size-dependent term in Eq. (7) diminishes in a tall flow; if neglected, the remaining terms give exactly the inertial rheology at steady flow. The transition between these extremes is marked by velocity profiles displaying a linear character. It can also be observed that, in agreement with the data and in contrast to the Bagnold profile, the model does not always require a zero strain-rate at the free surface.

## 5 Discussion of the Froude number

We now discuss the mean flow speed of granular media down an incline. Pouliquen ^{8}, and others thereafter ^{9, 10, 11}, have observed that the Froude number of the flow, , appears to be a relatively well-defined, one-to-one function of the relative height, , for all angles and heights as long as is not close to . In this range, they find

(13) |

where, for glass beads, . When results of the nonlocal theory are plotted in this fashion, Fig. 3(a), we do not observe a one-to-one correspondence between the Froude number and the relative height. We have two comments on this point, that may explain the discrepancy.

### 5.1 A more precise curve

Collapse of the Froude number is sensitively dependent on the fit one uses for . In the form obtained from our theory, we notice some deviations from the experiment which magnify, unsurprisingly, as one approaches the asymptote at (see Fig. 1). A different fit function^{10}, empirically matched to the data in Pouliquen ^{8}, takes the form

(14) |

where (see the caption to Fig. 8 of Forterre and Pouliquen ^{10}). If we use the above relation for , a rather strong collapse of the Froude number versus relative height emerges (see Fig. 3(b)). Thus, it may be that the lack of a collapse of the Froude number is attributable to deviation of the predicted curve.

### 5.2 Consistent local rheology

An alternative explanation could be found if we revisit the history of the local inertial rheology. The rheology in Eq. (3) was arrived at, coincidentally, by fitting data for glass beads flowing down inclines (see Appendix A of Jop et al. ^{17}). First, a form for was fitted empirically from experiments (i.e., Eq. (14)). Then was calculated in terms of an as-yet-unknown function (assuming fully local rheology and, consequently, the Bagnold velocity profile, Eq. (12)). This result was subsequently substituted, along with the empirical function, into Eq. (13) to solve for . The result is a local rheology having the form

(15) |

which generates Eq. (3), upon substituting Eq. (14) for . As was done in Jop et al. ^{17}, the term containing the square-root is replaced with a constant ( for glass beads) because the function does not vary much in the range of use. We see that, historically, the fit function chosen for gives rise to the functional form one obtains for the inertial rheology.

Let us repeat this logic in a thought experiment. Suppose, through Eq. (15), we replace the local rheology in our theory with the one corresponding to our fit for the function, Eq. (11). The result is

(16) |

With this, the nonlocal fluidity dynamics are expressible conveniently as

(17) |

Only the term has changed so the previous stability argument identically produces our function, and Eq. (16) is now obtained in the homogeneous flow limit as the new local response^{‡}^{‡}‡In fact, the steady-state-only system corresponding to Eq. (17) is given through Eq. (4) with local rheology Eq. (16) and cooperativity length Eq. (8). We note that the coefficient of in Eq. (17) is for the moment only defined for as needed for the Froude number analysis, but any monotonic extension of for may be assumed without affecting the analysis nor the corresponding steady-state-only relation..

We apply the above system to inclined plane flows. Letting and , writing , recalling the definition of , and allowing , Eq. (17) produces steady solutions given by

(18) |

The only varying parameter in the above is the ratio . Hence, all solutions have the form, . Likewise, the velocity field obeys

(19) |

Finally, we arrive at a relation for the Froude number:

(20) |

We conclude that when the same curve we derive theoretically is used to fit the local inertial relation, the emergent nonlocal theory requires a collapse between the Froude number and for all . This is an interesting result because while the local relation on its own would require such a collapse by definition, this is the full nonlocal solution. When is large enough, the second-derivative term in Eq. (18) is negligible (excepting a rapid variation near the lower boundary), and the solution for approaches that of a purely local material behavior. Likewise the function obeys for , matching Eq. (13). Numerical integration of Eq. (17) verifies these analytical properties (see Fig. 3(c)). We observe that , which is sensible as this corresponds to the stopping height, and our model does not include bistable terms. Although it is unclear in discrete simulations and experiments what the precise behavior of the Froude number is near , for round grains on a fully rough surface many data trends show the Froude number exceeds zero at this height ^{8} and other data suggests it equals zero ^{24}, as our model purports.

A short commentary is in order. The above result comes at the cost of adjusting the relation from the commonly used form of Eq. (7), to a form compatible with the theory’s own predicted function. The new relation, Eq. (16), is qualitatively similar to the previous, with increasing monotonically from to as increases. One difference is that the new relation has at whereas the former has a positive slope. There has been a debate recently over the behavior of near . Some experimental data suggests, in fact, a negative initial slope for ^{20, 25}, while other fits suggest a steeper-than-linear power-law behavior near ^{26}. It bears noting that our previous work on nonlocal fluidity has focused solely on quasi-static flows, in which almost everywhere and the only important aspect of the local rheology is the value of .

## 6 Conclusion

We have applied the theory of nonlocal granular fluidity to the canonical problem of size-dependence in granular inclined plane flows. The theory, as calibrated to glass beads based on prior data, predicts an curve that matches experimental data rather well. Further, the theory predicts flow profiles that vary in shape in a fashion consistent with existing discrete simulation data, marked by an upward curvature for layer thicknesses close to and transitioning to the Bagnold profile for large . While the Froude number does not exhibit a direct collapse against we have identified two possible explanations for this and corresponding ways in which the collapse can indeed be obtained.

Perhaps the most compelling result herein is the prediction of the function, which suggests the nonlocal fluidity concept could be used to model other size-sensitive flow stoppage phenomena, such as silo jamming. The famous Beverloo correlation ^{27}, an empirical functional form that gives silo flow rate in terms of aperture opening size, indicates a critical opening size at which flow stops. This effect may be roughly analogous to the stopping height observed for inclined flows, and it would be useful to apply the nonlocal fluidity model in this geometry to determine if the Beverloo correlation can be obtained from the nonlocal model.

## Acknowledgements

KK acknowledges funds from NSF-CBET-1253228 and the MIT Department of Mechanical Engineering, and DLH acknowledges funds from the Brown University School of Engineering.

## References

- da Cruz et al. 2005 F. da Cruz, S. Emam, M. Prochnow, J.-N. Roux and F. Chevoir, Phys. Rev. E., 2005, 72, 021309.
- Jop 2008 P. Jop, Phys. Rev. E, 2008, 77, 032301.
- Koval et al. 2009 G. Koval, J.-N. Roux, A. Corfdir and F. Chevoir, Phys. Rev. E, 2009, 79, 021306.
- Kamrin 2010 K. Kamrin, Int. J. Plasticity, 2010, 26, 167–188.
- Nichol et al. 2010 K. Nichol, A. Zanin, R. Bastien, E. Wandersman and M. van Hecke, Phys. Rev. Lett., 2010, 104, 078302.
- Reddy et al. 2011 K. Reddy, Y. Forterre and O. Pouliquen, Phys. Rev. Lett., 2011, 106, 108301.
- Wandersman and van Hecke 2014 E. Wandersman and M. van Hecke, Europhys. Lett., 2014, 105, 24002.
- Pouliquen 1999 O. Pouliquen, Phys. Fluids, 1999, 11, 542–548.
- Silbert et al. 2003 L. E. Silbert, J. W. Landry and G. S. Grest, Phys. Fluids, 2003, 15, 1–10.
- Forterre and Pouliquen 2003 Y. Forterre and O. Pouliquen, J. Fluid Mech., 2003, 486, 21–50.
- MiDi 2004 G. D. R. MiDi, Eur. Phys. J. E, 2004, 14, 341–365.
- Kamrin and Koval 2012 K. Kamrin and G. Koval, Phys. Rev. Lett., 2012, 108, 178301.
- Henann and Kamrin 2013 D. L. Henann and K. Kamrin, P. Natl. Acad. Sci. USA, 2013, 110, 6730–6735.
- Kamrin and Koval 2014 K. Kamrin and G. Koval, Comput. Part. Mech., 2014, 1, 169–176.
- Henann and Kamrin 2014 D. L. Henann and K. Kamrin, Int. J. Plasticity, 2014, 60, 145–162.
- Henann and Kamrin 2014 D. L. Henann and K. Kamrin, arXiv:1408.3884, 2014.
- Jop et al. 2005 P. Jop, Y. Forterre and O. Pouliquen, J. Fluid Mech., 2005, 541, 21–50.
- Jop et al. 2006 P. Jop, Y. Forterre and O. Pouliquen, Nature, 2006, 441, 727–730.
- Bocquet et al. 2009 L. Bocquet, A. Colin and A. Ajdari, Phys. Rev. Lett., 2009, 103, 036001.
- Dijksman et al. 2011 J. A. Dijksman, G. H. Wortel, L. T. van Dellen, O. Dauchot and M. van Hecke, Phys. Rev. Lett., 2011, 107, 108303.
- Schoefield and Wroth 1968 A. Schoefield and P. Wroth, Critical State Soil Mechanics, McGraw-Hill, 1968.
- Aranson and Tsimring 2002 I. S. Aranson and L. S. Tsimring, Phys. Rev. E, 2002, 65, 061303.
- Pouliquen 2004 O. Pouliquen, Phys. Rev. Lett., 2004, 93, 248001.
- Weinhart et al. 2012 T. Weinhart, A. R. Thornton, S. Luding and O. Bokhove, Granul. matter, 2012, 14, 531–552.
- Kuwano et al. 2013 O. Kuwano, R. Ando and T. Hatano, Geophys. Res. Lett., 2013, 40, 1295–1299.
- Peyneau and Roux 2008 P.-E. Peyneau and J.-N. Roux, Phys. Rev. E, 2008, 78, 011307.
- Beverloo et al. 1961 W. Beverloo, H. Leniger and J. Van de Velde, Chem. Eng. Sci., 1961, 15, 260–269.