Discrete Element Method simulations of standing jumps in granular flows down inclines

This paper describes a numerical set-up which uses Discrete Element Method to produce standing jumps in flows of dry granular materials down a slope in two dimensions. The grain-scale force interactions are modeled by a visco-elastic normal force and an elastic tangential force with a Coulomb threshold. We will show how it is possible to reproduce all the shapes of the jumps observed in a previous laboratory study: diffuse versus steep jumps and compressible versus incompressible jumps. Moreover, we will discuss the additional measurements that can be done thanks to discrete element modelling.


Introduction
A standing jump is a rapid change in height and velocity in a free-surface flow, demarcating supercritical from subcritical flows.The phenomenon is well-known for hydraulic jumps on a smooth horizontal bottom, where the jumps obey Bélanger equation: h * /h = ( √ 1 + 8Fr 2 −1)/2, where h and h * hold respectively for the height before and after the jump, and Fr = u/ gh is the Froude number of the flow before the jump, with u the depth-averaged velocity and g the acceleration of gravity.However, this equation is not suitable for a granular flow which can occur only down an incline because of the frictional nature of granular materials, and which may be compressible [1].A recent study [2] confirmed the deviation from Bélanger equation for grains, and even for unfrictional incompressible flows like water when the bottom is rough or inclined.Some applications like the design of avalanche protection dams need an accurate prediction of the geometry of the jumps formed in a flow of dry granular materials whatever the incoming regime [1].As highlighted in [2], the jump height ratio does not depend on Fr only but is rather a function of a number of parameters: where L is the jump length, ζ the incline slope, ρ and ρ * hold for the density before and after the jump, respectively, and μ e is the effective friction (see [2]).The coefficients β, β * , k and k * will be defined in Sec.2.3.In particular, Eq.1 suggests that h * /h can be predicted for any incoming flow and boundary conditions provided L/h is known.
e-mail: segolene.mejean@irstea.frOur study aims to reproduce numerically a standing jump formed in a flow of granular materials thanks to the Discrete Element Method (DEM), in order to measure precisely the macroscopic parameters needed for Eq.1, and to decipher the internal structure of the jump.In granular media, the jump is defined as the part of the flow between the incoming flow, where the free-surface is parallel to the incline, and the outgoing flow, where the slope of the freesurface is constant and equal to the critical slope ζ 0 below which no flow can occur, as proposed in [1].For the largest ζ, the flow is slightly accelerating (non-uniform, but still at steady state).In those cases, the slope angle of the freesurface is however very small compared to the changes in the jump, and is not considered here (like in [1]).The main variables defining a jump are summarized in Fig. 1.The depth-averaged velocity after the jump is noted ū * , and the mean diameter of the grains is d.
Section 2 describes the DEM contact-laws, the numerical set-up used to create a jump, and the techniques to measure all the variables of interest.Section 3 shows that this method allows to create a wide variety of jumps, and presents some of the first results.Finally, a short conclusion is given on the main challenging issues this granular jump numerical set-up will allow to investigate.

DEM simulations of standing jumps
This section describes how granular flows are simulated thanks to DEM, using YADE open-source software [3].

Microscopic contact laws
In order to model the grains interacting each other, we chose the classical following viscoelastic law: In this equation, F n and F s are the normal and tangential forces.F n is the sum of a linear spring (of stiffness k n , and proportional to the normal overlap δ n ) and a dashpot of damping coefficient c n that depends on the restitution coefficient e. F s is incremented at each time step dt as a linear spring (of stiffness k s proportional to the derivative of the tangential overlap δs ) restricted to a Coulomb threshold force defined by the friction coefficients μ and μ b for grain-grain and grain-wall interactions, respectively.In this study, the tests are carried out in two dimensions (2D), with only one particle across the width.This is made by blocking one degree of freedom in translation and two degrees of freedom in rotation for each particle.

DEM set-up to produce the jumps
The numerical set-up is shown in Fig. 2. A reservoir permanentely filled up with grains feeds an incline.At the exit, a gate initially retains the grains.When there are enough grains in the incline, the gate moves up.Every 150 time steps, the outflow discharge is calculated and compared with the inflow.The gate then slightly moves, up or down to adjust the discharge and guarantee a steady state.This is a slight difference from the laboratory model where the gate was adjusted by the operator and fixed [1].The strong advantage of DEM is that many control variables can be changed to see their effect on the jump: the slope angle ζ, the height H of the reservoir exit (that controls the inflow), but also the grain diameter, the interparticle friction and the restitution coefficient.All parameters have a default value, and each parameter was varied-the other ones being fixed at the default value (see Tab.1).Table 1: Default values (first line) and ranges (second line) of the varying parameters for the 53 preliminary tests.

Geometry of the jump
When a simulation runs, the first step is to reach a steady state.Then, the free-surface is calculated by discretising the incline, and identifying the highest grain in each cell.This is made a hundred times at several time steps to obtain a smooth time-averaged free-surface at continuum scale.The latter allows to identify the beginning of the jump, which is the location where the free-surface is not parallel to the bottom anymore-when the derivative of the function describing it becomes non-zero (Fig. 2).We also identify the end of the jump, where the slope of the freesurface becomes constant-meaning that the derivative is constant.Once a steady state is reached and both the beginning and the end of the jump are well identified, a new simulation is done to record the variables of interest.

Velocity profiles
An example of velocity profiles recorded before and after the jump is shown in Fig. 3.The raw DEM data sets exhibit a noticeable scattering caused by the (discrete) fluctuating nature of the flowing grains, but the mean velocities are well-described at continuum scale and they can be fitted by a Bagnold profile with basal slip, as already suggested in [1,4], which is given by the following equation: where u b denotes the sliding velocity of the basal layer, which is taken at z = z 0 = 1.5d (following [1,4]) and A is the Bagnold pre-factor, which was let free in the fits of Fig. 3.The hypothesis of the Bagnold profile seems very appropriate for the velocity of the incoming flow.But Fig. 3 shows that it is less satisfying for the outgoing flow, because of a change in concavity next to the free-surface that will need further investigation in the future.
Measuring accurately the velocities thanks to DEM makes it possible to calculate precisely the Boussinesq coefficient, β before the jump (or β * after the jump), which is defined by the relation β = ū2 /ū 2 .The variations of β and β * with Fr and Fr * respectively are given in Fig. 4(a).It clearly shows that the usual assumption that β = 1 is very accurate for median values of the Froude number that we find most of the time in the incoming flows, while it may become an underestimate for very low values of the Froude number as found in the downstream flow.

Earth pressure coefficients k,k *
The earth pressure coefficients, noted k and k * before and after the jump respectively, relate the normal stresses through the relation σ xx = kσ zz .In this study, k and k * were not directly measured but derived from soil mechan- ics concepts, as proposed by Savage and Hutter [5]: where ϕ s = tan Because the inter-particle friction coefficient used in DEM is known to be different from the macroscopic (internal) friction coefficient ( [6]), our assumption for Equation 4remains questionable.Future work will consider sheardependent friction coefficients, such as developed by [7].

A rich variety of granular jumps
The method presented in Sec.2.2 allows to create a wide range of jumps, and confirms the results obtained in the laboratory [1] as shown in Fig. 5.We could observe steep or diffuse jumps, compressible or incompressible jumps, and identify the presence or absence of a recirculation zone.Figure 6 displays pictures of jumps-with the internal streamlines drawn, obtained for different ζ and H/d.

Jump steepness
The two characteristic length-scales of a jump are its height ratio h * /h and its relative length L/h (if not neglected like in Bélanger equation).Combining the two leads to the steepness coefficient S = (h * − h)/L with respect to the slope of the bottom.S can grow by increasing ζ, and weakly depends on H/d (Fig. 7).It is also correlated with the presence of recirculation as a very steep jump will be instable, thus promoting the recirculation.A larger opening of the reservoir leads to a higher discharge and a denser jump, more stable.As such, the recirculation appears for higher ζ and S. It is worthwhile to note that jumps with recirculation were also observed for very low μ, typically smaller than μ b .This result-likely to unravel the key role played by friction in the birth of jump recirculation, will need further investigation in the future.

Jump compressibility
The granular jump compressibility is defined as the ratio φ * /φ (where ρ = ρ P φ) between the volume fractions before and after the jump.φ * was systematically found to be close to 0.81 (whatever the incoming flow and boundary conditions), which is nearly equal to the 2D granular random close packing.Then, once the incoming flow is dilute (φ < φ * ), a compressible jump is obtained.Our simulations showed that it is possible to tend toward an incompressible jump (φ * /φ → 1) by different ways: increasing H/d, decreasing ζ, decreasing d or decreasing μ.

Conclusion
The current paper presented preliminary DEM simulations of standing jumps formed in flows of granular materials down an incline.The first results are consistent with the recent findings in laboratory tests [1].This novel numerical set-up now offers the possibility to study in detail how the jump geometry and its internal structure are influenced by a number of parameters, and will allow to investigate the energy dissipation inside the jump, as proposed in [2].

Figure 1 :
Figure 1: Sketch of a granular jump and notations k n and k s are calculated from the values of Young's modulus E-which was taken equal to 1 × 10 6 to respect the limit of rigid grains condition, and Poisson coefficient ν taken equal to 0.3: k n = Ed/2 and k s = νk n , considering one diameter d for all particles.A polydispersity of 15% around the chosen diameter d = 4 cm is considered to avoid crystallization effects.The bed-friction coefficient was μ b = 0.25.The grain density ρ P = 2500 kg m 3 to mimick the density of glass beads usually used in laboratory models.The time-step was dt = 3 × 10 −4 s.

Figure 5 :
Figure 5: h * /h versus Fr obtained with the present numerical experiments (red circles) compared to the laboratory experiments conducted by [1] (green stars).
−1 (μ b ) is the bed friction angle and ϕ e = tan −1 (μ) the internal friction angle.When μ becomes lower than μ b , we can consider k = k * = 1.The variations of k and k * with μ according to Eq.4 are plotted in Fig.4(b).

Figure 6 :Figure 7 :
Figure 6: (H/d, ζ) phase diagram for granular jumps.In each picture of the jumps, the streamlines are drawn, thus allowing us to identify the presence or absence of a recirculation zone