Simulation of saltation motion using LBE based methods

The numerical model of the motion of the circular particle close to the bed in an open channel with a rugged bed based on the lattice Boltzmann method (LBM) is presented. The LBM is used as a DNS approach in which hydrodynamic forces are expressed as sum of contributions from fluid elements interacting with the moving particle. The corresponding numerical simulation for the saltation motion which represents a dominant mode of the bed load transport is developed. Flow is driven by the logarithmic velocity profile at the inlet of a two dimensional channel with a bed formed by semi-circles of variable radii in a bed of particles. Translational and rotational movements of the particle are induced by gravitational force on one hand, and by hydrodynamic forces on the other hand. The LBE (lattice Boltzmann equation) based simulation provides the opportunity to study the behavior of saltation motion in the moderate and high Reynolds number regimes. Most of the input parameters, including boundary conditions or flow conditions, are adjustable within a range of values. Stability issues of the simulation are considered and a resolution using a combination of different LBE models and the extension of computational resources is proposed. Finally, an enhancement of the simulation for more complex processes is suggested.


Introduction
Particle movement influences the flow of the fluid and its trajectory is affected by the flow in return. The influences can be classified within three types of interactions: particles with the fluid, mutual (particle-particle) collisions and particle collisions with the bed (particle-bed). Simulations of particle motions and their interactions close to the bed represent a significant focus recently in hydrodynamics.
Bed load transport is defined as the motion of particles close to the bed which occurs in three different modes: rolling, saltation, and suspension. These modes merge into themselves mutually either in the mentioned or the reversed order depending on flow conditions (shear velocity, fluid density and viscosity), particle attributes (size, shape and density) and bed characteristics (roughness and ruggedness expressed in terms of bed particles size and their respective packing density). As particles within a single cluster differ in shape, size and density, they can move in various modes in the same time.
It has been shown that saltation represents the dominant mode of the bed load transport. Saltation is a complicated motion of particles which consists of free motion in the flow, mutual collisions of particles and collisions with the bed. The characteristic pattern of the saltation mode is determined mainly by particle-bed collisions.
The lattice Boltzman method (LBM) is a two decade old developed numerical approach originating from the lattice gas automata methods used for the simulation of complex fluid flows, e.g., [1][2][3]. LBM represents a second order, efficient computational scheme due to its inherent locality and explicitness. Moreover, the possibility of straightforward parallelization yields another considerable advantage of the traditional numerical approaches to fluid flow problems. a e-mail: dolansky@ih.cas.cz Traditionally, models of saltation are composed from four ingredients: equations for the fluid flow (e.g., Navier-Stokes), equations for motion for particles in the flow (derived from Newton laws), impulse equations describing particle-bed collisions and a particle-particle collision model. These equations are then discretized by finite difference or element methods and standard numerical methods are applied to the resulting algebraic equations.
In contrast, LBE methods offer quite a general and robust tool to solve the above mentioned stages of the saltation motion within a unified frame when it is used as a direct numerical simulation (DNS) technique. Particle-fluid systems represent a significant category of applications of the LBE methods that have recently been developed and examined.
The aim of this work is to develop and test a LBE based simulation of the motion of the particle in the flow founded on the mesoscopic -scale between molecular and continuous (macroscopic) domain -model of the fluid. The mathematical model is designed for the saltation motion in an open channel and implemented employing LBM. It is shown that LBM can be equally used to model saltation motion and the resulting simulation produces likely and expectable results in comparison with the traditional approaches mentioned above without a loss of efficiency of the computations.
From the reason of the initial testing of the LBE based simulation, it is supposed that particle to particle collisions occur quite rarely, and hence it is sufficient to consider the motion of the single particle in the flow. The flow is driven by the logarithmic velocity profile at the inlet of a two dimensional channel with a bed formed by semi-circles of variable radii in a bed of particles. The saltating particle has a circular shape and it is not deformable. Translational and rotational movements of the particle are induced by gravitational force, and by hydrodynamic forces. The latter forces are modeled on the mesoscopic level as interaction between the fluid elements and the particle, see Sect. 3.
The LBE based simulation -after some modifications, see Sect. 5.2 -enables the chance to study the behavior of saltation motion in the moderate and high Reynolds number regimes. Most of the input parameters including, boundary conditions or flow conditions, are adjustable within a range of values.
This paper proceeds as follows. In the next section, the numerical model for flow and particle motion and its collisions with the bed are described at the macroscopic level. In Sect. 3 essentials of LBM are introduced. In Sect. 4 the numerical model of the motion of the particle in the flow is implemented in the mesoscpic level, i.e., using the LBE approach. Results and consequent observations are presented in Sect. 5 which consists of graphs demonstrating the outputs of the simulation and stability considerations concerning the simulation. Finally, the paper is summarized in Sect. 6 where also prospective plans are proposed.

Flow field
We assume that the flow is determined by the incompressible Navier-Stokes equations with constant kinematic viscosity ν. We suppose that the mean velocity of the flow is determined by the logarithmic law of the wall close to the bed as it can describe flows at higher Reynolds numbers. k = 0.4 stands for the Karman constant, u * for shear velocity and y 0 for the distance from the bed at which the velocity u goes to zero. The domain of the process is confined within the boundaries of two types: open boundaries (inlet, outflow and free water surface), and solid boundary (bed). Each boundary is represented by a different boundary condition.
Flow velocity is driven by the logarithmic velocity profile at the inlet of the channel, cf. Eq. (1). The velocity profile is induced long enough to achieve a steady state.
At the outflow, the Neumann free flow -the so called "do nothing" -boundary conditions are imposed which corresponds to the normal gradients of the velocity set to zero (n represents the outer normal to the boundary), e.g., [4]. These conditions are often used if any boundary conditions for the velocity at the outflow are not prescribed and represent a popular choice for open boundaries for the incompressible Navier-Stokes equations.
As it is supposed that the process takes place in an open channel; the water level is described by a free surface boundary condition which is characterized by constant pressure The bed of the channel is composed of half-circles with radii equal to the radius of the saltating particle. The noslip boundary condition -identification of the fluid velocity with the velocity of the corresponding surface -is supposed on the bed u(x bed ) = 0, as well as on the particle surface: where the subscript "bn" denotes the boundary nodes of the moving particle.

Particle motion
The motion of a particle in the flow is determined by actions of both body and hydrodynamic forces. They are summed into a resultant net force F net which occurs on the righthand side of the Newton equations and thus determines the particle motion (m = ρ p V is mass of the particle). The net force results from the summation of body and surface forces where F g stands for gravitational force and hydrodynamic forces are represented by drag force F d , force due to added mass F m , the Magnus force F M and the Basset force F B . The net force develops a torque on the saltating particle which determines its surface angular velocity ω from relation where J is the particle momentum of inertia and r is the vector pointing from the center of the particle to a boundary point. It has been shown that rotation influences both motion and collisions significantly, and especially for higher Reynolds particle number, rotation has to be taken into account as it essentially affects particle trajectories [5].

Particle-bed collision model
Descriptions of particle-bed collisions use to be based on impulse equations which result into relations between velocities before and after a collision. If the approach by Niño and Garcia [6] is adopted then where θ in represents the incidence angle, θ b is the angle between the tangent line passing through the contact point and the bed line and θ r is computed from where e and f represents restitution and friction coefficients respectively. Since the bed is irregular -in this case made of circles -a random number generator is used to pick out one location from all possible contact points and thus to determine the angle θ b [7]. Rotation has to be taken into account and the particle angular velocity is expressed appropriately. It is supposed that the particle rebounds without the slip in the current phase of development of the simulation. In other words, the bouncing particle does not slip along the surface of the bed during the collision which allows to express the angular velocity after collision (also solution of the impulse equations) as where m tot represent the sum of the particle mass m and the added mass, and v t corresponds to the tangent velocity component just before the collision, cf. [7]. It should be mentioned that although collisions with the bed determine the characteristic pattern of the saltation motion; particle-particle collisions affect the characteristic parameters of the saltating motion as height and length of jumps or the vertical particle distributions.
However, as in this work the saltation motion of a single particle is considered, it is not necessary to define the particle-particle collision model though as it will be included into the developed simulation.

LBM intro
The LBE approach is designed to view flow processes from a mesoscopic (situated in between microscopic and macroscopic levels) view. Formulation of a numerical model based on the LBE method is determined by postulating mesoscopic structure determined by choice of: (a) a lattice grid specified by space dimension and nodes layout, e.g., quadratic, hexagonal, cubic, (b) a set of discrete velocities c i connecting the nodes, (c) a collision operator Γ i for interactions of fluid elements at the nodes and (d) implementation of boundary conditions, cf. figure 1.
Fluid elements that can move only along discrete velocities c i in a selected lattice are represented by distribution which gives the probability of finding a fluid element in location x with a velocity of c i in time t. Collision rules that guarantee conservation of mass and momentum are applied on the fluid elements at grid nodes. The collision and propagation process is represented by the lattice Boltzmann equation Due to collision; the distribution function f i , aligned in the direction c i at position x and time t, is updated with Γ i after which propagates to the position x + c i Δt at the next time step t + Δt. The collision operator expresses the propensity to the local equilibrium f eq and it is often approximated by a simpler Bhatnagar-Gross-Krook operator where change of the distribution function f i is proportional to its difference from equilibrium f eq . The parameter τ expresses the rate of tendency to local equilibrium and it must be defined so as to imply the desired macroscopic dynamics. The equilibrium distribution function f eq corresponds to the Maxwell distribution of molecule velocities which is expanded to the second order in the Mach number Ma = u/c s where parameters w i represent weights in directions c i normalized as i w i = 1. Hence, the so called collision step represents the collision of fluid elements in a node, or tendency of distribution function to the equilibrium distribution (8).
Macroscopic quantities (density ρ and velocity u) of the flow are obtained as moments over particle distributions It can be shown that such a simple collision-propagation scheme leads back -by the method of the asymptotic analysis -to the macroscopic Navier-Stokes equations in the incompressible limit. From the relations for collision operator (7) and the equilibrium distribution (8) it is evident that the calculation of the post-collision distributions f c i is localized, i.e., it depends only on the distributions in a node and its direct neighbors. This fact determines the parallel nature of the LBE based methods. Moreover, in the stream step distributions f i spread along straight lines with constant speeds c i (in a chosen lattice). Thus, collision step (9) along with stream step (11) represent a two-step algorithmization of the lattice Boltzmann equation (6).

Boundary conditions
The implementation of different boundaries types is more or less straightforward within the LBE frame. However, any implementation scheme has to allow passage -by the asymptotic analysis -from a mesoscopic formulation of boundary conditions to their macroscopic counterparts (within certain accuracy limits). Four boundary conditions mentioned above -open boundaries of the inlet, outflow and free surface side, and the solid boundary of the bed -require the usage of different implementation techniques. The outline of the process domain along with the boundaries and their implementation is depicted in figure 1.
The no-slip boundary condition on the solid surface of the bed made of semi-circles -as well as the particle solid surface -is implemented by the so called bounceback scheme which consists in the inversion of distribu- along the directions of incident to the boundary nodes. In the case of moving a boundary, the bounced distributions must be modified (12) due to the exchange of momentum between fluid elements and the moving particle, cf. Eq. (14), where v tot = v + r × ω denotes the total velocity at the boundary nodes. The implementations can differ in treatment of the object interior, e.g., [8,9], or its surface curvature [10]. The specified velocity profile at the inlet is implemented by the Zou-He boundary scheme [11] which allows to define velocity u(x in , t) or pressure p(x in , t) on the boundary. In this case, only the non-equilibrium part of distributions are bounced-back in the normal direction with respect to the boundary where index i runs over only unknown distributions. This assumption helps to express the distributions along with imposing macroscopic velocity/pressure conditions. The same scheme is used to implement the free surface boundary with constant pressure which represents the water level. The corresponding LBE schema to the "do-nothing" condition given in terms of the hydrodynamical variables, see Eq. (2), is derived with the help of asymptotic analysis and implemented by two relations for the direction opposite to the outer normal, c i = −n, and for the other directions c i −n, e.g., [12]. The equilibrium distribution f eq i , see Eq. (8), is defined for constant density ρ = 1 instead of fluctuating density given by Eq. (10). The lattice nodes on which the outflow is located are denoted x out and ν * is the lattice viscosity.

Particle motion
In the case that LBM is used as a DNS technique, movements of an object in the flow is the effect of interaction of fluid elements (represented by distributions f i ) with the solid surface of the object. The interaction consists of the action of fluid elements on the particle motion, and the action of the particle on the fluid elements, i.e., on the flow. The first action moves the particle and is determined by the net force F net of the fluid elements, Eq. (13). The second interaction is implemented as a special -moving -case of bounce-back boundary conditions, cf. Eq. (12).
The particle motion is caused by the momentum transfer Δp from the incident fluid elements. Time rate of this momentum transfer defines the force by which the flow acts on the particle. The net force F net which enters the Newton equation (3) is the sum of the force due to the momentum transfer between the particle and fluid elements and the forces due to the nodes covered/uncovered by the particle within time step Δ t cf. figure 2. As in the macroscopic case, the net force is the sum (4) of particular hydrodynamic forces (drag, added mass, Magnus, Basset etc.). The momentum contribution Δp forms an impulse force on the particle given by with Δt = 1, cf. [9].
The impulse force F cov exerted on the particle due to covered nodes is given by the momentum sum from the previous time step Total momentum in uncovered nodes forms the negative impulse force increment where ρ(x un ) and u(x un ) are obtained as averages of these quantities over neighboring nodes, cf. figure 2(b).

EPJ Web of Conferences
x y u n c o v e r e d a r e a c o v e r e d a r e a Figure 2. The net force which determines particle motion is composed of three particular forces. (a) force due to momentum exchange in the boundary nodes of the particle. (b) forces due to the momentum of fluid elements in nodes covered by the particle and the momentum of uncovered particle nodes.

Integration of equations of motion
The hydrodynamic forces on the right-hand side of the Newton equation are the effect of interaction of fluid elements with the moving boundary. To update the particle position (x(t), y(t)) and velocity (v x (t), v y (t)) the left-hand side of the Newton equations has to be evaluated every time step. In other words, a computational procedure for integration of the Newton equations (3) is required. For this purpose, the leap-frog algorithm is chosen as it is simple, possesses second order accuracy, is invariant under time reversal and has other global properties, e.g., it is an area preserving map in phase space when integrating the Newton equation over a finite time, e.g., [13]. The equation of the motion of the particle can expressed as which are approximated by two (leap-frog) iteration formulas executed every half-time step.

Conversion into lattice units
To complete implementation of the numerical model, all the physical input quantities describing the flow and the particle have to be converted into their lattice counterpart. It is achieved by setting the conversion factors from the system of physical units (meter and second) into the system of lattice units (lattice unit and time step). These factors result from the relations between physical time and space elements, Δx and Δt, and lattice time and space elements which are set to equal to unit, Δx * = Δt * = 1.
If the physical length [L] = m is discretized by a lattice with N nodes then for the space element. The time element Δt should not be too small, otherwise it would result in a too small value of the lattice viscosity Consequently, and therefore the LBE method becomes unstable because of incapability to dissipate the energy due to very large velocity gradients. The time step Δ t can be determined by various ways depending on the specific model. Namely, it can be defined from ratio of physical and lattice sound velocity or by setting the relaxation parameter τ > 0.5 ensuring stability from or by assigning a maximal lattice velocity of the flow from relation with respect to a characteristic macroscopic velocity u char . Additionally, a requirement on incompressibility Ma = u * max /c s 1 must be satisfied for all cases. It is possible to choose any of the above mentioned conversion method in the simulation.

Results and observations
The developed simulation which models the particle motion and the velocity field of the fluid is initially expected to reproduce correct physical behavior to demonstrate its capacity to compete with traditional CFD methods.
Most of the input parameters are adjustable within a range of values, e.g., maximal inlet velocity u char , shear velocity u * , radius of both the moving particle and the bed particles, density of both the fluid and the particle or physical viscosity ν. Such a variability in the physical input parameters leads to good flexibility in modeling various physical conditions of the particle-fluid system.
For the resulting illustrations bellow it is supposed that the process is performed in an open channel of the length L = 1 m and height L/10 with bed particles of radii r b = L/100. The flow is driven by the logarithmic velocity profile with maximal value u(0, L/10) = 0.4 m.s −1 at the inlet where the fluid density ρ = 1000 kg.m −3 . The kinematic viscosity of water is ν = 10 −6 m 2 .s −1 , but for the stability reasons, see Sect. 5.2, fluid of viscosity ν = 10 −4 m 2 .s −1 is examined. Thus the Reynolds number reduces from Re = u char L/ν ≈ 10 5 (the particle Reynolds number Re p ≈ 10 4 ) to Re ≈ 10 3 (Re p ≈ 10 2 ).
The radius of the moving particle is assumed to be equal to the radii of the bed particles, i.e., r = r b = L/100. The saltating particle is released in the steady state flow (generated by the inlet logarithmic velocity profile flow) from x 0 = L/5 and y 0 = L/3 with zero translational and rotational initial velocity (v x0 , v y0 ) = (0, 0) and ω = 0. The process is examined with the sand-like particle of density  Figure 3. Trajectory of the particle versus length and height of the channel for the rotating and non-rotating-particle. The particle trajectories are depicted with respect to the length and height of the channel expressed in lattice space units. It can be seen that rotation around the particle center influence both length and height of jumps, and also their number within the examined domain. Specifically, the rotating particle undergoes less of larger jumps in the observed part of the channel.

Comparison of rotating and non-rotating particle
From comparison of results for rotating and non-rotating particle can be seen how this assumption affects behavior of the saltating particle. Thus, figures below serves as illustrations of the basic visual outputs of the simulation which are used to characterize the particle motion and the fluid flow.
In the first figure 3 two particle trajectories (x(t), y(t)) are compared corresponding to the rotating and non-rotating particle. The trajectories are depicted with respect to the length and height of the channel expressed in lattice space units. It can be seen that rotation around the particle center influences both length and height of jumps, and also number of jumps within the examined domain. Specifically, the rotating particle undergoes less of larger jumps in the observed part of the channel.
The trajectory shapes considered above implies that also translational velocity components (v x (t), v y (t)) differ for both cases. Specially, the rotating particle moves faster in both components of velocity vector v which also results in different time spent within the observed domain. Hence, the rotating particle spends in the process domain about one fifth less of lattice time steps with respect to the non-rotating one, cf. figure 4.
For illustration also the angular velocity with respect to the number of lattice time steps is depicted where the collision events can be detected as discontinuous changes in the quantity development. The angular velocity is described by a graph with few jump contributions -determined by change in angular velocity within the collision, cf. Eq. (5) -compensated by successive increase (in absolute value) following after the collision, cf. figure 5.
The velocity field (u x (t), u y (t)) of the flow is driven by the logarithmic velocity profile at the inlet on one side, and by the interaction of the flow with the moving particle on the other side.  , v y (t)) are compared for a rotating or non-rotating particle. The characteristic alternating pattern can be seen for the y-component as the particle undergoes collisions with the bed. The translational velocity components (v x (t), v y (t)) differ for both cases: The rotating particle moves faster in both components of velocity vector v which also results in different time spent within the observed domain. The rotating particle spends in the process domain about one fifth less of lattice time steps with respect to the non-rotating one.
Before the particle is released it behaves as a stationary solid circular obstacle with two asymmetric eddies (determined by the asymmetric initial position of the particle with respect to the height of the channel) separated from the surface of the particle. The eddies maintain themselves till the particle starts to move. Then the eddies start to circle around the particle and deform dramatically.
In the figures below, sequence of four detailed snapshots illustrates the velocity field around the saltating particle in free motion, before and after the collision for either rotating or non-rotating case. It can be seen that in the vicinity of the rotating particle the flow velocity is larger and fluctuates more intensely due to the greater total veloc- ity on the particle surface than in the non-rotating process, cf. figure 6.
Also the coarse-grained subdivision (discussed afterward in Sect. 5.2) of the lattice can be detected. Obviously, the proposed refinement of the lattice should introduce more precision into the simulation of the saltation motion and thus enable its significant improvement.

Stability issues
The very low viscosity of water results in high Reynolds number Re = uL/ν regime of the flow. This means that the relaxation parameter τ ≈ 0.5, cf. Eq. (15), and that the simulation becomes potentially unstable. The instability issues are cured in various ways.
Better stability or an increase of the relaxation parameter could be achieved by refining the lattice grid. However, the simulation operates currently on a rather coarse-grain lattice with respect to the physical dimensions of the experiment because of various reasons. First, the rather coarse grid allows immediate observations of potentially unstable regions in the process domain. Secondly, grid refinement yields significantly increase in demands on computational resources. This is will be achieved by using the CUDA GPU computing technology transformation which can be applied on the Matlab code which takes advantage of the Parallel Computing Toolbox. Due to the inherently parallel nature of the lattice Boltzmann it is expected that it will result in a satisfactory acceleration of the calculation. Consequently, it enables to refine the lattice grid significantly.
Except for the refinement of the grid, the stability properties of the simulation can be improved by modifications of the LBE model. It is worth to mention the best known LBE based methods for which stability reasons have been developed: regularized, entropic and multi-relaxation time LBM.
The regularized LBM is based on a replacement of the distribution functions by their regularized counterparts which guarantees symmetries resulting from the Chapman-Enskog expansion, e.g., [14,15]. The entropic LBM is based on the second thermodynamic law and uses the Boltzmann H-theorem [16,17]. The multi-relaxation time LBM, e.g., [18], replaces the LBGK collision model Eq. (7) with a multi-relaxation-time collision model which allows adjustment relaxation times independently for each of moment of distribution, cf. Eqs. (10).
One of the mentioned stabilization techniques, specifically entropic LBM, is employed for the region located in the surrounding of the moving particle where the largest flow velocity fluctuations occur. To avoid overly large velocity gradients corrections to distributions values f i of the ring nodes are made around the object by employing the entropic lattice Boltzmann model. Except for the above considered methods, another possible way to cure instabilities of this type could be to use adaptive mesh methods in neighboring nodes around the particle, e.g, [19].

Concluding remarks and prospective plans
This work presents a test of the capabilities of LBE based methods for simulating the saltation motion of the particle. Although the LBE method is subject to particle-fluid system studies, e.g., [20] it has not yet been used as a DNS strategy simulation tool for the saltation process.
Simulation is developed in the Matlab programming environment. Matlab represents a suitable initial programming tool as it enables to reflect the structure of the LBM, e.g., support for matrix algebra. Modular structure of the code results in enhanced flexibility and allows for additions of specific modules, e.g., for various boundary and initial conditions, various particle-bed collision models, incorporation of a particle-particle collision models etc.
Plans for the future are to develop the recent simulation tool to significant flexibility which would allow to test LBE based methods on all the types of particle motions (e.g., rolling, saltation, suspension) for cluster of particles in high Reynolds number regime flows. Further, results of the simulation should be comparable with outputs of the experiments carried out in the Institute.
Obviously, in such a case, enhancement and improvements of the current simulation will be necessary. Enhancement of the simulation consists in extension to the almost arbitrary surface of the bed (or the particle) which is possible due to the quite straightforward implementation even very complex boundaries, introducing many mutually colliding particles (computational complexity grows linearly with number of moving particles), ensuring stability for very large Reynolds numbers and finally a transition into the third dimensional process.
The first and second extension mainly represent the choice of corresponding models (e.g., particle-particle collision), and elaboration and modification of LBM (e.g., treatment of curved boundaries). The third requirement can be partly solved by the enrichment of LBM with methods mentioned in Sect. 5.2, eventually with their combinations. However, both stability requirements and transition to 3D mean significant increase in demands on computational resources. Extension of computational capacities will be achieved with use of the CUDA GPU computing technology which can take advantage of the parallel nature of LBM. 02020-p.7 Figure 6. Velocity field (u x (t), u y (t)) resulting from interaction of the particle and the fluid for the rotating (upper row) and non-rotatingparticle (lower roe). Sequence of four detailed snapshots illustrates the velocity field around the saltating particle in free motion, before and after the collision for either rotating or non-rotating case. It can be seen that in the vicinity of the rotating particle the flow velocity is larger and fluctuates more intensely due to the greater total velocity on the particle surface than in the non-rotating process.
Among the prospective benefits of the simulation could be possibility to simulate input parameters from wide ranges and examine initial and boundary conditions which may be difficult to prepare experimentally. Thus it could also serve as a guide or help for design of prospective experiments.