Model of incompressible Navier-Stokes with interface-capturing equation

Bulk and local densities

We consider two fluids \(A\) and \(B\) of constant densities \(\rho_A\) and \(\rho_B\) and constant viscosities \(\nu_A\) and \(\nu_B\). Those properties are called bulk densities and bulk viscosities. We introduce a phase index \(\phi(\boldsymbol{x},t)\) of value 0 in the first phase \(A\) and +1 in the second phase \(B\). With the phase-field \(\phi(\boldsymbol{x},t)\) we can define local densities depending on position and time \(\varrho_A(\boldsymbol{x},t)\) and \(\varrho_B(\boldsymbol{x},t)\) by:

(523)\[\varrho_{A}(\boldsymbol{x},t)=\rho_{A}(1-\phi(\boldsymbol{x},t))\]

and

(524)\[\varrho_{B}(\boldsymbol{x},t)=\rho_{B}\phi(\boldsymbol{x},t)\]

The total local density is

(525)\[\varrho(\boldsymbol{x},t)=\rho_{B}\phi(\boldsymbol{x},\,t)+\rho_{A}(1-\phi(\boldsymbol{x},t))\]

With that definition, if \(\phi(\boldsymbol{x},t)=0\), then \(\varrho(\boldsymbol{x},t)=\rho_A\) and if \(\phi(\boldsymbol{x},t)=1\), then \(\varrho(\boldsymbol{x},t)=\rho_B\).

Mass balance

A mass balance for each local density yields

(530)\[\begin{split}\begin{eqnarray} \frac{\partial\varrho_{A}}{\partial t}+\boldsymbol{\nabla}\cdot(\varrho_{A}\boldsymbol{u}+\rho_{A}\boldsymbol{j}_{A})& = & -\dot{m}'''\\ \frac{\partial\varrho_{B}}{\partial t}+\boldsymbol{\nabla}\cdot(\varrho_{B}\boldsymbol{u}+\rho_{B}\boldsymbol{j}_{B})& = & +\dot{m}''' \end{eqnarray}\end{split}\]

where \(\varrho_{B}\boldsymbol{u}\) and \(\varrho_{A}\boldsymbol{u}\) are two advective fluxes and \(\rho_{A}\boldsymbol{j}_{A}\) and \(\rho_{B}\boldsymbol{j}_{B}\) are two diffusive fluxes. Here we assume that \(\boldsymbol{j}_{A}\) and \(\boldsymbol{j}_{B}\) are equal and opposite:

(527)\[\boldsymbol{j}=\boldsymbol{j}_{B}=-\boldsymbol{j}_{A}\]

On the right-hand side, \(\dot{m}'''\) is a production source term: if it is added in one equation, then that quantity is substracted from the other equation. The three primes \('''\) means that the production is volumic. Its physical dimension is \([\text{M}]/([\text{L}]^3.[\text{T}])\). Its value is non zero when phase change occurs. For two immiscible fluids we can consider that \(\dot{m}'''=0\).

Expressed with \(\phi\), those two equations write

(528)\[\begin{split}\begin{eqnarray} \frac{\partial\phi}{\partial t}+\boldsymbol{\nabla}\cdot(\boldsymbol{u}\phi+\boldsymbol{j})& = & +\frac{\dot{m}'''}{\rho_{B}}\\ \frac{\partial(1-\phi)}{\partial t}+\boldsymbol{\nabla}\cdot(\boldsymbol{u}(1-\phi)-\boldsymbol{j}) & = & -\frac{\dot{m}'''}{\rho_{A}} \end{eqnarray}\end{split}\]

By summing those two equations, we obtain:

(529)\[\boldsymbol{\nabla}\cdot\boldsymbol{u}=\dot{m}'''\left(\frac{1}{\rho_{B}}-\frac{1}{\rho_{A}}\right)\]

Incompressible two-phase flows without phase change

Without phase change \(\dot{m}'''=0\) and we retrieve the classical mass balance:

(530)\[\boldsymbol{\nabla}\cdot\boldsymbol{u}=0\]

For flux \(\boldsymbol{j}\), two hypotheses are used.

  • For the first one, the flux is given by the gradient of chemical potential \(\mu_{\phi}\):

(531)\[\begin{split}\begin{eqnarray} \boldsymbol{j} & = & -\mathcal{M}_{\phi}\boldsymbol{\nabla}\mu_{\phi}\\ & = & -M_{\phi}\boldsymbol{\nabla}\left[2\phi(1-\phi)(1-2\phi)-\frac{W^2}{8}\boldsymbol{\nabla}^{2}\phi \right] \end{eqnarray}\end{split}\]

and Eq. (528) becomes the CH model:

(532)\[\frac{\partial\phi(\boldsymbol{x},t)}{\partial t}=\boldsymbol{\nabla}\cdot\left\{M_{\phi}\boldsymbol{\nabla}\left[2\phi(1-\phi)(1-2\phi)-\frac{W^2}{8}\boldsymbol{\nabla}^{2}\phi \right]\right\}\]

  • For the second one, the flux is chosen as standard diffusive flux \(\boldsymbol{j}_{diff}\) with counter term flux \(\boldsymbol{j}_{CT}\) (see Conservative Allen-Cahn (CAC) model for origin and interpretation of counter term):

(533)\[\boldsymbol{j}=-M_{\phi}\boldsymbol{\nabla}\phi+\frac{4}{W}\phi(1-\phi)\boldsymbol{n}_{\phi}\]

and Eq. (528) becomes the CAC model:

(534)\[\frac{\partial\phi}{\partial t}+\boldsymbol{\nabla}\cdot(\boldsymbol{u}\phi)=\boldsymbol{\nabla}\cdot\left[M_{\phi}\left(\boldsymbol{\nabla}\phi-\frac{4}{W}\phi(1-\phi)\boldsymbol{n}_{\phi}\right)\right]\]

Model of Navier-Stokes with interface-capturing equation

The model of two-phase flows is simply composed of incompressible Navier-Stokes equations which hold in both fluids \(A\) and \(B\). An interpolation with \(\phi\) is performed for local densities and local viscosities, and an additional force term is added in the impulsion balance equation: the capillary force representative of surface tension \(\sigma\) between both fluids. The model writes:

(535)\[\begin{split}\begin{eqnarray} \boldsymbol{\nabla}\cdot\boldsymbol{u} & = & 0\\ \varrho(\phi)\underbrace{\left[\frac{\partial\boldsymbol{u}}{\partial t}+\boldsymbol{\nabla}\cdot(\boldsymbol{u}\boldsymbol{u})\right]}_{\text{Acceleration}}& = & \underbrace{-\boldsymbol{\nabla}p_{h}}_{\text{Pressure force}}+\underbrace{\boldsymbol{\nabla}\cdot\left[\varrho(\phi)\vartheta(\phi)(\boldsymbol{\nabla}\boldsymbol{u}+\boldsymbol{\nabla}\boldsymbol{u}^{T})\right]}_{\text{Viscous force}}+\underbrace{\mu_{\phi}\boldsymbol{\nabla}\phi}_{\text{Capillary force }\boldsymbol{F}_{c}}+\underbrace{\varrho(\phi)\boldsymbol{g}}_{\text{Buoyancy}} \end{eqnarray}\end{split}\]

where \(\boldsymbol{u}\) is the mean velocity, \(p_h\) is the hydrodynamic pressure, \(\boldsymbol{g}\) is the gravity. The local properties of fluids are noted \(\varrho(\phi)\) for density and \(\vartheta(\phi)\) for kinematic viscosity. Their expressions will will summarized below. The physical dimensions of each term must be \([\text{F}]/[\text{L}]^3\) where \([\text{F}]\) is used for force: \([\text{F}]=[\text{M.L}]/[\text{T}]^2\). For example for the buoyancy force term

(536)\[[\varrho(\phi)\boldsymbol{g}]=\frac{[\text{M}]}{[\text{L}]^3} \frac{[\text{L}]}{[\text{T}]^2}=\frac{[\text{F}]}{[\text{L}]^3}\]

The capillary force \(\boldsymbol{F}_c=\mu_{\phi}\boldsymbol{\nabla}\phi\) has been formulated in [1] and [2]. That form is equivalent to the surface tension force \(-\delta_d \sigma \kappa \boldsymbol{n}\) where \(\kappa\) is the curvature, \(\sigma\) is the surface tension, \(\boldsymbol{n}\) is the normal vector at the interface and \(\delta_d\) is the Kronecker’s symbol but expressed in the phase-field framework, i.e. with a diffuse interface:

(537)\[\delta_{d}=\frac{3}{2}W\left|\boldsymbol{\nabla}\phi\right|^{2}\]

That term is homogeneous to an inverse of length: \([\delta_d]=1/[\text{L}]\). We can prove that

(538)\[\begin{split}\begin{eqnarray} \boldsymbol{F}_{c} & = & \mu_{\phi}\boldsymbol{\nabla}\phi\\ & = & \Bigl[4H\phi(\phi-1)(\phi-1/2)-\zeta\boldsymbol{\nabla}^{2}\phi\Bigr]\boldsymbol{\nabla}\phi\\ & = & -\frac{3}{2}W\sigma\Bigl[\Delta\phi-\frac{16}{W^{2}}\phi(1-\phi)(1-2\phi)\Bigr]\boldsymbol{\nabla}\phi\\ & = & -\frac{3}{2}W\sigma\kappa\left|\boldsymbol{\nabla}\phi\right|\boldsymbol{\nabla}\phi\\ & = & -\delta_{d}\sigma\kappa\boldsymbol{n} \end{eqnarray}\end{split}\]

where we have used

  • Eq. (429) for \(\mu_{\phi}\) in the second line.

  • Replace \(H\) and \(\zeta\) by \(H=12\sigma/W\) and \(\zeta=(3/2)W\sigma\) in the third line.

  • Eq. (521) for \(\kappa\left|\boldsymbol{\nabla}\phi\right|\) in the fourth line.

  • Multiply and divide by \(\left|\boldsymbol{\nabla}\phi\right|\) to make appear \(\boldsymbol{n}\) in the last line.

We can check that the capillary force is homogeneous to a volumic force:

(539)\[[\delta_{d}\sigma\kappa\boldsymbol{n}]=\frac{1}{[\text{L}]}\times \frac{[\text{F}]}{[\text{L}]}\times \frac{1}{[\text{L}]}=\frac{[\text{F}]}{[\text{L}]^3}\]

When expressed by its potential form:

(540)\[[\mu_{\phi}\boldsymbol{\nabla}\phi]=\frac{\text{E}}{[\text{L}]^3}\times \frac{1}{[\text{L}]}=\frac{[\text{F.L}]}{[\text{L}]^3}\times \frac{1}{[\text{L}]}=\frac{[\text{F}]}{[\text{L}]^3}\]

Two-phase incompressible Navier-Stokes

The two-phase flows model is composed of incompressible Navier-Stokes equations. The first one is the mass balance equation

(541)\[\boldsymbol{\nabla}\cdot\boldsymbol{u}=0\]

The second one is the impulsion balance equation:

(542)\[\varrho(\phi)\left[\frac{\partial\boldsymbol{u}}{\partial t}+\boldsymbol{\nabla}\cdot(\boldsymbol{u}\boldsymbol{u})\right]=-\boldsymbol{\nabla}p_{h}+\boldsymbol{\nabla}\cdot\left[\varrho(\phi)\vartheta(\phi)\left(\boldsymbol{\nabla}\boldsymbol{u}+\boldsymbol{\nabla}\boldsymbol{u}^{T}\right)\right]+\mu_{\phi}\boldsymbol{\nabla}\phi+\varrho(\phi)\boldsymbol{g}\]

where the local density \(\varrho(\phi)\) is defined by Eq. (525):

(543)\[\varrho(\boldsymbol{x},t)=\rho_{B}\phi(\boldsymbol{x},\,t)+\rho_{A}(1-\phi(\boldsymbol{x},t))\]

and the kinematic viscosity \(\vartheta(\phi)\) is interpolated with the harmonic mean:

(544)\[\frac{1}{\vartheta(\phi)}=\frac{\phi}{\nu_{B}}+\frac{1-\phi}{\nu_{A}}\]

In capillary force of Eq. (542) the chemical potential \(\mu_{\phi}\) is defined by

(545)\[\mu_{\phi}=\frac{3}{2}\sigma W\left[\frac{16}{W^2}\phi(1-\phi)(1-2\phi)-\boldsymbol{\nabla}^{2}\phi\right]\]

The phase-field \(\phi\) follows one of the two phase-field equation. The Conservative Allen-Cahn equation writes

Interface-capturing model 1: Conservative Allen-Cahn equation

(546)\[\frac{\partial\phi}{\partial t}+\boldsymbol{\nabla}\cdot(\boldsymbol{u}\phi)=\boldsymbol{\nabla}\cdot\left[M_{\phi}\left(\boldsymbol{\nabla}\phi-\frac{4}{W}\phi(1-\phi)\boldsymbol{n}_{\phi}\right)\right]\]

where \(M_{\phi}\) is the mobility coefficient of interface, \(W\) is the interface width and the unit normal vector \(\boldsymbol{n}_{\phi}\) is defined by

(547)\[\boldsymbol{n}_{\phi}=\frac{\boldsymbol{\nabla}\phi}{|\boldsymbol{\nabla}\phi|}\]

Interface-capturing model 2: Cahn-Hilliard equation

Whereas the Cahn-Hilliard equation writes

(548)\[\frac{\partial\phi(\boldsymbol{x},t)}{\partial t}=\boldsymbol{\nabla}\cdot\left\{M_{\phi}\boldsymbol{\nabla}\left[2\phi(1-\phi)(1-2\phi)-\frac{W^2}{8}\boldsymbol{\nabla}^{2}\phi \right]\right\}\]

Bibliography

Section author: Alain Cartalade