Governing equations

The mathematical formulation of the Eulerian two-phase flow model sedFoam is obtained by averaging local and instantaneous mass and momentum conservation equations over fluid and dispersed particles. Different averaging operators can be used, ensemble averaging (Drew, 1983) or spatial averaging (Jackson, 2000), and provided that the mathematical derivation is done properly, the different approaches should lead to the same conservation equations (Zhang and Prosperetti, 1997; Jackson, 1997). The resulting governing equations can be considered as the counterpart of the clear fluid Navier-Stokes equations for single phase flow. In order to apply these equations to turbulent flow, in which turbulent emeotions are generated by flow shear much larger than the grain scale, additional turbulence averaging or filtering is required. In the present model, turbulence-averaged Eulerian two-phase flow equations are derived by following a similar procedure presented in Hsu et al. (2003); Hsu and Liu (2004).

Turbulence-averaged two-phase flow governing equations

The mass conservation equations for the particle phase and fluid phase are written as:

\[ \frac{\partial {\alpha}}{\partial {t}} + \frac{\partial{{\alpha}u^a_i}}{\partial{x_i}}=0, \tag{1} \]
\[ \frac{\partial \beta}{\partial {t}} + \frac{\partial{\beta u^b_i}}{\partial{x_i}}=0, \tag{2} \]

where $\alpha$ and $\beta=1-\alpha$ are the particle and fluid volume concentrations, $u^a_i, u^b_i$ are the sediment and fluid phase velocities, and $i=1, 2, 3$ represents streamwise, spanwise and vertical component, respectively. The momentum equations for fluid and particle phases can be written as:

\[ \frac{\partial {\rho^a}{\alpha}u^a_i}{\partial {t}} + \frac{\partial{{\rho^a}{\alpha}u^a_iu^a_j}}{\partial{x_j}}=-{\alpha}\frac{\partial{p}}{\partial{x_i}}+{\alpha}f_i-\frac{\partial{\tilde{p}^a}}{\partial{x_i}}+\frac{\partial{\tau^a_{ij}}}{\partial{x_j}} + {\alpha}{\rho^a}g_i  + {\alpha}f_i + {\alpha}{\beta}K(u^b_i-u^a_i)- S_{US} \ {\beta}K\nu_t^{b}\frac{\partial{\alpha}}{\partial{x_i}}, \tag{3} \]

\[ \frac{\partial {\rho^b}{\beta}u^b_i}{\partial {t}} + \frac{\partial{{\rho^b}{\beta}u^b_iu^b_j}}{\partial{x_j}}=-\beta\frac{\partial{p}}{\partial{x_i}}+\beta f_i+\frac{\partial{\tau^b_{ij}}}{\partial{x_j}} + {\beta}{\rho^b}g_i +{\beta}f_i-{\alpha}{\beta}K(u^b_i-u^a_i)+ S_{US} \  {\beta}K{\nu_t^{b}}\frac{\partial{\alpha}}{\partial{x_i}}, \tag{4} \]

where $\rho^a, \rho^b$ are the particle and fluid density, respectively, $g_i$ is the gravitational acceleration, $f_i$ is an external volume force and $p$ is the fluid pressure. $f_i$ is the external force that drives the flow. The fluid stress $\tau^b_{ij}$ includes fluid grain-scale (viscous) stress and fluid Reynolds stresses (see section Fluid phase shear stress) and   $\tilde{p}^a, \tau^a_{ij}$ are the particle normal stress and shear stress (see section Particle phase stress). The last two terms on the right-hand-side (RHS) of  eq3  and eq4 are the momentum coupling between the fluid phase and particle phase through the drag force, where $K$ is the drag parameter. In particular the second to the last term represents the averaged drag force due to mean relative velocity between fluid and particle phases, while the last term represents the fluid turbulent suspension term, also called drift velocity by Simonin (1991). This term is due to the correlation of sediment concentration and fluid velocity fluctuations and the gradient transport assumption is adopted here for its closure. Hence, $\nu_t^{b}$ is the turbulent viscosity to be calculated using a turbulence closure, and $S_{US} = 1/\sigma_c$ is the inverse of the the Schmidt number. This term is equivalent to the turbulent suspension flux of the Rouse profile in the two-phase flow formalism (see  Chauchat (2018) appendix 1 for a detailed demonstration). Other forces such as the lift force or the added mass force could play a role in sediment transport, according to Jha and Bombardelli (2010), the lift force in dilute suspended sediment transport only represent 4% of the drag force and the added mass force can be of order of 10%  in the near bed region.

The drag parameter $K$ , is modeled following Schiller and Naumann (1933):

\[ K=0.75 C_d \frac{\rho^b}{d_{eff}} \parallel \mathbf{u^b}-\mathbf{u^a}\parallel\beta^{-h_{Exp}} \tag{5} \]

where $d_{eff}= {\psi}d$ is the effective sediment diameter, in which $\psi$ is the shape factor and $d$ is the particle diameter. The hindrance function $\beta^{-h_{Exp}}$ represents the drag increase when the particle volume concentration increases. $h_{Exp}$ is the hindrance exponent that depends on the particulate Reynolds number (Di Felice (1994)). For simplicity, the value of $h_{Exp}$ is assumed to be constant (default value is 2.65), and its value can be specified from the constant/transportProperties file in SedFoam-2312. This hypothesis is valid for particulate Reynolds numbers lower than unity or larger than 300, within this range the exponent value decreases down to $h_{Exp}\approx2$ . The drag coefficient $C_d$ is calculated as:

\[ C_d=\left\{ \begin{array}{r} \frac{24}{Re_p}(1+0.15 Re_p^{0.687}), \quad \quad Re_p \leq 1000 \\ 0.44, \quad  \quad Re_p>1000 \end{array}, \right. \tag{6} \]

in which, the particulate Reynolds number $Re_p$ is defined as: $Re_p=\beta\parallel\mathbf{u^b-u^a}\parallel d_{eff}/\nu^b$ , where $\nu^b$ stands for the fluid kinematic viscosity. This drag model can be chosen by the keyword `‘GidaspowSchillerNaumann’' in the file constant/interfacialProperties and it is especially well adapted for dealing with suspended particles. For situations in which the fluid flow inside the porous sediment bed has to be accurately solved, other drag models are available in SedFoam-2312. The major issue in developing an Eulerian two-phase flow model is to provide closure laws for turbulence closures and granular stress models. This will be extensively discussed in the following subsections, particularly different modeling options available in SedFoam-2312  are presented.

Fluid phase shear stress

Because the present model equations are obtained by averaging over turbulence, the fluid stresses consist of a large-scale component $R_{ij}^{bt}$ (i.e., Reynolds stress) and a grain-scale stress $r_{ij}^b$ , which includes the viscous stress and an additional effect due to fluid-particle interaction at grain scale. The total fluid stress is written as:

\[ \tau^b_{ij}=R_{ij}^{bt} +r_{ij}^b=\rho^b\beta\Big[2\nu_{Eff}^{b} \ S_{ij}^b-\frac{2}{3}k\delta_{ij}\Big], \tag{7} \]

where, $\nu_{Eff}^{b}=\nu_t^{b}+\nu^{mix}$ is the fluid phase effective viscosity with $\nu_t^{b}$ being the eddy viscosity, and $\nu^{mix}$ is the mixture viscosity. $S_{ij}^b$ is the deviatoric part of the fluid phase strain rate tensor defined as:

\[ S_{ij}^b = \frac{1}{2}\big(\frac{\partial{u^b_i}}{\partial{x_j}} + \frac{\partial{u^b_j}}{\partial{x_i}}\big)-\frac{1}{3}\frac{\partial{u^b_k}}{\partial{x_k}}\delta_{ij}, \tag{8} \]

The Reynolds stress tensor $R_{ij}^{bt}$ is modeled as:

\[ R_{ij}^{bt} =\rho^b\beta\Big[2 \nu_t^{b} \ S_{ij}^b  -\frac{2}{3}k\delta_{ij}\Big], \tag{9} \]

and the viscous stress tensor is modeled as:

\[ r_{ij}^{b} =2\rho^b\beta\nu^{mix} \ S_{ij}^b, \tag{10} \]

In SedFoam-2312, several different viscosity or turbulence closures are implemented, and these models can be selected according to specific flow conditions ranging from laminar to turbulent flows. The mixture viscosity can be selected in combination with a granular rheology model for the granular stresses (see section Dense granular rheology). For the turbulent eddy viscosity, modified turbulence closures are implemented for sediment transport applications.

Particle phase stress

In sediment transport applications, the particle stresses are important mechanisms to support the particle immersed weight in concentrated regions of sediment transport  (Hsu et al., 2004; Cheng et al., 2017). In these regions, the momentum exchanges due to particle collisions and/or enduring contacts exert dispersive stresses on a collection of particles. The particle phase stress tensor can be split into the normal and off-diagonal components that correspond to the particle pressure $\tilde{p}^a$ and the particle shear stress $\tilde{\tau}^a_{ij}$ . Though, details vary with the model selection, the particle normal stresses (or pressure) can be generally classified into two contributions: a shear induced or collisional component (super-script ‘a’) and a permanent contact component (super-script ‘ff’) (Johnson and Jackson, 1987):

\[ \tilde{p}^a=p^{ff}+p^{a}, \tag{11} \]

The first term of the particle pressure is due to enduring contact in highly concentrated region, where the sediment bed is quasi-static/immobile. This normal pressure increases rapidly when the sediment concentration is close to the maximum packing limit, and prevents unphysical sediment concentration in the sediment bed. Thus, this element is important to model a full transport profile including the quasi-static sediment bed. The permanent contact component $p^{ff}$ is calculated as:

\[ p^{ff}=\left\{ \begin{array}{r} 0, \alpha < \alpha_{min}^{Fric} \\ Fr\frac{(\alpha-\alpha_{min}^{Fric})^{\eta_0}}{(\alpha_{max}-\alpha)^{\eta_1}}, \alpha \geq \alpha_{min}^{Fric}, \tag{12} \end{array} \right\} \]

where $\alpha_{min}^{Fric}=0.57$ , $\alpha_{max}=0.635$ for spheres and $Fr$ , $\eta_0$ and $\eta_1$ are empirical coefficients. Following Cheng et al.(2017), the values are set to: $Fr=0.05, \eta_0=3$ and $\eta_1=5$ . Again, these coefficients can be reset in file constant/kineticTheoryProperties or  *constant/granularRheologyProperties* by specifying the keywords and values.

In modern sediment transport modeling framework, the two major threads of modeling approach for shear-induced/collisional particle normal stress and shear stress are the kinetic theory of granular flows and the dense granular flow rheology. They are implemented in this version of SedFoam-2312 (see subsections Kinetic theory and Dense granular rheology), respectively. As a result of different closures for particle normal stresses $p^a$ , the closure for the total particle shear stress $\tilde{\tau}^a_{ij}$ also varies, and they are described in the next two subsections.

 Dense granular rheology

Concerning the shear-induced contribution to the particle pressure, two different expressions can be adopted depending on the flow regime. Depending on the local Stokes and particulate Reynolds numbers, the regime of the granular flow rheology can change from free fall or grain inertia regime to viscous and turbulent regimes (Andreotti et al., 2013), and the definition of the controlling parameter ‘ $I$ ’ are different accordingly. In SedFoam-2312, the grain inertia and viscous regimes have been implemented and can be selected using the keywords ‘MuI’ and ‘MuIv’ respectively. Following  the work of Boyer et al. (2011), the shear-induced pressure under the viscous regime can be written as:

\[ p^{a}= \Big(\frac{B_{\phi} \alpha}{\alpha_{max}-\alpha}\Big)^2 \rho^a d^2 \parallel\mathbf{\sf S^a}\parallel, \tag{13} \]

where $B_{\phi}$    is   a   parameter   of   the   dilatancy   law (Boyer et al., 2011).

While in the grain inertia regime:

\[ p^{a}= \Big(\frac{B_{\phi} \alpha}{\alpha_{max}-\alpha}\Big)^2 \nu^a  \parallel\mathbf{\sf S^a}\parallel^2, \tag{14} \]

On the other hand, the particle phase shear stress is written as:

\[ \tilde{\tau}^a_{ij}=\mu({\rm I}) \ \tilde{p}^a \ \frac{S^a_{ij}}{\sqrt{2 \ S^a_{ij}{\cdot}S^a_{ij}}}, \tag{15} \]

where $S_{ij}^a$ is the deviatoric part of sediment phase strain rate tensor:

\[ S_{ij}^{a}=\frac{1}{2}\Big(\frac{\partial{u^a_i}}{\partial{x_j}}+\frac{\partial{u^a_j}}{\partial{x_i}}\Big)-\frac{1}{3}\dfrac{\partial{u^a_k}}{\partial{x_k}}\delta_{ij}, \tag{16} \]

The dynamic friction coefficient $\mu$ depends on the dimensionless controlling number ‘ $I$ ’. To be consistent with the definitions used in the kinetic theory, we still introduce the particle shear viscosity $\mu^a$ and frictional shear viscosity $\nu_{Fr}^a$ . However, we simply set $\mu^a=0$ , and  define the frictional shear viscosity only as  (Chauchat and Médale, 2014):

\[ \nu_{Fr}^a=\frac{\mu(I) \ \tilde{p}^a}{\rho^a \ \left(\parallel\mathbf{\sf S^a}\parallel^2+D_{small}^2\right)^{1/2}}, \tag{17} \]

where $\parallel\mathbf{\sf S^a}\parallel$ is the norm of the shear rate tensor (eq16 ). $D_{small}$ is a regularization parameter, which is introduced to avoid singularity. It is set to be $D_{small}=10^{-6} s^{-1}$ for all the simulations of granular rheology presented herein except stated otherwise.

Additionally,   dilatancy effects and pore pressure feedback on the granular medium are taken into account (if the user turns on this feature) modifying the elastic pressure (eq12) such as:

\[ p^{ff}=\left\{ \begin{array}{r} 0, \alpha < \alpha_{pl} \\ Fr\frac{(\alpha-\alpha_{pl})^{\eta_0}}{(\alpha_{max}-\alpha)^{\eta_1}}, \alpha \geq \alpha_{pl}, \tag{18} \end{array} \right\} \]

where the evolution of the plastic volumetric fraction ( $\alpha_{pl}$ ) is governed by the following equation:

\[ \frac{\partial \alpha_{pl} }{\partial t} + \mathbf{u^a} \cdot \nabla \alpha_{pl} = - \alpha_{pl} \delta \parallel\mathbf{\sf S^a}\parallel, \tag{19} \]

More details on the dilatancy model are described in Montellà et al. (2021) .

 Kinetic theory

For the kinetic theory, following Gidaspow (1994), the particle collisional stress is calculated as:

\[ \tau_{ij}^{a}=2\mu^a \ S_{ij}^{a}+ \lambda \frac{\partial{u^a_k}}{\partial{x_k}}\delta_{ij}.  \tag{20} \]

For more details, the interested reader is referred to Chauchat et al. (2017) .