A TURBULENCE-RESOLVING EULERIAN TWO-PHASE MODEL FOR SEDIMENT TRANSPORT

A turbulence-resolving two-phase Eulerian model for sediment transport is developed by extending the new solver twoPhaseEulerSedFoam. The development and validation of a turbulence-averaged multi-dimensional two-phase model for sediment transport, twoPhaseEulerSedFoam version 1.0, was recently completed and the model was disseminated to the research community via the Community Surface Dynamics Modeling System (CSDMS) model repository. In order to resolve flow turbulence and turbulence-sediment interaction, large eddy simulation (LES) with subgrid turbulence closure is further implemented and simulation is carried out in three-dimension. Closures on particle stresses are based on kinetic theory of granular flow for binary collision and phenomenological closure for stresses of enduring contact. One of the main concerns in a turbulence-resolving simulation is regarding the domain size and finest grid resolution so that the flow domain is sufficiently large to contain the largest eddy and meanwhile, a significant amount of the turbulence energy is resolved appropriately. This paper discusses these crucial issues for typical sheet flow condition in oscillating water tunnels. Preliminary results also show that the model is able to capture the sediment burst events during the flow reversal.


INTRODUCTION
Multi-dimensional sediment transport occurs in many coastal applications such as scour (Amoudry and Liu, 2009), momentary bed failure (Madsen, 1974;Sleath, 1999;Foster et al., 2006) and bedforms (e.g.Penko et al., 2010).In the past decade, several Eulerian two-phase numerical models have been developed (Dong andZhang, 1999, 2002;Hsu et al., 2004;Li et al., 2008;Amoudry and Liu, 2009;Bakhtyar et al., 2010).These two-phase models can resolve the full profiles of sediment transport without the need to divide the transport into bedload and suspended load components.Hence, they can be used to better understand concentrated region of transport and provide insight on sediment pickup and deposition driven by unsteady flows.However, most of the existing two-phase models are based on the turbulence-averaged approach and many of them are further simplified into one-dimensional-vertical (1DV) formulation.Hence, they are restricted to turbulence-averaged sheet flow applications.In the 1DV model of Kranenburg et al. (2014), they showed that model results are sensitive to the parameterization of turbulence-sediment interaction, which is highly empirical.Several laboratory measurements of sheet flow in oscillatory water tunnel also reveal that burst events of sediment concentration during the flow reversal are commonly observed for fine sand (O'Donoghue and Wright, 2004;Dohmen-Janssen, 1999).Because such burst events occur near flow reversal where flow may not be fully turbulent (Jensen et al., 1989), a turbulence-averaged model may not be able to capture such process due to several assumptions adopted in typical turbulence-averaged closure models.Hence, to improve our existing sediment transport modeling capability, a turbulence-resolving two-phase model is needed.
Recently, a multi-dimensional two-phase Eulerian model is developed, and disseminated to the research community through the Community Surface Dynamics Modeling System (CSDMS) model repository.This new model (twoPhaseEulerSedFoam Version 1.0) has the capability of capturing inhomogeneous flow features, and bed instabilities under intense oscillatory flow (Cheng and Hsu, 2014).Although the numerical model is developed in three spatial dimensions using the open-source CFD library of solvers, OpenFOAM (Rusche, 2002), the model is only validated for turbulence-averaged formulation with a kclosure for sheet flow in steady channel flow (Sumer et al., 1996) and oscillatory flow (O' Donoghue and Wright, 2004).Hence, we are motivated in this study to develop a turbulence-resolving two-phase Eulerian model for sediment transport by extending twoPhaseEulerSedFoam.
Large eddy simulation (LES) has been applied successfully to various multiphase flow problems, such as bubbly turbulent flow (Deen et al., 2001;Lakehal et al., 2002).In LES modeling strategy, only large eddies and their subsequent cascade with length scales greater than the grid size are resolved directly, while sufficiently small eddies are modeled using the resolved fields.The principle behind LES is that large eddies carry most of the turbulent energy, and it is justified to assume that their interactions with dispersed phases are the dominant mechanism driving sediment transport.The small motions, including the interaction of dispersed particles with surrounding turbulence, have minor influence on the mean flow and transport, thus they can be modeled more easily.As a result of the scale separation, the computational cost is significantly reduced compared with direct numerical simulation, in which all turbulence scales have to be resolved.On the other hand, LES approach is expected to be more accurate than the turbulence-averaged approach because the later parameterize all the scales of turbulence.
The aforementioned LES modeling strategy can be considered as the guidelines, however, the numerical implementation is not trivial.In this study, we present our first step in developing and carrying out LES modeling for two-phase sediment transport.The model formulation is first discussed and followed by a brief summary of the numerical scheme.Using typical flow condition for sheet flow sediment transport in oscillating water tunnel, we critically examine the criteria for (minimum) domain size and (maximum) mesh size in three spatial dimensions in order to carry out an appropriate turbulence-resolving two-phase sediment transport simulation that follows the LES modeling guidelines.

MODEL FORMULATIONS
The two-phase Eulerian formulations can be obtained by averaging over carrier fluid and dispersed particles (Drew, 1983).As shown in Fig. 1, with the four-way coupled two-phase flow formulation along with appropriate closure models, the resulting model can resolve the full dynamics of sediment transport from the (porous) immobile bed, to the highly concentrated regions of transport dominated by enduring contact forces, to less dense region dominated by particle collision and turbulent suspension, and to the dilute suspended load region driven solely by flow turbulence.Particularly, the concentrated regions of sediment transport can be resolved by including closures of particle stresses and fluid-particle interactions in the governing equations.Hence, using the two-phase modeling approach for sediment transport, the resulting model does not require bedload/suspended load assumptions, which are commonly adopted by the single-phase flow approach.More details of this modeling approach are discussed next.In the present multi-phase Eulerian model, the continuum assumption is applied to both fluid phase and sediment phase.The LES methodology is based on the separation of large scale and small scale in a turbulent flow.To achieve this separation of scales, low-pass filtering operation is applied to the Eulerian two-phase equations for a fluid-particle system.Due to the filtering operation, the velocity field is decomposed into resolved and unresolved (subgrid scale) parts: where u f i , u s i are the resolved velocities, ũ f i , ũs i are the instantaneous velocities, and u f i , u s i are the unresolved velocities.The 3D LES formulation for the carrier fluid phase is well established (Deen et al., 2001;Yuu et al., 2001;Zhou et al., 2004).For the dispersed sediment phase, a similar procedure can be applied.However, for simplicity, the subgrid motion of sediment phase is established only by applying kinetic theory of granular flows to model particle collisions and frictional stress in the enduring contact region.In other words, the effect of subgrid turbulent motion on particle transport (or conversely the effect of particle on sub-grid turbulence) is neglected at this point.Assuming that there is no mass transfer between the two phases, the mass conservation equations for fluid phase and sediment phase are written as: where φ is sediment volumetric concentration, u f i , u s i are ith component of fluid and sediment velocities, respectively, and i = 1, 2, 3 represents streamwise, spanwise and vertical components.
The momentum equations for fluid phase and sediment phase are written as: where ρ f , ρ s are fluid density and sediment density, respectively, g i is the gravitational acceleration, p f is the fluid pressure and τ f i j is the fluid shear stress, which includes fluid viscous stress and stresses associated with unresolved turbulence.Particle pressure p s and particle stress τ s i j are calculated from kinetic theory of granular flow and phenomenological closure of contact stresses.M f s i and M s f i represent the inter-phase momentum transfer between fluid phase and particle phase, and M f s i = −M s f i .Closures of the momentum transfer term and the stress terms for sediment transport are discussed next.

Inter-phase Momentum Transfer
In this framework of Eulerian two-phase flow formulation, both fluid and sediment phases are considered as continuum, and they are immiscible.The momentum of these two phases are coupled through the inter-phase momentum transfer terms.In general, the interaction between fluid phase and particle phase includes the drag force, the added mass force, the Basset force, the lift force (Maxey and Riley, 1983) and the effect of turbulence fluctuations on the effective momentum transfer.Typically, the drag force dominates in many sediment transport applications, and hence for simplicity we neglect the other terms such as lift force, added mass force and Basset force.The momentum exchange has the following form: where the last term on the right-hand-side (RHS) of Eq. ( 7) is the inter-phase pressure correction term, and the first term on the RHS of Eq. ( 7) represents drag force due to relative velocity between fluid and particle phases.Though in grain scale, processes such as flow around particle also contributes to the inter-phase momentum exchange, its complexity is beyond the scope of this study.
For the closure of β, we adopt that suggested by Ding and Gidaspow (1990), who combined Ergun (1952) for dense sediment concentration (φ ≥ 0.2) and Wen and Yu (1966) for lower sediment concentration (φ < 0.2): where d is the grain diameter, and C d is expressed as: in which, Re p = (1 − φ)|u f − u s |d/ν f is the particle Reynolds number, and ν f is the fluid molecular viscosity.

Subgrid Turbulence Closure
The filtering of nonlinear convection term on the left-hand-side of Eq. ( 5) leads to an additional stresslike term that requires closure.This subgrid stress represents the effect of the unresolved scales on the resolved scales.In the present fluid-particle system, the unresolved subgrid stress should also include additional effects due to fluid-particle interaction at grain scale but they are neglected in this paper.Thus, the total fluid stress is written as: where ν f sgs is the subgrid eddy viscosity, and it is modeled by the standard Smagorinsky model: where C s = 0.167 is the model constant, and S f is the deviatoric part of the resolved fluid strain rate tensor: and the filter length ∆ is defined as: where ∆ x,y,z is the grid size in x, y, and z directions, respectively.

Closures on Particle Stresses
To resolve the full dynamics from dilute suspension region to concentrated regions of moderate to high sediment concentration, the closure of particle stresses adopted in this study includes two parts.For moderate sediment concentration, it is assumed that binary collisions dominate inter-granular interactions and a closure based on the kinetic theory is adopted.For large sediment concentration, binary collision eventually become invalid and inter-granular interaction is dominated by enduring contact/frictional forces among particles.Hence, the closures of particle pressure and particle stress both consist of a collisionalkinetic component and a quasi-static component (Johnson and Jackson, 1987;Hsu et al., 2004): The collisional component is first discussed.In the kinetic theory, particle stress and particle pressure are quantified by particle velocity fluctuations due to binary collisions.To quantify the strength of particle velocity fluctuation, the concept of granular temperature Θ is introduced (Jenkins and Savage, 1983) for dry granular flow consists of smooth, slightly inelastic, spherical particles.In the present two-phase flow condition, we adopted the balance equation for granular temperature suggested by Ding and Gidaspow (1990): where κ is the conductivity of granular temperature, γ s is the energy dissipation rate due to inelastic collision and J int is the production (or dissipation) due to the interaction with the carrier fluid phase, the closures are summarized in Table 1.Collisional normal stress (Ding and Gidaspow, 1990): Fluid-particle interaction: Radial distribution function (Carnahan and Starling, 1969): Dissipation rate due to inelastic collision (Ding and Gidaspow, 1990): Collisional shear stress (Gidaspow, 1994): Frictional normal stress (Johnson and Jackson, 1987): (Gidaspow, 1994): Frictional tangential stress: Particle bulk viscosity (Gidaspow, 1994): Frictional viscosity (Srivastava and Sundaresan, 2003): Granular temperature conductivity (Gidaspow, 1994): Sediment shear rate: When the volumetric concentration of particles becomes close to random loose packing (φ ≈ 0.57), particles are constantly in contact with each other, and particulate energy can be dissipated by friction between sliding particles (Tardos, 1997).Thus, when the sediment concentration exceeds random loose packing concentration φ f , frictional stress model need to be adopted.In the present model, we adopt the simple model of Johnson and Jackson (1987) for particle pressure, and the tangential stress is calculated by the model of Srivastava and Sundaresan (2003), which is able to capture the transition of solid-like feature to fluid-like feature of the sediment bed.More details are summarized in Table 1.
In Table 1, e = 0.8 is the restitution coefficient during the collision, θ f is the angle of repose and is taken to be θ f = 28 • for sand, the random loose-packing concentration is set to be φ f = 0.57, and the random close-packing limit is set to be φ max = 0.635.In sediment transport, the quasi-static component of particle pressure and particle stress play a definite role to ensure the existence of an immobile sediment bed and a low mobility layer of enduring contact (Hsu et al., 2004).Hence, the empirical coefficients presented here are calibrated to ensure that a stable sediment bed can be established below the mobile transport region.The following values are adopted in the present LES study: F = 0.05, m = 3, n = 5.

MODEL RESULTS
The present LES model is developed from the turbulence-averaged model, twoPhaseEulerSedFoam Version 1.0 (Cheng and Hsu, 2014).The turbulence-averaged model has been calibrated/validated with laboratory flume data of Sumer et al. (1996) for sand transport in steady channel flow under sheet flow condition and medium sand transport in oscillatory flow of O' Donoghue and Wright (2004).The model validation can be found in Cheng and Hsu (2014).By modifying the fluid turbulence model to standard Smagorinsky model, and ignoring the effect of fluid turbulence on sediment granular temperature, the modified model is directly applied to investigate the necessary (minimum) domain size and (maximum) grid size in order to carry out an appropriate turbulence-resolving simulation for typical sheet flow condition in oscillating water tunnel.We like to point out that for a turbulence-resolving approach, this verification procedure is necessary before a rigorous calibration of the model coefficients (model validation) can be further carried out.
O'Donoghue and Wright ( 2004) measured sediment concentration in both the bedload and suspended load regions under oscillatory sheet flow for a range of flow condition and three different grain diameters.The oscillatory flow is driven by a mean streamwise pressure gradient.Outside the wave boundary layer, the shear stress vanishes, and the momentum equation reduces to: where p f is the mean pressure, U 0 is the free stream velocity.Here, we consider a representative case of sinusoidal flow and U 0 = U m sin(ωt), where U m is the free-stream velocity magnitude, ω = 2π/T is the wave frequency, and T is the wave period.Relevant flow parameters are summarized in table 2.
According to the flow condition shown in Table 2, the Stokes boundary layer thickness is δ = 2ν f /ω ≈ 1.26 × 10 −3 m, and the Stokes Reynolds number is 1890.After conducting several numerical experiments, the domain size in streamwise and spanwise is chosen to be 160δ and 80δ, respectively.In vertical direction, the initial bed depth is set to be 40δ, and the water depth is 100δ.This domain size is discretized into 128 × 64 × 180 grids in x, y and z directions, respectively.Uniform grids are used in streamwise and spanwise directions with ∆ x = ∆ y = 1.56 × 10 −3 m, and non-uniform grids are used in vertical direction with minimum grid size ∆ zmin = 4.5 × 10 −4 m near the mobile bed layer (z = 0.03 ∼ 0.08 m), and the grid size increases gradually away from this layer.Periodic boundary condition is imposed in both streamwise and spanwise directions, no-slip wall boundary is used on the bottom boundary, and symmetric boundary is imposed on the top boundary.In this preliminary simulation, fine sand of grain diameter 0.15 mm is specified for sediment phase.The basic concept of LES requires that the grid resolution is fine enough to resolve the inertial subrange, and the domain size is large enough to resolve the largest eddies.This domain size requirement is first verified by autocorrelation of streamwise, spanwise and vertical velocity components in streamwise and spanwise directions (see Fig. 2).The flow field is subtracted from the plane in the homogeneous directions (normal to the vertical direction), and the vertical location evaluated in the figure is at z = 0.055 m.It is confirmed that the two-point correlation of each velocity component vanishes when their distance is larger than 0.4L x in x direction or 0.5L y in y direction, i.e., about half of the domain size.Thus, the domain size implemented here is large enough to capture the most energetic eddies.Meanwhile, to evaluate grid resolution, the energy spectrum in wave number k x (streamwise) and k y (spanwise) fields are shown in Fig. 3.It is evident that the slope of −5/3 (the slope of inertia subrange) is resolved for two (one) order of magnitude of energy cascade in the streamwise (spanwise) direction, which indicates that the grid resolution is sufficiently fine to appropriately resolve significant turbulence energy.
It has been reported that short duration of concentration peaks occur during flow reversal above the pick-up layer, and these peaks are much more pronounced for fine sand (Dohmen-Janssen, 1999; O'Donoghue and Wright, 2004).Since a grain diameter of 0.15 mm is used in the present simulation, simulation results are examined between flow peak and flow reversal.In Fig. 4, the iso-surfaces of dilute concentration at φ = 0.001 (green) and high concentration at φ = 0.3 (blue) are shown during flow peak and flow reversal, respectively.It can be observed that during flow peak, the averaged distance between iso-surface of φ = 0.001 and iso-surface of φ = 0.3 is larger, suggesting a larger mobile bed layer thickness.The spatial variability due to resolved turbulence is present but rather mild.On the other hand, although the averaged distance between iso-surface of φ = 0.001 and iso-surface of φ = 0.3 is smaller during flow reversal, we observe large spatial variability, especially for dilute concentration.This indicates that for fine sand, local upward ejection of fine sand can be observed during flow reversal, and this feature is absent during flow peak.This upward sediment ejection is similar to the sediment burst events observed during flow reversal in the literature.Sediment bursts observed in the present simulation is further illustrated by showing the side view of sediment concentration (Fig. 5) at the center plane (in the spanwise direction).To obtain a better contrast, the contour level is selected to range from 0 to 0.015.and the red and green finger-like structures in the side view is associated with the upward ejection discussed in Fig. 4.Moreover, we can see that sediment concentration in the core of the bursts can exceed about 1%.

CONCLUSION
In this paper, we report the progress of the development of a turbulence-resolving two-phase Eulerian model for sediment transport.The turbulence-resolving model has been developed based on the recently published multi-dimensional turbulence-averaged open-source code, twoPhaseEulerSedFoam (Version 1.0).The mathematical formulation is similar to twoPhaseEulerSedFoam, except that the fluid turbulence model is replaced with the standard Smagorinsky model, and the kinetic theory for granular flow is adopted without consideration of effect of fluid turbulence.The more crucial task is the three-dimensional implementation, including domain size and grid size in order to resolve turbulence and the chaotic nature of the nonlinear system.As a first step, the LES model is directly applied to oscillatory flow condition of O' Donoghue and Wright (2004).Appropriate domain size and minimum grid size for the given flow condition are demonstrated.Moreover, the preliminary results also suggest that the turbulence-resolving model is able to capture the sediment burst event for fine sand during flow reversal.
However, future work is required to complete the development of turbulence-resolving simulation tool.For example, the subgrid stress closure for flow turbulence in the present model is based on standard Smagorinsky model.However, several previous work has shown that this closure model is too simple for boundary layer flow and it is relatively dissipative for transient turbulent flow.Meanwhile, the effect of fluid turbulence on the sediment velocity fluctuations and subgrid turbulent suspension are simply neglected in the present model.More comprehensive model validations with available laboratory data will be carried out and these neglected terms may be recovered and calibrated in the near future.

Figure 1 :
Figure 1: Illustration of model domain

Figure 3 :
Figure 3: Energy spectrum in k x (panel (a)) and k y (panel (b)) fields

Figure 5 :
Figure 5: Side (xz-plane) view of sediment concentration contour at the middle plane in the spanwise direction.

Table 1 :
Closures for particle stresses