Modeling of Fluid-Solid Two-Phase Geophysical Flows

Abstract

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.

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=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]:

$$\frac{\partial {\rho }_k}{\partial t}+\nabla \cdot {\mathbf{u}}_k=0,$$ E1

and

$$\frac{\partial \rho {\mathbf{u}}_k}{\partial t}+\nabla \cdot \left({\rho }_k{\mathbf{u}}_k{\mathbf{u}}_k\right)=\nabla \cdot {\mathbf{T}}_k+\rho \mathbf{g},$$ E2

where ${\rho }_k$ is the density, ${\mathbf{u}}_k$ is the velocity, and $\mathbf{g}$ is the acceleration due to gravity. The stress tensor ${\mathbf{T}}_k$ includes two components:

$${\mathbf{T}}_k=-p_k\mathbf{I}+{\tau }_k$$ E3

where $p^k$ is the microscopic pressure and ${\tau }^k$ is the microscopic stress tensor.

Because the fluid phase and the solid phase are immiscible, at any time $t$ , a point in space $\mathbf{x}$ can be occupied only by one phase, not both. This fact can be described mathematically by the following phase function $c_k\left(\mathbf{x},t\right)$ for phase $k$ :

$$c_k\left(\mathbf{x},t\right)=\left\{\begin{array}{ll}1,& \text{if}\hspace{0.17em}\text{the}\hspace{0.17em}\text{point}\mathbf{x}\text{is}\hspace{0.17em}\text{occupied}\mathrm{by}\text{phase}k\\ 0,& \text{if}\hspace{0.17em}\text{the}\hspace{0.17em}\text{point}\mathbf{x}\text{is}\hspace{0.17em}\text{not}\hspace{0.17em}\text{occupied}\mathrm{by}\text{phase}k\end{array}.\right.$$ E4

The volumetric concentration of phase $k$ is directly related to the probability of occurrence of phase $k$ at a given location $\mathbf{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⟩$ . The operator $⟨\cdots ⟩$ 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

$$\frac{\partial c_k}{\partial t}+{\mathbf{u}}^i\cdot \nabla c_k=0,$$ E5

where ${\mathbf{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 $\nabla c_k$ is zero except at the interface between two phases where $\nabla c_k$ behaves like a delta-function [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:

$$c_k\frac{\partial {\varphi }_k}{\partial t}=\frac{\partial c_k{\varphi }_k}{\partial t}-{\varphi }_k\frac{\partial c_k}{\partial t}=\frac{\partial c_k{\varphi }_k}{\partial t}+{\varphi }_k{\mathbf{u}}^i\cdot \nabla c_k,$$ E6

and

$$c_k\nabla \cdot \left({\rho }_k{\mathbf{u}}_k\right)=\nabla \cdot \left(c_k{\rho }_k{\mathbf{u}}_k\right)-\left({\rho }_k{\mathbf{u}}_k\right)\cdot \nabla c_k.$$ E7

The resulting equations governing the ensemble average motion of phase $k$ are [8]

$$\frac{\partial ⟨c_k{\rho }_k⟩}{\partial t}+\nabla \cdot ⟨c_k{\rho }_k{\mathbf{u}}_k⟩=⟨{\rho }_k\left({\mathbf{u}}_k-{\mathbf{u}}^i\right)\cdot \nabla c_k⟩,$$ E8

and

$$\frac{\partial ⟨c_k{\rho }_k{\mathbf{u}}_k⟩}{\partial t}+\nabla \cdot ⟨c_k{\rho }_k{\mathbf{u}}_k{\mathbf{u}}_k⟩=\nabla \cdot ⟨c_k{\mathbf{T}}_k⟩+⟨c_k{\rho }_k\mathbf{g}⟩+{\tilde{\mathbf{m}}}_k$$ E9

with

$${\tilde{\mathbf{m}}}_k=⟨{\rho }_k{\mathbf{u}}_k\left({\mathbf{u}}_k-{\mathbf{u}}^i\right)-{\mathbf{T}}_k\cdot \nabla c_k⟩,$$ E10

Note that $\nabla 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, ${\mathbf{u}}_k-{\mathbf{u}}^i=0$ . As a result, the right-hand side of Eq. (8) is zero and

$${\tilde{\mathbf{m}}}_k=-⟨{\mathbf{T}}_k\cdot \nabla c_k⟩,$$ E11

which is the density of the interfacial force [8]. Physically, ${\mathbf{T}}^k\cdot \nabla c_k$ is the microscopic density of the force acting on a surface whose normal direction is defined by $\nabla c_k$ .

After using Eq. (3) for ${\mathbf{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

$$\frac{\partial \tilde{c}{\tilde{\rho }}_k}{\partial t}+\nabla \cdot \left[\tilde{c}{\tilde{\rho }}_k{\widehat{\mathbf{u}}}_k\right]=0$$ E12

and

$$\frac{\partial {\tilde{c}}_k{\tilde{\rho }}_k{\widehat{\mathbf{u}}}_k}{\partial t}+\nabla \cdot \left[\tilde{c}{\tilde{\rho }}_k{\widehat{\mathbf{u}}}_k{\widehat{\mathbf{u}}}_k\right]=\tilde{c}{\tilde{\rho }}_k\mathbf{g}+\nabla \cdot \left[\tilde{c}\left({\tilde{p}}_k\mathbf{I}+{\tilde{\tau }}_k\right)\right]+\nabla \cdot \left[\tilde{c}{\tilde{\tau }}_k^{\prime }\right]+{\tilde{\mathbf{m}}}_k,$$ E13

where ${\tilde{c}}_k=⟨c_k⟩$ 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 ${\tilde{\rho }}_k$ , pressure ${\tilde{p}}_k$ , stress tensor ${\tilde{\tau }}_k$ , and velocity ${\widehat{\mathbf{u}}}_k$ , defined by

$${\tilde{\rho }}_k=\frac{⟨c_k{\rho }_k⟩}{⟨c_k⟩},{\tilde{p}}_k=\frac{⟨c_k{p}_k⟩}{⟨c_k⟩},{\tilde{\tau }}_k=\frac{⟨c_k{\tau }_k⟩}{⟨c_k⟩},{\widehat{\mathbf{u}}}_k=\frac{⟨c_k{\rho }_k{\mathbf{u}}_k⟩}{⟨c_k{\rho }_k⟩}$$ E14

and ${\tilde{\mathit{t}}}_k^{\prime }$ represents the $c$ -weighted ensemble average of microscopic momentum flux associated with the fluctuation of the velocity ${\mathbf{u}}_k$ around the ensemble averaged velocity ${\widehat{\mathbf{u}}}_k$

$${\tilde{\tau }}_k^{\prime }=-\frac{⟨c_k{\rho }_k{\mathbf{u}}_k^{\prime }{\mathbf{u}}_k^{\prime }⟩}{⟨c_k⟩},{\mathbf{u}}_k^{\prime }={\mathbf{u}}_k-{\widehat{\mathbf{u}}}_k$$ E15

For compressible materials ${\tilde{\rho }}_k$ is not a constant. However, for incompressible materials

$${\tilde{\rho }}_k=\frac{⟨c_k{\rho }_k⟩}{⟨c_k⟩}\equiv {\rho }_k,{\widehat{\mathbf{u}}}_k=\frac{⟨c_k{\rho }_k{\mathbf{u}}_k⟩}{⟨c_k{\rho }_k⟩}=\frac{⟨c_k{\mathbf{u}}_k⟩}{⟨c_k⟩}$$ E16

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, and ${\tilde{\mathbf{m}}}_f+{\tilde{\mathbf{m}}}_s=0$ , the governing equations reduce to

$$0=\left(1-{\tilde{c}}_s\right){\tilde{\rho }}_f\mathbf{g}-\nabla \left[\left(1-{\tilde{c}}_s\right){\tilde{p}}_f\right]-{\tilde{\mathbf{m}}}_s,$$ E17

for the fluid phase, and

$$0={\tilde{c}}_s{\tilde{\rho }}_s\mathbf{g}-\nabla \left[{\tilde{c}}_s{\tilde{p}}_s\right]+{\tilde{\mathbf{m}}}_s,$$ E18

for the solid phase.

Because ${\tilde{p}}_f$ is the hydrostatic pressure in this case, i.e., $\nabla {\tilde{p}}_f={\tilde{\rho }}_f\mathbf{g}$ , it then follows that

$${\tilde{\mathbf{m}}}_s={\tilde{p}}_f\nabla {\tilde{c}}_s$$ E19

which, physically, is the buoyancy acting on the solid phase. Now Eq. (18) becomes

$$0={\tilde{c}}_s{\tilde{\rho }}_s\mathbf{g}-\nabla \left[{\tilde{c}}_s{\tilde{p}}_s\right]+{\tilde{p}}_f\nabla {\tilde{c}}_s$$ E20

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 ${\tilde{p}}_s={\tilde{p}}_f+{\stackrel{⌣}{p}}_s$ , with ${\tilde{p}}_f$ being the total fluid pressure and ${\stackrel{⌣}{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

$$\frac{\partial \left(1-c\right){\rho }_f}{\partial t}+\nabla \cdot \left[\left(1-c\right){\rho }_f{\mathbf{u}}_f\right]=0,$$ E21

and

$$\begin{array}{ll}\frac{\partial \left(1-c\right){\rho }_f{\mathbf{u}}_f}{\partial t}& +\nabla \cdot \left[\left(1-c\right){\rho }_f{\mathbf{u}}_f{\mathbf{u}}_f\right]\\ & =\left(1-c\right){\rho }_f\mathbf{g}+\nabla \cdot \left[\left(1-c\right)\left(-p_f\mathbf{I}+{\tau }_f\right)\right]+\nabla \cdot \left[\left(1-c\right){\tau }_f^{\prime }\right]-\mathbf{m}.\end{array}$$ E22

The ensemble averaged equations governing the motion of the solid phase are

$$\frac{\partial c{\rho }_s}{\partial t}+\nabla \cdot \left(c{\rho }_s{\mathbf{u}}_s\right)=0,$$ E23

and

$$\frac{\partial c{\rho }_s{\mathbf{u}}_s}{\partial t}+\nabla \cdot \left[c{\rho }_f{\mathbf{u}}_s{\mathbf{u}}_s\right]={\rho }_{s}c\mathbf{g}+\nabla \cdot \left[c\left(-p_f\mathbf{I}-p_s\mathbf{I}+{\tau }_s\right)\right]+\nabla \cdot \left(c{\tau }_s^{\prime }\right)+\mathbf{m}.$$ E24

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 ${\tau }_s^{\prime }$ , ${\tau }_f^{\prime }$ , ${\tau }_s$ , ${\tau }_f$ , $p_s$ , and $\mathbf{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⟩$ and ${\mathbf{u}}_{ f}$ associated with the turbulent flow [9] and (ii) including a term in the model for $\mathbf{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 $\mathbf{m}$ , the interphase force $\mathbf{m}$ should include the so-called general buoyancy $p_f\nabla c$ and a component $\mathbf{f}$ which includes drag force, inertial force, and lift force

$$\mathbf{m}=\mathbf{f}+p_f\nabla c\equiv \mathbf{f}-c\nabla p_f+\nabla \left({\mathit{cp}}_f\right).$$ E25

This expression for $\mathbf{m}$ has been derived by [11] using a control volume/surface approach. For most fluid-solid two-phase geophysical flows, the drag force dominates $\mathbf{f}$ [9] and thus $\mathbf{f}$ can be modeled by

$$\mathbf{f}=c{\rho }_s\frac{{\mathbf{u}}_f-{\mathbf{u}}_s}{{\tau }_p},$$ E26

where ${\tau }_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.

2.2 Favre averaged equations

The volumetric concentration and the velocities can be written as

$$c=\overline{c}+c^{\prime \prime },p_f={\overline{p}}_f+p_f^{\prime \prime },{\mathbf{u}}_f={\overline{\mathbf{u}}}_f+{\mathbf{u}}_f^{\prime \prime },{\mathbf{u}}_s={\overline{\mathbf{u}}}_s+{\mathbf{u}}_s^{\prime \prime },$$ E27

where the Favre averages are defined as

$$\begin{array}{ll}{\overline{\rho }}_s=\frac{\overline{c{\rho }_s}}{\overline{c}},& {\overline{\rho }}_f=\frac{\overline{\left(1-c\right){\rho }_f}}{1-\overline{c}},{\overline{\mathbf{u}}}_s=\frac{\overline{c{\rho }_s{\mathbf{u}}_s}}{\overline{c{\rho }_s}},{\overline{\mathbf{u}}}_f=\frac{\overline{\left(1-c\right){\rho }_f{\mathbf{u}}_f}}{\overline{\left(1-c\right){\rho }_f}},\end{array}$$ E28

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:

$$\frac{\partial {\overline{\rho }}_f\left(1-\overline{c}\right)}{\partial t}+\nabla \cdot \left[{\overline{\rho }}_f\left(1-\overline{c}\right){\overline{\mathbf{u}}}_f\right]=0,$$ E29

$$\begin{array}{ll}\frac{\partial {\overline{\rho }}_f\left(1-\overline{c}\right){\overline{\mathbf{u}}}_f}{\partial t}& +\nabla \cdot \left[{\overline{\rho }}_f\left(1-\overline{c}\right){\overline{\mathbf{u}}}_f{\overline{\mathbf{u}}}_f\right]={\overline{\rho }}_f\left(1-\overline{c}\right)\mathbf{g}-\left(1-\overline{c}\right)\nabla {\overline{p}}_f+\overline{c^{\prime \prime }\nabla p_f^{\prime \prime }}\\ & +\nabla \cdot \overline{\left(1-c\right)\left({\tau }_f+{\tau }_f^{\prime }+{\tau }_f^{\prime \prime }\right)}-\overline{c}{\overline{\rho }}_s\frac{{\overline{\mathbf{u}}}_f-{\overline{\mathbf{u}}}_s}{{\tau }_p}-\frac{{\overline{\rho }}_s}{{\tau }_p}\left(\overline{c{\mathbf{u}}_f^{\prime \prime }}\right),\end{array}$$ E30

for the fluid phase, with ${\tau }_f^{\prime \prime }$ being defined by

$${\tau }_f^{\prime \prime }=-{\rho }_f{\mathbf{u}}_f^{\prime \prime }{\mathbf{u}}_f^{\prime \prime },$$ E31

and

$$\frac{\partial {\overline{\rho }}_s\overline{c}}{\partial t}+\nabla \cdot \left[{\overline{\rho }}_s\overline{c}{\overline{\mathbf{u}}}_s\right]=0,$$ E32

$$\begin{array}{ll}\frac{\partial {\overline{\rho }}_s\overline{c}{\overline{\mathbf{u}}}_s}{\partial t}& +\nabla \cdot \left[{\overline{\rho }}_s\overline{c}{\overline{\mathbf{u}}}_s{\overline{\mathbf{u}}}_s\right]={\overline{\rho }}_s\overline{c}\mathbf{g}-\overline{c}\nabla {\overline{p}}_f-\overline{c^{\prime \prime }\nabla p_f^{\prime \prime }}-\nabla \overline{{\mathit{cp}}_s}\\ & +\nabla \cdot \overline{c\left({\tau }_s+{\tau }_s^{\prime }+{\tau }_s^{\prime \prime }\right)}+\overline{c}{\overline{\rho }}_s\frac{{\overline{\mathbf{u}}}_f-{\overline{\mathbf{u}}}_s}{{\tau }_p}+\frac{{\overline{\rho }}_s}{{\tau }_p}\left(\overline{c{\mathbf{u}}_f^{\prime \prime }}\right),\end{array}$$ E33

for the solid phase, with ${\tau }_s^{\prime \prime }$ being defined by

$${\tau }_s^{\prime \prime }=-{\rho }_s{\mathbf{u}}_s^{\prime \prime }{\mathbf{u}}_s^{\prime \prime }$$ E34

It is remarked here that the terms $\left(1-\tilde{c}\right)\nabla {\tilde{p}}_f$ in Eq. (30) and $\tilde{c}\nabla {\tilde{p}}_f$ in Eq. (33) have been obtained by using the expression for $\mathbf{m}$ given in Eq. (25).

In order to close these averaged equations, closure models are required for the following terms: $\overline{c\left({\tau }_s+{\tau }_s^{\prime }+{\tau }_s^{\prime \prime }\right)}$ , $\overline{\left(1-c\right)\left({\tau }_f+{\tau }_f^{\prime }+{\tau }_f^{\prime \prime }\right)}$ , $\overline{c{\mathbf{u}}_f^{\prime \prime }}$ , and $\overline{c^{\prime \prime }\nabla p_f^{\prime \prime }}$ . The last term can be neglected based on an analysis of their orders of magnitude by Drew [12]. The term $\overline{c{\mathbf{u}}_f^{\prime \prime }}$ is approximated by the following gradient transport hypotheses:

$$\overline{c{\mathbf{u}}_f^{\prime \prime }}=-\frac{{\nu }_{\mathit{\text{ft}}}}{{\sigma }_c}\nabla c$$ E35

where ${\nu }_{\mathit{\text{ft}}}$ is the eddy viscosity and ${\sigma }_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:

$$\overline{c\left({\tau }_s+{\tau }_s^{\prime }+{\tau }_s^{\prime \prime }\right)}=\overline{c}{\overline{\tau }}_s,\overline{\left(1-c\right)\left({\tau }_f+{\tau }_f^{\prime }+{\tau }_f^{\prime \prime }\right)}=\left(1-\overline{c}\right){\overline{\tau }}_f,\overline{{\mathit{cp}}_s}=\overline{c}{\overline{p}}_s$$ E36

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

$$\frac{\partial {\rho }_f\left(1-c\right)}{\partial t}+\nabla \cdot {\rho }_f\left(1-c\right){\mathbf{u}}_f=0,$$ E37

$$\begin{array}{ll}\frac{\partial {\rho }_f\left(1-c\right){\mathbf{u}}_f}{\partial t}& +\nabla \cdot \left[{\rho }_f\left(1-c\right){\mathbf{u}}_f{\mathbf{u}}_f\right]\\ & ={\rho }_f\left(1-c\right)\mathbf{g}-\left(1-c\right)\nabla p_f+\nabla \cdot \left(1-c\right){\tau }_f\\ & -\left\{c{\rho }_s\frac{{\mathbf{u}}_f-{\mathbf{u}}_s}{{\tau }_p}+\frac{{\rho }_s}{{\tau }_p}\frac{{\nu }_{\mathit{\text{ft}}}}{{\sigma }_c}\nabla c\right\},\end{array}$$ E38

for the fluid phase and

$$\frac{\partial {\rho }_{s}c}{\partial t}+\nabla \cdot {\rho }_{s}c{\mathbf{u}}_s=0,$$ E39

$$\begin{array}{ll}\frac{\partial {\rho }_{s}c{\mathbf{u}}_s}{\partial t}& +\nabla \cdot \left({\rho }_{s}c{\mathbf{u}}_s{\mathbf{u}}_s\right)\\ & ={\rho }_s\mathit{cg}-c\nabla p_f-\nabla \left({\mathit{cp}}_s\right)+\nabla \cdot c{\tau }_s\\ & +\left\{c{\rho }_s\frac{\left({\mathbf{u}}_f-{\mathbf{u}}_s\right)}{{\tau }_p}-\frac{{\rho }_s}{{\tau }_p}\frac{{\nu }_{\mathit{\text{ft}}}}{{\sigma }_c}\nabla c\right\},\end{array}$$ E40

for the solid phase.

3. Closure models

3.1 Stresses for the fluid phase

The stress tensor for the fluid phase ${\tau }_f$ includes two parts: a part for the viscous stress, ${\tau }_f^v$ , and the other part for the turbulent Reynolds stress, ${\tau }_f^t$

$${\tau }_f={\tau }_f^v+{\tau }_f^t$$ E41

The viscous stress tensor ${\tau }_f^v$ is usually computed by

$${\tau }_f^v=-{\rho }_f\left(\frac{2}{3}{\nu }_f\nabla \cdot {\mathbf{u}}_f\right)\mathbf{I}+2{\rho }_f{\nu }_f{\mathbf{D}}_f$$ E42

where ${\nu }_f$ is the kinematic viscosity of the fluid phase and ${\mathbf{D}}_f=\left[\nabla {\mathbf{u}}_f+{\left(\nabla {\mathbf{u}}_f\right)}^T\right]/2$ , where the superscript $T$ denotes a transpose. Some studies [13] suggested modifying ${\nu }_f$ to consider the effect of the solid phase; other studies [14], however, obtained satisfactory results even without considering this effect.

The stress tensor ${\tau }_f^t$ is related to the turbulent characteristics, which need to be provided by solving a turbulent closure model such as $k-\mathrm{ϵ}$ or $k-\omega$ model. For a $k-\mathrm{ϵ}$ model with low-Reynolds-number correction [15], ${\tau }_f^t$ can be computed by

$${\tau }_f^t=-{\rho }_f\left(\frac{2}{3}k+\frac{2}{3}{\nu }_f^t\nabla \cdot {\mathbf{u}}_f\right)\mathbf{I}+2{\rho }_f{\nu }_f^t{\mathbf{D}}_f$$ E43

where $k$ is the turbulence kinetic energy and ${\nu }_f^t$ is the eddy viscosity of the fluid phase, given by

$${\nu }_f^t=f_{\mu }{C}_{\mu }{k}^2/\mathrm{ϵ}$$ E44

with $\mathrm{ϵ}$ being the turbulent dissipation of the fluid phase to be provide by solving the $k-\mathrm{ϵ}$ equation. The coefficient $f_{\mu }=\text{exp}\left[-3.4/{\left(1+{\text{Re}}_t/50\right)}^2\right]$ represents the low-Reynolds-number correction with ${\text{Re}}_t=k^2/{\nu }_f\mathrm{ϵ}$ . The coefficient $C_{\mu }$ is usually assumed to be a constant.

The equations governing $k$ and $\mathrm{ϵ}$ are similar to those for clear water [15]

$$\begin{array}{ll}& \frac{\partial {\rho }_f\left(1-c\right)k}{\partial t}+\nabla \cdot \left[{\rho }_f\left(1-c\right){\mathbf{u}}_{f}k\right]=\left(1-c\right){\mathit{t}}_f:\nabla {\mathbf{u}}_f-{\rho }_f\left(1-c\right)\mathrm{ϵ}\\ & +\nabla \cdot \left[{\rho }_f\frac{{\nu }_f^t}{{\sigma }_c}\left(1-c\right)k\right]-\left\{\left({\rho }_s-{\rho }_f\right)\frac{{\nu }_f^t}{{\sigma }_c}\nabla c\cdot \mathbf{g}+\frac{2{\rho }_{s}c\left(1-\alpha \right)k}{{\tau }_p}\right\},\end{array}$$ E45

and

$$\begin{array}{ll}& \frac{\partial {\rho }_f\left(1-c\right)\mathrm{ϵ}}{\partial t}+\nabla \left[{\rho }_f\left(1-c\right){\mathbf{u}}_f\mathrm{ϵ}\right]=\frac{\mathrm{ϵ}}{k}\left[C_{\mathrm{ϵ}1} f_1\left(1-c\right){\tau }_f:\nabla {\mathbf{u}}_f-C_{\mathrm{ϵ}2} f_2{\rho }_f\left(1-c\right)\mathrm{ϵ}\right]\\ & +\nabla \cdot \left[{\rho }_f\frac{{\nu }_f^t}{{\sigma }_{\mathrm{ϵ}}}\left(1-c\right)\mathrm{ϵ}\right]-\frac{\mathrm{ϵ}}{k}{C}_{\mathrm{ϵ}3}\left\{\left({\rho }_s-{\rho }_f\right)\frac{{\nu }_f^f}{{\sigma }_c}\nabla c\cdot \mathbf{g}+\frac{2{\rho }_{s}c\left(1-\alpha \right)k}{{\tau }_p}\right\},\end{array}$$ E46

where coefficients $C_{\mathrm{ϵ}1}$ , $C_{\mathrm{ϵ}2}$ , ${\sigma }_{\mathrm{ϵ}}$ , ${\sigma }_k$ , and $f_2$ are model parameters, whose values can be taken the same as those in the $k-\mathrm{ϵ}$ 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_{\mathrm{ϵ}3}=1$ is usually adopted in the literature [28]; however, it is remarked that the value of $C_{\mathrm{ϵ}3}$ is not well understood at the present and a sensitivity test to understand how the value of this $C_{\mathrm{ϵ}3}$ on the simulation results is recommended. The parameter $\alpha$ reflects the correlation between the solid-phase and fluid-phase turbulent motions and is given by

$$\alpha ={\left(1+\frac{{\tau }_p}{\text{min}\left({\tau }_l,{\tau }_c\right)}\right)}^{-1},$$ E47

where ${\tau }_l=0.165k/\mathrm{ϵ}$ is a time scale for the turbulent flow and ${\tau }_c$ is a time scale for particle collisions given by [16]

$${\tau }_c=\left[{\left(\frac{c_{\mathrm{rcp}}}{c}\right)}^{\frac{1}{3}}-1\right]d{\left(\frac{{\rho }_s}{p_s}\right)}^{1/2}$$ E48

with $c_{\mathrm{rcp}}$ being the random close packing fraction and $d$ being the particle diameter. $c_{\mathrm{rcp}}$ is 0.634 for spheres [17]. The term ${\left(c_{\mathrm{rcp}}/c\right)}^{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-\mathrm{ϵ}$ 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 $p_s$ and ${\tau }_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$ :

$$p_s=p_s^t+p_s^r+p_s^e,$$ E49

where $p_s^t$ accounts for the turbulent motion of solid particles (important for dilute flows); $p_s^r$ reflects the rheological characteristics of dense flows and includes the effects such as fluid viscosity, enduring contact, and particle inertial; $p_s^e$ 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_s^e$

$$\begin{array}{l}{p}_s^e=K{\left[\text{max}\left(c-c_o,0\right)\right]}^{\chi }\left\{1+\text{sin}\left[\text{max}\left(\frac{c-c_o}{c_{\mathrm{rcp}}-c_o},0\right)\pi -\frac{\pi }{2}\right]\right\},\end{array}$$ E50

where $c_o$ is random loose packing fraction and coefficients $K$ and $\chi$ 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_s^r$ and $p_s^t$ 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

$${\mathit{t}}_s=-\left(\frac{2}{3}{\rho }_s{\nu }_s\nabla \cdot {\mathbf{u}}_s\right)+2{\rho }_s{\nu }_s{\mathbf{D}}_s,$$ E51

The kinematic viscosity of the solid phase ${\nu }_s$ is computed by the sum of two terms:

$${\nu }_s={\nu }_s^v+{\nu }_s^t,$$ E52

where ${\nu }_s^v$ and ${\nu }_s^t$ represent the visco-plastic and turbulence effects, respectively. This model for ${\nu }_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_s^t$ and ${\nu }_s^t$ can be computed by

$$p_s^t=\frac{2}{3}{\rho }_s\mathit{\alpha k},$$ E53

and

$${\nu }_s^t={\mathit{\alpha \nu }}_f^t.$$ E54

where the coefficient $\alpha$ 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+{\mathit{aI}}_i^2$ , 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{\rho }_f{\nu }_f{D}^s/{\mathit{cp}}_s$ where ${\nu }_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 $I_i=2{\mathit{dD}}^s/\sqrt{{\mathit{cp}}_s/{\rho }_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 ${\mathrm{st}}_v=I_i^2/I_{\nu }$ . Some formulas have been proposed in the literature to describe $c-I$ and $\eta -I$ relationships, where $\eta =T^s/p_s$ with $T^s$ being the second invariant of ${\tau }_s$ .

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

$$c=\frac{c_c}{1+{\mathit{bI}}^{1/2}}$$ E55

where $c_c$ is a critical concentration and $b$ is a model parameter. Trulsson et al. [21] proposed

$$\eta ={\eta }_1+\frac{{\eta }_2-{\eta }_1}{1+I_o/I^{1/2}},$$ E56

where ${\eta }_1=\text{tan}{\theta }_s$ with ${\theta }_s$ being the angle of repose and ${\eta }_2$ and $I_o$ are constants. Based on Eqs. (56) and (55), the following expressions for $p_s^r$ and ${\nu }_s^v$ can be derived [16]:

$${\nu }_s^v=\frac{\left(p_s^r+p_s^e\right)\eta }{2{\rho }_s{D}_s},$$ E57

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

$$p_s^r=\frac{2b^{2}c}{{\left(c_c-c\right)}^2}\left({\rho }_f{\nu }_f+2a{\rho }_s{d}^2{D}_s\right)D_s,$$ E58

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

3.3 Closure models for particle response time

The drag force between the two phases is modeled through the particle response time ${\tau }_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

$$-{\rho }_f\left(1-c\right)g-\left(1-c\right)\frac{\partial p_f}{\partial z}-\frac{c{\rho }_s\left(w_f-w_s\right)}{{\tau }_p}=0,$$ E59

and

$$-{\rho }_s\mathit{cg}-c\frac{\partial p_f}{\partial z}+\frac{c{\rho }_s\left(w_f-w_s\right)}{{\tau }_p}=0,$$ E60

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

$$\left(1-c\right)w_f+{\mathit{cw}}_s=0.$$ E61

Combining Eqs. (59) and (61) yields

$$-\frac{\partial p_f}{\partial z}=\frac{-c{\rho }_s{w}_s}{{\left(1-c\right)}^2{\tau }_p}+{\rho }_{f}g.$$ E62

Substituting Eqs. (61) and (62) into Eq. (60) leads to

$${\tau }_p=\frac{{\rho }_s{w}_s}{{\left(1-c\right)}^2\left({\rho }_s-{\rho }_f\right)g},$$ E63

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

$$\frac{w_s}{w_0}={\left(1-c\right)}^n,$$ E64

where the coefficient $n$ is related to the particle Reynolds number ${\mathit{\text{Re}}}_s=w_{0}d/{\nu }_f$

$$n=\left\{\begin{array}{ll}4.65,& {\mathit{\text{Re}}}_s<0.2\\ 4.4{\mathit{\text{Re}}}_s^{-0.33},& 0.2\le {\mathit{\text{Re}}}_s<1\\ 4.4{\mathit{\text{Re}}}_s^{-0.1},& 1\le {\mathit{\text{Re}}}_s<500\\ 2.4,& 500\le {\mathit{\text{Re}}}_s\end{array}\right..$$ E65

The terminal velocity of a single particle $w_0$ can be computed by

$$w_0=\sqrt{\frac{4\mathit{dg}}{3C_d}\frac{{\rho }_s-{\rho }_f}{{\rho }_f}},$$ E66

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:

$$C_d=\frac{24}{{\mathit{\text{Re}}}_p}+\frac{6}{1+\sqrt{{\mathit{\text{Re}}}_p}}+0.4,$$ E67

where ${\mathit{Re}}_p=\mid {\mathbf{u}}_f-{\mathbf{u}}_s\mid d/{\nu }_f$ . Combing Eqs. (63)–(67) yields

$${\tau }_p=\frac{{\rho }_s}{{\rho }_f}\frac{d^2}{{\nu }_f}\frac{{\left(1-c\right)}^{n-2}}{18+\left(4.5/\left(1+\sqrt{{\mathit{\text{Re}}}_p}\right)+0.3\right){\mathit{\text{Re}}}_p}.$$ E68

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.

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

$$-\frac{\partial p_f}{\partial x}=\frac{c{\rho }_s{u}_f}{\left(1-c\right){\tau }_p},$$ E69

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

$$-\frac{\partial p_f}{\partial x}=a_F{\rho }_f\left(1-c\right)u_f+b_F{\rho }_f{\left(1-c\right)}^2{u}_f^2,$$ E70

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:

$$a_F=\frac{a_E{c}^3{\nu }_f}{{\left(1-c\right)}^2{d}^2},b_F=\frac{b_{E}c}{g{\left(1-c\right)}^{3}d},$$ E71

Comparing Eqs. (69) and (70) and using Eq.(71) give

$${\tau }_p=\frac{{\rho }_s{d}^2}{{\rho }_f{\nu }_f}\frac{1}{a_E{c}^2+b_E{\mathit{Re}}_p},$$ E72

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\approx 2\times {10}^{-4}$ m, $k_p\approx {10}^{-10}\sim {10}^{-11}{\mathrm{m}}^2$ [30], which gives $a_E\approx 1.6\times {10}^3\sim 1.6\times {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 ${\kappa }_p$ . According to Darcy’s law for seepage [29], the pressure gradient can also be written as

$$-\frac{\partial p_f}{\partial x}=\frac{{\rho }_f{\nu }_f\left(1-c\right)u_f}{{\kappa }_p},$$ E73

where ${\kappa }_p$ is the permeability. Combining Eqs. (69) and (73) gives

$${\tau }_p=\frac{c{\rho }_s{\kappa }_p}{{\left(1-c\right)}^2{\rho }_f{\nu }_f}$$ E74

When the flow is very slow, Eqs. (70), (71), and (73) suggest that

$$a_E=\frac{d^2}{k_p{\left(1-c\right)}^2},$$ E75

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

$$\frac{w}{w_s}={\left(1-c\right)}^{n-1}{\left[\text{max}\left(1-c/c_m,0\right)\right]}^{c_m},$$ E76

where $c_m$ is the maximum concentration at which $w=0$ . In this study, $c_m=c_o$ is adopted because when $c\ge c_o$ , contact networks can form in the granular material.

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

$${\tau }_p=\frac{{\rho }_s}{{\rho }_f}\frac{d^2}{{\nu }_f}\frac{{\left(1-c\right)}^{n-3}{\left[\text{max}\left(1-c/c_m,0\right)\right]}^{c_m}}{18+\left(4.5/\left(1+\sqrt{{\mathit{\text{Re}}}_p}\right)+0.3\right){\mathit{\text{Re}}}_p}.$$ E77

We stress that $c=c_m$ will lead to ${\tau }_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:

$${\tau }_p=\left\{\begin{array}{ll}\frac{{\rho }_s}{{\rho }_f}\frac{d^2}{{\nu }_f}\frac{{\left(1-c\right)}^{n-3}{\left[\text{max}\left(1-c/c_m,0\right)\right]}^{c_m}}{18+\left(4.5/\left(1+\sqrt{{\mathit{\text{Re}}}_p}\right)+0.3\right){\mathit{\text{Re}}}_p},& \text{for}c<c_r\\ \frac{{\rho }_s{d}^2}{{\rho }_f{\nu }_f}\frac{1}{a_E{c}^2+b_E{\mathit{\text{Re}}}_p},& \text{for}c\ge c_r\end{array}\right.$$ E78

where $c_r$ 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 $c=c_r$ . The concentration at the point joining the two models ( $c_r$ ) is problem-dependent and can be found in principle by solving the following equation:

$$\begin{array}{l}\frac{{\left(1-c_r\right)}^{n-3}{\left[\text{max}\left(1-c_r/c_m,0\right)\right]}^{c_m}}{18+\left(4.5/\left(1+\sqrt{{\mathit{\text{Re}}}_p}\right)+0.3\right){\mathit{\text{Re}}}_p}=\frac{1}{a_E{c}_r^2+b_E{\mathit{\text{Re}}}_p}.\end{array}$$ E79

For given values of $a_E$ and $b_E$ , Eq. (79) implicitly defines $c_r$ as a function of ${\mathit{Re}}_p$ .

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 $c\to 0$ , Rusche [34] suggests that the momentum equations (Eqs. (38) and (40)) should be converted into the following “phase-intensive” form by dividing ${\rho }_f\left(1-c\right)$ and ${\rho }_{s}c$ :

$$\begin{array}{ll}\frac{\partial {\mathbf{u}}_f}{\partial t}+\nabla \cdot \left({\mathbf{u}}_f{\mathbf{u}}_f\right)& -\left(\nabla \cdot {\mathbf{u}}_f\right){\mathbf{u}}_f=\mathbf{g}-\frac{1}{{\rho }_f}\nabla p_f+\frac{1}{{\rho }_f}\nabla \cdot {\tau }_f-\frac{{\tau }_f\cdot \nabla c}{{\rho }_f\left(1-c\right)}\\ & -\frac{c{\rho }_s}{{\rho }_f\left(1-c\right)}\frac{\left({\mathbf{u}}_f-{\mathbf{u}}_s\right)}{{\tau }_p}+\frac{{\rho }_s}{{\rho }_f\left(1-c\right){\tau }_p}\frac{{\nu }_f^t}{{\sigma }_c}\nabla c\end{array}$$ E80

and

$$\begin{array}{ll}\frac{\partial {\mathbf{u}}_s}{\partial t}+\nabla \cdot \left({\mathbf{u}}_s{\mathbf{u}}_s\right)& -\left(\nabla \cdot {\mathbf{u}}_s\right){\mathbf{u}}_s=\mathbf{g}-\frac{1}{{\rho }_s}\nabla p_f-\frac{1}{{\rho }_{s}c}\nabla {\mathit{cp}}_s\\ & +\frac{1}{{\rho }_s}\nabla \cdot {\tau }_s+\frac{{\tau }_s\cdot \nabla c}{{\rho }_{s}c}+\frac{\left({\mathbf{u}}_f-{\mathbf{u}}_s\right)}{{\tau }_p}-\frac{1}{c{\tau }_p}\frac{{\nu }_f^t}{{\sigma }_c}\nabla c\end{array}$$ E81

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

$${\mathbf{u}}_f=\frac{{\mathbf{A}}_H^f}{{\mathbf{A}}_D^f}+\frac{\mathbf{g}}{A_D^f}-\frac{\nabla p_f}{{\rho }_f{\mathbf{A}}_D^f}+\frac{{\rho }_{s}c{\mathbf{u}}^s}{{\rho }_f{\mathbf{A}}_D^f\left(1-c\right){\tau }_p}+\frac{{\rho }_s}{{\rho }_f{\mathbf{A}}_D^f\left(1-c\right){\tau }_p}\frac{{\nu }_f^t}{{\sigma }_c}\nabla c$$ E82

$${\mathbf{u}}_s=\frac{{\mathbf{A}}_H^s}{{\mathbf{A}}_D^s}+\frac{\mathbf{g}}{{\mathbf{A}}_D^s}-\frac{\nabla p_f}{{\rho }_s{A}_D^s}-\frac{\nabla p_s}{{\rho }_s{\mathbf{A}}_D^s}-\frac{p_s\nabla c}{{\rho }_s{\mathbf{A}}_D^{s}c}+\frac{{\rho }_s{\mathbf{u}}_f}{{\mathbf{A}}_D^s{\tau }_p}-\frac{1}{{\mathbf{A}}_D^{s}c{\tau }_p}\frac{{\nu }_f^t}{{\sigma }_c}\nabla c$$ E83

where ${\mathbf{A}}^{\beta }$ ( $\beta$ = $s$ or $f$ ) denotes the systems of linear algebraic equations arising from the discretization of either Eqs. (82) or (83). The matrix ${\mathbf{A}}^{\beta }$ is decomposed into a diagonal matrix, ${\mathbf{A}}_D^{\beta }$ , and an off-diagonal matrix, ${\mathbf{A}}_O^w$ . Also, ${\mathbf{A}}_H^w={\mathbf{b}}^w-{\mathbf{A}}_O^{\beta }{\mathbf{u}}^{\beta }$ with ${\mathbf{b}}^{\beta }$ 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 ${\mathbf{A}}_D^{\beta }$ and ${\mathbf{A}}_H^{\beta }$ , 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 ${\mathbf{u}}_s$ and Eq. (39) to calculate $c$ , then $c$ may increase rapidly toward $c_c$ , leading to an infinite $p_s$ for large $c$ . This can be avoided by using a prediction-correction method to compute ${\mathbf{u}}_f$ and ${\mathbf{u}}_s$ . This is achieved by splitting Eq. (83) into a predictor ${\mathbf{u}}_s^{\ast }$ and a corrector. The predictor is

$${\mathbf{u}}_s^{\ast }=\frac{{\mathbf{A}}_H^s}{{\mathbf{A}}_D^s}+\frac{\mathbf{g}}{{\mathbf{A}}_D^s}-\frac{\nabla p_f}{{\rho }_s{A}_D^s}+\frac{{\rho }_s{\mathbf{u}}_f}{{\mathbf{A}}_D^s{\tau }_p}$$ E84

which is corrected by the following corrector

$${\mathbf{u}}_s={\mathbf{u}}_s^{\ast }-\left(\frac{p_s\nabla c}{{\rho }_s{\mathbf{A}}_D^{s}c}-\frac{1}{{\mathbf{A}}_D^{s}c{\tau }_p}\frac{{\nu }_f^t}{{\sigma }_c}\nabla c\right)$$ E85

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$ :

$$\frac{\partial c}{\partial t}+\nabla \cdot \left(c{\mathbf{u}}_s^{\prime }\right)=\nabla \cdot \left(\frac{p_s}{{\rho }_s{\mathbf{A}}_D^s}+\frac{1}{{\mathbf{A}}_D^s{\tau }_p}\frac{{\nu }_f^t}{{\sigma }_c}\right)\nabla c$$ E86

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

For the velocity-pressure coupling, Eq. (82) is similarly solved using a predictor ${\mathbf{u}}_f^{\ast }$ and a corrector. The predictor is

$${\mathbf{u}}_f^{\ast }=\frac{{\mathbf{A}}_H^f}{{\mathbf{A}}_D^f}+\frac{\mathbf{g}}{A_D^f}+\frac{{\rho }_{s}c{\mathbf{u}}_s}{{\rho }_f{\mathbf{A}}_D^f\left(1-c\right){\tau }_p}+\frac{{\rho }_s}{{\rho }_f{\mathbf{A}}_D^f\left(1-c\right){\tau }_p}\frac{{\nu }_f^t}{{\sigma }_c}\nabla c$$ E87

which is corrected by the following corrector

$${\mathbf{u}}_f={\mathbf{u}}_f^{\ast }-\frac{\nabla p_f}{{\rho }_f{\mathbf{A}}_D^f}$$ E88

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):

$$\nabla \cdot \left[\left(1-c\right){\mathbf{u}}_f+c{\mathbf{u}}_s\right]=0$$ E89

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

$${\widehat{\mathbf{u}}}_s=\frac{{\mathbf{A}}_H^s}{{\mathbf{A}}_D^s}+\frac{\mathbf{g}}{{\mathbf{A}}_D^s}-\frac{\nabla p_s}{{\rho }_s{\mathbf{A}}_D^s}-\frac{p_s\nabla c}{{\rho }_s{\mathbf{A}}_D^{s}c}+\frac{{\rho }_s{\mathbf{u}}_f}{{\mathbf{A}}_D^s{\tau }_p}-\frac{1}{{\mathbf{A}}_D^{s}c{\tau }_p}\frac{{\nu }_f^t}{{\sigma }_c}\nabla c$$ E90

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

$$\nabla \cdot \left[\left(1-c\right){\widehat{\mathbf{u}}}_f+c{\widehat{\mathbf{u}}}_s\right]=\nabla \cdot \left[\frac{1-c}{{\rho }_f{\mathbf{A}}_D^f}+\frac{c}{{\rho }_s{\mathbf{A}}_D^s}\right]\nabla p_f$$ E91

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-\mathrm{ϵ}$ 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 $c\to 0$ , Eq. (83) becomes singular. To avoid this, $1/c$ is replaced by $1/\left(c+\mathit{\delta c}\right)$ in numerical computations, where $\mathit{\delta c}$ is a very small number, say ${10}^{-6}$ . When $c\le \mathit{\delta c}$ , only a very small amount of solid particles are moving with the fluid; replacing $1/c$ by $1/\left(c+\mathit{\delta c}\right)$ may introduce error in computing ${\mathbf{u}}^s$ ; to avoid this error, we can set ${\mathbf{u}}_s={\mathbf{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 $c{\mathbf{u}}_s$ is very small when $c\le {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$ , ${\mathbf{u}}_f$ , ${\widehat{\mathbf{u}}}_s$ , and $p_f$ obtained at the previous time step, and it is outlined below:

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

2. Compute ${\mathbf{u}}_s^{\ast }$ from Eq. (84).

3. Solve Eq. (86) for $c$ .

4. Compute ${\mathbf{u}}_s$ from Eq. (85).

5. Compute ${\widehat{\mathbf{u}}}_s$ from Eq. (90).

6. Compute ${\mathbf{u}}_f^{\ast }$ from Eq. (87).

7. Solve Eq. (91) for $p_f$ .

8. Repeat Eqs. (5)–(7) for $n$ times (say $n$ = 1).

9. Compute ${\mathbf{u}}_f$ from Eq. (88).

10. Set ${\mathbf{u}}_s={\mathbf{u}}_f$ for very dilute region, specifically $c\le {10}^{-6}$ .

11. Repeat Eqs. (1)–(10) with the updated $c$ , ${\mathbf{u}}_f$ , ${\widehat{\mathbf{u}}}_s$ , and $p_f$ until the residuals of Eqs. (80), (86), and (91) are smaller than the tolerance (say ${10}^{-5}$ ).

12. Solve Eqs. (45) and (46) for $k$ and $\mathrm{ϵ}$ , 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 ${\mathbf{u}}_f$ , Eq. (81) for ${\mathbf{u}}_s$ , Eq. (86) for $c$ , Eq. (91) for $p_f$ , Eq. (45) for $k$ , and (46) for $\mathrm{ϵ}$ ; the convergence criteria are set at the residuals not exceeding ${10}^{-8}$ . 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 ${\mathbf{u}}_s={\mathbf{u}}_f$ 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 $\mathrm{\Delta }t$ and the size of the cell where such particle is located, is defined as

$$\mathrm{CFL}=\sum \frac{\mathrm{abs}\left({\mathbf{u}}^{ }\cdot {\mathbf{S}}^{ }\right)}{2V}\mathrm{\Delta }t,$$ E92

where in ${\mathbf{u}}^{ }=\left(1-c\right){\mathbf{u}}_f^j+c{\mathbf{u}}_s^j$ , the subscript “ $j$ ” represents the $j^{\mathit{th}}$ face of the cell, ${\mathbf{S}}^{ }$ is a unit normal vector, $V$ is the volume of the cell, and $\mathrm{\Delta }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 $\text{max}\left({\left.\mathrm{CFL}\right|}_{c>c_o}\right)$ < 0.005. The time step is recommended to be in the range of ${10}^{-5}$ and ${10}^{-4}$ s.

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 ${\rho }_s=2650$ kg/m³ 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²/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].

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

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/m³, respectively. The density and the dynamic viscosity of the liquid are 1010 kg/m³ 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-\mathrm{ϵ}$ 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.

Acknowledgments

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.