Algorithm

The numerical implementation of the present Eulerian two-phase flow sediment transport model is based on an open-source finite volume CFD library called OpenFOAM. Taking advantage of the numerical discretization schemes and framework of finite volume solvers in OpenFOAM, the two-phase flow governing equations are implemented by modifying the solver twoPhaseEulerFoam (Rusche, 2002; Weller, 2002; Peltola, 2009). OpenFOAM uses the Finite Volume Method (FVM) over a collocated grid arrangement. The collocated arrangement stores all dependent variables at the cell center and the same Control Volume (CV) is used for all variables to minimize the computational effort. The advantage of the FVM is that the system of partial differential equations can be discretized on arbitrary three-dimensional structured or unstructured meshes. Thus, complex geometries can be easily handled. The Gauss theorem is applied to the convection and diffusion terms leading to conservative schemes.

To illustrate the numerical discretization, the fluid phase momentum equation are taken as an example. Rearranging the fluid phase momentum equation by dividing $\beta\rho^b$ , the resulting equation can be written as

\[ \dfrac{\partial{\mathbf{u}^b}}{\partial{t}}+\nabla\cdot(\mathbf{u^b}\mathbf{u^b})-(\nabla\cdot\mathbf{u^b})\mathbf{u^b}=-\dfrac{1}{\rho^b}\nabla{p^*} - \dfrac{\alpha K}{\rho^b}(\mathbf{u^b}-\mathbf{u^a}) + \dfrac{K}{\rho^b} \dfrac{1}{\sigma_c} \ \nu^{bt}\nabla{\alpha}+ \dfrac{\mathbf{f}}{\rho^b}+ \dfrac{1}{\beta}\nabla\cdot\mathbf{\tau^b} \tag{1} \]

The last term in the above equation, the gradient of fluid phase shear stress, can be written and expanded as follows:

\[ \dfrac{1}{\beta}\nabla\cdot\mathbf{\tau^b}= \nabla\cdot\left( \nu^{b}_{Eff} \nabla\mathbf{u}^{b}\right) + \nu^{b}_{Eff} \dfrac{\nabla\beta}{\beta} \nabla\mathbf{u}^{b} + \dfrac{1}{\beta}\nabla\cdot\left\{\beta \ \nu^{b}_{Eff} \left[ \nabla\mathbf{u}^{b T} -\dfrac{2}{3} \nabla\cdot \mathbf{u}^b \right]\right\} \tag{2} \]

In the above equation, the first two terms on the RHS are treated implicitly while the last two terms are treated explicitly. By substituting the expanded shear stress formulation in the momentum equation, the following equation is obtained:

\[\dfrac{1}{\beta}\nabla\cdot\mathbf{\tau^b}= \nabla\cdot\left( \nu^{b}_{Eff} \nabla\mathbf{u}^{b}\right) + \nu^{b}_{Eff} \dfrac{\nabla\beta}{\beta} \nabla\mathbf{u}^{b} + \dfrac{1}{\beta}\nabla\cdot\left\{\beta \ \nu^{b}_{Eff} \left[ \nabla\mathbf{u}^{b T} -\dfrac{2}{3} \nabla\cdot \mathbf{u}^b \right]\right\} \tag{3} \]

\[ \dfrac{\partial{\mathbf{u}^b}}{\partial{t}}+\nabla\cdot(\mathbf{u^b}\mathbf{u^b}) -(\nabla\cdot\mathbf{u^b})\mathbf{u^b} - \nabla\cdot\left( \nu^{b}_{Eff} \nabla\mathbf{u}^{b}\right) - \nu^{b}_{Eff} \dfrac{\nabla\beta}{\beta} \nabla\mathbf{u}^{b} + \dfrac{\alpha K}{\rho^b} \mathbf{u^b} = -\dfrac{1}{\rho^b}\nabla{p^*} + \dfrac{\alpha K}{\rho^b}\mathbf{u^a} + \dfrac{1}{\sigma_c} \ \dfrac{K \nu^{bt} }{\rho^b}\nabla{\alpha} + \dfrac{\mathbf{f}}{\rho^b} + \dfrac{1}{\beta}\nabla\cdot\left\{\beta \ \nu^{b}_{Eff} \left[ \nabla\mathbf{u}^{b T} -\dfrac{2}{3} \nabla\cdot \mathbf{u}^b \right]\right\} \tag{4} \]

It is more convenient to rewrite the above equation into a matrix form:

\[ \big[\mathbf{\sf A^b}\big] \cdot\mathbf{u^b}= \mathbf{H^b}+\mathbf{R^b}+ \dfrac{\alpha K}{\rho^b}\mathbf{u^a} -\dfrac{1}{\rho^b}\nabla{p^*}\tag{5}\]

The matrix $\big[\mathbf{\sf A^b}\big]$ is composed of the diagonal terms of the algebraic system associated with Eq5 , whereas $\mathbf{H^b}$ includes the off-diagonal terms and the source terms. $\mathbf{R^b} $ is composed of the explicit drag term, the turbulent suspension term, the gravity term and the explicit diffusion terms:

\[ \mathbf{R^b}= \dfrac{\mathbf{f}}{\rho^b} + \dfrac{1}{\beta}\nabla\cdot\left\{\beta \ \nu^{b}_{Eff} \left[ \nabla\mathbf{u}^{b T} -\dfrac{2}{3} \nabla\cdot \mathbf{u}^b \right]\right\} + \dfrac{1}{\sigma_c} \dfrac{K \ \nu^{bt}}{\rho^b} \nabla{\alpha}\label{eq_mom_fluid_7} \tag{6} \]

The same process can be carried out for the solid phase momentum equation that leads to:

\[ \dfrac{\partial{\mathbf{u}^a}}{\partial{t}} +\nabla\cdot(\mathbf{u^a}\mathbf{u^a}) -(\nabla\cdot\mathbf{u^a})\mathbf{u^a} - \dfrac{1}{\tilde{\alpha}} \nabla\cdot\left( \nu^{a}_{Fr} \nabla\mathbf{u}^{a}\right) - \nabla\cdot\left( \nu^{a}_{Eff} \nabla\mathbf{u}^{a}\right) - \nu^{a}_{Eff} \dfrac{\nabla\alpha}{\tilde{\alpha}} \nabla\mathbf{u}^{a} + \dfrac{\beta K}{\rho^a} \mathbf{u^a} = - \dfrac{1}{\tilde{\alpha} \rho^a} \nabla{p^{ff}} \displaystyle -\dfrac{1}{\rho^a}\nabla{p^*} + \dfrac{\beta K}{\rho^a}\mathbf{u^b} -\dfrac{1}{\sigma_c} \dfrac{\beta K \ \nu^{bt}}{\tilde{\alpha} \rho^a} \nabla{\alpha}+ \mathbf{g} \left(1-\dfrac{\rho^b}{\rho^a}\right)+ \dfrac{\mathbf{f}}{\rho^a} - \dfrac{1}{\tilde{\alpha} \rho^a} \nabla{p^{a}} + \dfrac{1}{\tilde{\alpha}} \left\{\nabla\cdot\left[(\alpha \nu^{a}_{Eff}+\nu^a_{Fr}) \nabla\mathbf{u}^{a\ T}\right]\right\} +\dfrac{1}{\tilde{\alpha}} \left\{\nabla\left[\left(\lambda-\dfrac{2}{3}( \alpha \nu^{a}_{Eff}+\nu^a_{Fr} ) \right)\nabla\cdot{\mathbf{u}^a}\right]\right\} \tag{7} \]

where $\alpha$ at the denominator is substituted by $\tilde{\alpha}=\alpha+\alpha_{Small}$ to avoid dividing by zero when the solid phase volume concentration vanishes. This partial differential system of equations can also be written as a matrix equation as follows:

\[ \big[\mathbf{\sf A^a}\big]\cdot\mathbf{u^a}= \mathbf{H^a}+\mathbf{R^a} + \dfrac{\beta K}{\rho^a}\mathbf{u^b}-\dfrac{1}{\rho^a}\nabla{p^*} \tag{8} \]

where the source term $\mathbf{R^a}$ contains the following terms:

\[ \mathbf{R^a}= -\dfrac{1}{\sigma_c} \dfrac{\beta K \ \nu^{bt}}{\tilde{\alpha} \rho^a} \nabla{\alpha} + \mathbf{g} \left(1-\dfrac{\rho^b}{\rho^a}\right)+ \dfrac{\mathbf{f}}{\rho^a}- \dfrac{1}{\tilde{\alpha} \rho^a} \nabla{p^{a}} + \dfrac{1}{\tilde{\alpha}}\left\{\nabla\cdot\left[(\alpha \nu^{a}_{Eff}+\nu^a_{Fr}) \nabla\mathbf{u}^{a\ T}\right]\right\} +\dfrac{1}{\tilde{\alpha}}\left\{\nabla\left[\left(\lambda-\dfrac{2}{3}( \alpha \nu^{a}_{Eff}+\nu^a_{Fr} ) \right)\nabla\cdot{\mathbf{u}^a}\right]\right\}\label{eq_mom_sed5} \tag{9} \]

Following Rusche (2002) the terms involving the ratio of particle phase volume gradient to the volume concentration are treated at the cell face level in the predictor-corrector algorithm. However, we noted the exception of the particle phase normal stress $p^{ff}$ gradient for which a reconstruction of the surface normal gradient at the cell centre allows to get more stable solutions.

In sedFoam, two options are proposed to treat the explicit shear stress terms (i) they can be solved at the cell faces and added in pEqn using the keyword 'faceMomentum' = true (the solver is using the code given in the pUf folder) or (ii) they can be solved at the cell centers and added in UEqn using the keyword 'faceMomentum' = false (the solver is using the code given in the pU folder). By default, the explicit shear stress terms are solved at the cell centers and the keyword 'faceMomentum' = false.

The advantage of separating the right hand side of the momentum equations as the sum of two terms: $\mathbf{R}$ and $\mathbf{H}$ , is for writing the pressure-velocity algorithm. A similar method as Rhie and Chow (1983) one can be applied for the gradient terms. The details of the velocity-pressure algorithm are presented in the next subsection.

Velocity-pressure algorithm

The pressure-velocity coupling and the consequent oscillations in the pressure fields are resolved using the Rhie and Chow method. The PISO algorithm is used to solve fluid and particle velocities (Rusche2002,Weller2002,Peltola2009).

Step 1: Velocity predictor at the cell centres

First, the intermediate velocities $\mathbf{u^{a*}}$ , $\mathbf{u^{b*}}$ (denoted as HabyA and HbbyA) are computed using the corresponding momentum equations (equations without the pressure gradient term and the explicit part of the drag term):

\[ \mathbf{u^{a*}}= \big[\mathbf{\sf A^a}\big]^{-1}\mathbf{H^a} = HabyA \tag{10} \]
\[ \mathbf{u^{b*}}= \big[\mathbf{\sf A^b}\big]^{-1}\mathbf{H^b} = HbbyA \tag{11} \]

where $\big[\mathbf{ A^a}\big]^{-1}=rUAa$ and $\big[\mathbf{ A^b}\big]^{-1}=rUAb$ represent the inverse matrices of $\big[\mathbf{A^a}\big]$ and $\big[\mathbf{A^b}\big]$ , respectively. These intermediate velocities do not satisfy the mass conservation equations.

Step 2: Velocity fluxes at the cell faces

The right hand side of equations (Eq5 - Eq8 ) are computed at the cell faces:

\[ \Phi_R^a = \Big[A^a_f\Big]^{-1}\mathbf{R^a_f}\quad \text{and} \quad \Phi_R^b = \Big[A^b_f\Big]^{-1} \mathbf{R^b_f} \tag{12} \]

and the velocity flux associated with the predictor step is interpolated at the cell faces:

\[ \Phi_U^{a} = \left(\Big[A^a\Big]^{-1}\mathbf{H^a}\right)_f \quad \text{and} \quad \Phi_U^{b} = \left(\Big[A^b\Big]^{-1} \mathbf{H^b}\right)_f \tag{13} \]

Introducing the notation:

\[ AK_a = \left( \Big[A^a\Big]^{-1} \beta K / \rho_a \right)_f \quad \text{and} \quad AK_b = \left( \Big[A^b\Big]^{-1} \alpha K / \rho_b \right)_f \tag{14} \]

one can write the volume averaged velocity flux at the cell faces as:

\[ \Phi^{*} = \alpha_f \left[\Phi_U^{a} +\Phi_R^a + AK_a \Phi^b\right] +\beta_f \left[ \Phi_U^{b} +\Phi_R^b + AK_b \Phi^a \right] \tag{15} \]

where $\Phi^a=\mathbf{u^a}\vert_f.\mathbf{n}\vert_f S_f$ and $\Phi^b=\mathbf{u^b}\vert_f.\mathbf{n}\vert_f S_f$ denote the fluid and particle phases velocity fluxes at the previous iteration or time step and at the cell faces, respectively, and $S_f$ is the cell face area associated with face $f$ . At this stage, the drag force is partially treated explicitly for the mixture flux $\Phi^*$

Step 3: Velocity corrector at the cell faces

Writing the semi-discrete velocity equation including the pressure gradient term leads to the following equations:

\[ \Phi^{a**} = \Phi_U^{a} +\Phi_R^{a}+\left(\dfrac{\beta K}{\rho^a \ \big[A^a\big]}\right)_f \Phi^{b**} - \dfrac{\mathbf{\nabla^{\perp}}p^*\vert_f}{\rho^a \big[A^a\big]_f} \tag{16} \]

\[ \Phi^{b**} = \Phi_U^{b} +\Phi_R^{b} +\left(\dfrac{\alpha K}{\rho^b \ \big[A^b\big]}\right)_f \Phi^{a**} - \dfrac{\mathbf{\nabla^{\perp}}p^*\vert_f}{\rho^b \big[A^b\big]_f} \tag{17} \]

Taking the divergence of the volume averaged mixture velocity given by the velocity correction equation and imposing the incompressibility constraint, $\nabla\cdot\mathbf{U^{**}}=\nabla\cdot(\alpha \mathbf{u^{a**}} + \beta \mathbf{u^{b**}})=0$ , one can built the pressure equation as a function of the predicted velocity or predicted face fluxes:

\[ \int_{Vp} \nabla\cdot \left[\left(\dfrac{\alpha}{\rho^a \big[A^a\big]} + \dfrac{\beta}{\rho^b \big[A^b\big]}\right)\nabla{p^*}\right] dV= \int_{Vp}\nabla\cdot\mathbf{U^{*}} dV \quad \\ \text{or} \quad \oint_{S} \left(\dfrac{\alpha_f}{\rho^a \big[A^a\big]_f} + \dfrac{\beta_f}{\rho^b \big[A^b\big]_f}\right)\nabla^\perp{p^*}\vert_f \ \mathbf{n}\vert_f \ dS= \oint_S \mathbf{U^{*}}\vert_f\ \mathbf{n}\vert_f \ dS \tag{18}\]

\[ \int_{Vp} \nabla\cdot\left[\alpha\mathbf{u^{a**}} + \beta\mathbf{u^{b**}}\right] dV= \oint_S \left[\alpha_f\mathbf{u^{a**}}\vert_f + \beta_f\mathbf{u^{b**}}\vert_f\right]\cdot \mathbf{n}\ dS =0 \tag{19} \]

where the subscript $f$ denotes variables interpolated at the cell faces. The two expressions shown above are equivalent by using the Gauss theorem. At the discrete level, this equation is written as:

\[ \sum_f\big[\alpha_f\Phi^{a**}+\beta_f{\Phi^{b**}}\big] = 0 \tag{20} \]

substituting the velocity corrector equations ( Eq16 - Eq17 ) into the previous equation, the Poisson equation for the pressure reads:

\[ \sum_f \left(\dfrac{\alpha_f}{\rho^a \big[A^a\big]_f} + \dfrac{\beta_f}{\rho^b \big[A^b\big]_f}\right) \nabla^\perp{p^*}\vert_f \ \mathbf{n}\vert_f \ S_f= \sum_f \Phi^{*} \tag{21} \]

This equation leads to a matrix system written at the cell faces. The resulting algebraic system is usually solved using an multigrid solver (GAMG). The resulting pressure field $p^*$ is used for the correction step in which the mixture velocity face flux is corrected using equations ( Eq16 - Eq17 ):

\[ \Phi^{**} = \Phi^* - \sum_f \mathbf{\nabla^{\perp}}p^*\vert_f S_f \tag{22} \]

And the fluid and particle phase face fluxes are corrected accordingly

\[ \Phi_s^a = \Phi_U^a+\Phi_R^a - \dfrac{\mathbf{\nabla^{\perp}}p^*\vert_f}{\rho^a \big[A^a\big]_f} \tag{23} \]

\[ \Phi_s^b = \Phi_U^b+\Phi_R^b - \dfrac{\mathbf{\nabla^{\perp}}p^*\vert_f}{\rho^b \big[A^b\big]_f}\tag{24} \]

in which the explicit drag contributions coming from the other phase are still missing. We can rewrite the corrected fluxes for each phase as:

\[ \Phi^{a**} = \Phi_s^a +AK_a \Phi^{b**} \tag{25} \]
\[ \Phi^{b**} = \Phi_s^b +AK_b \Phi^{a**} \tag{26} \]

and the face flux associated with the relative velocity as:

\[ \Phi^{r**} =\Phi^{a**} -\Phi^{b**}= \Phi_s^{a} + AK_a (\Phi_s^b + AK_b \Phi^{a**}) - \left[ \Phi_s^{b} + AK_b (\Phi_s^a + AK_a \Phi^{b**}) \right] \tag{27} \]
\[ \Phi^{r**} = \Phi_s^{a} + AK_a \Phi_s^b - \left[ \Phi_s^{b} + AK_b \Phi_s^a \right] + AK_a AK_b \left[ \Phi^{a**} - \Phi^{b**} \right]\tag{28} \]
\[ (1-AK_a AK_b) \Phi^{r**} = \Phi_s^{a} + AK_a \Phi_s^b - \left[ \Phi_s^{b} + AK_b \Phi_s^a \right] \tag{29} \]
\[ \Phi^{r**} = \dfrac{1}{1-AK_a AK_b} \left\{\Phi_s^{a} + AK_a \Phi_s^b - \left[ \Phi_s^{b} + AK_b \Phi_s^a \right] \right\}\tag{30} \]

This new expression for the relative flux allows to treat the drag force almost implicitly. This is supposed to stabilize the coupling between the two momentum equations and allows for higher time steps with regard to the particle response time.

The predicted fluxes for each phase is then obtained from the mixture flux plus a correction coming from the relative flux. Let's write the mixture velocity and relative velocity expression:

\[ \Phi^{**}=\alpha_f \Phi^{a**} + \beta_f \Phi^{b**} \quad \text{and} \quad \Phi^{r**} = \Phi^{a**} - \Phi^{b**} \tag{31} \]

from which we can deduce the following relationships:

\[ \Phi^{a}= \Phi^{**} + \beta_f \Phi^{r**} \quad \text{and} \quad \Phi^{b}= \Phi^{**} - \alpha_f \Phi^{r**} \tag{32} \]

The same type of procedure is applied for the reconstruction of the velocities at the cell centres: $U,U^a,U^b$ .

In order to ensure the mass conservation an iterative procedure of N cycles is sometimes required. From our experience, three iterations (N=3) is usually enough for a convergence. The finite volume discretisation of the equations have not been shown here but all the details can be found in Jasak (1996) and Rusche (2002).

Summary of the solution procedure

The numerical solution procedure for the proposed two-phase flow model is outlined as follow:

  • Solve for sediment concentration $\alpha$ ;
  • Update the volume concentration of fluid: $\beta=1-\alpha$ ;
  • Update the drag parameter $K$ in the drag term;
  • Solve for the fluid turbulence closure, update $k$ , $\varepsilon$ or $\omega$ (depends on the turbulence closure $k-\varepsilon$ or $k-\omega$ ), and then calculate the eddy viscosity and effective fluid total viscosity;
  • Solve for the particle phase stress (kinetic theory model or the dense granular rheology);

  • PISO-loop, solving velocity-pressure coupling for $N$ loops:

    a. Construct the coefficient matrix $ \big[\mathbf{\sf A^a}\big]$ and $\big[\mathbf{\sf A^b}\big]$ and explicit array $\mathbf{H^a}$ and $\mathbf{H^b}$ using Eq5 and Eq8 .

    b. update the other explicit source terms $\mathbf{R^a}$ and $\mathbf{R^b}$ , Eq6 and Eq9 ;

    c. Calculate $\mathbf{u^{a*}}$ , $\mathbf{u^{b*}}$ using Eq. ustara and ustarb without fluid pressure gradient term;

    d. Construct and solve the pressure Eq. Eq21 ;

    e. Correct fluid and particle velocities after solving pressure and update fluxes ( Eq23 - Eq24 );

    f. go to (a-e) if the number of loops is smaller than $N$ (no tolerance criteria).

  • Advance to the next time step

In the above solution procedure, the velocity-pressure coupling steps are looped for $N$ times. The advantage of this loop is to avoid velocity-pressure decoupling caused by the direct solving method. From our numerical practices, the loop number $N=1$ to 3 is usually enough to give reasonably accurate results, and shows good convergence, especially for steady flows.

The time step, $\Delta t$ , can be set to a constant value or adjusted automatically based on two Courant numbers, one related to the local flow velocity and the local grid size $Co=1/2 \sum_f \Phi_f \Delta t / V_p$ (the same as for single phase problems) and one related to the relative velocity $Co_r=1/2 \sum_f \left|\Phi^a_f-\Phi^b_f\right| \Delta t / V_p$ which is specific to the coupling of the fluid and sediment phase momentum equations in the two-phase flow model. The most limiting time step is used as the criterion for setting the adjustable time step. Our practice is to set these two Courant numbers to 0.3 and 0.1, respectively.