Application of a Eulerian two-phase flow model to scour processes

In this paper, the application of a two-phase flow model to scour processes is presented. The model is first calibrated against experimental data of unidirectional sheet-flow (one-dimensional configuration). The model is then applied to multi-dimensional configurations for the scour under a submarine pipeline and around a vertical pile. The results show that quantitative results can be obtained at the upstream sides of structures, the lee-wake erosion driven by the vortex shedding deserves further research.


Introduction
Over the last decade, significant research efforts have been devoted to modeling sediment transport using two-phase flow models either with Eulerian-Lagrangian or Eulerian-Eulerian approaches.Contrary to classical sediment transport models relying on empirical formula for the bed-load and pickup flux for suspended-load [1], these new models are based on physical grounds allowing for in-depth investigation of sediment transport mechanisms at the intermediate scale, i.e. ranging between the particle scale and the mean flow length scales.Recently, we have developed a community multi-dimensional Eulerian two-phase flow sediment transport model, sedFOAM, under the platform OpenFOAM (https://github.com/sedfoam/sedfoam)[2,3].The aim of this contribution is to present some recent developments and the application of the model to scour around structures.
As a first step the model is tested against experimental data for unidirectional sheet flows using one-dimensional vertical simulations.In a second part, the application of the Eulerian two-phase model to the scour phenomenon below a pipeline (2D configuration) and around a cylindrical pile (3D configuration) are presented.For these latter applications, the k-omega turbulence model from Wilcox (2008) [4] has been implemented and adapted for two-phase flow models to accurately predict the flow dynamics in presence of adverse pressure gradient.

Two-phase flow model
The two-phase flow model presented in this paper is the model introduced by Chauchat et al. (2017) [3].In the two-phase flow formulation, both solid and fluid phases are described by Eulerian equations.The fluid and solid mass conservation equations reads: where φ is the sediment phase concentration, u s i and u f i are the sediment and fluid phase velocities respectively.The momentum equations for the solid and fluid phases reads: where ρ s and ρ f are the solid and fluid density, g i is the acceleration of gravity, p is the fluid pressure, ps is the solid phase normal stress, τ s i j and τ f i j are the fluid and solid phase shear stresses.The fluid phase shear stress can be written as: with k the turbulent kinetic energy and ν E f f the effective velocity defined by ν E f f = ν f t + ν mix with ν f t the turbulent viscosity calculated by a turbulence closure model and ν mix the mixture viscosity following the model proposed by Boyer et al. (2011) [5]: where φ max = 0.635 is the maximum value for the solid phase concentration.
The last two terms of the right hand side of both momentum equations represent the drag force coupling the two phases with σ c the Schmidt number and K the drag parameter modeled according to: where d is the particles diameter, h Exp is the hindrance exponent controlling the drag increase with increasing solid concentration and C d the drag coefficient calculated by the empirical formula given by Schiller and Naumann (1933) [6]: where Re p is the particulate Reynolds number defined by : Re p = (1 − φ) u f − u s d/ν f with ν f the fluid kinematic viscosity.
The solid phase stress tensor is decomposed into the normal stress components and the off diagonal components namely the particle pressure ps and particle shear stress τ s i j .The particle pressure is the sum of the pressure induced by collision p s and the pressure induced by permanent contacts between the particles p f f (see [3] for more details).
The particle phase shear stress is modeled using a frictional rheology [7] and implemented using a regularization technique [8]: where S s i j is the velocity shear rate for the solid phase, D small is a regularisation parameter and µ(I) is the friction coefficient that depends on the inertial number I = ∇u s d ρ s / ps calculated as: where µ s is the static friction coefficient and µ 2 and I 0 are empirical coefficients (for sand : µ s = 0.63, µ 2 = 1.13 and I 0 = 0.6).
Several turbulence models are available in sedFOAM to calculate the eddy viscosity ν f t : a mixing length turbulence model [9], a k − ε model [10,11] and a k − ω model [3] (see [3] for more details).In this contribution the revisited version of the k − ω [4] (referred as k − ω2006 in this paper) has been implemented and adapted to the two-phase flow model.
The mixing length model is a turbulence model for 1D configuraions only.The eddyviscosity is calculated as ν t f = l 2 m ∇u f using only one algebraic equation: The two phase version of the two-equation turbulence models differs from classical ones for single-phase flows by the presence of additional terms in the TKE and turbulence dissipation (ε or ω) equations that account for the sediment induced turbulence damping through turbulent drag work and density stratification.Only the TKE equation is recalled here to illustrates these two terms: where P is the TKE production by the mean flow, D is the TKE diffusion term and is the rate of TKE dissipation.The above k-equation is similar to the clear fluid one, except for the last two terms on the Right Hand Side (RHS) that represents drag dissipation and density stratification, respectively.In the drag-induced damping term, the parameter α is introduced to characterize the degree of correlation between particles and fluid velocity fluctuations that can be characterized by the Stokes number S t = t p /t l where t p = ρ s /((1 − φ) K) is the particle response time and t l = k/(6ε) is the eddy turn over time.Cheng et al. [2] proposed to use the following function for α: α = e −B•S t where B is an empirical coefficient.The ε and ω equations are empirically constructed by analogy with the k equation in a similar way as for single fluid models.Empirical coefficients are introduced in front of all the terms on the RHS.The difference between the k − ω2006 and the standard k − ω model are the presence of an additional cross diffusion term in the specific dissipation equation and the presence of a stress limiter in the eddy viscosity definition.The k − ω2006 is very similar to the k − ω SST model from Menter (1994) [12] and is the best RANS model to deal with adverse pressure gradient situations [4].

One-dimensional configurations
Before applying the model to multi-dimensional configurations, a calibration on a onedimensional configuration for which experimental data are available is required.The same configuration as the one presented in Chauchat et al. (2017) [3] is used here-in to specifically test the k − ω2006 model.The configuration of Revil-Baudard et al. (2015) [13] on unidirectional sheet flow is used.The particles are made of non-spherical lightweight PMMA particles with density ρ s = 1190 kg.m −3 and mean grain size diameter d = 3 ± 0.5 mm.The fluid is water with density ρ f = 1000 kg.m −3 , kinematic viscosity ν f = 10 −6 m 2 .s−1 and the bed friction velocity is u * = 0.05 m.s −1 .
The computational domain is 0.17m high discretized using a uniform grid composed of 400 cells.The flow is driven by a constant pressure gradient and the model is run until steady state is reached.For PMMA particles, the rheological model parameters are set to µ s = 0.52, µ 2 = 0.96 and I 0 = 0.6.Four simulations are carried-out using the four turbulence models presented in section 2. Figure 1 shows the comparison of these simulations with experimental data [13] in terms of averaged mixture velocity profile (U = u s × φ + u f × (1 − φ)), volume fraction profile (φ) and shear stress profile (τ xz ).Compared with the results presented in Chauchat et al. (2017) [3], the B value has been calibrated to B = 1.Using this value, the results in terms of velocity and concentration profiles are improved.In particular, the strong turbulence damping observed in the experiments with a von Karman constant measured at κ = 0.23 can be reproduced by all four turbulence models.The k − ω2006 gives very similar results than the one obtained with the standard k − ω model.The calibration of the B value confirms the findings of Cheng et al. (2018) [14] that the turbulence damping under sheet flow conditions is due to the turbulent work of the drag force.In Cheng et al. [14] the two-phase flow model was based on kinetic theory of granular flows for the granular stress and Large Eddy Simulation for turbulence modeling.The authors found a value of B = 0.25 from the above mentioned turbulence resolving simulations.The fact that the B value is different when using the dense granular flow rheology illustrates the non-linear coupling between granular stress and fluid turbulence.This point requires further investigation that lies beyond the scope of the present study.

2D scour under a submarine pipeline
In this section the two-dimensional configuration of scour under a submarine pipeline is investigated numerically using the two-phase flow model described in section 2. Following Lee et al. [15] the configuration of Mao [16] is used as a benchmark.The pipeline has a diameter D = 5 cm placed just above a sediment bed made of medium sand with median diameter d 50 = 360 µm and density ρ s = 2600 kg.m −3 .The incoming current has a mean velocity U = 0.87 m.s −1 corresponding to an undisturbed Shields number θ = u 2 * /(s − 1)gd 50 = 0.33 where s = ρ s /ρ f is the density ratio.Initially, the pipeline is just laid on the sediment bed with no embedment.Mao [16] measured the sediment bed profile at different times until scour equilibrium was reached.Three different phases can be distinguished: onset, tunneling and lee-wake erosion.The computational domain (see figure 2) is 1.75m long and 0.305m high.An unstructured mesh has been generated using snappyHexMesh with 124 000 cells having grid resolution ranging from ∆x ∈ [7.5×10 −4 , 3×10 −3 ]m and ∆y ∈ [7.5×10 −4 , 10 −2 ]m refined around the cylinder and in the bed interface region.The outlet boundary condition for the reduced pressure is a homogeneous Dirichlet condition (p − ρgy = 0).For both the sediment and the fluid phase velocities, a homogeneous Neumann boundary condition is used for outgoing flows and a homogeneous Dirichlet boundary condition is used otherwise.This allows to avoid a non-physical incoming flow at the outlet when vortex structures are advected through the outlet.A logarithmic velocity profile (law of the wall) and constant turbulent quantities are imposed at the inlet in the same way as in Lee et al. [15].
Figure 3 shows the comparison of sediment bed profiles at t = 11, 18, and 25 seconds obtained with the two-phase flow model presented in section 2 with Mao experimental results [16] and Lee et al. [15] numerical results.Lee et al. [15] used a two-phase flow model having a granular rheology similar to the one used herein and a low Reynolds k − ε turbulence model from Launder and Sharma.The present two-phase numerical simulations were carried out using a standard k − ε and the k − ω2006 turbulence models.
The numerical results obtained with the above mentioned model and the k − turbulence model are very similar to the one obtained by Lee et al. [15].At the early stages, the comparison with experimental data from Mao [16] is very good showing that the two-phase flow model is able to reproduce the tunneling stage.However, for t ≥ 18 s the height of the sand dune at the downstream side of the pipe is overestimated.In the experiments, after the scour Fig. 3. Comparison of bed profiles at t = 11, 18, and 25 obtained using the two-phase flow model with standard k − ε turbulence model, experimental data from Mao [16] and numerical results from Lee et al. [15].hole depth reaches 0.3D, vortices are shed in the wake of the cylinder.These organized flow structures are responsible for the erosion and the resulting bed profile downstream of the pipeline, this process is known as the lee-wake erosion stage.As demonstrated in Lee et al. [15] the k − ε turbulence model is not able to reproduce the vortex shedding and is, therefore, not able to reproduce the lee-wake erosion stage.For simulations with the k − ω2006, the depth of the scour hole is largely under-estimated.The ongoing work is to understand and improve the behaviour of the k − ω2006 in the sediment bed for this configuration.The k − ω2006 is able to reproduce the vortex shedding phenomenon and should improve the numerical prediction of the lee-wake erosion.

3D scour around a vertical cylinder
The 3D configuration of the scour around a vertical cylindrical pile reported in Roulund et al. [17] is used as a benchmark for the two-phase flow model.A vertical cylinder of diameter D=0.1m is placed in a erodible sediment bed made of medium sand, with density ρ s =2650 kg.m −3 , mean grain size diameter d 50 =260 µm.The incoming current has mean velocity U=0.326 m.s −1 leading to a Reynolds number based on the pile diamater of Re D about 4.6 × 10 4 .For such value of the Reynolds number vortex shedding is expected at the lee side of the cylinder.Hydrodynamic simulations showed that vortex-shedding was correctly predicted A typical snapshot of the live-bed simulation at t=600s (10 minutes) of dynamics is presented in figure 4.a showing the vertical elevation of the iso-surface of sediment concentration φ= 0.57, a proxy for the bed elevation.A comparison with figure 33 in Roulund et al. [17] shows that the two-phase flow model is qualitatively able to reproduce the semi-circular shaped scour mark at the upstream side of the pile.The downstream bed morphology differs from the one reported in Roulund et al. [17] and this may be attributed to the interaction between vortices and sediment dynamics.A finer mesh might be required to better capture these interactions.depth S/D at the upstream side of the pile.The two-phase flow numerical results are compared with experimental and numerical results from Roulund et al. [17].The two-phase flow model is able to reproduce quantitatively the upstream scour depth evolution up to 300s.For times greater than 300s the maximum dimensionless scour depth tends to saturate at a value of S/D = 0.63, whereas the experimental data shows a continuous increase of the scour depth with a value of S/D = 0.8 at t = 600s.The numerical results from Roulund et al. [17] obtained using a classic sediment transport model (bed-load only) are not better than the present two-phase simulations.To the best of our knowledge, the present result is the first attempt to apply a two-phase flow model to a 3D scour configuration.Even if the computational time is significant (480 hours corresponding to 20 days on 224 CPUs) the good agreement obtained for the scour evolution at the upstream side of the pile proves that two-phase flow modeling is a useful tool to investigate sediment transport processes in 'complex' flow configurations.

Conclusion
In this contribution the application of a multi-dimensional two-phase flow model to 2D and 3D scour processes has been presented.For both configurations, scour below a pipeline and scour around a cylindrical pile, the two-phase flow model is able to quantitatively predict the scour erosion at the upstream side of the structures.Using the revisited k − ω from Wilcox (2008) [4] vortex shedding can be predicted at the lee side of structures.But further work is needed to improve the modeling of sediment-vortex interactions as illustrated in the 3D configuration.The results presented herein is a proof of concept that two-phase flow numerical simulations are now possible for scour processes in multi-dimensional configurations.With the increase of computers power, it could become a relevant engineering tool in hydraulic in the future.In terms of research perspectives, this new generation of models will be used to provide insight into fine scale sediment transport mechanisms involved in the scour process.

Fig. 1 .
Fig. 1.Comparison of two-phase numerical results with experiments of Revil-Baudard et al. (2015) [13] (black symbols) in terms of volume averaged velocity profiles in the left panel, sediment concentration in the center panel and Reynolds shear stress (blue lines) and granular stress (red lines) in the right panel, using the dense granular flow rheology (µ(I)) with the four turbulence models: Mixing Length ML, k − ε, k − ω and k − ω2006.

Fig. 2 .
Fig.2.Sketch of the geometry and the boundary conditions used for the computational domain.

Figure 4 .
b shows the time evolution of the maximum dimensionless scour

Fig. 4 .
Fig.4.Bed elevation after 600s of dynamics (a) and time evolution of the dimensionless scour depth at the upstream and at the downstream edge of pile.