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.
- granular flows
- submarine landslides
- sediment transport
- continuum model
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 . 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 , and it is not practical to use this method to model large-scale geophysical flow problems in the foreseeable future . 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 , 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 ; 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.
2. 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 two-phase 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 . 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].
2.1 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 (for the fluid phase and for the solid phase ) at the microscopic scale are the following equations for the conservation of mass and momentum [8, 10]:
where is the density, is the velocity, and is the acceleration due to gravity. The stress tensor includes two components:
where is the microscopic pressure and is the microscopic stress tensor.
Because the fluid phase and the solid phase are immiscible, at any time , a point in space can be occupied only by one phase, not both. This fact can be described mathematically by the following phase function for phase :
The volumetric concentration of phase is directly related to the probability of occurrence of phase at a given location at the time and can be obtained by taking ensemble average of . Using the phase function given in Eq. (4), the volumetric concentration of phase is obtained by taking the ensemble average of , denoted by . The operator means taking an ensemble average of its argument.
There are several methods to derive the ensemble averaged equations governing the motion of phase . This chapter treats the phase function as a general function and uses it to define the derivatives of the phase function with respect to time and space and the equation governing the evolution of . As stated in Drew , the phase function 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 is
where 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 is zero except at the interface between two phases where behaves like a delta-function .
The ensemble averaged equations governing the motion of phase are obtained by multiplying Eqs. (1) and (2) with 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:
The resulting equations governing the ensemble average motion of phase are 
Note that is not zero only on the interface of the region occupied by phase (grain boundary). For the fluid-solid two-phase flows, the interface of phase must satisfy the no-slip and no-flux conditions; therefore, . As a result, the right-hand side of Eq. (8) is zero and
which is the density of the interfacial force . Physically, is the microscopic density of the force acting on a surface whose normal direction is defined by .
where is the volumetric concentration of phase . Other ensemble averaged quantities used in Eqs. (12) and (13) to describe the motion of phase at the macroscopic scale are density , pressure , stress tensor , and velocity , defined by
and represents the -weighted ensemble average of microscopic momentum flux associated with the fluctuation of the velocity around the ensemble averaged velocity
For compressible materials 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 , both phases are not moving, and , the governing equations reduce to
for the fluid phase, and
for the solid phase.
Because is the hydrostatic pressure in this case, i.e., , it then follows that
which, physically, is the buoyancy acting on the solid phase. Now Eq. (18) becomes
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 as , with being the total fluid pressure and 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 by as well by and drop the symbols representing the ensemble averages hereinafter. The ensemble averaged equations governing the motion of the fluid phase are
The ensemble averaged equations governing the motion of the solid phase are
where 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 , , , , , and .
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 and associated with the turbulent flow  and (ii) including a term in the model for to account for the turbulent dispersion . 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 , the interphase force should include the so-called general buoyancy and a component which includes drag force, inertial force, and lift force
where 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.
2.2 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 being defined by
for the solid phase, with being defined by
In order to close these averaged equations, closure models are required for the following terms: , , , and . The last term can be neglected based on an analysis of their orders of magnitude by Drew . The term is approximated by the following gradient transport hypotheses:
where is the eddy viscosity and 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.
3. Closure models
3.1 Stresses for the fluid phase
The stress tensor for the fluid phase includes two parts: a part for the viscous stress, , and the other part for the turbulent Reynolds stress,
The viscous stress tensor is usually computed by
where is the kinematic viscosity of the fluid phase and , where the superscript denotes a transpose. Some studies  suggested modifying to consider the effect of the solid phase; other studies , however, obtained satisfactory results even without considering this effect.
The stress tensor is related to the turbulent characteristics, which need to be provided by solving a turbulent closure model such as or model. For a model with low-Reynolds-number correction , can be computed by
where is the turbulence kinetic energy and 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 equation. The coefficient represents the low-Reynolds-number correction with . The coefficient is usually assumed to be a constant.
The equations governing and are similar to those for clear water 
where coefficients , , , , and are model parameters, whose values can be taken the same as those in the model for clear fluid under low-Reynolds-number conditions . 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. is usually adopted in the literature ; however, it is remarked that the value of is not well understood at the present and a sensitivity test to understand how the value of this 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 is a time scale for the turbulent flow and is a time scale for particle collisions given by 
with being the random close packing fraction and being the particle diameter. is 0.634 for spheres . The term 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 . The 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.
3.2 Stresses for the solid phase
The closure models for and used in Lee et al.  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.  suggested the following model for :
where accounts for the turbulent motion of solid particles (important for dilute flows); reflects the rheological characteristics of dense flows and includes the effects such as fluid viscosity, enduring contact, and particle inertial; 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.  can be used to compute
where is random loose packing fraction and coefficients and are model parameters. For spheres, ranges from 0.54 to 0.634, depending on the friction . The coefficient is associated with the Young’s modulus of the compact bed, and the other terms are related to material deformation.
The closure models for and 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 is computed by the sum of two terms:
where and represent the visco-plastic and turbulence effects, respectively. This model for 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  suggests that and can be computed by
For dense fluid-solid two-phase flows, the visco-plastic rheological characteristics depend on a dimensionless parameter , where is the viscous number, is the inertial number, and is a constant . The viscous number is defined by where is the kinematic viscosity of the fluid and 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 . Some formulas have been proposed in the literature to describe and relationships, where with being the second invariant of .
where is a critical concentration and is a model parameter. Trulsson et al.  proposed
which considers the solid phase in its static state as a very viscous fluid and
where is a constant. In Lee et al. , and were taken.
3.3 Closure models for particle response time
The drag force between the two phases is modeled through the particle response time . Three representative models for particle response time are introduced in this section.
3.3.1 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 . In this case, Eqs. (38) and (40) reduce to
where and are the vertical velocities of the fluid and solid phases, respectively.
Because net volume flux through any horizontal plane must be zero, we have
where the solid-phase velocity is also called the hindered settling velocity . The hindered velocity is smaller than the terminal velocity of a single particle, , due to the influence of volumetric concentration (including many-body hydrodynamic interactions). Richardson and Zaki  suggested
where the coefficient is related to the particle Reynolds number
The terminal velocity of a single particle can be computed by
It is remarked that Eq. (64) is validated only for . When the concentration is so high that contact networks form among particles, , becomes zero; when this happens, Eq. (64) is no longer valid any more.
3.3.2 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 points in the direction of the flow and is the velocity component in -direction.
For this problem, Forchheimer  suggested
where and are two model parameters. Several formulas for computing and can be found in previous studies. The following two expressions for and suggested by Engelund  are recommended for the applications presented at the end of this chapter:
where and are two model parameters depending on the composition of the solid phase. The parameter is associated with as will be shown later. For m, , which gives for . The parameter 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 . According to Darcy’s law for seepage , the pressure gradient can also be written as
which means that the particle response time can be related to the permeability.
3.3.3 A hybrid model
where is the maximum concentration at which . In this study, is adopted because when , contact networks can form in the granular material.
We stress that will lead to and thus an infinite drag force. Physically, when the volumetric concentration is greater than some critical value, say , 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:
where is the concentration at the intercept point of Eq. (72) and Eq. (77). The transition from Eq. (77) to Eq. (72) is continuous at the intercept point where . The concentration at the point joining the two models () is problem-dependent and can be found in principle by solving the following equation:
For given values of and , Eq. (79) implicitly defines as a function of .
4. Numerical implementation with OpenFOAM
4.1 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.
4.2 Semidiscretized forms of the governing equations
where (= or ) denotes the systems of linear algebraic equations arising from the discretization of either Eqs. (82) or (83). The matrix is decomposed into a diagonal matrix, , and an off-diagonal matrix, . Also, with relating to the second to final terms on the right-hand side of either Eqs. (82) or (83). OpenFOAM® built-in functions are used to compute and , which depend on the discretization schemes. For example, Lee et al.  and Lee and Huang  used a second-order time-implicit scheme and a limited linear interpolation scheme for all variables except for velocity. To interpolate velocities, the total-variation-diminishing (TVD) limited linear interpolation scheme is adopted for velocity.
4.3 A prediction-correction method
If Eq. (83) is directly used to calculate and Eq. (39) to calculate , then may increase rapidly toward , leading to an infinite for large . This can be avoided by using a prediction-correction method to compute and . This is achieved by splitting Eq. (83) into a predictor and a corrector. The predictor is
which is corrected by the following corrector
This predictor-corrector scheme can improve the numerical stability by introducing a numerical diffusion term. To see this, we combine Eqs. (39) and (85) to obtain the following equation describing the evolution of :
The right-hand side of Eq. (86) now has a diffusive term introduced by the numerical scheme. High sediment concentration and large increase the numerical diffusion (the right-hand side of Eq. (86)) and thus can avoid a rapid increase of and the numerical instability due to high sediment concentration.
For the velocity-pressure coupling, Eq. (82) is similarly solved using a predictor and a corrector. The predictor is
which is corrected by the following corrector
Substituting Eq. (88) into Eq. (37) gives a pressure equation. However, when using this pressure equation to simulate air-water flows, numerical experiments have shown that the lighter material is poorly conserved . The poor conservation of lighter material can be avoided by combining Eqs. (37) and (39) into the following Eq. (37):
The numerical diffusion term on the right-hand side of Eq. (91) can help improve the numerical stability.
The prediction-correction method presented here deals with velocity-pressure coupling and avoids the numerical instability caused by high concentration. The turbulence closure model is also solved in “phase-intensive” forms. For other details relating to the numerical treatments, the reader is referred to “twoPhaseEulerFoam,” a two-phase solver provided by OpenFOAM®.
4.3.1 Outline of the solution procedure
When , Eq. (83) becomes singular. To avoid this, is replaced by in numerical computations, where is a very small number, say . When , only a very small amount of solid particles are moving with the fluid; replacing by may introduce error in computing ; to avoid this error, we can set , 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 is very small when . Because the maximum value of 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 ,, , and obtained at the previous time step, and it is outlined below:
Compute from Eq. (84).
Solve Eq. (86) for .
Compute from Eq. (85).
Compute from Eq. (90).
Compute from Eq. (87).
Solve Eq. (91) for .
Compute from Eq. (88).
Set for very dilute region, specifically .
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 Eq. (80) for , Eq. (81) for , Eq. (86) for , Eq. (91) for , Eq. (45) for , and (46) for ; the convergence criteria are set at the residuals not exceeding . Because Eqs. (80), (81), (86), and (87) are coupled, additional residual checks need to be performed at step 11; however, the residual for Eq. (81) is not checked because is enforced in step 10.
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 and the size of the cell where such particle is located, is defined as
where in , the subscript “” represents the face of the cell, is a unit normal vector, is the volume of the cell, and 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 can be avoided. Therefore, it is recommended that < 0.005. The time step is recommended to be in the range of and s.
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  and Lee et al. .
5.1 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. ). 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.  are briefly described. The experimental setup of Chatterjee et al.  is shown in Figure 2. To numerically simulate the experiment of , we use the same sand and dimensions to set up the numerical simulations: quartz sand with kg/m3 and 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 m; the sediment reservoir length is m; the overflow weir on the right end has a height of m; the upstream inflow discharge rate at the sluice opening is m2/s, which translates into an average horizontal flow velocity m/s under the sluice gate. As an example, the computed development of scour depth is shown in Figure 3 together with the measurement of Chatterjee et al. .
The problem involves also an air-water surface, which can be tracked using a modified volume-of-fluid method introduced in . 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 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.
5.2 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. . A 1:1 scale two-phase flow simulation was performed by Lee and Huang  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/m3, respectively. The density and the dynamic viscosity of the liquid are 1010 kg/m3 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. . 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  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.  is much viscous than water.
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 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 two-phase flow models presented in this chapter.
This material presented here is partially based upon work supported by the National Science Foundation under Grant No. 1706938 and the Ministry of Science and Technology, Taiwan [MOST 107-2221-E-032-018-MY3]. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation.