Onset of sediment transport in mono- and bidisperse beds under turbulent shear flow

In this paper, we present a study on the onset of sediment motion of particles under turbulent shear flow. We use a fully resolved coupled lattice Boltzmann–discrete element solver to study the behavior of monodisperse and bidisperse beds of spheres. For monodisperse beds, we find that the onset of sediment motion is in agreement with the predictions with models from the literature. For bidisperse beds, segregation phenomena at a short timescale are found that lead to increased entrainment of larger particles. This shows that, in addition to long-term segregation phenomena such as riverbed armoring, segregation also occurs at much shorter timescales.


Introduction
The study of sediment transport and its onset is important in many fields of science and engineering, ranging from coastal engineering to soil erosion. In all such cases, understanding the mechanisms that govern the entrainment of sediment particles is crucial. The most important early study on the onset of sedimentation was conducted by Shields [30]. His work consisted of analytic considerations about relevant parameters and of an exhaustive set of measurements. In order to identify criteria for sediment transport, Shields used the friction Reynolds number Re f = u f d p /ν and the dimensionless shear stress = τ b / g(ρ p − ρ f )d p , where τ b is the shear stress at the bed surface, u f = √ τ b ν the friction velocity, d p the particle diameter, ν the kinematic viscosity of the fluid, ρ p,f the particle and fluid densities, and g = 9.81 m/s 2 the acceleration due to gravity. By thorough experiments, he found a curve c (Re f ) above which particle motion can be expected. This curve is known as the Shields curve and is shown in Fig. 1. A huge body of work on sediment motion under shear flow followed, which is summarized in the review of Buffington and Montgomery [7] who performed a thorough analysis of over 600 data sets on incipient motion. Additionally, Dey and Papanicolau [12] summarized theoretical developments on the matter. Most of this work is focused on natural sediments with a relatively uniform size distribution and a high angle of repose. Dey [11] developed a semi-analytic model in which the latter is used as an additional parameter. This approach is able to accurately reproduce sedimentation onset for a wide range of materials. However, little is found on the threshold for sediments consisting of more than one size fraction. Wilcock [37] compiled results of experiments on bimodal sediments and constructed empirical models in which the critical shear stress for each particles size fraction is related to the degree of bimodality of the mixture. Other authors [36,38] adapted and refined these models, although no final consensus has been found.
Starting from the perspective of a single particle, the local bed morphology clearly plays a crucial role. The probability of a particle to be picked up by the flow significantly decreases when other particles are in close vicinity [1,10]. Dancey et al. [10] found that the value of c is highly sensitive to the planar packing density λ (area covered by particles over total area). Already a single near-by particle decreases, the chance of entrainment [1]. This underlines the necessity of detailed information on the flow field on the sub-particle scale to understand the relevant mechanisms.
Resolved numerical modeling of particle beds under shear flow is computationally demanding, especially in the case of turbulent flow. Therefore, not many studies in that field exist. Pan and Banerjee [25] and Shao et al. [29] investigated channel flow with buoyant particles, but put more emphasis on the turbulence structure in the fluid than on particle dynamics. Kidanemariam and Uhlmann performed simulations of particle beds under laminar [21] and turbulent [20] flow and focused on the time evolution of the beds. All articles mentioned above addressed monodisperse systems, and to the best of the authors' knowledge, no detailed simulations on polydisperse systems are found in the literature.
In this paper, we present fully resolved simulations of mono-and bidisperse beds of spheres under turbulent shear flow. After explanation of the numerical methods and a validation example, we analyze the behavior of the bed with respect to motion threshold and identify a fast mechanism for segregation that could be of relevance in situations where the flow velocity varies frequently. To conclude, we give an outlook how to integrate such small-scale simulations with high resolution into less detailed ones on larger scales via an LBDEM magnification lens.

Computational methods
In this section, the simulation methods employed in the present study are explained. All algorithms were implemented by creating a coupling framework between the lattice Boltzmann code Palabos [15] and the discrete element method code LIGGGHTS [22]. This allows us to take advantage of robust and well-tested software platforms. The coupling code, LBDEMcoupling 1 [27,28] was released under a free license.

The lattice Boltzmann method
The lattice Boltzmann method (LBM) [17] is based on kinetic theory. Technical details on the theoretical background and implementation issues are available in the literature [32,39], and most developments in LBM can be found in the reviews by Chen and Doleen [8] and Aidun and Clausen [3]. Instead of directly manipulating macroscopically observable quantities such as pressure or velocity, the fluid is described by a set of particle populations { f i } traveling in directions {c i }. These populations live on a regular grid, hence the name lattice Boltzmann. The macroscopic quantities can be obtained by The dynamics of the { f i } are governed by the LBGK (lattice BGK, using the BGK approximation [6]) equation where F i is a forcing term to be discussed below, and 1 https://github.com/ParticulateFlow/LBDEMcoupling-public.
is the equilibrium distribution for a fluid with velocity u and density ρ. Using a Chapman-Enskog expansion, the Navier-Stokes equation can be derived from Eq. 2 and 3. Details can be found for example in the book by Wolf-Gladrow [39]. The constants w i and c s depend on the discretization scheme chosen. τ is the relaxation time and represents the timescale at which the populations relax toward their equilibrium values. It is connected to the kinematic viscosity ν via In the present work, the D3Q19 discretization scheme was used. It allows particle populations to travel to the next and second-next neighbors on the grid with lattice vectors The velocity weights are given by and the speed of sound is c s = 1 √ 3 . A lattice Boltzmann simulation is performed in repeating two steps: First, a new set of populations { f * i } is computed using the right-hand side of Eq. 2 (collision step). Then, the populations are redistributed according to the velocity vectors c i : Each population f * i (x) is shifted to the position x + c i . This is given by the left-hand side of Eq. 2 (streaming step). The simulation then continues with a new set of off-equilibrium populations at each grid node.
To describe sub-grid turbulence, a Smagorinsky model [18,31] was used. The sub-grid turbulence is accounted for by an additional turbulent viscosity where |S| is the norm of the strain rate tensor, and C S = 0.1 . . . 0.2 is the Smagorinsky "constant." In the case of LBM, the filter width = 1. The turbulent viscosity is computed at each cell prior to the collision step, and the relaxation time τ is recomputed according to Eq. 4 using a total viscosity ν total = ν + ν t .
Body forces are applied as follows [16]: Given an acceleration G, the macroscopic velocity is computed as and a forcing term is added in Eq. 2.

The discrete element method
To model the behavior of solid particles, the discrete element method (DEM) was used. It was first applied to granular assemblies interacting via pairwise forces by Cundall and Strack [9]. The trajectory of a particle i at position x i with mass m i can be computed by integrating Newton's second law where F

(int)
i j is a pairwise force due to particle-particle interaction and F (ext) is some external force such as gravity or drag. Granular assemblies are usually modeled by rigid spheres, for which the contact detection is very simple: If the overlap in terms of the spheres' radii r i, j is positive, then they are in contact. The overlap given in Eq. 11 is also used to compute the force between the two spheres, which can be split up in a normal and a tangential component where v ⊥ i j and v i j are the normal and tangential relative velocities of the collision partners. The normal and tangential forces can be computed from material and collision parameters as follows [4,34]: with being the effective radius and mass, and the effective Young's modulus and shear modulus of the collision partners. ν i, j is Poisson's ratio for particles i and j, respectively. Further, can be computed from the coefficient of restitution α i j [4]. For a contact starting at t 0 , the tangential overlap is calculated as following Cundall and Strack [9]. The tangential force is limited to μ f F ⊥ to account for sliding friction between the particles, where μ f is the coefficient of sliding friction. Rolling friction was included using the directional constant torque model by Zhou et al. [2,40], in which an additional torque contribution proportional to the coefficient of rolling friction μ r is added to both particles. ω i j is the relative rotational velocity of the particles.

Coupling methodology
To couple the fluid and the solid phase, the immersed moving boundary method of Noble and Torczinsky [24] was used. The presence of a moving solid is modeled by blending between the BGK collision operator Ω BGK i in Eq. 2 and an additional collision operator Ω s i . The modified LBGK equation now reads with where u s is the velocity of the particle and the blending function B s is given by The solid fraction f s at each grid cell is used to blend between fluid collision and solid bounceback. The force and torque on a particle with center of mass at x 0 are then evaluated as

Corona boundary condition
To drive the flow, a force-based corona boundary condition was used. An acceleration is applied at each lattice site x in a region at the boundary of the simulation domain. This is actually a simple proportional controller with proportional gain a o u max where the corona base strength a 0 is a free parameter and u max = max u ext . While this acceleration is applied, the domain is still modeled as periodic in streamwise direction. This allows to impose a large-scale flow field while retaining small-scale turbulence. To achieve this, the value a 0 must be chosen adequately: If a 0 is too small, the flow will not be driven properly. If it is too large, small-scale structures will be damped out in the corona. Rudimentary correctness of this approach was checked by looking at the turbulent energy density spectrum in the bulk flow (not shown). Currently, we determine the correct value by a trial-and-error approach, but further studies will be focused on a more thorough way of finding a 0 . One reason for choosing such a boundary condition is that the present work should be understood within the bigger picture of an LBDEM magnification lens. This concept is further elaborated in the outlook (Sect. 6).

Validation
The following validation case has already been published in [27,28], but is repeated here for completeness. Feng et al. Fig. 2 Comparison of measurements by Ten Cate et al. [33] for a position and b velocity with our simulations. Background images with symbols indicating measurement data were reproduced from ref. 33 [13] showed that the combination of methods outlined above is suitable for simulations of turbulent flow with immersed particles. Therefore, this validation case is meant to show the correctness of the implementation. For this purpose, we used results from Ten Cate et al. [33]. They performed experiments on a sedimenting sphere in a tank of fluid with a setup designed to allow for easy simulation. The tank measured 100×100×160mm and the sphere diameter was d s = 15mm. The tank was filled with silicon oil of varying density and viscosity. Initially, the gap between the sphere and the tank bottom was 8d s . While the sphere was allowed to settle under gravity, its position and velocity were recorded using a PIV system. The Reynolds number of this system can be defined as where u ∞ is the theoretical sedimentation velocity in an infinite medium, and ν is the kinematic viscosity of the fluid. Experiments and our simulations were performed at Re = 1.5, 4.1, 11.6, 32.2. Simulations were conducted with N = 8 grid points along the sphere diameter, resulting in a grid spacing of x = 1.875mm. A rather low maximum LB velocity u lb,max = 0.002, and thus low Mach number, was chosen to avoid any compressibility effects. Figure 2a , b shows comparisons of position and velocity over time between simulation and experiment. We see very good agreement for all Reynolds numbers except for the lowest one (Re = 1.5), where the final stage is not depicted accurately. It is known that the LBM cannot properly compute the forces between moving objects when their distance becomes smaller than the grid spacing, so that lubrication corrections-the research topic of Ten Cate et al. [33]might be necessary. Since in our simulations, no lubrication forces were considered, the slowest approach was not modeled accurately. The acceleration and maximum velocity, however, are correct. We therefore conclude that our imple-

Setup
A box with edge length l = 2 cm as depicted in Fig. 3 was chosen as simulation domain. In streamwise (x) and lateral (y) directions, periodic boundary conditions were used. On the top, we imposed a free-slip boundary condition while on the bottom, a no-slip condition was used. Up to a height of h bed = 5mm, the domain was filled with a particle bed. To drive the flow, we used the force-based corona boundary condition as described in Sect. 2.4 to enforce a generic shear flow profile following Schlichting [26]. The corona acceleration (Eq. 25) was only imposed near the inlet, at 0 ≤ x ≤ 2mm (gray  The domain was resolved with N = 180 grid cells along the edge lengths, resulting in 9 cells per small particle diameter. This resolution corresponds to y + ≈ 1m and was found to give reasonable accuracy in the validation tests. Further, the resolution is comparable to those used in other studies [13,20]. Of course, a finer grid would have been desirable to fully resolve the boundary layer of the particles, but this was computationally prohibitively expensive. With respect to the particle diameter, the domain size was 20 3 d p . This is comparable to the boxes used in some other studies [21,25,29]. However, Kidanemariam and Uhlmann [20] used a domain with streamwise extent of ≈ 300d p and still reported influence of the box size. Therefore, we assume that our box was large enough to obtain qualitative findings, but cannot exclude box size influence on quantitative results.

Results and discussion
In the following section, the results of our simulations are presented and discussed. To that end, a few quantities need to be defined. The first question to be answered is when a particle is considered "moving." There are different definitions and approaches, such as the zero flux expansion by the United States Waterways Experiment Station (USWES) [35] and the classification by visual observation by Kramer [23]. Both methods suffer from a degree of arbitrariness. In the present work, a particle was considered "moving" when |v part | > 0.01 m/s. This empirical threshold was chosen to exclude most particles only wiggling due to fluctuations but to include those that actually move streamwise. Similar to the proposed methods described above, this definition is problem-dependent and arbitrary to a certain degree. Further, the fraction of moving particles was defined as where n m is the number of moving particles and n total is the total number of particles in the domain. The fraction of small particles can be defined as and an analogous definition is made for f m,large .

Monodisperse beds
In Fig. 4a, the fraction of moving particles over time is shown for the simulations with d p = 1mm. Once transport was established, the fraction did not vary much anymore. Only the smallest considered fluid velocity u top = 0.4 m/s could not sustain particle motion. Figure 5 shows the locations of the four cases in the Shields plot together with a case with d p = 2mm, u top = 1 m/s, relative to the Shields curve. While the small grains were set in motion for all except the lowest fluid velocity, the large particles could not be moved in any of the investigated cases.
Interestingly, the parameter setting corresponding to u top = 0.6 m/s lies below the Shields curve, but nonetheless According to Zhou et al. [41], spheres with rolling and sliding friction as used in the simulation typically show an angle of repose between 23 • (d s = 1mm) and 26 • (d s = 2mm). To take this parameter into account, we employed a semi-analytic model proposed by Dey [11] and added the corresponding curves for φ = 20 • and φ = 25 • in Fig. 5. The case with u top = 0.6 m/s lies between them, while the case with d p = 2mm is very close to the φ = 20 • curve. Dey's model predicts no movement for the large grains and movement for u top = 0.6 m/s, d s = 1mm for the small particles, just as found in the simulations.

Bidisperse bed
In contrast to the investigated monodisperse beds, the fraction of moving particles in the bidisperse bed varied strongly over time, especially for the two higher velocities. It can be seen in Fig. 4b that initially f m went up to values comparable to those of monodisperse beds. After about 1 s (u top = 1 m/s) and 0.7 s (u top = 0.8, 0.6 m/s), respectively, a notable drop in f m is visible. For the lowest velocity u top = 0.4 m/s, again no bed movement other than some initial reordering was found.
Closer inspection revealed that the most pronounced drop corresponding to u top = 1.0 m/s was accompanied by a change in particle motion. Figure 6 shows the fractions of moving small, large, and all particles over time. For u top = 0.8 m/s, the large particles showed almost no movement during the whole simulation, while f m,small , and with it f m , decreased over time from 0.04 to 0.02. For u top = 1.0 m/s, however, a different situation was found: After approximately 0.5 s, the large particles started moving. Once their motion was fully established, that of the of small particles drastically decreased until f m,small ≈ f m,large ≈ f m . At that point, particles of both sizes were in continuous motion.
The fact that there was motion of large particles far below their critical shear stress from the monodisperse case is consistent with the findings of Wilcock [37] who found that for beds of strongly bimodal sediments, the friction Reynolds number and the critical stress need to be computed using the average particle diameter D m = d p,small n small + d p,large n large n total .
(31) Table 2 shows the values computed with D m = 1.3 mm corresponding to the investigated mixture for u top = 0.8 m/s and u top = 1.0 m/s. Only for the higher of the two velocities, the shear stress exceeded the critical value c . According  to Wilcock [37], however, c is the critical shear stress for motion of all particle fractions; thus, no movement of small particles should be expected at u top = 0.8 m/s. We ascribe our deviating results to the low coefficient of rolling friction used in our simulations as discussed in the previous section.
To understand the influence of small particles on close-by large particles and their ability to be taken up by the flow, we analyzed the planar packing density of the top layer over time, which was defined as all particles with center at z > 0.0045 m that were not in motion. Its planar packing density is given by [10] and represents the fraction of the area covered by the projection of the particles onto the bed area. λ is not to be confused with the volumetric packing density. Figure 8 illustrates the concept by showing the top layer of the case with u top = 1 m/s at different times. Dancey et al. [10] studied the onset of motion in regular packings of glass beads for 0.026 ≤ λ ≤ 0.91 and found that the critical dimensionless stress c showed a strong dependence on λ (see Table  3). Its evolution in our simulations is displayed in Fig. 7.
As the small particles were set into motion, λ decreased with increasing u top . Of particular interest is the case with u top = 1 m/s. When the large particles started moving at approximately t = 2.5 s, the planar packing density had dropped from 0.42 to 0.26. According to Table 3, only 29.2% of the critical shear stress at λ = 0.45 were then necessary to initiate particle motion. This criterion was ful- filled ( / c = 53.9%; monodisperse for d p = 2 mm) once a sufficient amount of the small particles had been washed away. The situation was different at u top = 0.8 m/s. Here, λ = 0.29, and from the data of Dancey et al. [10], it can be estimated that 56% of the critical shear stress would have been needed to move the large particles, but / c = 48.91%. This agrees with the observation of some wiggling but no real transport. Thus, the drop of λ was the reason for the onset of the large particles' movement. This shows that not only the flow but also the local bed morphology has a huge impact on sediment motion.
Our findings are consistent with observed segregation phenomena in bidisperse beds known as riverbed armoring: a slow process in which the bed segregates in a top region with large particles and a region below with smaller particles [5,19]. Recently, it has been found that such granular creep is very similar to kinetic sieving processes in dry granular flows [14]. The segregation phenomena we observed The last entry in the table was found by linear interpolation yield a similar outcome, but occur at much shorter time scales.

Conclusions and outlook
We presented simulations of monodisperse and bidisperse beds under turbulent shear flow and compared the (non) appearance of mass transport with empirical correlations. Given some additional considerations, our results showed the correctness of our model: The study of monodisperse particles implies that material parameters, in particular the angle of repose, need to be taken into account in the correlations. Polydisperse systems may be regard as mixtures and phenomenological formulae then have to be evaluated in terms of the average particle diameter. If one wants to study segregation occurring on the bed surface, we recommend to investigate each size fraction separately and include the influence of other fractions reflecting in the local bed morphology in the calculations. The transition from the presented, small test cases to real-sized riverbeds is not straight forward. Fully resolved LBDEM simulations of processes spanning multiple time and length scales are prohibitively demanding for even the largest supercomputers. Therefore, we are aiming at magnification lens concepts: A full-scale simulation using continuum models is conducted. In this coarse simulation, a fully resolved, small-scale simulation is embedded and receives flow events, for example, via a corona boundary condition similar to the one presented in Sect. 2.4. The results of the small-scale simulation can then either be post-processed individually or fed back to trigger events in the large-scale model. To achieve this, appropriate boundary conditions need to be developed and validated.