Air Entrainment Modeling in the SPH Method: A Two-Phase Mixture Formulation with Open Boundaries

Abstract

Air entrainment within water is a common feature of flows over hydraulic works – spill over a dam, wave breaking on a dike, etc. – and its accurate modeling is a key to better design such structures. The Smoothed Particle Hydrodynamics (SPH) method appears as a natural way to model such highly distorted flows. To avoid computationally prohibitive costs related to the full discretization of bubbles or drops, a mixture model for high density ratio flows relying on a volume-based formulation with relative velocity between phases has first been developed and validated in Fonty et al. (Proceedings of 13th international SPHERIC workshop. Galway, Ireland, 2018; Int J Multiph Flow 111:158–174, 2019. https://doi.org/10.1016/j.ijmultiphaseflow.2018.11.007). Instead of having a once and for all assigned phase as in multifluid SPH, each particle now carries both phases through their respective volume fractions. In the present work, in order to handle practical air entrainment application cases, the open boundary formulation described in Ferrand et al. (Comput Phys Commun 210:29–44, 2017. https://doi.org/10.1016/j.cpc.2016.09.009) is adapted to this mixture model. Then, after introducing turbulence through a $k{-}\epsilon$ model, a specific closure on the air bubbles relative velocity is proposed including a Stokesian drag term and turbulent diffusion. This model is then applied to two cases of air entrainment: a stepped spillway for interfacial aeration and a plunging jet for local aeration. Finally a 3D industrial test case of discharge-control structure at the La Coche power plant (France) is considered. While valuable insights are obtained for the volume fraction field, further investigations are required to improve the modeling of the flow dynamics.

Appendices

Appendices

1.1 Time Integration of the Continuity Equation

The following derivation follows faithfully (Ferrand et al. 2017), but replacing the density $\rho$ by the inverse of the volume noted $\sigma$, making some adjustments when needed. From the continuity Eq. (5), one has:

$$\begin{aligned} \displaystyle {\frac{d \sigma _{a}}{d t}}=-\,\sigma _{a}D_{a}^{\gamma }\{\varvec{j}\} \end{aligned}$$

(54)

with:

$$\begin{aligned} D_{a}^{\gamma }\{\varvec{j}\}=-\,\frac{1}{\gamma _{a}\sigma _{a}}\sum _{b\in \left( \mathcal {F}\cup \mathcal {V}\right) }\theta _{b}\left( \varvec{j}_{a}-\varvec{j}_{b}\right) \cdot \varvec{\nabla }w_{ab}+\frac{1}{\gamma _{a}\sigma _{a}}\sum _{s\in \mathcal {S}}\sigma _{s}\left( \varvec{j}_{a}-\varvec{j}_{s}\right) \cdot \varvec{\nabla }\gamma _{as} \end{aligned}$$

(55)

Let us make the distinction between the Eulerian fluid velocity $\varvec{J}_{b}$ and the Lagrangian particle velocity $\varvec{j}_{b}$. Additional terms then appear in the divergence operator computation and lead to:

$$\begin{aligned} \displaystyle {\frac{d \sigma _{a}}{d t}}=\displaystyle {\frac{1}{\gamma _{a}}}\sum _{b\in \left( \mathcal {F}\cup \mathcal {V}\right) }\theta _{b}\left( \varvec{j}_{a}-\varvec{j}_{b}\right) \cdot \varvec{\nabla }w_{ab}-\delta \sigma _{a}^{i/o}-\displaystyle {\frac{1}{\gamma _{a}}}\sum _{s}\sigma _{s}\left( \varvec{j}_{a}-\varvec{j}_{s}\right) \cdot \varvec{\nabla }\gamma _{as}+\frac{\sigma _{a}}{\gamma _{a}}\delta \gamma _{a}^{i/o} \end{aligned}$$

(56)

With the definitions:

$$\begin{aligned} \delta \sigma _{a}^{i/o}&= \frac{1}{\gamma _{a}}\sum _{v\in \mathcal {V}^{i/o}}\theta _{v}\left( \varvec{J}_{v}-\varvec{j}_{v}\right) \cdot \varvec{\nabla }w_{av} \end{aligned}$$

(57)

$$\begin{aligned} \delta \gamma _{a}^{i/o}&= \sum _{s\in \mathcal {S}^{i/o}}\frac{\sigma _{s}}{\sigma _{a}}\left( \varvec{J}_{s}-\varvec{j}_{s}\right) \cdot \varvec{\nabla }\gamma _{as} \end{aligned}$$

(58)

where $\mathcal {V}^{i/o}$ and $\mathcal {S}^{i/o}$ are respectively the sets of vertex particles and segments belonging to the open boundaries. In a Lagrangian frame, (29) leads to:

$$\begin{aligned} \displaystyle {\frac{d w_{ab}}{d t}}&= \left( \varvec{j}_{a}-\varvec{j}_{b}\right) \cdot \varvec{\nabla }w_{ab} \end{aligned}$$

(59)

$$\begin{aligned} \displaystyle {\frac{d \gamma _{a}}{d t}}&= \sum _{s\in \mathcal {S}}\left( \varvec{j}_{a}-\varvec{j}_{s}\right) \cdot \varvec{\nabla }\gamma _{as} \end{aligned}$$

(60)

If one makes the approximation:

$$\begin{aligned} \displaystyle {\frac{d \gamma _{a}}{d t}}\approx \sum _{s\in \mathcal {S}}\frac{\sigma _{s}}{\sigma _{a}}\left( \varvec{j}_{a}-\varvec{j}_{s}\right) \cdot \varvec{\nabla }\gamma _{as} \end{aligned}$$

(61)

as volume variations remain limited, (56) now writes:

$$\begin{aligned} \displaystyle {\frac{d \sigma _{a}}{d t}}=\displaystyle {\frac{1}{\gamma _{a}}}\displaystyle {\frac{d }{d t}}\left( \sum _{b\in \left( \mathcal {F}\cup \mathcal {V}\right) }\theta _{b}w_{ab}\right) -\delta \sigma _{a}^{i/o}-\displaystyle {\frac{\sigma _{a}}{\gamma _{a}}}\displaystyle {\frac{d \gamma _{a}}{d t}}+\frac{\sigma _{a}}{\gamma _{a}}\delta \gamma _{a}^{i/o} \end{aligned}$$

(62)

Hence:

$$\begin{aligned} \displaystyle {\frac{d }{d t}}\left( \gamma _{a}\sigma _{a}\right) =\gamma _{a}\displaystyle {\frac{d \sigma _{a}}{d t}}+\sigma _{a}\displaystyle {\frac{d \gamma _{a}}{d t}}=\displaystyle {\displaystyle {\frac{d }{d t}}}\left( \sum _{b\in \left( \mathcal {F}\cup \mathcal {V}\right) }\theta _{b}w_{ab}\right) -\gamma _{a}\delta \sigma _{a}^{i/o}+\sigma _{a}\delta \gamma _{a}^{i/o} \end{aligned}$$

(63)

The temporal integration of the continuity equation between $t^{n}$ and $t^{n+1}$ leads to:

$$\begin{aligned} \left( \gamma _{a}\sigma _{a}\right) ^{n+1}-\left( \gamma _{a}\sigma _{a}\right) ^{n}=\sum _{b\in \left( \mathcal {F}\cup \mathcal {V}\right) }\left( \theta _{b}^{n+1}w_{ab}^{n+1}-\theta _{b}^{n}w_{ab}^{n}\right) -\int _{t^{n}}^{t^{n+1}}\gamma _{a}\delta \sigma _{a}^{i/o}+\int _{t^{n}}^{t^{n+1}}\sigma _{a}\delta \gamma _{a}^{i/o} \end{aligned}$$

(64)

With the virtual displacement $\delta r_{a}^{i/o}=\delta t\left( \varvec{J}_{a}^{n}-\varvec{j}_{a}^{n}\right)$, the virtual variations terms write:

$$\begin{aligned} \int _{t^{n}}^{t^{n+1}}\gamma _{a}\delta \sigma _{a}^{i/o}&= \sum _{v\in \mathcal {V}^{I/O}}\theta _{v}^{n}\left( w\left( \varvec{r}_{av}^{n}+\delta \varvec{r}_{v}^{i/o}\right) -w\left( \varvec{r}_{av}^{n}\right) \right) \end{aligned}$$

(65)

$$\begin{aligned} \int _{t^{n}}^{t^{n+1}}\sigma _{a}\delta \gamma _{a}^{i/o}&= \frac{1}{2}\sum _{s\in \mathcal {S}^{I/O}}\sigma _{s}^{n}\left( \varvec{\nabla }\gamma _{as}\left( \varvec{r}_{as}^{n}+\delta \varvec{r}_{s}^{i/o}\right) +\varvec{\nabla }\gamma _{as}\left( \varvec{r}_{as}^{n}\right) \right) \cdot \delta \varvec{r}_{s}^{i/o} \end{aligned}$$

(66)

where $\delta r_{s}^{i/o}=\delta t\left( \varvec{J}_{s}^{n}-\varvec{j}_{s}^{n}\right)$ and $\delta r_{v}^{i/o}=\delta t\left( \varvec{J}_{v}^{n}-\varvec{j}_{v}^{n}\right)$. In the present work, the factor $\sigma _{s}^{n}$ was approximated by $\sigma _{a}^{n}$, consistently with the approximation (61). For an analytical computation of $\gamma _{a}$:

$$\begin{aligned} \int _{t^{n}}^{t^{n+1}}\sigma _{a}\delta \gamma _{a}^{i/o}=\sum _{s\in \mathcal {S}^{I/O}}\sigma _{s}^{n}\left( \gamma _{as}\left( \varvec{r}_{as}^{n}+\delta \varvec{r}_{s}^{i/o}\right) -\gamma _{as}\left( \varvec{r}_{as}^{n}\right) \right) \end{aligned}$$

(67)

1.2 Riemann Problem Formulation

We follow here a similar reasoning to Ferrand et al. (2017). Neglecting the relative velocity and the viscous efforts at the boundary, the homogeneous (i.e. no relative velocity) mixture model projected along the wall normal (n index refers to the normal component and $\tau$ index to the tangential component) writes in non-conservative form:

$$\begin{aligned} \displaystyle {\frac{\partial }{\partial t}}\varvec{W}+\varvec{B}\left( \varvec{W}\right) \displaystyle {\frac{\partial \varvec{W}}{\partial \varvec{n}}}=0 \end{aligned}$$

(68)

where

$$\begin{aligned} \varvec{W}=\left( \begin{array}{c} \sigma \\ \alpha \\ j_{n} \\ j_{\tau } \end{array}\right) \text { and } \varvec{B}\left( \varvec{W}\right) =\left( \begin{array}{cccc} j_{n} &{} 0 &{} \sigma &{} 0 \\ 0 &{} j_{n} &{} 0 &{} 0 \\ \displaystyle {}\frac{c^{2}}{\sigma } &{} D &{} j_{n} &{} 0 \\ 0 &{} 0 &{} 0 &{} j_{n} \end{array}\right) \end{aligned}$$

where $D=\frac{1}{\rho }\left( \rho ^{\alpha }\left( c^{\alpha }\right) ^{2}-\rho ^{\beta }\left( c^{\beta }\right) ^{2}\right) \left( \frac{\sigma }{\sigma _{0}}-1\right)$ is deduced from the linearized state law (28). The eigenvalues of $\varvec{B}$ are therefore $j_{n}$ (of multiplicity 2) and $j_{n}\pm c$ where the sound speed writes:

$$\begin{aligned} c=\sqrt{\frac{\alpha \rho ^{\alpha }\left( c^{\alpha }\right) ^{2}+\beta \rho ^{\beta }\left( c^{\beta }\right) ^{2}}{\alpha \rho ^{\alpha }+\beta \rho ^{\beta }}} \end{aligned}$$

(69)

Following Smoller (1994), one can compute the associated k-Riemann invariants that are detailed in Table 10 where

$$\begin{aligned} \psi \left( \alpha ,\sigma \right) =c\left( \alpha \right) \ln \left( \frac{\sigma }{\sigma _{0}}\right) \end{aligned}$$

(70)

As illustrated on Fig. 25, depending on the sign of $\lambda _{0}$, the state of the boundary segment is determined by the first or second state. As flows considered are subsonic, $\lambda _{-1}<0$ and $\lambda _{+1}>0$. $\lambda _{-1}$ is considered as a ghost wave so that data at exterior state and state 1 are considered to be equal. $\lambda _{0}$ is a contact discontinuity so that the associated k-Riemann invariants hold and simplify in the relations $j_{n,1}=j_{n,2}$ and $c\left( \alpha _{1}\right) \ln \left( \frac{\sigma _{1}}{\sigma _{0}}\right) =c\left( \alpha _{2}\right) \ln \left( \frac{\sigma _{2}}{\sigma _{0}}\right)$. To relate the fluid velocity along the normal of the segment and the pressure, one needs to find the relation between state 2 and the interior state (determined through a SPH interpolation with a renormalization by a Shepard filter) across the wave $\lambda _{+1}$. The characteristic wave $\lambda _{+1}$ can correspond to two possible discontinuities:

• Expansion wave: characteristics are diverging, a smooth transition connects both states and:

  $$\begin{aligned} \lambda _{+1,2}<\lambda _{+1,int}\rightarrow \text { Riemann Invariant } \end{aligned}$$

  (71)

• Shock wave: characteristics are diverging, a smooth transition connects both states and:

  $$\begin{aligned} \lambda _{+1,2}>\lambda _{+1,int}\rightarrow \text { Rankine-Hugoniot } \end{aligned}$$

  (72)

For the shock wave, the Rankine-Hugoniot jump conditions expressed here in terms of conservative variables write:

$$\begin{aligned} S\left( \rho _{2}-\rho _{int}\right)&= \rho _{2}j_{n,2}-\rho _{int}j_{n,int} \end{aligned}$$

(73)

$$\begin{aligned} S\left( \rho _{2}\alpha _{2}-\rho _{int}\alpha _{int}\right)&= \rho _{2}\alpha _{2}j_{n,2}-\rho _{int}\alpha _{int}j_{n,int} \end{aligned}$$

(74)

$$\begin{aligned} S\left( \rho _{2}j_{n,2}-\rho _{int}j_{n,int}\right)&= p_{2}+\rho _{2}j_{n,2}^{2}-p_{int}-\rho _{int}j_{n,int}^{2} \end{aligned}$$

(75)

$$\begin{aligned} S\left( \rho _{2}j_{\tau ,2}-\rho _{int}j_{\tau ,int}\right)&= \rho _{2}j_{n,2}j_{\tau ,2}-\rho _{int}j_{n,int}j_{\tau ,int} \end{aligned}$$

(76)

where S is the speed of the shock. Unknown pressure terms are computed thanks to the state Eq. (28). For the density, we use:

$$\begin{aligned} \frac{\rho }{\alpha \rho ^{\alpha }+\beta \rho ^{\beta }}=\frac{\sigma }{\sigma _{0}} \end{aligned}$$

(77)

Combining (73) and (76) with (74), one gets $\alpha _{2}=\alpha _{int}$ and $j_{\tau ,2}=j_{\tau ,int}$. Equations (73) and (75) can be combined to determine $\rho _{2}$ and then deduce $u_{n,2}$:

$$\begin{aligned} \rho _{2}\rho _{int}\left( j_{n,2}-j_{n,int}\right) ^{2}=\left( \rho _{2}-\rho _{int}\right) \left( p_{2}-p_{int}\right) \end{aligned}$$

(78)

This relation is explicit if $\rho _{2}$ is known but implicit if $j_{n,2}$ is known. One should then iterate. Assuming that $\rho$ variations are small, we deduced a first guess from the linearization of the relation that proved to be sufficiently accurate. For the other quantities to compute, at an inlet (i.e. $\lambda _{0}>0$), tangential velocity and volume fraction are defined by the user. They are determined from the interior state at an outlet (i.e. $\lambda _{0}<0$). In pseudo-code, it leads to the schemes:

Fig. 25

Riemann problems configurations (Ferrand et al. 2017). The state imposed at the boundary is underlined

Table 10 k-Riemann invariants for each eigenvalue

Analytical Solution for the Two-Phase Mixture Poiseuille Flow

At steady state, the volume fraction equation of (6) becomes:

$$\begin{aligned} \varvec{\nabla }\cdot \left( \alpha \beta \varvec{v}^{r}\right) =0 \end{aligned}$$

(79)

Under the longitudinal periodicity condition, it writes:

$$\begin{aligned} \displaystyle {\frac{d }{d z}}\left( \alpha \beta \varvec{v}^{r}\cdot \varvec{e}_{z}\right) =0 \end{aligned}$$

(80)

The no-flux conditions at the upper and lower walls imply that at steady state the Eq. (80) becomes after integration:

$$\begin{aligned} \varvec{v}^{r}\cdot \varvec{e}_{z}=0 \end{aligned}$$

(81)

The volume fraction equation will therefore depend on the chosen closure on the relative velocity. Starting with the closure $\varvec{v}^{r}=\varvec{v}_{0}-K\varvec{\nabla }\alpha /\alpha$, the momentum equation of (3) becomes in this framework:

$$\begin{aligned} \left\{ \begin{array}{l} \displaystyle {\frac{d p}{d z}}=-\,\rho g \\ \displaystyle {\frac{d }{d z}}\left( \rho \nu \displaystyle {\frac{d }{d z}}\varvec{j}\cdot \varvec{e}_{x}\right) +\rho \varvec{F}\cdot \varvec{e}_{x}=0 \end{array}\right. \end{aligned}$$

(82)

In the simplified momentum Eq. (82), the dynamic viscosity is variable and depends on the volume fraction solution of Eq. (81). Let us nondimensionalize the system using $z_{\star }=z/e$, $j_{\star }=\varvec{j}\cdot \varvec{e}_{x}/U_{m0}$ with $U_{m0}=\rho ^{\beta }Fe^{2}/\left( 3\mu ^{\beta }\right)$ (discharge for the usual single-fluid Poiseuille flow), $p_{\star }=p/\left( \rho ^{\beta }ge\right)$ and introduce the Péclet number $\textsf {Pe}=e|\varvec{v}_{0}\cdot \varvec{e}_{z}|/K$ as the ratio of convective and diffusive transports. Noting the density ratio $R_{\rho }=\left( \rho ^{\alpha }-\rho ^{\beta }\right) /\rho ^{\beta }$ and viscosity ratio $R_{\mu }=\left( \mu ^{\alpha }-\mu ^{\beta }\right) /\mu ^{\beta }$, the system becomes:

$$\begin{aligned} \left\{ \begin{array}{l} \displaystyle {\frac{d \alpha }{d z_{\star }}}=\textsf {Pe}\,\alpha \\ \displaystyle {\frac{d p_{\star }}{d z_{\star }}}=-\,\left( 1+R_{\rho }\alpha \right) \\ \displaystyle {\frac{d }{d z_{\star }}}\left[ \left( 1+R_{\mu }\alpha \right) \displaystyle {\frac{d j_{\star }}{d z_{\star }}}\right] =-\,3\left( 1+R_{\rho }\alpha \right) \end{array}\right. \end{aligned}$$

(83)

Constant dynamic viscosity The nondimensionalized solution for a constant dynamic viscosity $\mu =\mu ^{\alpha }=\mu ^{\beta }$ writes:

• Volume fraction

  $$\begin{aligned} \alpha \left( z_{\star }\right) =\alpha _{1}\exp \left( \textsf {Pe}\,z_{\star }\right) \end{aligned}$$

  (84)

• Pressure profile

  $$\begin{aligned} p_{\star }\left( z_{\star }\right) =p_{B\star }+1-z_{\star }+\alpha _{1}R_{\rho }\frac{\exp \left( \textsf {Pe}\right) }{\textsf {Pe}}\left[ 1-\exp \left( -\textsf {Pe}\left( 1-z_{\star }\right) \right) \right] \end{aligned}$$

  (85)

  where $p_{B\star }$ is the background pressure nondimensionalized as the total pressure.

• Longitudinal velocity

  $$\begin{aligned} j_{\star }(z_{\star })&= \displaystyle {}\frac{3}{2}\left( 1-z_{\star }^{2}\right) +\frac{3}{\textsf {Pe}^{2}}\left[ \textsf {Li}_{2,R}\left( z_{\star }\right) +\textsf {Pe}\,z_{\star }\ln _{R}\left( z_{\star }\right) \right. \nonumber \\&\displaystyle {}\left. -\frac{r}{R}\ln _{R}\left( z_{\star }\right) +C_{1}\left( \ln _{R}\left( z_{\star }\right) -\textsf {Pe}\,z_{\star }\right) -C_{2}\right] \end{aligned}$$

  (86)

  where $\textsf {Li}_{2}$ is the dilogarithm function that writes $\textsf {Li}_{2}(x)=-\,\int _{0}^{x}\ln (1-t)/t\,dt$Abramowitz and Stegun (1964). If $|x|\le 1$, one can write a series expression $\textsf {Li}_{2}(x)=\sum _{n}x^{n}/n^{2}$. We introduced the notations $\ln _{R}\left( z_{\star }\right) =\ln \left( 1+\alpha _{1}R_{\mu }\exp \left( \textsf {Pe}\,z_{\star }\right) \right)$ and $\textsf {Li}_{2,R}\left( z_{\star }\right) =\textsf {Li}_{2}\left( -\alpha _{1}R_{\mu }\exp \left( \textsf {Pe}\,z_{\star }\right) \right)$. $C_{1}$ and $C_{2}$ are deduced from the no-slip condition at walls:

  $$\begin{aligned} C_{1}&= \displaystyle {}\frac{\textsf {Li}_{2,R}(1)-\textsf {Li}_{2,R}(-1)+\textsf {Pe}\left[ \ln _{R}\left( 1\right) +\ln _{R}\left( -1\right) \right] -R_{\rho }/R_{\mu }\left[ \ln _{R}\left( 1\right) -\ln _{R}\left( -1\right) \right] }{2\textsf {Pe}-\ln _{R}\left( 1\right) +\ln _{R}\left( -1\right) } \end{aligned}$$

  (87)
  $$\begin{aligned} C_{2}&= \textsf {Li}_{2,R}(1)+\textsf {Pe}\ln _{R}(1)-\frac{R_{\rho }}{R_{\mu }}\ln _{R}(1)+C_{1}\left( \ln _{R}(1)-\textsf {Pe}\right) \end{aligned}$$

  (88)

$\alpha _{1}$ is computed thanks to the conservation of volume (integrating over the height of the channel) for a given initial uniform profile of $\alpha \left( z\right) =\alpha _{0}$:

$$\begin{aligned} \frac{\alpha _{1}}{\alpha _{0}}=\frac{\textsf {Pe}}{\sinh \left( \textsf {Pe}\right) } \end{aligned}$$

(89)

To avoid complete separation of phases, (84) gives a condition on $\alpha _{1}$:

$$\begin{aligned} 0\le \alpha _{1}\le \exp \left( -\textsf {Pe}\right) \end{aligned}$$

(90)

And therefore a condition on the initial uniform volume fraction $\alpha _{0}$ using (89):

$$\begin{aligned} 0\le \alpha _{0}\le \frac{1-\exp \left( -2\textsf {Pe}\right) }{2\textsf {Pe}} \end{aligned}$$

(91)