Modeling of Fluid-Solid Two-Phase Geophysical Flows

Fluid-solid two-phase flows are frequently encountered in geophysical flow problems such as sediment transport and submarine landslides. It is still a challenge to the current experiment techniques to provide information such as detailed flow and pressure fields of each phase, which however is easily obtainable through numerical simulations using fluid-solid two-phase flow models. This chapter focuses on the Eulerian-Eulerian approach to two-phase geophysical flows. Brief derivations of the governing equations and some closure models are provided, and the numerical implementation in the finite-volume framework of OpenFOAM® is described. Two applications in sediment transport and submarine landslides are also included at the end of the chapter.


Introduction
Fluid-solid two-phase flows are important in many geophysical problems such as sediment erosion, transport and deposition in rivers or coastal environment, debris flows, scour at river or marine structures, and submarine landslides. Behaviors of fluid-solid two-phase flows are very different from those of liquid-gas two-phase flows where bubbles are dispersed in the liquid or droplets dispersed in the gas. Vast numbers of experiments on various scales have been carried out for different applications of fluid-solid two-phase flows; these experiments have advanced our understanding of bulk behaviors of some important flow characteristics. However, development of measurement techniques suitable for collecting data that contribute to understanding important physics involved in fluid-solid two-phase flows is a still-evolving science. With the modern computer technology, many data that are not obtainable currently in the experiment can be easily produced by performing time-dependent, multidimensional numerical simulations. Of course, empirical closure models required to close the governing equations still need high-quality experimental data for model validation.
Numerical approaches to two-phase flows include Eulerian-Eulerian approach, direct numerical simulations (DNS) based on Eulerian-Lagrangian formulations (Lagrangian point-particle approach), and fully resolved DNS approach [1]. Fully resolved DNS can resolve all important scales of the fluid and particles, but these simulations are currently limited to about 10 k uniform-size spheres on a Cray XE6 with 2048 cores [2], and it is not practical to use this method to model large-scale geophysical flow problems in the foreseeable future [1]. Lagrangian point-particle approach uses Eulerian formulation for the fluid phase and Lagrangian formulation for tracking the instantaneous positions of the particles. Lagrangian point-particle simulations make use of semiempirical relationships to provide both hydrodynamic force and torque acting on each particle and thus avoid modeling processes on scales smaller than Kolmogorov scale [1], making it possible to include more particles and run in a domain larger than that for fully resolved DNS. The application of Lagrangian point-particle approach is crucially dependent on the availability and accuracy of such semiempirical relationships. A recent study shows that good results can be obtained for about 100k uniform-size spherical particles in a vertical channel flow [3]; however, using this approach to investigate large-scale two-phase flow problems is still beyond the current computing capacity. Two-phase Eulerian-Eulerian approach treats both the fluid and particle phases as continuum media and is suitable for solving large-scale two-phase flow problems.
Eulerian-Eulerian two-phase flow models based on large-eddy-simulations solve a separate set of equations describing conservation of mass, momentum, and kinetic energy for each phase [4][5][6][7] and thus have the potential to consider all important processes involved in the interactions between the two phases through parameterization of particle-scale processes. This chapter introduces the basics of Eulerian-Eulerian two-phase flow modeling, its implementation in the finite-volume framework of OpenFOAM®, and two applications in geophysical flow problems.

Governing equations for fluid-solid two-phase flows
Let us consider a mixture of fluid and solid particles. Fluid can be gas, water, or a mixture of water and gas. In DNS and Lagrangian point-particle approaches to twophase flows, the flow field is solved by solving the Navier-Stokes equations, and the motion of each particle is determined by the Newton's equation of motion. In Eulerian-Eulerian two-phase flow approaches, however, the motions of individual particles are not of the interest, and the focus is on the macroscopic motion of the fluid and solid particles instead. For this purpose, the solid particles are modeled as a continuum mass through an ensemble averaging operation, which is based on the existence of possible equivalent realizations. After taking ensemble average, the mixture of fluid and particles consists of two continuous phases: the fluid (water, gas, or a mixture of water and gas) is the fluid phase, and the solid particle is the solid phase. Both phases are incompressible. The motions of the fluid and solid phases are governed by their own equations, which are obtained by taking ensemble average of the microscopic governing equations for each phase [8]. Even though some aspects of fluid-solid interaction can be considered through the ensemble average, the ensemble averaging operation itself, however, does not explicitly introduce any turbulent dispersion in the resulting equations. To consider the turbulent dispersion in the Eulerian-Eulerian description of the fluid-solid two-phase flows, another averaging operation (usually a Favre average) is needed to consider the correlations of turbulent components [5,9].

Ensemble averaged equations
At the microscopic scale, the fluid-solid mixture is a discrete system. The purpose of performing an ensemble averaging operation is to derive a set of equations describing this discrete system as a continuous system at the macroscopic scale, where the typical length scale should be much larger than one particle diameter.
In the Eulerian-Eulerian approach to two-phase flows, it is assumed that the equations governing the motion of phase k (for the fluid phase k ¼ f and for the solid phase k ¼ s) at the microscopic scale are the following equations for the conservation of mass and momentum [8,10]: and where ρ k is the density, u k is the velocity, and g is the acceleration due to gravity. The stress tensor T k includes two components: where p k is the microscopic pressure and τ k is the microscopic stress tensor. Because the fluid phase and the solid phase are immiscible, at any time t, a point in space x can be occupied only by one phase, not both. This fact can be described mathematically by the following phase function c k x; t ð Þ for phase k: The volumetric concentration of phase k is directly related to the probability of occurrence of phase k at a given location x at the time t and can be obtained by taking ensemble average of c k . Using the phase function given in Eq. (4), the volumetric concentration of phase k is obtained by taking the ensemble average of c k , denoted by c k h i. The operator ⋯ h i means taking an ensemble average of its argument.
There are several methods to derive the ensemble averaged equations governing the motion of phase k. This chapter treats the phase function as a general function and uses it to define the derivatives of the phase function c k with respect to time and space and the equation governing the evolution of c k . As stated in Drew [8], the phase function c k can be treated as a generalized function whose derivative can be defined in terms of a set of test functions. These test functions must be sufficiently smooth and have compact support so that the integration of a derivative of the phase function, weighed with the test function, is finite. The equation describing the evolution of c k is where u i is the velocity of the interface between the region occupied by the fluid phase and the region occupied by the solid phase. It is stressed here that ∇c k is zero except at the interface between two phases where ∇c k behaves like a deltafunction [8].
The ensemble averaged equations governing the motion of phase k are obtained by multiplying Eqs. (1) and (2) with c k and performing an ensemble average operation on every term in the resulting equations. When performing ensemble average operations, Reynolds' rules for algebraic operations, Leibniz' rule for time derivatives, Gauss' rule for spatial derivatives, and the following two identities are used: and The resulting equations governing the ensemble average motion of phase k are [8] ∂ and Note that ∇c k is not zero only on the interface of the region occupied by phase k (grain boundary). For the fluid-solid two-phase flows, the interface of phase k must satisfy the no-slip and no-flux conditions; therefore, u k À u i ¼ 0. As a result, the right-hand side of Eq. (8) is zero and which is the density of the interfacial force [8]. Physically, T k Á ∇c k is the microscopic density of the force acting on a surface whose normal direction is defined by ∇c k . After using Eq. (3) for T k in Eq. (9), the ensemble averaged equations can be further written in terms of the ensemble averaged qualities describing the motion of phase k as and wherec k ¼ c k h i is the volumetric concentration of phase k. Other ensemble averaged quantities used in Eqs. (12) and (13) to describe the motion of phase k at the macroscopic scale are densityρ k , pressurep k , stress tensorτ k , and velocityû k , defined byρ andt 0 k represents the c-weighted ensemble average of microscopic momentum flux associated with the fluctuation of the velocity u k around the ensemble averaged velocityû k For compressible materialsρ k is not a constant. However, for incompressible materialsρ Now we examine the limiting case where the fluid-solid system is at its static state. Because the phase functions for the two phases satisfy c f þ c s ¼ 1, both phases are not moving, andm f þm s ¼ 0, the governing equations reduce to for the fluid phase, and 0 ¼c sρs g À ∇c sps þm s , Â for the solid phase. Becausep f is the hydrostatic pressure in this case, i.e., ∇p f ¼ρ f g, it then follows thatm s ¼p f ∇c s (19) which, physically, is the buoyancy acting on the solid phase. Now Eq. (18) becomes 0 ¼c sρs g À ∇c sps Â Ã þp f ∇c s (20) which states that the weight of the solid particles is supported by the buoyancy and the interparticle forces. Therefore, the ensemble pressure of the solid phase can be written asp s ¼p f þ p s , withp f being the total fluid pressure and p s accounting for the contributions from other factors such as collision and enduring contact to the ensemble averaged pressure.
For brevity of the presentation, we shall denote simply c s by c as well c f by 1 À c and drop the symbols representing the ensemble averages hereinafter. The ensemble averaged equations governing the motion of the fluid phase are and The ensemble averaged equations governing the motion of the solid phase are and where p s denotes the contributions from interparticle interactions such as collision and enduring contact to the ensemble averaged pressure of the solid phase.
To close the equations for the fluid and solid phases, closure models are needed for τ 0 s , τ 0 f , τ s , τ f , p s , and m. It is remarked here that the definitions of the ensemble averages given in Eq. (14) do not consider the contribution from the correlations between the fluctuations of the velocities and the fluctuations of phase functions at microscopic scale; therefore, the effects of turbulent dispersion are not directly included in the ensemble averaged equations describing the motion of the each phase. In the literature, two approaches have been used to consider the turbulent dispersion: (i) considering the correlation between the fluctuations of c k h i and u f associated with the turbulent flow [9] and (ii) including a term in the model for m to account for the turbulent dispersion [8]. This chapter considers the turbulent dispersion using the first approach in the next section by taking another Favre averaging operation.
In the absence of the turbulent dispersion from m, the interphase force m should include the so-called general buoyancy p f ∇c and a component f which includes drag force, inertial force, and lift force This expression for m has been derived by [11] using a control volume/surface approach. For most fluid-solid two-phase geophysical flows, the drag force dominates f [9] and thus f can be modeled by where τ p is the so-called particle response time (i.e., a relaxation time of the particle to respond the surrounding flow). As expected, the particle response time should be related to drag coefficient and grain Reynolds number.

Favre averaged equations
The volumetric concentration and the velocities can be written as where the Favre averages are defined as and the overline stands for an integration with respect to time over a time scale longer than small-scale turbulent fluctuations but shorter than the variation of the mean flow field.
The averaged equations for the mean flow fields of the two phases are obtained by taking the following steps: (i) substituting Eq. (25) with Eq. (26) in Eqs. (22) and (24), (ii) substituting Eq. (27) in the equations obtained at step (i), and (iii) taking average of the equations obtained at step (ii) to obtain the following equations: for the fluid phase, with τ 00 f being defined by and for the solid phase, with τ 00 s being defined by It is remarked here that the terms 1 Àc ð Þ∇p f in Eq. (30) andc∇p f in Eq. (33) have been obtained by using the expression for m given in Eq. (25). In order to close these averaged equations, closure models are required for the following terms: c τ s þ τ 0 s þ τ 00 , cu 00 f , and c 00 ∇p 00 f . The last term can be neglected based on an analysis of their orders of magnitude by Drew [12].
The term cu 00 f is approximated by the following gradient transport hypotheses: where ν ft is the eddy viscosity and σ c is the Schmidt number, which represents the ratio of the eddy viscosity of the fluid phase to the eddy diffusivity of the solid phase. Furthermore, the following approximations are introduced: For brevity of the presentation, the symbols representing Favre averages are dropped hereinafter, and the final equations governing the conservation of mass and momentum of each phase are for the fluid phase and for the solid phase.

Closure models 3.1 Stresses for the fluid phase
The stress tensor for the fluid phase τ f includes two parts: a part for the viscous stress, τ v f , and the other part for the turbulent Reynolds stress, τ t The viscous stress tensor τ v f is usually computed by where ν f is the kinematic viscosity of the fluid phase and where the superscript T denotes a transpose. Some studies [13] suggested modifying ν f to consider the effect of the solid phase; other studies [14], however, obtained satisfactory results even without considering this effect. The stress tensor τ t f is related to the turbulent characteristics, which need to be provided by solving a turbulent closure model such as k À ϵ or k À ω model. For a k À ϵ model with low-Reynolds-number correction [15], τ t f can be computed by where k is the turbulence kinetic energy and ν t f is the eddy viscosity of the fluid phase, given by with ϵ being the turbulent dissipation of the fluid phase to be provide by solving the represents the low-Reynolds-number correction with Re t ¼ k 2 =ν f ϵ. The coefficient C μ is usually assumed to be a constant. The equations governing k and ϵ are similar to those for clear water [15] ∂ρ and where coefficients C ϵ1 , C ϵ2 , σ ϵ , σ k , and f 2 are model parameters, whose values can be taken the same as those in the k À ϵ model for clear fluid under low-Reynoldsnumber conditions [15]. There are two terms inside the curly brackets, and both terms account for the turbulence modulation by the presence of particles: the first term is associated with the general buoyancy, and the second term is due to the correlation of the fluctuating velocities of solid and fluid phases. C ϵ3 ¼ 1 is usually adopted in the literature [28]; however, it is remarked that the value of C ϵ3 is not well understood at the present and a sensitivity test to understand how the value of this C ϵ3 on the simulation results is recommended. The parameter α reflects the correlation between the solid-phase and fluid-phase turbulent motions and is given by where τ l ¼ 0:165k=ϵ is a time scale for the turbulent flow and τ c is a time scale for particle collisions given by [16] τ c ¼ c rcp c with c rcp being the random close packing fraction and d being the particle diameter.
c rcp is 0.634 for spheres [17]. The term c rcp =c À Á 1=3 À 1 is related to the ratio of the mean free dispersion distance to the diameter of the solid particle. It is remarked here that the presence of solid particles in the turbulent flow may either enhance (for large particles) or reduce (for small particles) the turbulence [18]. The k À ϵ model given here can only reflect the reduction of turbulence and thus is not suitable for problems with large particles. Other turbulence models [7,18] include a term describing the enhancement of turbulence; however, including that term in the present model may induce numerical instability in some cases.

Stresses for the solid phase
The closure models for p s and τ s used in Lee et al. [16] will be described here. In order to cover flow regimes with different solid-phase concentrations (dilute flows, dense flows, and compact beds), Lee et al. [16] suggested the following model for p s : where p t s accounts for the turbulent motion of solid particles (important for dilute flows); p r s reflects the rheological characteristics of dense flows and includes the effects such as fluid viscosity, enduring contact, and particle inertial; p e s accounts for the elastic effect, which is important when the particles are in their static state in a compact bed.
For solid particles in a compact bed, the formula proposed by Hsu et al. [19] can be used to compute p e s p e where c o is random loose packing fraction and coefficients K and χ are model parameters. For spheres, c o ranges from 0.54 to 0.634, depending on the friction [17]. The coefficient K is associated with the Young's modulus of the compact bed, and the other terms are related to material deformation. The closure models for p r s and p t s are closely related to the stress tensor and the visco-plastic rheological characteristics for the solid phase. The stress tensor for the solid phase can be computed by The kinematic viscosity of the solid phase ν s is computed by the sum of two terms: where ν v s and ν t s represent the visco-plastic and turbulence effects, respectively. This model for ν s can consider both the turbulence behavior (for dilute flows) and the visco-plastic behavior (for dense flows and compact beds).
Based on an analysis of heavy and small particles in homogeneous steady turbulent flows, Hinze [20] suggests that p t s and ν t s can be computed by and where the coefficient α is the same as that in Eqs.(45) and (46). For dense fluid-solid two-phase flows, the visco-plastic rheological characteristics depend on a dimensionless parameter I ¼ I v þ aI 2 i , where I v is the viscous number, I i is the inertial number, and a is a constant [21]. The viscous number is defined by I v ¼ 2ρ f ν f D s =cp s where ν f is the kinematic viscosity of the fluid and D s is the second invariant of the strain rate. Physically, the viscous number describes the ratio of the viscous stress to the quasi-static shear stress associated with the weight (resulting from the enduring contact). The inertial number is defined by , which describes the ratio of the inertial stress to the quasi-static stress. The relative importance of the inertial number to the viscous number can be measured by the Stokes number st v ¼ I 2 i =I ν . Some formulas have been proposed in the literature to describe c À I and η À I relationships, where η ¼ T s =p s with T s being the second invariant of τ s .
Following the work of Boyer et al. [22], Lee et al. [16] assumed where c c is a critical concentration and b is a model parameter. Trulsson et al. [21] proposed where η 1 ¼ tan θ s with θ s being the angle of repose and η 2 and I o are constants. Based on Eqs. (56) and (55), the following expressions for p r s and ν v s can be derived [16]: which considers the solid phase in its static state as a very viscous fluid and where b is a constant. In Lee et al. [7], a ¼ 0:11 and b ¼ 1 were taken.

Closure models for particle response time
The drag force between the two phases is modeled through the particle response time τ p . Three representative models for particle response time are introduced in this section.

A model based on the particle sedimentation in still water
The first model is based on particle sedimentation in still water, which can be simplified as a one-dimensional problem, where the steady sedimentation assures that there are no stresses in both the solid and fluid phases in the vertical direction z. In this case, Eqs. (38) and (40) reduce to and where w f and w s are the vertical velocities of the fluid and solid phases, respectively. Because net volume flux through any horizontal plane must be zero, we have Combining Eqs. (59) and (61) yields Substituting Eqs. (61) and (62) into Eq. (60) leads to where the solid-phase velocity w s is also called the hindered settling velocity [23]. The hindered velocity is smaller than the terminal velocity of a single particle, w 0 , due to the influence of volumetric concentration (including many-body hydrodynamic interactions). Richardson and Zaki [24] suggested where the coefficient n is related to the particle Reynolds number Re s ¼ w 0 d=ν f n ¼ The terminal velocity of a single particle w 0 can be computed by where C d is the drag coefficient for steady flows passing a single particle [25,26]. For spheres, the following formula of White [27] can be used: where Re p ¼ |u f À u s |d=ν f . Combing Eqs. (63)-(67) yields It is remarked that Eq. (64) is validated only for c , 0:4 [28]. When the concentration c is so high that contact networks form among particles, w s , becomes zero; when this happens, Eq. (64) is no longer valid any more.

A model based on the pressure drop in steady flows through a homogeneous porous media
Another model for particle response time can be derived by examining the pressure drop in the steady flow through a porous media. For a one-dimensional problem of a horizontal, steady flow through porous media, the terms containing the stresses of the fluid phase disappear, and Eq. (38) reduces to where the horizontal coordinate x points in the direction of the flow and u is the velocity component in x-direction.
For this problem, Forchheimer [29] suggested where a F and b F are two model parameters. Several formulas for computing a F and b F can be found in previous studies. The following two expressions for a F and b F suggested by Engelund [25] are recommended for the applications presented at the end of this chapter: Comparing Eqs. (69) and (70) and using Eq.(71) give where a E and b E are two model parameters depending on the composition of the solid phase. The parameter a E is associated with k p as will be shown later. For d ≈ 2 Â 10 À4 m, k p ≈ 10 À10 $ 10 À11 m 2 [30], which gives a E ≈ 1:6 Â 10 3 $ 1:6 Â 10 4 for c ¼ 0:5. The parameter b E varies from 1.8 to 3.6 or more [28,31,32]. For flow in a porous media, the particle response time can also be related to its permeability κ p . According to Darcy's law for seepage [29], the pressure gradient can also be written as where κ p is the permeability. Combining Eqs. (69) and (73) gives When the flow is very slow, Eqs. (70), (71), and (73) suggest that which means that the particle response time can be related to the permeability.

A hybrid model
Equation (64) is validated only for c , 0:4 [28]. To extend Eq. (64) to high concentration regions, Camenen [33] modified Eq. (64) to where c m is the maximum concentration at which w ¼ 0. In this study, c m ¼ c o is adopted because when c ≥ c o , contact networks can form in the granular material. Combining Eqs. (63), (76), and (66)-(67) gives We stress that c ¼ c m will lead to τ p ¼ 0 and thus an infinite drag force. Physically, when the volumetric concentration is greater than some critical value, say c r , Eq. (63) ceases to be valid, and Eq. (72) should be used. To avoid unnaturally large drag force between the two phases, we propose the following model for particle response time: For given values of a E and b E , Eq. (79) implicitly defines c r as a function of Re p .

Introduction to OpenFOAM
This section introduces how to use OpenFOAM® to solve the governing equations with the closure models presented in the previous section. OpenFOAM® is a C++ toolbox developed based on the finite-volume method; it allows CFD code developers to sidestep the discretization of derivative terms on unstructured grids.

Semidiscretized forms of the governing equations
To avoid numerical noises occurring when c ! 0, Rusche [34] suggests that the momentum equations (Eqs. (38) and (40)) should be converted into the following "phase-intensive" form by dividing ρ f 1 À c ð Þand ρ s c: ∂u f we can set u s ¼ u f , which means the solid particles completely follow the water particles; this does not affect the computations of other variables because the momentum of the solid phase cu s is very small when c ≤ 10 À6 . Because the maximum value of c is always smaller than 1, there is no singularity issue with Eq. (82). An iteration procedure is needed to solve the governing equations at each time step for the values of c,u f ,û s , and p f obtained at the previous time step, and it is outlined below: Figure 1 is a flowchart showing these 12 solution steps. In the absence of the solid phase, the numerical scheme outlined here reduces to the "PIMPLE" scheme, which is a combination of the "pressure implicit with splitting of operator" (PISO) scheme and the "semi-implicit method for pressure-linked equations" (SIMPLE) scheme. Iterations need to be done separately to solve To ensure the stability of the overall numerical scheme, the Courant-Friedrichs-Lewy (CFL) condition must be satisfied for each cell. The local Courant number for each cell, which is related to the ratio between the distance of a particle moving within Δt and the size of the cell where such particle is located, is defined as where in u j ¼ 1 À c ð Þu j f þ cu j s , the subscript "j" represents the j th face of the cell, S j is a unit normal vector, V is the volume of the cell, and Δt is the time step. The Courant number must be less than 1 to avoid numerical instability. Generally, max (CFL) <0.1 is suggested. The values of CFL for high concentration regions should be much smaller than those for low concentration regions so that rapid changes of c can be avoided. Therefore, it is recommended that max CFLj c>c o < 0.005. The time step is recommended to be in the range of 10 À5 and 10 À4 s.

Applications
This section briefly describes two examples that have been studied using the two-phase flow models described. The problem descriptions and numerical setups for these two problems are included here; for other relevant information, the reader is referred to Lee and Huang [35] and Lee et al. [38].

Scour downstream of a sluice gate
A sluice gate is a hydraulic structure used to control the flow in a water channel. Sluice gate structures usually have a rigid floor followed by an erodible bed. The scour downstream of a sluice gate is caused by the horizontal submerged water jet issuing from the sluice gate. It is of practical importance to understand the maximum scour depth for the safety of a sluice gate structure. Many experimental studies have been done to investigate the maximum scour depth and the evolution of scour profile (e.g., Chatterjee et al. [39]). For numerical simulations, this problem includes water (fluid phase) and sediment (solid phase) and is best modeled by a liquid-solid two-phase flow approach. In the following, the numerical setup and main conclusions used in Lee et al. [38] are briefly described. The experimental setup of Chatterjee et al. [39] is shown in Figure 2. To numerically simulate the experiment of [8], we use the same sand and dimensions to set up the numerical simulations: quartz sand with ρ s ¼ 2650 kg/m 3 and d ¼ 0:76 mm is placed in the sediment reservoir, with its top surface being on the same level as the top surface of the apron; the sluice gate opening is 2 cm; the length of apron is 0:66 m; the sediment reservoir length is 2:1 m; the overflow weir on the right end has a height of 0:239 m; the upstream inflow discharge rate at the sluice opening is 0:204 m 2 /s, which translates into an average horizontal flow velocity V ¼ 1:02 m/s under the sluice gate. As an example, the computed development of scour depth d s is shown in Figure 3 together with the measurement of Chatterjee et al. [39].
The problem involves also an air-water surface, which can be tracked using a modified volume-of-fluid method introduced in [38]. A nonuniform mesh is used in the two-phase flow simulation because of the air-water interface, the interfacial momentum transfer at the bed, and the large velocity variation due to the water jet. The finest mesh with a vertical mesh resolution of 2d is used in the vicinity of the sediment-fluid interface; this fine mesh covers the dynamic sediment-fluid  interface during the entire simulation. In regions away from the sediment-fluid interface or regions where the scouring is predicted to be negligible (e.g., further downstream the scour hole), the mesh sizes with a vertical resolution ranging from 3 to 5 mm are used. The aspect ratio of the mesh outside the wall jet region is less than 3.0. Since in the wall jet, horizontal velocity is significantly larger than the vertical velocity, the aspect ratio of the local mesh in the wall jet region is less than 5.0.
The scour process is sensitive to the model for particle response time used in the simulation. Because Eq. (72) can provide a better prediction of sediment transport rate for small values of Shields parameter, it is recommended for this problem. The two-phase flow model can reproduce well the measured scour depth and the location of sand dune downstream of the scour hole.

Collapse of a deeply submerged granular column
Another application of the fluid-solid two-phase flow simulation is the simulation of the collapse of a deeply submerged granular column. The problem is best described as a granular flow problem, which involves sediment (a solid phase) and water (fluid phase). Many experimental studies have been reported in the literature on this topic. This section describes a numerical simulation using the fluid-solid two-phase flow model described in this chapter. Figure 4 shows the experimental setup of Rondon et al. [40]. A 1:1 scale twophase flow simulation was performed by Lee and Huang [35] using the fluid-solid two-phase flow model presented in this chapter. The diameter and the density of the sand grain are 0.225 mm and 2500 kg/m 3 , respectively. The density and the dynamic viscosity of the liquid are 1010 kg/m 3 and 12 mPa s, respectively. Note that the viscosity of the liquid in the experiment is ten times larger than that for water at room temperature. For this problem, using a mesh of 1.0 Â 1.0 mm and the particle response model given by Eq. (78), the fluid-solid two-phase flow model presented in this chapter can reproduce well the collapse process reported in Rondon et al. [40]. Figure 5 shows the simulated collapsing processes compared with the measurement for two initial packing conditions: initially loosely packed condition and initially densely packed condition.
The two-phase model and closure models presented in this chapter are able to deal with both initially loose packing and initially dense packing conditions and reveal the roles played by the contractancy inside the granular column with a loose packing and dilatancy inside a granular column with a dense packing. One of the conclusions of Lee and Huang [35] is that the collapse process of a densely packed granular column is more sensitive to the model used for particle response time than that of a loosely packed granular column. The particle response model given by Eq. (78) performs better than other models; this is possibly because the liquid used in Rondon et al. [40] is much viscous than water.

Summary
This chapter presented a brief introduction to the equations and closure models suitable for fluid-solid two-phase flow problems such as sediment transport, submarine landslides, and scour at hydraulic structures. Two averaging operations were performed to derive the governing equations so that the turbulent dispersion, important for geophysical flow problems, can be considered. A new model for the rheological characteristics of sediment phase was used when computing the stresses of the solid phase. The k À ϵ model was used to determine the Reynolds stresses. A hybrid model to compute the particle response time was introduced, and the numerical implementation in the framework of OpenFOAM® was discussed. A numerical scheme was introduced to avoid numerical instability when the concentration is high. Two applications were describe to show the capacity of the twophase flow models presented in this chapter.