EULER-LAGRANGE SIMULATIONS OF SEDIMENT TRANSPORT IN OSCILLATORY BOUNDARY LAYERS WITH BEDFORMS

Accurate prediction of sediment transport in the presence of bedforms such as sand ripples requires an advanced understanding of how dynamic sediment beds interact with turbulent oscillatory flows. In this paper we propose a new approach for simulating these interactions, based on a fixed grid multiphase Euler-Lagrange simulation, that fully couples dynamic bed evolution to the motion of a sub-grid scale Lagrangian sediment phase. The sediment phase is evolved by computing hydrodynamic and inter-particle forces and torques acting on individual particles, and is coupled to the fluid phase through the volume-filtered Navier-Stokes equations. We validate the approach for sediment transport applications using hindered settling velocity tests, and show very good agreement with the experimental data of Baldock et al. (2004). We then apply the approach to simulate sediment transport and ripple bed morphology in oscillatory flow conditions corresponding to the experimental studies of Van der Werf et al. (2007). During the simulation, particles near ripple surface are isolated from immobile ones below allowing the computation to devote resources only to particles that may be become mobilized. Although preliminary in nature, the simulation results demonstrate that that the model can correctly capture the near bed velocities, suspended sediment concentrations, and pick-up of sediment by key vortical structures.


INTRODUCTION
In a coastal environment, sediment particles on the seabed surface will start to move once the shear force exerted by fluid exceeds the critical values.Such motion often leads to complex fluid-particle interactions and results in various bed formations, such as vortex ripples, that are frequently found on the near-shore beach (Van Rijn et al., 1993).These bed features play an important role in affecting surface wave propagation, sediment entrainment from the bed surface, and mixing within the wave bottom boundary layer (WBBL) and higher above.Typical approaches to simulate the bedform dynamics and sediment suspension are based on the shear force concept which was initially derived for sphere particle sitting on a plane bed under steady current.However, both experimental and theoretical studies have revealed that within the wave boundary layer over a rippled bed, strong vortex shedding dominants the particle entrainment processes, which is very different from that in a plane bed conditions (e.g.Andersen and Faraci (2003), Li and O'Connor (2007), van der A et al. (2010)).Consequently, the considerable uncertainties in the predictions that rely on such an simple concept are often difficult to quantify.In addition, existing turbulence closures also fail to represent coexistence of the low Reynolds number flows in the ripple trough and the high turbulence region above the ripple crest (Barr et al., 2004).More importantly, the migration of bedforms require complex deforming meshes and transport closures to deal with the large bed slope or severe scour.Although noticeable progresses have been made in recent years (Marieu et al., 2008), these difficulties still remains as major challenge.
Increasingly, alternative multi-phase approaches are attracting more attention because the bed deformation can be dealt with naturally in such a model as part of the solutions.Different modelling techniques have been proposed in the literature, including the Euler-Euler approach or mixture theory approaches, which in essence treat the sediment as an extra fluid phase.Similar to the model in Giri and Shimizu (2006), Penko et al. (2013) used a mixture model to simulate ripple dynamics under oscillatory flows with noticeable success.The model results highlighted the 3D nature of the vortex which needs to be taken in account for the sediment suspension.However, models of this type often have limited inclusion of fluid-particle and interparticle interactions due to the assumption that the sediment behaves as a continuum.In high concentration regions, especially below the bed surface where the concentration is approaching packed bed limit, these models require many empirical approximations that potentially hinder their generic application to much wider range of conditions.Instead, in the present study, we propose a new approach, based on a fixed grid multiphase Euler-Lagrange (E-L) simulation, which fully couples dynamic bed evolution to the motion of a sub-grid scale Lagrangian sediment phase.The motion of sediment particles are therefore represented naturally and the fluid-particle interactions can be simulated on a physical sound basis.
It should be noted that the simulation of bedforms using our E-L approach presents unique challenges and opportunities.The largest challenge is probably the extremely large separation of space and time-scales.Collision durations between sand particles a few hundred microns in size last on the order of micro-seconds, and must be resolved by the particle time step, dt p .However, full scale sand ripples are 10s of centimetres in length and appear in flows with oscillatory periods of 5-10 seconds or more (e.g.Nielsen (1992)).However, only a small portion of the grains making up a sand ripple need to be simulated with such a small timestep.These are the grains near the surface of the ripple, which are subject to the erosion/mobilization processes driven by the flow.In this work, we present a novel approach to isolate and simulate only these grains as the ripple evolves, reducing the computational cost by roughly an order of magnitude for a full scale ripple simulation.

NUMERICAL MODEL
In this paper, we adopt a fixed grid, multiphase Euler-Lagrange simulation that fully couples dynamic bed evolution to the motion of a Lagrangian sediment phase.Below we outline the governing equations for both the sediment and fluid phase.

Sediment Motion
Motion of the sediment phase is computed by evaluating the sum of hydrodynamic and inter-particle forces and torques acting on individual Lagrangian particles.Particle positions, x p , velocities, u p , and angular velocities, ω p , are integrated in time according to the following equations of motion: Here, m p = π 6 d 3 p and i p = m p d 2 p 10 are the particle mass and moment of inertia respectively.The particles are assumed to be spherical, and their diameter, d p , is considered to be on the order of, or smaller than, the Eulerian grid spacing.In this work, we consider forces resulting from drag, lift, added mass, pressure, gravity, and collisions.Torques, and thus particle rotation, are considered to result from both inter-particle collisions and local rotation of the flow.Due to the unresolved nature of the fluid-particle interface, expressions for the hydrodynamic and inter-particle forces forces and torques acting on the particles employ closures developed from theory, experiments and fully resolved simulation (DNS).The expressions used to compute each force in the RHS sums of equations 2 and 3 are tabulated in Table 1.Where appropriate, any closure model, its origin, and user user defined parameters are noted.The gravity force, F g is simply the weight of the particle.When combined the pressure force, F p , which involves the gradient of the total fluid pressure, the particles respond to both buoyancy and local variations in dynamic pressure.The drag force, F d is computed based on the difference between particle velocity and the fluid velocity evaluated at the particle, u f |p .The influence of solid particle concentration on the effective drag coefficient of solid suspensions has been known for some time (Richardson and Zaki (1954)).Experiments and fully resolved simulations have allowed for parametrization of the effective drag coefficient, C d , as a function of Reynolds number, Re p = ρ p d p (u p − u f |p )/µ f , and solid fraction, θ p .In this work, the correlation of Beetstra et al. (2007) is chosen because it is valid for the entire range of solid fractions (0 < θ p < 0.6) and particle Reynolds numbers encountered in the WBBL.We consider only shear based contributions to the lift force, F l , and we use the closure of Saffman (1965) for the lift coefficient, C l .The expression for added mass force, F am , is based on an isolated rigid sphere in non-uniform flow (Crowe et al., 2011), with added mass coefficient, C am = 0.5.
In the high concentration sediment bed, particle motion is dominated by collisions and enduring contacts between particles.The collision force, F c , is based on a soft-sphere discrete element model.When two particles, denoted i and j and illustrated in Figure 1, come into contact, equal and opposite repulsive forces based on the spring-damper model of Cundall and Strack (1979) are generated.At any instant the  2007) 5.32 37.2 0 20 ≤ Re r < 50 6.44 32.2 0 50 ≤ Re r < 100 6.45 32.1 0 particles are separated by distance, x i − x j .The total relative velocity of the two spheres at the contact point can be written as, where n ij is the unit normal vector pointing from particle i to j.This can be decomposed into the normal and tangential components, Upon collision, a repulsive force is generated in the normal direction based on the overlap, δ i j = x i − x j − d e f f , and relative normal velocity, u n i j , of the two spheres.The spring constant, k c , force range, α = 0.075, coefficient of restitution, e c and coefficient of friction, ϑ are chosen on a case by case basis.These parameters set the collision duration, ), which must be resolved by the particle timestep, dt p .
Torques acting on the individual particles are computed using two simple models.Due to oblique collisions, inter-particle contacts will result in rotation.We model this using the simple Coulomb friction law, which as been demonstrated to provide reasonable accuracy and good speed for densely packed systems (Capecelatro and Desjardins, 2012).A hydrodynamic couple, T h is computed based on relative rate of particle and fluid rotation,

Fluid Motion
The sediment motion is coupled to the Eulerian fluid motion by solving the volume filtered Navier-Stokes equations (Anderson and Jackson, 1967) In this form, the conservation of mass and momentum equations are modified to account for the volume of fluid which is locally displaced by the sediment.For dynamic sediment beds, temporal and spatial gradients of the sediment fraction can be very high, and we expect that this attribute of the volume filtered equations may be important for correctly capturing bed morphology.In addition to the volume exclusion effects, Eqn.6 also contains an inter-phase momentum transfer term, f p→ f .This term includes the equal and opposite reaction from the particle surface onto the flow.It can be expressed as, Here, F x cv , x p is the filter kernel which transfers a property from the Lagrangian sediment phase to the fixed Eulerian grid.Note that body forces (gravity) and collision forces are excluded from this reaction.The fluid volume fraction is computed using the same filter operator to transfer the Lagrangian particle volume, V p = πd p 3 6 , to the fixed grid.
In this work, a Gaussian filter function is used, with standard deviation, σ = 2d p .In order to account for the increased viscosity of non-dilute regions of the flow, a model for the effective mixture viscosity is sometimes adopted in E-L simulations using the volume filtered equations (Patankar and Joseph, 2001;Capecelatro and Desjardins, 2012).An Eiler's equation (Eilers, 1941) model for the mixture viscosity of a particle suspension has been applied successfully within a mixture theory framework (Penko et al., 2009(Penko et al., , 2013)), and it is adopted for use in the present model.The viscosity of the fluid phase is, Here, µ 0 is the nominal fluid viscosity in the absence of sediment, and θ CP is the close packed sediment fraction, taken to be 0.615.Additional contributions to the effective viscosity due to unresolved turbulent stresses can be modelled using LES, which will be an area of future effort.In this work, the fluid motion is resolved to near particle scale, and it is assumed that such contributions will be small.The continuous fluid phase and Lagrangian sediment phase are evolved by solving equations 1-3, 6 and 6 using a new structured grid implementation of the algorithm presented by Shams et al. (2011).Extensive validation of the single phase solver can be found in Skitka (2013), andin Finn et al. (2014) for the Lagrangian sediment solver.In general the particles do not all settle at the same velocity because they induce an unsteady flow and form clusters of fast and slow moving particles.This is illustrated in Figure 2a for the case θ p = 0.3.Each particle is coloured by its normalized fall velocity v p /v t , where v t is the unhindered terminal velocity of a single particle (v t = 0.048ms −1 in this case).The normalized particle fall velocity (in a frame which the mean fluid velocity is zero) is averaged over all grains and plotted vs θ p alongside the experimental data in Figure 2b.We also plot the empirical relationship of Richardson and Zaki (1954) for these conditions (Re t = 16.8, n = 3.37) .Overall, the agreement with experiments is excellent over almost the entire range of solid fractions.At the highest solid fractions, θ p ≥ 0.55, the numerical simulations over-predict the experimental results somewhat.This is likely due to micro-structural phenomena, occurring at scales smaller than the particle size, which cannot be resolved in the present approach.

SIMULATION OF FULL SCALE RIPPLE DYNAMICS
The model is now applied to oscillatory flow over a single, full scale, mobile sand ripple, in conditions corresponding to the experiments of Van der Werf et al. (2007).Their dataset is unique in that it provides detailed velocity and concentration data for full-scale, well characterized ripples obtained in a closed loop oscillatory flow tunnel.Their two dimensional sand ripples were generated from an initially flat bed by a velocity skewed free-stream (see inset of Figure 3).In this paper, we consider case MR5b63 from the experimental database, where the equilibrium ripple wavelength and height were measured to be, λ = 41cm and η = 7.6cm respectively.The skewness of the flow resulted in asymmetric, forward leaning ripples, a net transport in the "offshore" (negative x) direction and a ripple migration "onshore" (positive x).Reproducing time dependent velocity and concentration fields in these conditions represents a significant challenge for any model, and we use this difficult test to show the capability of the coupled E-L approach.Figure 3 illustrates several aspects of the simulation setup.Our domain dimensions, shown in Figure 3a, are chosen to be L x = λ in the stream-wise direction, L y = 4η in the wall normal direction, and L z = 30δ x in the span-wise direction, where δ s is the stokes layer thickness based on the free stream velocity.The domain is periodic in the stream-wise and span-wise directions.A rough wall is imposed at the bottom wall (y=-9.6cm)by fixing any particles coming into contact with this boundary.At the top boundary (y=21cm), a zero stress condition is imposed.The sand used in the experiments was well sorted with d 10 , d 30 , d 50 , d 70 , d 90 equal to 0.25mm, 0.35mm, 0.44mm, 0.53mm, and 0.66mm respectively.In the simulations, only a single sand size of d p = d 50 = 0.44mm is used, although our future simulations will investigate the effects of the particle diameter distribution on transport characteristics.The initial ripple shape is created in a two step procedure.First, 19 million particles are released from random positions above the rough wall, and allowed to settle into a random close-pack configuration shown in Figure 3a.We then use the ripple profile generated by Van der Werf et al. (2008) as a periodic fit to the experimental measurements, to create the ripple shape shown in Figure 3b.After trimming away particles, roughly 12.5 million particles remain in the domain.The fluid domain is discretized with a uniform 768 × 576 × 64 grid, which corresponds to a cubic spacing just larger than the particle diameter, ∆ = 0.5mm = 1.13d p .Particle collision parameters are chosen to be k c = 1000, e = 0.5, ϑ = 0.1.
The flow is started from rest and a time dependent body force is applied to both the flow and particles (see Calantoni and Puleo (2006) for the importance of the latter) so that the free stream velocity, U(t) matches the second order Stokes free stream velocity, Here, ω = 2π/T is the angular frequency, T = 5s is the period, and u 1 and u 2 are 0.54 and 0.095 m/s respectively.The flow is run for a total of four periods which, in this preliminary study, seems at least qualitatively adequate for any initial transients to dissipate.

Isolation of mobile sediment
Simulating particle-particle and particle-fluid interactions for 12.5 million sand particles for 20 seconds represents a significant computational overhead.To minimize the amount of effort dedicated to computing the soft-sphere particle interactions, the bottleneck of the simulation, we employ a novel approach to isolate and simulate only the grains near the surface of the ripple.At regular intervals, roughly every 10ms, a mobility flag is computed for each particle, based on its depth, D p , below the ripple surface.The particle flag takes one of the three following values: • Mobile (D p < 9d p ): The particle is mobile, and will be advanced according to equations 1-3, using the time-step dt p .dt p is set to be roughly 0.1τ c • Fixed (9d p < D p < 18d p ): The particle position is fixed, but the particle interacts with mobile particles through collisions.It will also interact with the flow via the local volume fraction, θ p and the inter-phase coupling term, f 2w .Because collisions with mobile particles must be computed, a time-step of dt p is used.
• Dormant (18d p < D p ): The particle position is fixed, and the particle interacts with the flow only.
No collision computations are required.These forces acting on these particles are computed using the (larger) fluid time-step, dt f .In Figure 3b, each particle in the initial configuration is coloured corresponding to its mobility flag.The value of D p is computed by first integrating θ p in the y direction, taking advantage of the structured Cartesian grid, and then interpolating the result to the individual particles.The ripple surface is determined as the height where the θ p > 0.6.As the ripple evolves, so will D p and the mobility flag should be recomputed.This dynamic procedure is shown in Figure 4 for several instances in the 4th cycle, illustrating how the ripple shape deforms significantly because of the flow.The particles are sorted and stored in memory based on this flag, allowing for efficient access of the mobile particles as a subset of all particles in the domain.The overhead associated with this computing the mobility flag, and rearranging the particles in memory is negligible, and results in a roughly order of magnitude speed-up for this simulation compared to the case where all particles are considered to be mobile.

Experimental comparison
Because of the limited amount of flow periods simulated, ensemble averaging of the velocity and concentration fields is not possible.Instead, we take advantage of the two-dimensional nature of the ripple, and perform span-wise averaging of the results over the 4th cycle.By averaging over the 64 control volumes in the z direction, span-wise averaging produces the two dimensional fields u(x, y, t) and θ p (x, y, t) that we will interrogate below and compare with experimental data.In some instances there is still a significant amount of noise in the span-wise averaged variables, and additional simulations will be performed in the near future to facilitate ensemble averaging over several periods.
A comparison of the near bed velocity at 8 points along the ripple with the experimental PIV measurements is made in Figure 5.The probe locations were chosen to match the experimental measurement Turning our attention to suspended sediment concentrations, we make two types of comparisons with the experimental data.In Figure 6a, we compare the cycle average suspended sediment concentration as a function of height above the ripple crest with the experimentally measured (ABS) concentration.Within 1.5η of the crest there is excellent agreement between the measurements and simulation.Above this level, a strange departure of the two curves is observed.We believe that this is due to contamination by the upper (zero stress) boundary at y/η ≈ 3.This boundary dampens any turbulence and suppresses vertical transport of sediment higher up into the water column.The result is the observed bulge in the concentration curve at y/η ≈ 2. Future simulations will extend the simulation by several ripple heights to hopefully remove this effect.In Figure 6b we compare the time-series of suspended sediment concentration at three points above the ripple crest with the ensemble averaged ABS data.It is difficult to make precise comparisons because of the lack of temporal resolution and noise in the simulation data.Maximum and minimum concentration magnitudes predicted by the simulation do seem to be in line with the experimental measurements at each location.Future simulations with ensemble averaging over several cycles should facilitate a better comparison of this result.
Finally, we examine the interaction between the vorticity and suspended sediment during a single period.Despite a two dimensional ripple geometry in this case, the turbulent flow near the bed does have a three dimensional structure.We leave a three-dimensional analysis of this flow for later work, and instead focus on how the E-L results can be used to investigate the interactions of the bed with the large scale 2D flow features with large span-wise (z) vorticity component.Figure 7 shows contours of the z vorticity component in the lower part of the boundary layer, along with a contour (black line) of high suspended sediment concentration (13 g/L).The letters in the individual frames correspond to to the phase indicated by Figure 4.The flow accelerates to the right during the first two frames (Figure 7a-b), causing significant shear over the ripple crest.The flow decelerates and reverses in the trough ahead of the free stream creating a large CW vortex in the lee side trough (Figure 7c).The stress on the ripple due to the creation of this lee vortex is large enough to significantly deform the ripple (Figure 7c-d).As the free stream flow reverses, the lee vortex is ejected and caries with it a large amount of sediment which was picked up from the lee side slope (Figure 7d-e).The flow accelerates less quickly to the left and creates a similar but weaker vortex in the stoss side trough (Figure 7f-g).This vortex is less effective at entraining and ejecting sediment from the stoss side trough, and as a result there is a net transport to the left (offshore), even while the ripple migrates to the right (onshore).

CONCLUSION
In this paper, we have demonstrated the capabilities of a coupled Euler-Lagrange model for simulating the dynamics of full scale sand ripples in oscillatory flow.In this approach, the sediment phase is evolved by computing hydrodynamic and inter-particle forces on each individual grain, while the fluid phase is governed by the volume-filtered Navier-Stokes equations.The model is first validated by comparing the hindered settling velocity of a suspension of glass beads with the experimental data of Baldock et al. (2004), showing very good agreement over a wide range of solid fractions.It is then extended to the case of oscillatory flow over a full scale sand ripple corresponding to the experiments of Van der Werf et al. (2007).The computational demand for involved in such a simulation is extremely high, mostly due to the separation of collisional and hydrodynamic space and time-scales.To make the simulations tractable, mobile particles near the surface of the ripple are isolated from the remaining particles, and only these particles are used in the expensive inter-particle collision computations.Despite the preliminary nature of the results, they demonstrate that the E-L model is a very capable tool for studying dynamics of sediment-turbulence interaction in the WBBL.Our future goals involve incorporating poly-disperse sediments and a LES turbulence closure in the model, and obtaining longer simulation times so that ensemble averaging of the results can be performed.
with the torque coefficients suggested by Pan et al. (2001), valid for rotational Reynolds numbers, Re R = d 2 p |ω rel | 4ν f , up to 100.
VALIDATION: HINDERED SETTLING OF A PARTICLE SUSPENSIONIn the WBBL, sediment concentrations range from very dilute to the random close packed limit.A successful simulation of sediment dynamics in this environment must faithfully capture the sediment-fluid interactions over the entire range of solid fractions.In order to validate the behaviour of the coupled Euler-Lagrange model for WBBL flows, we consider the hindered settling of solid particles in an unbounded domain.A triply periodic domain with L x = L y = L z = 64d p is chosen for the simulations and discretized with 64 3 control volumes with spacing ∆ = d p .Identical, non-overlapping particles are seeded at random in the box to obtain target solid fractions in the range 0.05 ≤ θ p ≤ 0.6.The fluid and particles have the properties corresponding to the fluidization experiments ofBaldock et al. (2004) using uniform glass beads in water (d p = 0.35mm, ρ p = 2500kg/m 3 ρ f = 1000kg/m 3 , µ 0 = 0.001kg(ms) −1 ).Particle collision properties are chosen to be k c = 1000, e = 0.65, ϑ = 0.15.The particles are initially at rest and settle under the action of gravity (g y = −9.81ms−2 ), while a uniform body force is applied to the fluid in the positive y direction to balance the weight of the particles.The simulations last for a total time, t = 30 ρ p d 2 p 18µ , which was determined to be sufficiently long for the (average) settling velocity to achieve a stationary value.

Figure 2 :
Figure 2: Hindered settling of a particle suspension under the action of gravity.(a) shows particles coloured by instantaneous fall velocity for mean solid fraction, θ p = 0.3.(b) shows fall velocity normalized by the single particle terminal velocity, v t , as a function of θ p .Comparisons are made with the Richardson and Zaki (1954) correlation, and the experimental measurements of Baldock et al. (2004).

Figure 3 :
Figure 3: Configuration of the dynamic ripple simulation

Figure 4 :
Figure 4: Dynamic isolation of mobile particles near ripple surface.

Figure 7 :
Figure 7: 2D, span-wise averaged contours of z vorticity component (red = ccw rotation, blue = cw rotation).The black line shows the region of high suspended sediment concentration (13 g/L).

Table 1 :
Expressions, parameters, and closure relations governing Lagrangian particle motion