Granular front for flow down a rough incline : about the value of the shape factor in depths averaged models

The thin layer (or Boundary Layer scaling) approximation of the Navier Stokes equation with μ(I) rheology for dry granular flows over an inclined plate is presented. It is called "Granular RNS/P" (Reduced Navier Stokes/ Prandtl). Integrated over the depth, it gives the classical depth averaged equations (Shallow Water type equations, popularized by Savage Hutter). But the process of depth averaging needs hypothesis on the velocity shape profile trough a "shape factor coefficient". In the case of the displacement of the front of an avalanche down a rough incline, we present comparisons of "Granular RNS/P" simulations versus depth averaged model (which compares well with experimental data). One conclusion is that the shape factor goes from 5/4, the Bagnold value, to a unit value, corresponding to a flat velocity profile, in a region very close to the front of the avalanche.


Introduction
Considering dry avalanches, landslides, rockfalls, snow avalanches, debris flow, or general geophysical flow as depth averaged thin layer equations is a very common and powerful point of view.These simplified continuous models contain the dominant physics: gravity driven and Coulomb friction dominated flow.Therefore the system is smaller and faster to solve numerically.Indeed, the case of dry avalanches (the one we will look at here) is up to now too demanding for complete continuous flows such as in [10].The obvious reason being the small aspect ratio of the avalanches, they are thin and long.This induces numerical difficulties.
As a practical example of such flows reproduced in the laboratory, we have revisited Pouliquen [12,13] seminal work (small-scale model of a granular layer of glass beads flowing on a rough incline, see figure 1).In [17] the experiment set up is reproduced.In this paper, larger Froude numbers are tested.It is shown that the solution for the global shape of the front of the avalanche proposed by [7,12,13] is no more valid for Froude numbers of order one.The analysis is based on depth averaged equations in which several hypothesis are changed.One is about the "shape factor" (defined as α latter, see Eq. 3).The usual approximation is to take a shape factor of value one (supposing a "plug flow").In [17] it is shown that depth averaged equations may be solved with an exact implicit solution with a any shape factor α value.From comparisons with experiments, the value α = 5/4 (corresponding to a Bagnold flow) is found to be the best one.The effect is clear for Froude numbers of order one.For small e-mail: pyl@ccr.jussieu.frFroude, the value of α plays no role.This good comparison fails in the vicinity of the font it self.If one uses the value α = 5/4, the very front of the avalanche is not sharp.It presents an unphysical exponential decay.If α = 1, the very front of the avalanche is sharp.It presents a finite slope, which is more physical.This is one of the reasons for the use of α = 1 in the literature.
In this work, we will present a set of equations "Granular Reduced Navier Stokes/ Prandtl" using the frictional rheology μ(I) ( [6,8]).These equations have been proposed by [4].This set comes when reducing the number of terms in Navier Stokes thanks to thin layer approximation (Prandtl Boundary layer scaling).We present a numerical method to solve this system.An integration across the depth gives then the classical depth averaged model.We will then numerically solve these "Granular RNS/P" equations and compare to depth averaged solutions.The outcome of the simulation will be the "shape factor", it will be a result rather than an hypothesis.We will see that it changes from α = 5/4 to α = 1 in the close neighborhood of the front.This explains why the results of [17] are good for α = 5/4, and removes the unphysical exponential decay as α → 1 at the front.

Granular RNS/P equations
First we present the equations.We start from the continuous Navier Stokes incompressible model (as in [10]) with the μ(I) rheology.We write the thin layer system of equations associated as introduced by [4] recently.It means that some derivatives are dropped in the equations leading A rough sketch of the experimental setup of ref [17] describing a granular avalanche of hight h(x, t) of discharge Q(x, t) along a plate at angle θ.The position of the front is x f .to a kind of boundary layer system of equations.It results in the following system with dimensions, we call it "Granular RNS/P" (Granular Reduced Navier Stokes/ Prandtl): Z b (x) is the topography on which flows the granular fluid.
The constitutive law for the stress is here τ xy = ν eq ∂u ∂y with ν eq = μ(I)p ∂u/∂y , it depends on the inertial number ∂u ∂y (where d is the diameter of the grains) through a given μ(I) function: , Jop et al. [8] for details on μ(I)rheology.The boundary conditions are no slip at the wall u = v = 0 in y = 0, and no stress at interface p = 0 and ∂ y u = 0 in y = h(x, t), plus any initial velocity field.It is clearly a set of equations deduced from Navier Stokes equations in which we have incompressibility, the convective term, the hydrostatic pressure term and the dominant part of the derivative of the stress which follows a Druker-Prager law.This system is reminiscent to the RNS/P system ( [3,9]).This terminology comes from the efforts to reduce the set of Navier Stokes in adding only the transverse derivative of the viscous forces in the equations (see Rubin et al. [15,16]).In simple words, these equations are no more than the Prandtl Boundary Layer equations (see the fundamental book of Schlichting [19]), with the gravity terms as sources and different boundary conditions.

Numerical discretisation
System (1) is a kind of Boundary Layer system.Numerical resolution of steady boundary layer equations presents several traps, for instance it is well know that there is a singularity when steady equations are solved with a given pressure (Goldstein singularity, [19]).Furthermore, the unsteady boundary layer equations are either instable or may even lead to a finite time singularity if pressure is given.Fortunately, Audusse et al. [1] proposed recently a resolution of system (1) extended by [4].The idea of resolution is in fact a clever transverse discretization of (1) in several layers, each layer being a Saint-Venant one (so that in this paper, there is not any reference to boundary layer theory but it is focussed on Shallow Water literature).Those layers are solved by finite volumes.Each layer interacting one with the two others: the upper and the lower, through the viscous force (viscosity ν eq ), and by cross exchange of mass and momentum (in reference [1] this exchange is done by the functions G i which are the difference between the velocity of the moving layer and the real transverse velocity).This implementation has been shown to be very robust.

Depth averaged models (or Saint-Venant Savage-Hutter 1D models)
We have defined a full 2D problem (1) and we have proposed a numerical method to solve it.The usual model is the 1D depth averaged model.The 1D depth averaged model is a simplification of system (1).Indeed, integrating equation ( 1) from y = 0 to y = h(x, t), we have without new hypothesis: (2) At this point of the analysis comes the choice of the velocity profile to close the system to obtain the 1D model.In aerodynamics, the integral equations arising are the Von Kármán equations, and the closure maybe the Pohlhausen one (see [19]).In hydrodynamics, the integral equations arising are the Saint-Venant (Shallow Water) equations, and the closure involves plug flow and turbulent friction (or half a Poiseuille in viscous films).In granular flows, the integral equations arising are the Savage-Hutter equations.First the basal friction term τ xy | 0 has to be modeled, usually since [18] by a Coulomb friction: τ xy | 0 = μρgh (with a constant friction coefficient μ).Second, the shape of the profile has to be defined to link the flux of velocity (the volume flow rate) to the flux of kinetic energy, defining the shape factor (as a constant in those models): ( The usual approximation is to use a plug profile so that if we define the discharge as Q = h 0 udy then h 0 u 2 dy is approximated by Q 2 /h (i.e.α = 1 in the following equation, if the profile is supposed to be a Bagnold one, then α = 5/4, see [17]): These are the Savage-Hutter Saint-Venant equations for granular flows [18].This hierarchy of simplifications is justified by the fact that the avalanches are very thin compared to their longitudinal extent.Indeed, numerical simulation of Navier Stokes equations with μ(I) rheology (as in [10]) is very difficult, and numerical simulations of these thin layer approximations ( [4] or [18]) are more tractable.
The approximation α = 1 is used in nearly every paper as it helps the numerics and as it gives a sharp front (Poulinquen [12,13]).

Analytical solution of the front shape, comparison with experiments
It has been shown in [17] that an analytical solution of Savage-Hutter Saint-Venant system (4) may be obtained in the case of the flow over an inclined flat plate with Z b = − tan θ.It is an implicit solution were the position is function of the reduced hight: valid for any α, see reference [17] for details and for the exact expression of F which involves arc tangent and logarithm functions.This implicit analytical solution extends the Gray and Ancey's one [7] valid only for plug flows or for small Froude.
In [17] a flume of 2 meters long and of 40 cm of width, with a variable inclination angle starting by a silo has been build.Several avalanches at various angles, and with various initial height of aperture of the silo have been performed (see reference for details).With optical techniques (laser sheet), profiles of the front span wise and stream wise were measured.The aim was to explore larger Froude number values.The constant front velocity as in [12] was reobtained.Very good comparisons between measured fronts and analytical solution Eq. 5 for values of Froude of order one and α = 5/4 have been obtained.The effect of increase of inertia (by increasing the angle) is clearly reproduced by the analytical solution.
One of the drawbacks of the solution Eq. 5 is the existence of a precursor film with an exponential decay, which is not observed in the experiments.For α = 1, there is no such precursor film.

Discussion
In the depth averaged formulation of system (4), the shape factor is constant by construction.With this formulation an analytical solution may be obtained, this is Eq. 5.
It allows quantitative comparisons with experimental data as just explained and presented in [17].The sole problem being this spurious exponential decay, which is not observed in the experiments.Hence, now we turn to the numerical solution of (1) in the case of the front.We compute the shape factor as a result, it is no more an hypothesis.On figure 2 several computed shape factors are ploted for different values of Froude number.We observe that indeed the computed α from (1) is mostly equal to Bagnold's value 5/4.Nevertheless, close to the front it goes to one (indeed, longitudinal velocity at the very front is more like a plug flow instead of a Bagnold profile).
On each plot of figure 3, three curves are presented.First is the analytical solution with α = 1 or Fr = 0 (that is the same) corresponding to the Savage-Hutter Saint-Venant solution [7,13].This is for reference.Second is the present numerical solution of system (1).Third is the analytical solution with α = 5/4 and the corresponding Fr value from the second.As we can notice, the second and the third are very close.It is only very near the front that we see the differences.As the value of α is fixed to 5/4, the solution of equation ( 5) presents a slight difference when approaching the front.Furthermore, the solution of equation ( 5) presents the exponential unphysical precursor layer.The agreement between (system (1) and equation ( 5) with α = 5/4) and Savage-Hutter Saint-Venant (equation ( 5) with α = 1) is only good for small values of Froude number (Fr .5).Then, the higher the Froude the smaller the slope of the front and the larger the departure from system 1 to α = 1 models.

Conclusion and perspectives
This paper is devoted to the solution of a reduced system (1) issued from the non-Newtonian μ(I)-rheology of In red h(x, t = 150)/h max as function of (x − x f )/h max at same time (t = 150) for various increasing θ of various increasing initial thickness giving a Froude from 0.348 to 2.09 (from top to bottom).The analytical solution (equ.5) is plotted as well in blue for the corresponding values of θ and Fr (with α = 5/4).For reference the analytical solution (equ. 5 with α = 1 or F = 0) of the usual Savage-Hutter Saint-Venant resolution is plotted in black and labelled F = 0 (colors online).
the Navier-Stokes equations leading to kind of Prandtl equations.This system is showed to reproduce the expected behavior in the case of avalanche: in steady conditions far from entrance and far from the front we obtain the Bagnold solution.Our method allows to focus on the front, we notably compute the local value of the shape factor α. It starts from 1 in the very beginning of the front.This corresponds to a very first plug flow.It goes then to 5/4 the Bagnold expected value at constant height.The numerical resolution uses the robust multi layer approach in finite volumes as proposed by [1] and [4].The numerical resolution is shown to be very close to the analytical solution established in [17], this analytical solution was very close to experimental data.The conclusion is that the usual Savage-Hutter Saint-Venant equations (with α = 1) are not a very good description of the avalanche front at Froude of order one.A better description is with α = 5/4, but there is a small unphysical exponential tail.A better description is this Granular RNS/P description.Note that reduced system (1) seems not to have the instability, or "ill posed" behavior found by [2] in the full system.
Finally, those equations are in between full resolution (Navier-Stokes) and a too much simplified one (depth averaged Savage-Hutter Saint-Venant).They create the missing link between a too complicated Navier-Stokes resolution (useful to focus on details) and the Savage-Hutter Saint-Venant.This latter depth averaged method is too crude maybe, but it is fast and efficient enough to compute a reasonable and realistic description of an avalanche on an actual topography [11].

Figure 2 .
Figure 2. Numerical resolution of "Granular RNS/P".Plot of α the shape factor ratio of flux and kinetic flux (defined by equation 3) at a fixed time as a function of (x− x f )/h max the position before the front divided by the value of the thickness of the avalanche for three values of the Froude number (upper curve Fr = 0.43, middle curve Fr = 1.03, lower curve Fr = 2.09).Far from the front, the value starts from the Bagnold 5/4 value and goes towards 1 at the front on a scale of several thicknesses.The smaller the Froude, the better is the approximation α = 5/4, except close to the front.