Modeling flow arrest using a non-local rheology ?

Using an incompressible Navier-Stokes solver with a modified non-Newtonian viscosity, we implement the μ(I)-rheology and study the transient arrest for thick granular layers. In order to study the stopping height phenomenology Hstop, we implement the steady-state-only approximate version of the non-local model from [6, 7] and successfully incorporate this feature in the solver. We compare the velocity profiles for granular layer thicknesses close to Hstop given by the model with discrete simulations, and discuss the ability of the nonlocal model to describe the dynamics of the flow close to stopping height. We find that for flows with thickness close to Hstop, the velocity profiles given by the model do not reflect the behavior of discrete simulation: discrete velocity profiles remain self-similar when approaching Hstop, while the non-local rheology predicts a change of concavity near the bottom.


Introduction
The flow of granular matter on a rigid or erodible bed seems the simplest configuration one can think of, yet its intrinsic complexity has long been a challenge to understanding and modeling.As a matter of fact, surface flows and chute flows still form a much studied system.Those flows have a lot of industrial and geophysical applications.The existence of a finite range of slope values in which the flow is stationary and the evolution of the rheological properties of the flow within this range is well captured by the μ(I)-rheology [3].However, other features require the introduction of non-local effects in the rheological model, which allow for local information on the shear rate to propagate in the system: a distant shear, for instance, or a solid boundary.
In this contribution we are interested in the arrest of granular chute flows, either during the transient evolution towards arrest, or during steady flow close to the critical thicknesses H stop and H start [11].In particular, we try to address the relevance of non-local modeling to describe the arrest mechanisms.

Modeling non-local effects 2.1 The μ(I)-rheology and the non-local fluidity μ(I)-rheology
We consider the case of a dense granular flow in the xz plan, tilted at an angle θ and oriented in the x direction, which develops a steady velocity profile u(z) as illustrated in Fig. 1.As the flow develops, by analogy with the Coulombic friction, the local normal stress p and the local tangential stress τ are found in practice and assumed in the model to be proportional: μ being the analogue of a coefficient of friction, that may vary from point to point in a granular flow.The strength of the μ(I)-rheology lays in the fact that it relates the effective coefficient of friction μ to a dimensionless number reflecting the local state of the granular packing, the inertial number I [3], defined as: where |γ| = |du/dz| is the shear rate, d is the mean diameter of the grains and ρ their density.The normal stress p is set to zero at the surface (z = H).The following nonlinear expression on I is adopted to account for the shape of the μ(I) dependence [5]: where the values of the coefficients μ s , μ 2 and I 0 are material-dependent.By analogy with Newtonian fluids, we define an effective granular viscosity η eff proportional to the shear rate, such that from which we deduce the expression for the non-Newtonian viscosity used in the model: Therefore, the granular viscosity depends at the same time on the material characteristics, on the inertial number, on the local normal stress and on the shear rate itself.The previous equations have been written for a pure shear flow, see [5] for a 3D generalized formulation.

Non-local model
The non-local model implemented in this work, from [6,7], consists in introducing a granular fluidity defined by: and solving for each time step the following additional set of equations with values obtained from the previous local formulation, in order to obtain a global (non-local) value for the granular viscosity (written here for a pure shear flow: z dependence only): The field g is a state parameter and ξ is the plastic cooperativity length, proportional to the grain diameter d and calibrated with a constant A: In planar shear, d 2 g/dz 2 = 0 and the above reduces appropriately to the local rheology (g = g loc ), but in the presence of gradients, the second-derivative term spreads fluidity based on ξ.
An important characteristic of the model is the imposition of zero-shear on the wall, which corresponds g(z = 0) = 0, condition that will be propagated by the model.See [6,7] for more details.We implemented the steady-state-only approximate version of the model in Basilisk 1 , and solved the momentum balance along the x and z axis, together with 3. Full details on the resolution will be found in [8].Velocity profiles for steady flows with various values of H/d are then compared.

Asymptotic analysis for a chute flow
The previous system can be solved numerically.Fig. 2 shows the variation of the re-scaled g profiles ḡ(z) = g(z)/ max(g bag ), z = z/H, with the layer thickness H/d.g bag stands for the fluidity profile associated to the Bagnold solution, corresponding to the local formulation, for which u(z We observe that the profiles superimpose on the characteristic Bagnold profiles far from the wall, but close to the wall the condition ḡ(z = 0) = 0 creates a boundary layer whose relative width zbl = argmax z (ḡ) decreases with H/d.We want to use the fluidity profile to characterize the velocity profile u in this boundary layer.ḡ depends on two variables: the column height H/d and the reduced ordinate z, so ḡ = ḡ(H/d, z).From Fig. 2, we observe that in the boundary layer region ḡ can be well approximated by a straight: In addition, one can verify linear dependence between ∂ḡ/∂z on the wall and H/d: Thus the velocity profile u(H/d, z) on the boundary layer can finally be obtained from the definition of fluidity ∂u/∂z = μg, which gives Now, assuming μ nearly constant (it should actually equal tan(θ)) and using (10), we finally get 1 open source, available and documented on the website [ The model proposed in [6,7] to take into account nonlocal effects introduces a boundary layer close to the wall, whose relative width decreases with H/d and whose growth behaves like H/d • z2 .Hence, imposing a zero shear rate at the boundary is crucial to recover the H stop phenomenology.
Note that using the steady-state approximate version of the non-local model does not allow for a true arrest of the flow (by contrast with the full model, see [7]).Hence, a cut-off value for the velocity must be introduce to obtain an estimate of H stop .

Discrete simulations
A contact dynamics algorithm is applied to simulate the flow of perfectly rigid grains interacting through contact friction and non-elastic collisions [4,9].We consider 2D simple chute flow configurations.The grains are perfect disks of diameters distributed around the mean value d within an interval of ±0.4d.This slight size difference allows for the introduction of generic disorder without inducing significant size segregation at the time scales considered.They form a bed of height H and width L with periodic boundary conditions in the x direction (see illustration in Fig. 1), tilted at an angle θ.H is varied between 4d and 35d to recover the H start / H stop phenomenologies, as well as to evidence finite size effects in the shape of the velocity profiles.The numerical parameters setting the values of contact friction and restitution were kept constant and set to μ = 0.In the following, we perform series of discrete simulations to obtain insights on the structure of granular flows close to arrest and systematically compare them to the outcome of continuum modeling.

Comparing discrete and continuum flow behavior 4.1 Transient arrest
In a first configuration, we consider a discrete layer of grains with H = 35d flowing at an angle θ = 20 • .For this set of height and slope, the flow is steady and the velocity profiles is consistent with Bagnold's prediction [1].We then tilt instantaneously the layer at a lower angle θ = 15 • , for which the flow is no longer sustainable and comes to rest.During the deceleration phase, the velocity profile is plotted at consecutive instants averaged over a time interval T S = 0.15T D , where T D is the full duration of the deceleration phase (Fig. 3).We observe how the velocity profile evolves from a Bagnold-like shape to lower velocity nearly linear shape, and finally changes convexity before vanishing.The same evolution is simulated using the continuum μ(I)rheology described above and implemented in Basilisk [2]  with corresponding parameters, without introducing nonlocal effects.For the same time periods spanning the deceleration phase, the velocity profiles are plotted in Fig. 4.They compare well qualitatively with the discrete evolution shown in Fig. 3, and in particularly we observe that the change in the flow structure is reproduced.This suggests that the transient arrest of a granular flow does not require introducing non-local effects, as shown in [10]; the I dependence from the μ(I) is a sufficient physical ingredient.

Steady flow close to H start / H stop
In the previous configuration, the flow evolves towards arrest in time, and eventually comes to rest.We now investigate the case of a layer whose thickness is close to the critical case for which no flow occurs at all, but which is nevertheless able to flow in steady regime.Thereby we want to evidence the possible emergence of changes in the flow structure that could be seen as for-runners of the arrest condition.We consider discrete granular layers tilted at an angle θ = 20 , where u 0 is the bottom slip velocity (Fig. 5).We recover self-similar profiles resembling the Bagnold prediction, with no non-local effect noticeable on the velocity profile shapes.We now perform a similar experiment applying continuum modeling.Introducing non-local effects in the model is necessary to recover the H stop phenomenology.For the continuum flow of a granular material similar to the one from [7] at θ = 24 • , Fig. 6 shows the change on the velocity profile shape for values of H close to H stop .Profiles are re-scaled by the max velocity value at the top of the layer, u(H), and superimposed.As H → H stop , there is a change of concavity near the wall, as expected from the analysis in Section 2.2.We notice a qualitative difference of behavior between discrete results for H start and continuum results for H stop .

Conclusion
Using a continuum model, we were able to reproduce the behavior of granular chute flows down an inclined plane for the steady regime and for the transient arrest using the local μ(I)-rheology model.These results were tested against discrete contact dynamics simulations.The local model, however, is not enough to reproduce the stopping height phenomenology.Thus we implemented a non-local model based on the propagation of a physical quantity called granular fluidity [6,7], which allowed us to get the desired H stop for granular media spreading a condition of zero shearing on the contact with the wall.We studied the impact of the steady-state-only approximate version of the model from [6,7] on the velocity profile shapes as H → H stop .In order to question whether the observed change of concavity near the bottom is physical, we performed discrete numerical simulations using contact dynamics for H close to H start .The discrete simulations exhibit consistent self-similar Bagnold profiles, until a sharp bifurcation causing the grains not to flow anymore occurs.No non-local effect was perceptible on the shape of the velocity profiles.These preliminary results suggest that the flow arrest may not be fully understood within the framework proposed in [6,7].Additional work need be done, including solving the full time-dependent version of the model.It is worth stressing that close to H stop , only few grain diameters do compose the layer, situation that might be challenging for continuum modeling in general.

Figure 1 .
Figure 1.Diagram of the physical problem studied: dense granular flow down an inclined plane.Gray particles correspond to the periodicity of the domain.

Figure 2 .
Figure 2. Numerical resolution of the granular fluidity g shape varying the number of grains.In black, the fluidity without the non-local model.g(z = 0) = 0 is imposed when solving(7).
5 and e = 0.5.The bottom of the flow is made of fixed glued grains of diameter d, forming a bumpy bottom condition.The dimensional values used for the simulations are d = 5 × 10 −4 m, gravity | g | = 9.81 m/s 2 and grains density ρ = 0.1 kg/m 2 .

Figure 3 .Figure 4 .
Figure 3. Discrete simulation: evolution of the velocity profile for a thick-layer flow towards arrest.
Discrete simulation, normalized velocity profiles for flows with H close to H start at an angle θ = 20 • .Profiles observed are self-similar.Figure 6. Continuum simulation, normalized velocity profiles for a different material than the one from Fig. 5.At an angle θ = 24 • , behavior as H → H stop is characterized by a change of concavity near the wall.Here H stop ≈ 2d.other cases a steady flow develops, giving an estimate of H start ≈ 10d.Thus, we are studying configurations for which H can reach remarkably close to H start .For those cases in which the layer flows, namely H/d = 10 to 16, we superimpose the velocity profiles after re-scaling them as ũ • , and of different thicknesses H = 8d, 9d, ... , 16d.For H = 8d and 9d, no flow is observed, but for all