Derivation of incompressible Navier-Stokes with interface-capturing equation model

Mass balance & interface 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:

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

and

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

The total local density is

(465)\[\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 and interface equation

A mass balance for each local density yields

(470)\[\begin{split}\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{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:

(467)\[\boldsymbol{j}_{\phi}=\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

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

By summing those two equations, we obtain:

(469)\[\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:

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

For interface tracking equation, Eq. (468) becomes:

\[\frac{\partial\phi}{\partial t}+\boldsymbol{\nabla}\cdot(\boldsymbol{u}\phi)=-\boldsymbol{\nabla}\cdot\boldsymbol{j}_{\phi}\]

where \(\boldsymbol{j}_{\phi}\) has to be determined. With the material derivative:

(471)\[\frac{d}{dt}\,\hat{=}\,\frac{\partial}{\partial t}+\boldsymbol{u}\cdot\boldsymbol{\nabla}\]

the interface-tracking equation writes:

\[\begin{split}\frac{\partial\phi}{\partial t}+\boldsymbol{u}\cdot\boldsymbol{\nabla}\phi+\phi\boldsymbol{\nabla}\cdot\boldsymbol{u}&=-\boldsymbol{\nabla}\cdot\boldsymbol{j}_{\phi}\\ \frac{d\phi}{dt}+\phi\boldsymbol{\nabla}\cdot\boldsymbol{u}&=-\boldsymbol{\nabla}\cdot\boldsymbol{j}_{\phi}\end{split}\]

Derivation of constitutive laws

Summary of balance equations

Finally, we summarize the three balance equations with the material derivative. First, the mass balance equation writes:

(472)\[\frac{d\rho_{0}}{dt}=-\rho_{0}\boldsymbol{\nabla}\cdot\boldsymbol{u}\]

The impulsion balance equation is:

(473)\[\rho_{0}\frac{d\boldsymbol{u}}{dt}=\boldsymbol{\nabla}\cdot{\color{red}\overline{\overline{\boldsymbol{T}}}}\]

where the stress tensor \({\color{red}\overline{\overline{\boldsymbol{T}}}}\) has to be determined. At last, the balance of phase-field \(\phi\) writes:

(474)\[\frac{d\phi}{dt}+\phi\boldsymbol{\nabla}\cdot\boldsymbol{u}=-\boldsymbol{\nabla}\cdot{\color{red}\boldsymbol{j}_{\phi}}\]

where \({\color{red}\boldsymbol{j}_{\phi}}\) has also to be determined.

Objective

Objective

In Section Time-evolution: Cahn-Hilliard and Allen-Cahn models the diffusive flux \(\boldsymbol{j}_{diff}\) has been derived such as \(\partial\mathscr{F}[\phi]/\partial t < 0\). Here we consider the total energy \(\mathscr{E}_{tot}\) which is a functional depending on two functions: the phase-field \(\phi\) and the velocity \(\boldsymbol{u}\). It is composed of two parts:

(475)\[\mathscr{E}_{tot}[\phi,{\color{red}\boldsymbol{u}}]=\int_{V}\biggl[\underbrace{\frac{1}{2}\rho_{0}\bigl|{\color{red}\boldsymbol{u}}\bigr|^{2}}_{\text{kinetic energy}}+\underbrace{\mathcal{F}(\phi,\boldsymbol{\nabla}\phi)}_{\text{potential energy}}\biggr]dV\]

where \(\rho_{0}\bigl|\boldsymbol{u}\bigr|^{2}/2\) is the kinetic energy (\(\rho_0\) is the constant density) and the potential energy is given by the free energy density:

(476)\[\mathcal{F}(\phi,\boldsymbol{\nabla}\phi)=f_{dw}(\phi)+\frac{\zeta}{2}\bigl|\boldsymbol{\nabla}\phi\bigr|^{2}\]

In Eq. (476), \(f_{dw}(\phi)\) is the double-well free energy density, and \(\zeta\) is the capillary coefficient.

The objectice is determine the constitutive laws for the stress tensor \(\overline{\overline{\boldsymbol{T}}}\) in Eq. (473) and the flux \(\boldsymbol{j}_{\phi}\) in Eq. (474) such as

(477)\[\frac{d\mathscr{E}_{tot}}{dt}-\int_{V}\underbrace{\lambda\boldsymbol{\nabla}\cdot\boldsymbol{u}}_{\hat{=}\mathscr{L}}\leqslant0\]

where \(\lambda\) is the lagrange multiplier of constraint \(\boldsymbol{\nabla}\cdot\boldsymbol{u}=0\).

Method

Express

(478)\[\frac{d\mathscr{E}_{tot}}{dt} =\frac{d}{dt}\int_{V}\left[\mathcal{F}(\phi,\boldsymbol{\nabla}\phi)+\frac{1}{2}\rho_{0}\bigl|\boldsymbol{u}\bigr|^{2}\right]dV\]

on the form

(479)\[\frac{d\mathscr{E}_{tot}}{dt}=-\mathcal{D}(V){\color{gray}+\underbrace{\mathcal{W}(V)+\Phi(\partial V)}_{\text{neglected here}}}\leqslant0\]

where \(\mathcal{D}\) is the dissipation with \(\mathcal{D}\geqslant0\), \(\mathcal{W}\) is the work of external forces and \(\Phi\) is the flux through surface.

For that purpose we use the Reynolds transport theorem:

(480)\[\frac{d}{dt}\left[\int_{V}\Psi dV\right]\,\hat{=}\,\int_{V}\left[\frac{d\Psi}{dt}+\Psi\boldsymbol{\nabla}\cdot\boldsymbol{u}\right]dV\]

Example:

\[\begin{split}\frac{d\mathscr{E}_{tot}}{dt}&=\frac{d}{dt}\int_{V}\left[\mathcal{F}(\phi,\boldsymbol{\nabla}\phi)+\frac{1}{2}\rho_{0}\bigl|\boldsymbol{u}\bigr|^{2}\right]dV\\ &=\int_{V}\left\{ \frac{d}{dt}\left[\mathcal{F}+\frac{1}{2}\rho_{0}\bigl|\boldsymbol{u}\bigr|^{2}\right]+\left[\mathcal{F}+\frac{1}{2}\rho_{0}\bigl|\boldsymbol{u}\bigr|^{2}\right]\boldsymbol{\nabla}\cdot\boldsymbol{u}\right\} dV\\ &=\int_{V}\biggl\{\left[\frac{d\mathcal{F}}{dt}+\mathcal{F}\boldsymbol{\nabla}\cdot\boldsymbol{u}\right]+\biggl[\cancel{\frac{1}{2}\frac{d\rho_{0}}{dt}\bigl|\boldsymbol{u}\bigr|^{2}}+\frac{1}{2}\rho_{0}\frac{d\bigl|\boldsymbol{u}\bigr|^{2}}{dt}+\cancel{\frac{1}{2}\rho_{0}\bigl|\boldsymbol{u}\bigr|^{2}\boldsymbol{\nabla}\cdot\boldsymbol{u}}\biggr]\biggr\} dV\\ &=\int_{V}\biggl\{\underbrace{\left[\frac{d\mathcal{F}}{dt}+\mathcal{F}\boldsymbol{\nabla}\cdot\boldsymbol{u}\right]}_{\hat{=}\mathcal{I}}+\underbrace{\frac{1}{2}\rho_{0}\frac{d\bigl|\boldsymbol{u}\bigr|^{2}}{dt}}_{\hat{=}\mathcal{K}}\biggr\} dV\end{split}\]

The next stages require to make appear balance equations in \(\mathcal{I}\) and \(\mathcal{K}\) and evaluate the differentials \(d\mathcal{F}\) and \(d\bigl|\boldsymbol{u}\bigr|^{2}\).

The details of derivation are presented in Details for deriving constitutive laws. Finally the constitutive laws which satisfy the condition \(\frac{d\mathscr{E}_{tot}}{dt}=-\mathcal{D}(V)\leqslant0\) can be chosen such as:

Constitutive laws

Stress tensor \(\overline{\overline{\boldsymbol{T}}}\)
(481)\[\overline{\overline{\boldsymbol{T}}}=-\overline{\overline{\boldsymbol{P}}}+\eta(\boldsymbol{\nabla}\boldsymbol{u}+\boldsymbol{\nabla}\boldsymbol{u}^{T})\]

with the pressure tensor is defined by

(482)\[\overline{\overline{\boldsymbol{P}}}=\Bigl[(p_{h}-\mathcal{F})\overline{\overline{\boldsymbol{I}}}+\zeta\boldsymbol{\nabla}\phi\otimes\boldsymbol{\nabla}\phi\Bigr]\]
Flux \(\boldsymbol{j}_{\phi}\)
(483)\[\boldsymbol{j}_{\phi}=-\mathcal{M}_{\phi}\boldsymbol{\nabla}\mu_{\phi}\]

where the chemical potential is defined by

(484)\[\mu_{\phi}=f_{dw}^{\prime}(\phi)-\zeta\boldsymbol{\nabla}^{2}\phi\]

By replacing (481) in Eq. (473) we obtain two terms in the right-hand side: \(-\boldsymbol{\nabla}\cdot\overline{\overline{\boldsymbol{P}}}\) and \(\boldsymbol{\nabla}\cdot\eta(\boldsymbol{\nabla}\boldsymbol{u}+\boldsymbol{\nabla}\boldsymbol{u}^{T})\). The first one can be written with its potential form:

Potential form of pressure tensor
(485)\[-\boldsymbol{\nabla}\cdot\overline{\overline{\boldsymbol{P}}}=-\boldsymbol{\nabla}p_{h}+\mu_{\phi}\boldsymbol{\nabla}\phi\]

The interpretation of term \(\mu_{\phi}\boldsymbol{\nabla}\phi\) is the capillary force \(\boldsymbol{F}_{c}\) (see proof below), and \(p_h\) is the hydrodynamic pressure ensuring the continuity equation.

Proof of Eq. (485):

\[\begin{split}-\boldsymbol{\nabla}\cdot\overline{\overline{\boldsymbol{P}}}&=-\boldsymbol{\nabla}p_{h}+\boldsymbol{\nabla}\mathcal{F}-\zeta\boldsymbol{\nabla}\cdot(\boldsymbol{\nabla}\phi\otimes\boldsymbol{\nabla}\phi)\\&=-\boldsymbol{\nabla}p_{h}+\boldsymbol{\nabla}\Bigl[f_{dw}+\frac{\zeta}{2}\bigl|\boldsymbol{\nabla}\phi\bigr|^{2}\Bigr]-\zeta\boldsymbol{\nabla}\cdot(\boldsymbol{\nabla}\phi\otimes\boldsymbol{\nabla}\phi)\\&=-\boldsymbol{\nabla}p_{h}+\boldsymbol{\nabla}f_{dw}+\cancel{\frac{\zeta}{2}\boldsymbol{\nabla}(\bigl|\boldsymbol{\nabla}\phi\bigr|^{2})}-\cancel{\frac{\zeta}{2}\boldsymbol{\nabla}(\bigl|\boldsymbol{\nabla}\phi\bigr|^{2})}-\zeta(\boldsymbol{\nabla}^{2}\phi)\boldsymbol{\nabla}\phi\\&=-\boldsymbol{\nabla}p_{h}+\Bigl[\underbrace{f_{dw}^{\prime}-\zeta(\boldsymbol{\nabla}^{2}\phi)}_{\equiv\mu_{\phi}}\Bigr]\boldsymbol{\nabla}\phi\end{split}\]

Model of incompressible Navier-Stokes with interface-capturing equation

Incompressible Navier-Stokes with capillary force term

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:

(486)\[\begin{split}\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{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

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

Interpretation of term \(\mu_{\phi}\boldsymbol{\nabla}\phi\)

Interpretation of term \(\mu_{\phi}\boldsymbol{\nabla}\phi\)

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}_{\phi}\) where \(\kappa\) is the curvature, \(\sigma\) is the surface tension, \(\boldsymbol{n}_{\phi}\) 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:

(488)\[\boldsymbol{F}_{c}=\mu_{\phi}\boldsymbol{\nabla}\phi=-\delta_{d}\sigma\kappa\boldsymbol{n}_{\phi}\hspace{1cm}\text{with}\hspace{1cm}\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 check that the capillary force is homogeneous to a volumic force:

(489)\[[\delta_{d}\sigma\kappa\boldsymbol{n}_{\phi}]=\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:

(490)\[[\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}\]

Proof

\[\begin{split}\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{split}\]

Where we have used

  • Eq. (365) 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. (461) 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.

Summary

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

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

The second one is the impulsion balance equation:

(492)\[\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. (465):

(493)\[\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:

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

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

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

Interface-capturing model 1: Cahn-Hilliard equation

Once the diffusive flux (483) is replaced in Eq. (474), the advective Cahn-Hilliard equation (Eq. (445)) is obtained:

(496)\[\frac{\partial\phi(\boldsymbol{x},t)}{\partial t}+\boldsymbol{\nabla}\cdot(\boldsymbol{u}\phi)=\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\}\]

That model can advantageously be replaced by the Conservative Allen-Cahn equation (Eq. (458)) for immiscible two fluid flows:

(497)\[\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

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

Appendix: details for deriving constitutive laws

References

Section author: Alain Cartalade