Model of the saltation transport by Discrete Element Method coupled with wind interaction

We study the Aeolian saltation transport problem by analysing the collision of incident energetic beads with granular packing. We investigate the collision process for the case where the incident bead and those from the packing have identical mechanical properties. We analyse the features of the consecutive collision process. We used a molecular dynamics method known as DEM (soft Discrete Element Method) with 20000 particles (2D). The grains were displayed randomly in a box (250X60). A few incident disks are launched with a constant velocity and angle with high random position to initiate the flow. A wind velocity profile is applied on the flowing zone of the saltation. The velocity profile is obtained by the calculi of the counter-flow due to the local packing fraction induced by the granular flow. We analyse the evolution of the upper surface of the disk packing. In the beginning, the saltation process can be seen as the classical "splash function" in which one bead hits a fully static dense packing. Then, the quasi-fluidized upper layer of the packing creates a completely different behaviour of the "animated splash function". The dilation of the upper surface due to the previous collisions is responsible for a need of less input energy for launching new ejected disks. This phenomenon permits to maintain a constant granular flow with a "small" wind velocity on the surface of the disk bed.


Introduction
Drifting sand is a serious threat to the life of many people and may seriously disturb infrastructure and technical installations.The dunes invade houses, agricultural land and perturb the circulation on the roads.It is therefore very important to understand the mechanism of saltation which is the primary way, in which sand grains are propelled by the wind along the surface in short hops [1,2].The saltating grains are very energetic and when impact a sand surface, they rebound and consequently eject other particles from the sand bed.The ejected grains, called reptating grains, contribute to the augmentation of the sand flux.Some of them can be promoted to the saltation motion.
The main outcomes of our numerical studies follow.Firstly, we make a short overview of the previous experimental and numerical results concerning the successive collision between incident beads and a granular packing.Then, we describe briefly our mechanical model based on the Discrete Element Method and the special wind conditions used in this study.Finally, we present the main results available from this study.different authors.They studied the effect of the impact angle [3] and the incident speed [4,5] on the collision process.The authors have used steel beads (4 mm by Mitha et al [3]), PVC beads (6 mm by Rioual et al. [5]) and sand grains (850 μm by Werner [4]).The experimental results showed that the mean coefficient of restitution (COR) for the impacting bead (defined as the ratio between the rebound speed and the incident speed) depends only on the impact angle.The restitution coefficient decreases with increasing impact angle.The vertical one can be larger than 1 for grazing incident angle.Concerning the ejection process, Mitha et al. [3] showed that the number of ejected beads does not vary significantly when the impact angle increases from 17 • to 31 • and that the average of vertical speed of ejection is of order of 3 g.d(where g is the gravitational acceleration and d the bead diameter).However, the number of ejected grains increases with increasing incident speed [5].
In addition to the experiments [6,7], some simulations of the collision process have been performed [8][9][10].The authors made simulation using a discrete element method (DEM).The numerical simulations allowed to extract a law for the incident and ejected beads behaviours The numerical results confirm the experimental result, in particular: -The dependence of the mean normal restitution coefficient for impacting bead with impact angle.This coefficient is inversely proportional to the sinus of the impact angle.
-The mean number of ejected beads increases linearly with incident speed and impact angle.
-The distribution of the velocities for the ejected beads can be fitted by an exponential law.These simulations were limited to a single shot and a single grain size but the DEM simulations can be easily extended for multiple grains sizes or for sequential or simultaneous shots of particles on partially immobile packings.This second case is the goal of our present study.We study successive collisions of incident energetic beads with granular packing in the context of Aeolian saltation transport (i.e. with the presence of lateral wind).We have limited our simulation to the 2D case to minimize the sample size and also the global time calculation, but as seen in the 3D splash [10], the lateral dispersion remains symmetric which allows to neglect this direction if we don't take into account the bed packing shape (i.e no counter effect in the flow direction).

Description of the numerical simulation
The full description of the different steps of the numerical simulations to mimic the collision process between an incident bead and a packing of identical beads was already made by Oger et al. [10].So we will briefly present here the two main steps of the numerical simulations.Indeed, in order to save time of computation, we have simplified the technique to create the static packings of disks prior to the active part of the simulation; we perform two consecutive and complementary steps: a pure geometric piling step then a dynamic one.The second step can be separated by an initiation one and a flowing one according how the flow is maintained.

Initial conditions
Firstly, we generated a pure geometrical piling of disks with a Gaussian distribution of diameters to avoid crystallization in which the particles are just touching each other.This small polydispersity is a partial answer of the natural bidisperse natural sand gains in dunes.From a mechanical point of view, this stage corresponds to a dense packing of perfectly rigid, immobile grains (i.e.no overlap yet).The total number of disks used for the packing is between 20.000 and 30.000 (150 − 250 large by 60 − 100 height).Periodic boundary conditions are assigned to the lateral walls during both the piling and the saltating simulations.
DEM model: Soft model approach Our mechanical model is similar to the twodimensional formulation of Savage [11].The ith particle is characterized by its radius r i , the position of its centre (x i , y i ).A "soft-particle" approach is used, where each particle can have multiple contacts that can persist for extended durations.Both normal and tangential forces develop at the contact between two particles.The normal and tangential contact forces increase as the centres of the particles approach each other.
The normal force F n at the contact is modelled as viscoelastic.It consists of an elastic contribution (a linear spring K n ) and viscous damping (a linear dashpot coefficient b n ).The force in the tangential direction is modelled also as a viscoelastic one; a linear spring (with constant K t ) and a linear dashpot (with a constant b t ) are used to generate a tangential contact force.The tangential force is also limited to a maximum value which is chosen according to a Coulomb friction law when slipping can occur according to the inter-grain coefficient of friction.For this mechanical model, the time step calculation remains constant.
It is convenient to cast the governing equations in a non dimensional form and perform the computations based upon these dimensionless equations [11].Hence, all lengths are nondimensionalized by D, the diameter of the largest sized particles used in the computations.Time is nondimensionalized by dividing by √ M/K n , where M is the mass of the largest particle.Velocities are nondimensionalized by dividing by D/ √ M/K n .This dimensionless approach lets us model either large bead collisions such as in Rioual et al. [5] experiments (6 mm) or sand grains (200 μm) [12][13][14].For each particle, the forces acting on it are calculated; then the particle is displaced according to the resultant of all these forces.Local displacements of each particle are recorded.Wind interaction For simplicity, in the initiation regime, the wind profile is definedfrom the calculation made by Creyssels et al. (see figure 9 in [6].Then, but when enough grains are moving, it is possible to calculate the horizontal wind acceleration on the grains and to deduce the wind profile from this profile. Running conditions So we assign to each disk, physical parameters corresponding to those of classical sand grains.The disk diameter is taken equal to 200 μm and the density is 2500 kg/m 3 .The coefficient of restitution and the friction coefficient will be selected according to the desired studies (inside a range of 0.75 − 0.95 and 0.2 − 0.4 respectively).Generally the COR for sand can be smaller (0.6-.7) [15]  -A bead becomes an ejecta if its centre moves, from its initial coordinate, more than one diameter in the vertical upward direction.
-In order to limit the possible flights of moving disks, we have defined a horizontal upper limit like a perfect mirror for the moving disks (same horizontal velocity component and reverse vertical one after touching the mirror).These conditions are biased compared to natural parabolic flights but the number of grains involved in this situation is always very small and avoiding infinite vertical dimension helps a lot in limiting the global time calculation.following some referee comments, incorporation of the results of the cut part of the upper parabolic flight can be made to mimic more precisely the grains acceleration due to the wind.
The initial shooting sequence is activated and the collision process proceeds: Firstly a consecutive series of launched disks at fixed angles and velocities and random horizontal positions is made.This part is to initiate the general movement of the disk knowing that the wind interaction is already present to accelerate the moving disks.Secondly, only wind interactions with the beads are kept.So, after that, a stabilised regime can occurs if the full wind interactions (i.e.wind acceleration and fluidized zone deceleration for the grains) are not too strong.During this stabilisation period, we start the calculi of the horizontal wind acceleration on the disks (Figure 2) and by consequence the pressure applied on them (Figure 3).

Results
There are many adjustable parameters in this study that can be tuned such as the restitution and friction coefficients, the stiffness of the grains, the wind velocity at high altitude (i.e.independent of the grain interaction), and so on...In previous researches, we have seen that the amount of moving disks can increase according to the energy flux increase.So the initial amount of input energy has to be limited to control and reach the desired initial stability of the global process.As already described, the initial initiation configuration is not defined as a stable one for the grain flow according to the wind velocity but the system will be able to reach this equilibrium if the wind input energy can be exactly counter-balanced by the grain friction and the loss of energy during the splash collision and the flowing graingrain ones.Of course, this evolution is only possible for given experimental conditions: the wind velocity has to be high enough to compensate the decrease of energy of the moving grains resulting from collisions with other moving grains or with the quasi-static bed of grains, otherwise the system will slow down and stop.Figure 4 confirms this evolution to a stability regime by a decrease of the granular temperature from the regime initiated by the forced energy input (i.e.high granular temperature due to launched grains) and the consecutive regime in which only action and counter-reaction between grains and between grains and wind are taken into consideration.
To perform some series of analysis of the input parameters, we have observed the evolution of the results of the granular temperature by only changing the restitution coefficient of the disks.
Figure 5 confirms that the increase of the restitution coefficient can keep more disks in movement and that this higher amount of moving disks can be maintained even during the stabilisation regime: i.e the granular flux is higher when the restitution coefficient is higher and up to a maximal value in which too many grains are moving inside the dense flowing zone.In this case, a large number of collisions between grains can damp the granular temperature   and generates in consequence a kind of saturation value for the higher possible granular temperature.Figure 6 shows the evolution of the number of moving disks with time up to the apparition of the stabilisation regime obtained with three different values of the restitution coefficients keeping constant all the other parameters of the model.We can see that the range of parameters are correctly selected as they generated a perfect equilibrium regime: i.e. after the initiation and stabilization regime (up to t= 0.15s), the number of moving grains is kept almost constant and much lower than the total number of disks available in the bed. Figure 7 shows the evolution of the number of moving disks versus time with different wind velocities keeping constant all the other mechanical parameters.So we have partially defined the range of possible couples of input parameters in which we can observe the constant stability regime of the saltation.Right now, we have to continue the observation of the transition zone between the dense quasi-static packing and the fully developed saltation one when the system has reached the right saltation regime.

Conclusion
Previously, we have demonstrated that using the saltating splash function defined for a unique shoot of a grain on a static grain bed is not enough to handle the full complexity of the saltation process.In the same manner, the energy kept by the bed decreases with higher global energy input.The regime with a constant number of disks inside the flying zone requires an energy balance between grain-wind friction, grain collisions and wind input energy.

Figure 1 .
Figure 1.Schematic representation of the disk packing (2).The lower part (3) is the damping zone.The initial lateral wind profile is represented in the upper empty zone (1).

Figure 2 .
Figure 2. Analysis of the horizontal Velocities for beads going up and down at a given height inside the saltation zone.The difference is linked to the acceleration due to the wind during the parabolic movement of the disks

Figure 3 .
Figure 3.The calculated pressure obtained from the velocities difference according to the equations defined by Creyssels et al. [6].

Figure 4 .
Figure 4. Evolution of the granular temperature versus the height of the saltation zone for three different times of evolution for packing of disks with a restitution coefficient equal to 0.8 and a wind velocity equal to 40 m/s.

Figure 5 .
Figure 5. Evolution of the granular temperature versus the height of the saltation zone for three different restitution coefficients at the time t= 0.6s.

Figure 6 .
Figure 6.Influence of the restitution coefficient on the number of moving disks versus time with a true wind velocities = 20 m/s.

Figure 7 .
Figure 7. Influence of the velocities effect on the number of moving disks with restitution coefficient = 0.5 .