Open access peer-reviewed chapter - ONLINE FIRST

Modeling of Fluid-Solid Two-Phase Geophysical Flows

By Zhenhua Huang and Cheng-Hsien Lee

Submitted: May 10th 2018Reviewed: September 11th 2018Published: December 19th 2018

DOI: 10.5772/intechopen.81449

Downloaded: 238


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
  • scour
  • continuum model
  • OpenFOAM®

1. 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.

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 [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].

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 k(for the fluid phase k=fand for the solid phase k=s) at the microscopic scale are the following equations for the conservation of mass and momentum [8, 10]:




where ρkis the density, ukis the velocity, and gis the acceleration due to gravity. The stress tensor Tkincludes two components:


where pkis the microscopic pressure and τkis the microscopic stress tensor.

Because the fluid phase and the solid phase are immiscible, at any time t, a point in space xcan be occupied only by one phase, not both. This fact can be described mathematically by the following phase function ckxtfor phase k:


The volumetric concentration of phase kis directly related to the probability of occurrence of phase kat a given location xat the time tand can be obtained by taking ensemble average of ck. Using the phase function given in Eq. (4), the volumetric concentration of phase kis obtained by taking the ensemble average of ck, denoted by ck. 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 k. This chapter treats the phase function as a general function and uses it to define the derivatives of the phase function ckwith respect to time and space and the equation governing the evolution of ck. As stated in Drew [8], the phase function ckcan 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 ckis


where uiis 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 ckis zero except at the interface between two phases where ckbehaves like a delta-function [8].

The ensemble averaged equations governing the motion of phase kare obtained by multiplying Eqs. (1) and (2) with ckand 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 kare [8]






Note that ckis 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 kmust satisfy the no-slip and no-flux conditions; therefore, ukui=0. As a result, the right-hand side of Eq. (8) is zero and


which is the density of the interfacial force [8]. Physically, Tkckis the microscopic density of the force acting on a surface whose normal direction is defined by ck.

After using Eq. (3) for Tkin Eq. (9), the ensemble averaged equations can be further written in terms of the ensemble averaged qualities describing the motion of phase kas




where c˜k=ckis the volumetric concentration of phase k. Other ensemble averaged quantities used in Eqs. (12) and (13) to describe the motion of phase kat the macroscopic scale are density ρ˜k, pressure p˜k, stress tensor τ˜k, and velocity ûk, defined by


and t˜krepresents the c-weighted ensemble average of microscopic momentum flux associated with the fluctuation of the velocity ukaround the ensemble averaged velocity ûk


For compressible materials ρ˜kis 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 cf+cs=1, both phases are not moving, and m˜f+m˜s=0, the governing equations reduce to


for the fluid phase, and


for the solid phase.

Because p˜fis the hydrostatic pressure in this case, i.e., p˜f=ρ˜fg, 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 p˜s=p˜f+ps, with p˜fbeing the total fluid pressure and psaccounting 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 csby cas well cfby 1cand 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 psdenotes 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 τs, τf, τs, τf, ps, 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 ckand u fassociated with the turbulent flow [9] and (ii) including a term in the model for mto 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 mshould include the so-called general buoyancy pfcand a component fwhich includes drag force, inertial force, and lift force


This expression for mhas 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 fcan be modeled by


where τpis 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 τf′′being defined by




for the solid phase, with τs′′being defined by


It is remarked here that the terms 1c˜p˜fin Eq. (30) and c˜p˜fin Eq. (33) have been obtained by using the expression for mgiven in Eq. (25).

In order to close these averaged equations, closure models are required for the following terms: cτs+τs+τs′′¯, 1cτf+τf+τf′′¯, cuf′′¯, and c′′pf′′¯. The last term can be neglected based on an analysis of their orders of magnitude by Drew [12]. The term cuf′′¯is approximated by the following gradient transport hypotheses:


where νftis the eddy viscosity and σcis 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 τfincludes two parts: a part for the viscous stress, τfv, and the other part for the turbulent Reynolds stress, τft


The viscous stress tensor τfvis usually computed by


where νfis the kinematic viscosity of the fluid phase and Df=uf+ufT/2, where the superscript Tdenotes a transpose. Some studies [13] suggested modifying νfto consider the effect of the solid phase; other studies [14], however, obtained satisfactory results even without considering this effect.

The stress tensor τftis 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], τftcan be computed by


where kis the turbulence kinetic energy and νftis the eddy viscosity of the fluid phase, given by


with ϵbeing the turbulent dissipation of the fluid phase to be provide by solving the kϵequation. The coefficient fμ=exp3.4/1+Ret/502represents the low-Reynolds-number correction with Ret=k2/νfϵ. The coefficient Cμis usually assumed to be a constant.

The equations governing kand ϵare similar to those for clear water [15]



ρf1cϵt+ρf1cufϵ=ϵkCϵ1 f11cτf:ufCϵ2 f2ρf1cϵ+ρfνftσϵ1cϵϵkCϵ3ρsρfνffσccg+2ρsc1αkτp,E46

where coefficients Cϵ1, Cϵ2, σϵ, σk, and f2are model parameters, whose values can be taken the same as those in the kϵmodel for clear fluid under low-Reynolds-number 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=1is usually adopted in the literature [28]; however, it is remarked that the value of Cϵ3is not well understood at the present and a sensitivity test to understand how the value of this Cϵ3on 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 τcis a time scale for particle collisions given by [16]


with crcpbeing the random close packing fraction and dbeing the particle diameter. crcpis 0.634 for spheres [17]. The term crcp/c1/31is 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.

3.2 Stresses for the solid phase

The closure models for psand τsused 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 ps:


where pstaccounts for the turbulent motion of solid particles (important for dilute flows); psrreflects the rheological characteristics of dense flows and includes the effects such as fluid viscosity, enduring contact, and particle inertial; pseaccounts 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 pse


where cois random loose packing fraction and coefficients Kand χare model parameters. For spheres, coranges from 0.54 to 0.634, depending on the friction [17]. The coefficient Kis associated with the Young’s modulus of the compact bed, and the other terms are related to material deformation.

The closure models for psrand pstare 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 νsis computed by the sum of two terms:


where νsvand νstrepresent the visco-plastic and turbulence effects, respectively. This model for νscan 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 pstand νstcan be computed by




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=Iv+aIi2, where Ivis the viscous number, Iiis the inertial number, and ais a constant [21]. The viscous number is defined by Iv=2ρfνfDs/cpswhere νfis the kinematic viscosity of the fluid and Dsis 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 Ii=2dDs/cps/ρs, 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 stv=Ii2/Iν. Some formulas have been proposed in the literature to describe cIand ηIrelationships, where η=Ts/pswith Tsbeing the second invariant of τs.

Following the work of Boyer et al. [22], Lee et al. [16] assumed


where ccis a critical concentration and bis a model parameter. Trulsson et al. [21] proposed


where η1=tanθswith θsbeing the angle of repose and η2and Ioare constants. Based on Eqs. (56) and (55), the following expressions for psrand νsvcan be derived [16]:


which considers the solid phase in its static state as a very viscous fluid and


where bis a constant. In Lee et al. [7], a=0.11and b=1were taken.

3.3 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.

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 z. In this case, Eqs. (38) and (40) reduce to




where wfand wsare 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 wsis also called the hindered settling velocity [23]. The hindered velocity is smaller than the terminal velocity of a single particle, w0, due to the influence of volumetric concentration (including many-body hydrodynamic interactions). Richardson and Zaki [24] suggested


where the coefficient nis related to the particle Reynolds number Res=w0d/νf


The terminal velocity of a single particle w0can be computed by


where Cdis the drag coefficient for steady flows passing a single particle [25, 26]. For spheres, the following formula of White [27] can be used:


where Rep=ufusd/νf. Combing Eqs. (63)(67) yields


It is remarked that Eq. (64) is validated only for c<0.4[28]. When the concentration cis so high that contact networks form among particles, ws, 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 xpoints in the direction of the flow and uis the velocity component in x-direction.

For this problem, Forchheimer [29] suggested


where aFand bFare two model parameters. Several formulas for computing aFand bFcan be found in previous studies. The following two expressions for aFand bFsuggested 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 aEand bEare two model parameters depending on the composition of the solid phase. The parameter aEis associated with kpas will be shown later. For d2×104m, kp10101011m2[30], which gives aE1.6×1031.6×104for c=0.5. The parameter bEvaries 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 κpis 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.

3.3.3 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 cmis the maximum concentration at which w=0. In this study, cm=cois adopted because when cco, contact networks can form in the granular material.

Combining Eqs. (63), (76), and (66)(67) gives


We stress that c=cmwill lead to τp=0and thus an infinite drag force. Physically, when the volumetric concentration is greater than some critical value, say cr, 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 cris 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 c=cr. The concentration at the point joining the two models (cr) is problem-dependent and can be found in principle by solving the following equation:


For given values of aEand bE, Eq. (79) implicitly defines cras a function of Rep.

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

To avoid numerical noises occurring when c0, Rusche [34] suggests that the momentum equations (Eqs. (38) and (40)) should be converted into the following “phase-intensive” form by dividing ρf1cand ρsc:




The solutions of Eqs. (80) and (81) are expressed in the following semidiscretized forms:


where Aβ(β= sor f) denotes the systems of linear algebraic equations arising from the discretization of either Eqs. (82) or (83). The matrix Aβis decomposed into a diagonal matrix, ADβ, and an off-diagonal matrix, AOw. Also, AHw=bwAOβuβwith bβ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 ADβand AHβ, which depend on the discretization schemes. For example, Lee et al. [16] and Lee and Huang [35] 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 usand Eq. (39) to calculate c, then cmay increase rapidly toward cc, leading to an infinite psfor large c. This can be avoided by using a prediction-correction method to compute ufand us. This is achieved by splitting Eq. (83) into a predictor usand 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 c:


The right-hand side of Eq. (86) now has a diffusive term introduced by the numerical scheme. High sediment concentration and large psincrease the numerical diffusion (the right-hand side of Eq. (86)) and thus can avoid a rapid increase of cand the numerical instability due to high sediment concentration.

For the velocity-pressure coupling, Eq. (82) is similarly solved using a predictor ufand 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 [36]. The poor conservation of lighter material can be avoided by combining Eqs. (37) and (39) into the following Eq. (37):


and using Eq. (89) to correct pf. The method proposed [37] can help avoid the numerical instability. To show this, we follow Carver [37] and define


and combine Eqs. (83) and (88)(90) to obtain the following equation


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 kϵ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 c0, Eq. (83) becomes singular. To avoid this, 1/cis replaced by 1/c+δcin numerical computations, where δcis a very small number, say 106. When cδc, only a very small amount of solid particles are moving with the fluid; replacing 1/cby 1/c+δcmay introduce error in computing us; to avoid this error, we can set us=uf, 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 cusis very small when c106. Because the maximum value of cis 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,uf, ûs, and pfobtained at the previous time step, and it is outlined below:

  1. Solve Eqs. (80) and (81).

  2. Compute usfrom Eq. (84).

  3. Solve Eq. (86) for c.

  4. Compute usfrom Eq. (85).

  5. Compute ûsfrom Eq. (90).

  6. Compute uffrom Eq. (87).

  7. Solve Eq. (91) for pf.

  8. Repeat Eqs. (5)(7) for ntimes (say n= 1).

  9. Compute uffrom Eq. (88).

  10. Set us=uffor very dilute region, specifically c106.

  11. Repeat Eqs. (1)(10) with the updated c, uf, ûs, and pfuntil the residuals of Eqs. (80), (86), and (91) are smaller than the tolerance (say 105).

  12. Solve Eqs. (45) and (46) for kand ϵ, and compute the related coefficients.

Figure 1 is a flowchart showing these 12 solution steps.

Figure 1.

A flow chart showing the solution procedure using OpenFOAM®.

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 uf, Eq. (81) for us, Eq. (86) for c, Eq. (91) for pf, Eq. (45) for k, and (46) for ϵ; the convergence criteria are set at the residuals not exceeding 108. 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 us=ufis 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 Δtand the size of the cell where such particle is located, is defined as

CFL=absu jS j2VΔt,E92

where in u j=1cufj+cusj, the subscript “j” represents the jthface of the cell, S jis a unit normal vector, Vis the volume of the cell, and Δtis 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 ccan be avoided. Therefore, it is recommended that maxCFLc>co< 0.005. The time step is recommended to be in the range of 105and 104s.

5. 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].

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. [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=2650kg/m3 and d=0.76mm 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.66m; the sediment reservoir length is 2.1m; the overflow weir on the right end has a height of 0.239m; the upstream inflow discharge rate at the sluice opening is 0.204m2/s, which translates into an average horizontal flow velocity V=1.02m/s under the sluice gate. As an example, the computed development of scour depth dsis shown in Figure 3 together with the measurement of Chatterjee et al. [39].

Figure 2.

A sketch of the experimental setup for scour induced by a submerged water jet.

Figure 3.

Comparison of the computed scour depth with measurements 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 2dis 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. [40]. A 1:1 scale two-phase 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/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. [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.

Figure 4.

A sketch of the experimental setup for the collapse of a deeply submerged granular column.

Figure 5.

The simulated collapsing processes for the initially loose condition (a)–(d) and the initially dense condition (e)–(h). The lines represent contours of the computed concentrations, and the symbols were experimental data of Rondon et al. [40]. The figure is adapted from Lee and Huang [35].

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.

6. 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 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.


chapter PDF

© 2018 The Author(s). Licensee IntechOpen. This chapter is distributed under the terms of the Creative Commons Attribution 3.0 License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

How to cite and reference

Link to this chapter Copy to clipboard

Cite this chapter Copy to clipboard

Zhenhua Huang and Cheng-Hsien Lee (December 19th 2018). Modeling of Fluid-Solid Two-Phase Geophysical Flows [Online First], IntechOpen, DOI: 10.5772/intechopen.81449. Available from:

chapter statistics

238total chapter downloads

More statistics for editors and authors

Login to your personal dashboard for more detailed statistics on your publications.

Access personal reporting

We are IntechOpen, the world's leading publisher of Open Access books. Built by scientists, for scientists. Our readership spans scientists, professors, researchers, librarians, and students, as well as business professionals. We share our knowledge and peer-reveiwed research papers with libraries, scientific and engineering societies, and also work with corporate R&D departments and government entities.

More About Us