Ken
Kamrin
*^{a} and
David L.
Henann
^{b}
^{a}Department of Mechanical Engineering, MIT, Cambridge, MA, USA. E-mail: kkamrin@mit.edu
^{b}School of Engineering, Brown University, Providence, RI, USA. E-mail: david_henann@brown.edu

Received
19th August 2014
, Accepted 29th September 2014

First published on 6th October 2014

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.

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.^{15–19} 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 also captures the phenomenology observed in the inclined plane geometry.

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.

(1) |

μ_{loc} (I) = μ_{s} + bI, | (2) |

μ_{loc} (I) = μ_{s} + Δμ/(I_{0}/I + 1) | (3) |

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 I 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^{15} and experiments.^{16} It is the first continuum model to quantitatively describe all experimentally-obtained flow data in the complex split-bottom family of geometries.^{16} 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.^{19}

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

(4) |

(5) |

(6) |

(7) |

(8) |

The dynamics of g 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 eqn (2) and (3)) may exist giving rise to bistable flow hysteresis in certain circumstances.^{24} To the latter point, Reynolds dilation during shear initiation induces transient variations in material strength, which are commonly described using critical-state models^{25} 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.^{22} The dynamical form of the nonlocal fluidity model shows mathematical similarities to order-parameter-based rheological approaches, which also feature a diffusing state variable.^{12}

μ = tanθ, P = ϕρ_{s}Gzcosθ. | (9) |

Fig. 1 (a) Setup of the inclined plane flow geometry. (b) Theoretically predicted H_{stop} locus (−), eqn (11), as calibrated for glass beads on a fully rough base, compared to experimentally determined values (○) from Pouliquen.^{8} |

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–11} the angle at which an initially flowing layer of grains comes to a stop, θ_{stop}, depends sensitively on the thickness of the layer when H is small. Inverting, one can extract a function H_{stop} (θ) 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 g = g_{loc} + δ approximation necessary to reduce the mathematics to steady-state-only form; i.e., material stops (g = 0) although the inclination can be steep enough for g_{loc} to be notably greater than zero. To study inclined flows, we must revert to the primitive, dynamical form of the nonlocal fluidity model, eqn (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.^{20} and are μ_{s} = 0.382 = tan(20.9°), μ_{2} = 0.643 = tan(32.76°), b = Δμ/I_{0} = 0.938, and ρ_{s} = 2450 kg m^{−3}. The nonlocal amplitude for glass beads was calibrated in Henann and Kamrin^{16} to be A = 0.48. 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 g = 0 and ν = 0 along the rigid surface at z = H and take ∂g/∂z = 0 at the free surface, z = 0. 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 g = 0 solution to be stable, as this decides if a system perturbed to flow returns to a static state. In eqn (7), a perturbation about g = 0 renders the g^{2} term negligible. The remaining PDE is linear and, substituting μ = tanθ, it is solved by

(10) |

(11) |

For a direct comparison, in Fig. 1(b), the above predicted form for H_{stop}(θ), 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.^{16} 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 H_{start}(θ) – 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 H_{stop}, though it shares the same qualitative features (e.g., θ_{1} and θ_{2} values where the curve asymptotes and vanishes respectively). Our analysis above is tailored to H_{stop} since it reflects the limiting behavior of a flowing solutions as H (or θ) is decreased, however we note that without a bistable term in our dynamics, a distinct H_{start} curve does not arise theoretically.

ν_{Bag}(z) ∝ H^{3/2} − z^{3/2}. | (12) |

Fig. 2 (a) Discrete particle simulation data of Silbert et al.^{9} for θ = 24° and many layer thicknesses (increasing from top to bottom). Dashed line is the Bagnold profile. (b) Theoretical profiles for the same relative heights, H/H_{stop} ≅ 1.1, 1.9, 2.7, 4.1, 6.4, as well as H/H_{stop} = 20. |

The model shows that for H near H_{stop} 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 g^{2} term in eqn (7) when flows are slow, which, if neglected, gives a cosinusoidal solution for g and ν ∼ 1 − sin(zπ/2H). As H increases to be much greater than H_{stop}, 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 eqn (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. Because the value of A determines the degree of nonlocality, and hence the degree of deviation from the Bagnold profile, it is clear that as A increases, the range of heights over which the profile remains concave up will increase. However, since increasing A also increases H_{stop} proportionally, it is the value of H/H_{stop} that has the most direct influence over the shape of the flow profile.

(13) |

Fig. 3 Predicted Froude number, , at inclinations θ = 22.5(+), 24(□), 25.5(○), 27(Δ), 28.5(∇), and 30(*) degrees while varying H from 12–80d, plotted against the relative height, utilizing H_{stop}(θ) as shown in (a) eqn (11) and (b) eqn (14). (c) Theoretical result when the local rheology is chosen consistent with the theoretical H_{stop}, per eqn (15). In (a)–(c), dashed line is the experimental collapse line, eqn (13), with β = 0.136. |

(14) |

(15) |

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

(16) |

(17) |

We apply the above system to inclined plane flows. Letting = z/H and , writing P = ϕρ_{s}Gzcosθ, recalling the definition of B, and allowing θ_{1} ≤ θ ≤ θ_{2}, eqn (17) produces steady solutions given by

(18) |

(19) |

(20) |

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

Perhaps the most compelling result herein is the prediction of the H_{stop} 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,^{30} 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.

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

## Footnotes |

† In the numerics, we allow a small pressure on the top surface corresponding to the weight of a layer of (1/2)d thickness, i.e., P_{top} = (1/2)ϕρ_{s}Gdcosθ, to ensure that the coefficient of the g^{2} term in eqn (7) remains bounded. |

‡ In fact, the steady-state-only system corresponding to eqn (17) is given through eqn (4) with local rheology eqn (16) and cooperativity length eqn (8). We note that the coefficient of g^{2} in eqn (17) is for the moment only defined for μ > μ_{s} as needed for the Froude number analysis, but any monotonic extension of 1/H_{stop} for μ < μ_{s} may be assumed without affecting the analysis nor the corresponding steady-state-only relation. |

This journal is © The Royal Society of Chemistry 2015 |