Incompressible Navier-Stokes equations for two-phase flows
For low Mach Navier-Stokes (NS), the EqDF are presented in Continuous Boltzmann equation. Here are presented some modification for simulating incompressible single-phase and two-phase flows. First we remind the macroscopic PDEs we want to simulate:
(189)\[\boldsymbol{\nabla}\cdot\boldsymbol{u}=0\]
(190)\[\varrho(\phi)\left[\frac{\partial\boldsymbol{u}}{\partial t}+\boldsymbol{\nabla}\cdot(\boldsymbol{u}\boldsymbol{u})\right]=-\boldsymbol{\nabla}p_{h}+\boldsymbol{\nabla}\cdot\eta(\boldsymbol{x},t)\left[\boldsymbol{\nabla}\boldsymbol{u}+(\boldsymbol{\nabla}\boldsymbol{u})^{T}\right]\]
where \(p_h\) is the hydrodynamic pressure, \(\boldsymbol{u}\) is the fluid velocity , and \(\varrho(\phi)\) is the interpolation of constant bulk densities: each phase is characterized by constant densities, e.g. \(\rho_0\) and \(\rho_1\). The total density \(\varrho(\boldsymbol{x},t)\) is a function of position and time with the phase-field:
(191)\[\varrho(\phi)=\rho_0\phi(\boldsymbol{x},t)+\rho_1(1-\phi(\boldsymbol{x},t))\]
Standard \(f_i^{eq}\) for incompressible single-phase
For simulating incompressible NS equations for single-phase flows, the EqDF writes:
(192)\[f_{i}^{eq}(\boldsymbol{x},t)=w_{i}\left[p_{h}+\rho_{0}c_{s}^{2}\left(\frac{\boldsymbol{c}_{i}\cdot\boldsymbol{u}}{c_{s}^{2}}+\frac{(\boldsymbol{c}_{i}\cdot\boldsymbol{u})^{2}}{2c_{s}^{4}}-\frac{\boldsymbol{u}\cdot\boldsymbol{u}}{2c_{s}^{2}}\right)\right]\]
where the macroscopic physical variables are the hydrodynamic pressure \(p_h\), the fluid velocity \(\boldsymbol{u}\), and the constant bulk density \(\rho_0\). The lattice is defined by its weights \(w_i\) and its moving directions \(\boldsymbol{c}_i\) where \(i=0,...,N_{pop}\). \(c_s\) is the lattice sound speed defined by \(c/\sqrt{3}\) which depends the the time and space discretization. The lattice speed \(c\) is defined by \(\delta x/\delta t\).
After collision and streaming, the moments are updated with the following relationships. First, the moment of order zero is the hydrodynamic pressure:
(193)\[p_{h}(\boldsymbol{x},t)=\sum_{i}f_{i}(\boldsymbol{x},t)\]
Next, the moment of order one is the impulsion:
(194)\[\rho_{0}\boldsymbol{u}(\boldsymbol{x},t)=\frac{1}{c_{s}^{2}}\sum_{i}f_{i}(\boldsymbol{x},t)\boldsymbol{c}_{i}\]
Equivalent algorithm of artificial compressibility
We can prove that the Chapman-Enskog expansion of LBE with EqDF (192) yields to the equivalent system of equations:
(195)\[\frac{\partial{\color{red}p_h}}{\partial t}+{\color{red}\rho_{0}c_{s}^{2}}\boldsymbol{\nabla}\cdot\boldsymbol{u}=0\]
(196)\[{\color{red}\rho_{0}\cancel{c_{s}^{2}}}\left[\frac{\partial\boldsymbol{u}}{\partial t}+\boldsymbol{\nabla}\cdot(\boldsymbol{u}\boldsymbol{u})\right]=-\boldsymbol{\nabla}{\color{red}p_h}\cancel{c_{s}^{2}}+\boldsymbol{\nabla}\cdot\left[{\color{red}\rho_{0}\cancel{c_{s}^{2}}}\nu(\boldsymbol{\nabla}\boldsymbol{u}+(\boldsymbol{\nabla}\boldsymbol{u})^{T})\right]+\mathcal{O}(\rho u^{3})\]
The mass balance Eq. (195) is one particular method of classical CFD to simulate the incompressible Navier-Stokes equations, called the articial compressibility algorithm. In LBM the compressibility coefficient is \(\beta=\rho_0c_s^{2}\).
\(f_i^{eq}\) for incompressible two-phase flows
Now we consider two-phase flows by assuming that each phase is incompressible.
From the previous EqDF, two methods are implemented in LBM_Saclay to account for that variation (191). The first method is a direct extension of (192). The second one is a slight modification which introduces the dimensionless pressure \(p_h^{\star}\).
Version 1 for variable density
That first version is a direct adaptation of Eq. (192) where the variable density \(\varrho(\phi)\) simply replaces the constant density \(\rho_0\):
(197)\[f_{i}^{eq}(\boldsymbol{x},t)=w_{i}\left[p_{h}+\varrho(\phi)c_{s}^{2}\left(\frac{\boldsymbol{c}_{i}\cdot\boldsymbol{u}}{c_{s}^{2}}+\frac{(\boldsymbol{c}_{i}\cdot\boldsymbol{u})^{2}}{2c_{s}^{4}}-\frac{\boldsymbol{u}\cdot\boldsymbol{u}}{2c_{s}^{2}}\right)\right]\]
Let us notice that the physical dimension of terms inside the brackets is equivalent to a pressure. The hydrodynamic pressure \(p_h\) and the impulsion are respectively obtained by the moment of order zero:
(198)\[p_{h}(\boldsymbol{x},t)=\sum_{i}f_{i}(\boldsymbol{x},t)\]
and its moment of order 1:
(199)\[\varrho(\phi)\boldsymbol{u}(\boldsymbol{x},t)=\frac{1}{c_{s}^{2}}\sum_{i}f_{i}(\boldsymbol{x},t)\boldsymbol{c}_{i}\]
Version 1: equivalent macroscopic PDE
Once the Chapman-Enskog expansion is performed, the macroscopic equations which recovered are with that EqDF are:
(200)\[\frac{\partial p_{h}}{\partial t}+\boldsymbol{\nabla}\cdot\left[\varrho(\phi)c_{s}^{2}\boldsymbol{u}\right]=0\]
for the mass balance equation, and
(201)\[\varrho(\phi)\cancel{c_{s}^{2}}\left[\frac{\partial\boldsymbol{u}}{\partial t}+\boldsymbol{\nabla}\cdot(\boldsymbol{u}\boldsymbol{u})\right]=-\boldsymbol{\nabla}p_{h}\cancel{c_{s}^{2}}+\boldsymbol{\nabla}\cdot\left[\varrho(\phi)\cancel{c_{s}^{2}}\nu(\boldsymbol{\nabla}\boldsymbol{u}+(\boldsymbol{\nabla}\boldsymbol{u})^{T})\right]+\mathcal{O}(\rho u^{3})\]
for the impulsion balance, where the kinematic viscosity is given by
(202)\[\nu=\frac{1}{3}\left(\tau-\frac{1}{2}\right)\frac{\delta x^{2}}{\delta t}\]
By expanding the divergence term, the mass balance Eq. (200) can be rewritten:
(203)\[\frac{\partial p_{h}}{\partial t}+\varrho(\phi)c_{s}^{2}\boldsymbol{\nabla}\cdot\boldsymbol{u}=-\boldsymbol{u}\cdot\boldsymbol{\nabla}(\varrho(\phi)c_{s}^{2})\]
Version 1: source term to add
The mass balance recovered by the Chapman-Ensokg procedure Eq. (200) slighty differs from the classic artificial mass balance equation because the density is a function of position \(\varrho(\phi)\). To correct that, it is needed to add a source term:
(204)\[S_p=\boldsymbol{u}\cdot\boldsymbol{\nabla}[\varrho(\phi)c_s^2]\]
in the LBE.
Version 2 for variable density
That version is implemented in the kernel NSAC_Comp. The pressure that is used in the equilibrium distribution function is a dimensionless pressure defined by:
(205)\[p_h^{\star}=\frac{p_{h}}{\varrho(\phi)c_{s}^{2}}\]
The EqDF writes:
(206)\[f_{i}^{eq}(\boldsymbol{x},t)=w_{i}\left[p_h^{\star}+\left(\frac{\boldsymbol{c}_{i}\cdot\boldsymbol{u}}{c_{s}^{2}}+\frac{(\boldsymbol{c}_{i}\cdot\boldsymbol{u})^{2}}{2c_{s}^{4}}-\frac{\boldsymbol{u}\cdot\boldsymbol{u}}{2c_{s}^{2}}\right)\right]\]
and the moments are
(207)\[p_h^{\star}(\boldsymbol{x},t)=\sum_{i}f_{i}(\boldsymbol{x},t)\]
(208)\[\boldsymbol{u}(\boldsymbol{x},t)=\sum_{i}f_{i}(\boldsymbol{x},t)\boldsymbol{c}_{i}\]
Version 2: equivalent macroscopic PDE
Once the Chapman-Enskog expansion is performed, the macroscopic equations which recovered are with that EqDF are:
(209)\[\frac{\partial p^{\star}}{\partial t}+\boldsymbol{\nabla}\cdot\boldsymbol{u}=0\]
(210)\[\frac{\partial\boldsymbol{u}}{\partial t}+\boldsymbol{\nabla}\cdot(\boldsymbol{u}\boldsymbol{u})=-\boldsymbol{\nabla}(p^{\star}c_{s}^{2})+\boldsymbol{\nabla}\cdot\left[\nu(\boldsymbol{\nabla}\boldsymbol{u}+(\boldsymbol{\nabla}\boldsymbol{u})^{T})\right]\]
In Eq. (210), the pressure gradient can be expanded:
(211)\[\begin{split}-\frac{1}{\varrho}\boldsymbol{\nabla}p_{h}&=-\boldsymbol{\nabla}(p^{\star}c_{s}^{2})+\frac{1}{\varrho}\boldsymbol{F}_{p}\\{\color{red}\boldsymbol{F}_{p}}&=-\frac{p_{h}}{\varrho}\boldsymbol{\nabla}\varrho\end{split}\]
The pressure term must be corrrected with a new force term \(\boldsymbol{F}_p\) to match the gradient term of Eq. (190). In that equation, The viscous term can be expanded:
(212)\[\frac{1}{\varrho(\phi)}\boldsymbol{\nabla}\cdot\left[\varrho(\phi)\nu(\boldsymbol{\nabla}\boldsymbol{u}+(\boldsymbol{\nabla}\boldsymbol{u})^{T})\right]=\boldsymbol{\nabla}\cdot\left[\nu(\boldsymbol{\nabla}\boldsymbol{u}+(\boldsymbol{\nabla}\boldsymbol{u})^{T})\right]+\nu\left[\boldsymbol{\nabla}\boldsymbol{u}+(\boldsymbol{\nabla}\boldsymbol{u})^{T}\right]\cdot\frac{\boldsymbol{\nabla}\varrho(\phi)}{\varrho(\phi)}\]
The second term of the right-hand side corresponds to the viscous term of Eq. (210). A supplementary force \(\boldsymbol{F}_v\) has to be added to match Eq. (190).
Version 2: forces to add
To match (210) to Eq. (190), two force terms must be added in the microscopic forcing term of LBE. The first one is the pressure force pressure force \(\boldsymbol{F}_p\) defined by
(213)\[\boldsymbol{F}_{p}=-\frac{p_{h}}{\varrho}\boldsymbol{\nabla}\varrho\]
The second one is the viscosity force \(\boldsymbol{F}_v\) defined by
(214)\[\boldsymbol{F}_{v}=\nu\left[\boldsymbol{\nabla}\boldsymbol{u}+(\boldsymbol{\nabla}\boldsymbol{u})^{T}\right]\cdot\boldsymbol{\nabla}\varrho(\phi)\]
Transport equations
For standard Advection-Diffusion Equation (ADE), the EqDF are presented in Continuous Boltzmann equation. We present here the EqDF for several ADE-type equations such as the Conservative Allen-Cahn equation, the Cahn-Hilliard equation and phase-field equation for crystal growth. For those phase-field equations, a new distribution function \(g_i\) is introduced. The lattice Boltzmann equation for that \(g_i\) with the BGK collision operator writes:
(215)\[g_{i}(\boldsymbol{x}+\boldsymbol{c}_{i}\delta t,t+\delta t)=g_{i}(\boldsymbol{x},t)-\frac{1}{\tau_{g}+0.5}\left[g_{i}(\boldsymbol{x},t)-g_{i}^{eq,\iota}(\boldsymbol{x},t)\right]+\delta t\mathcal{G}_i\]
where \(\mathcal{G}_i\) is the microscopic source term. In \(g_i^{eq,\iota}\), the superscript \(\iota\) can be \(ADE,CH,CAC\). Those options are described below. In the rest of this section, the coefficient \(M_{\phi}\) involved in the transport equation is related to the collision rate \(\tau_g\) by:
(216)\[M_{\phi}=\tau_{g}c_{s}^{2}\delta t\]
where \(c_s=c/\sqrt{3}\) and \(c=\delta x/\delta t\).
\(g_i^{eq}\) for ADE with source term
Target macroscopic equation
The ADE with a source term writes:
(217)\[\partial_{t}\underbrace{\phi}_{\mathcal{M}_{0}}+\boldsymbol{\nabla}\cdot(\underbrace{\boldsymbol{u}\phi}_{\boldsymbol{\mathcal{M}}_{1}})=\boldsymbol{\nabla}\cdot\Bigl[M_{\phi}\boldsymbol{\nabla}\cdot(\underbrace{\phi\boldsymbol{I}}_{\boldsymbol{\mathcal{M}}_{2}})\Bigr]+\mathscr{S}_{\phi}\]
where the notation \(\phi\) is used here as the conserved quantity (e.g. concentration), \(\boldsymbol{u}\) is the velocity, the notation \(M_{\phi}\) represents the diffusion coefficient and \(\mathscr{S}_{\phi}\) is the macroscopic source term to be defined. Inside the bracket, \(\boldsymbol{I}\) is the identity tensor of second order. In that equation, the moments are noted:
\(\mathcal{M}_{0}\): is the moment of order 0 (scalar)
\(\boldsymbol{\mathcal{M}}_{1}\): is the moment of order 1 (vector)
\(\boldsymbol{\mathcal{M}}_{2}\): is the moment of order 2 (tensor of rank 2)
Source term
The microscopic source term is simply defined by
(218)\[\mathcal{G}_i=w_i\mathscr{S}_{\phi}\]
Its moment of order zero is \(\mathscr{S}_{\phi}\) because \(\sum_iw_i=1\).
Equilibrium with source term
The EqDF \(g_i^{eq,ADE}\) must be defined such as \(\mathcal{M}_{0}=\phi\), \(\boldsymbol{\mathcal{M}}_{1}=\boldsymbol{u}\phi\) and \(\boldsymbol{\mathcal{M}}_{2}=\phi\boldsymbol{I}\). One can prove that one possible solution is
(219)\[g_{i}^{eq,ADE}(\boldsymbol{x},t)=w_{i}\phi\left[1+\frac{\boldsymbol{u}\cdot\boldsymbol{c}_{i}}{c_{s}^{2}}\right]=w_{i}\underbrace{\phi}_{\mathcal{M}_{0}\text{ and }\boldsymbol{\mathcal{M}}_{2}}+w_{i}\frac{(\overbrace{\boldsymbol{u}\phi}^{\boldsymbol{\mathcal{M}}_{1}})\cdot\boldsymbol{c}_{i}}{c_{s}^{2}}\]
When using a source term in LBE (215), a second order of accuracy is obtained provided that the following variable change \(\overline{g}_{i}^{eq,ADE}\) is used:
(220)\[\overline{g}_{i}^{eq,ADE}(\boldsymbol{x},t)=g_{i}^{eq,ADE}-\frac{\delta t}{2}\mathcal{G}_{i}\]
Moment
After collision and streaming, the new \(\phi\) function is obtained by the moment of order zero corrected by the source term \(\mathcal{G}_{i}\delta t/2\):
(221)\[\phi=\sum_{i}g_{i}+\frac{\delta t}{2}\mathscr{S}_{\phi}\]
\(g_i^{eq}\) for Cahn-Hilliard equation
Target macroscopic equation
The Cahn-Hilliard equation writes:
(222)\[\partial_{t}\underbrace{\phi}_{\mathcal{M}_{0}}+\boldsymbol{\nabla}\cdot(\underbrace{\boldsymbol{u}\phi}_{\boldsymbol{\mathcal{M}}_{1}})=\boldsymbol{\nabla}\cdot\Bigl[M_{\phi}\boldsymbol{\nabla}\cdot(\underbrace{\mu_{\phi}\boldsymbol{I}}_{\boldsymbol{\mathcal{M}}_{2}})\Bigr]\]
where The chemical potential is a function of \(\phi\) defined by
(223)\[\mu_{\phi}=4H\phi(\phi-1)(\phi-1/2)-\zeta\underbrace{\boldsymbol{\nabla}^{2}\phi}_{\text{Laplacian}}\]
The first term of the right-hand side is the derivative of double-well and the second one involves the Laplacian of \(\phi\). Compared to the ADE (217), in Cahn-Hilliard equation (222), the source term is null \(\mathscr{S}_{\phi}=0\), and the chemical potential \(\mu_{\phi}\) replaces \(\phi\) for the moment of second order.
Equilibrium for CH equation
To design the EqDF, its moment of order zero ant its moment of first order must remain unchanged i.e. \(\mathcal{M}_{0}=\phi\) and \(\boldsymbol{\mathcal{M}}_{1}=\boldsymbol{u}\phi\). But now its moment of second order must be equal to \(\boldsymbol{\mathcal{M}}_{2}=\mu_{\phi}\boldsymbol{I}\). One solution writes:
(224)\[g_{i}^{eq,\,CH}=\mathcal{A}_{i}(\phi,\,\mu_{\phi})+w_{i}\frac{\boldsymbol{c}_{i}\cdot(\overbrace{\boldsymbol{u}\phi}^{\boldsymbol{\mathcal{M}}_{1}})}{c_{s}^{2}}\]
where
(225)\[\begin{split}\mathcal{A}_{i}(\phi,\,\mu_{\phi})=\begin{cases}
\phi-3\mu_{\phi}(1-w_{0}) & \text{if }i=0\quad\mathcal{M}_{0}\\
3w_{i}\mu_{\phi} & \text{if }i\neq0\quad\boldsymbol{\mathcal{M}}_{2}
\end{cases}\end{split}\]
The chemical potential \(\mu_{\phi}\) appears in (225). For the computation of the Laplacian which is involved in its definition (223), see Section Gradients and Laplacian.
Moment
After collision and streaming, the new \(\phi\) function is obtained by the moment of order zero
(226)\[\phi=\sum_{i}g_{i}\]
\(g_i^{eq}\) for Conservative Allen-Cahn equation
Target macroscopic equation
The Conservative Allen-Cahn equation writes:
(227)\[\partial_{t}\underbrace{\phi}_{\mathcal{M}_{0}}+\boldsymbol{\nabla}\cdot(\underbrace{\boldsymbol{u}\phi}_{\boldsymbol{\mathcal{M}}_{1}})+\boldsymbol{\nabla}\cdot\biggl[\underbrace{M_{\phi}\frac{4}{W}\phi(1-\phi)\boldsymbol{n}}_{\boldsymbol{j}_{ct}\equiv\boldsymbol{\mathcal{M}}_{1}}\biggr]=\boldsymbol{\nabla}\cdot\Bigl[M_{\phi}\boldsymbol{\nabla}\cdot(\underbrace{\phi\boldsymbol{I}}_{\boldsymbol{\mathcal{M}}_{2}})\Bigr]\]
Compared to ADE with source term, a new term appears: the last one in the left-hand side. That term involves a counter term noted \(\boldsymbol{j}_ct\) which cancels the interface movement due to the curvature. The EqDF is designed such as its moments of order 0 and 2 are unchanged compared to ADE. Its moment of order 1 must be shlightly modified to account for the new counter term. Two methods are possible for that purpose.
Method 1: equilibrium for CAC without source term (\(\mathcal{G}_i=0\))
The first method consists in modifying the EqDF by adding a new equilibrium for that counter term. That \(g_{i}^{eq,CAC}\) is the sum of the classical ADE equiblibrium plus an equilibrium designed for \(\boldsymbol{j}_{ct}\):
(228)\[\begin{split}g_{i}^{eq,CAC}=g_{i}^{eq,ADE}+g_{i}^{eq,ct}\qquad\text{with }\begin{cases}
g_{i}^{eq,ADE} & =w_{i}\phi+(\boldsymbol{u}\phi)\cdot\left(\frac{w_{i}\boldsymbol{c}_{i}}{c_{s}^{2}}\right)\\
g_{i}^{eq,ct} & =\boldsymbol{j}_{ct}\cdot\left(\frac{w_{i}\boldsymbol{c}_{i}}{c_{s}^{2}}\right)
\end{cases}\end{split}\]
where \(g_{i}^{eq,ct}\) is simply obtained by the scalar product of \(\boldsymbol{j}_{ct}\) with \(w_i\boldsymbol{c}_{i}/c_{s}^{2}\) like the advective term \(\boldsymbol{u}\phi\).
(229)\[g_{i}^{eq,CAC}(\boldsymbol{x},\,t)=\underbrace{w_{i}\phi\left[1+\frac{\boldsymbol{c}_{i}\cdot\boldsymbol{u}}{c_{s}^{2}}\right]}_{g_{i}^{eq,ADE}}+\underbrace{M_{\phi}\left[\frac{4}{W}\phi(1-\phi)\right]\boldsymbol{n}}_{\boldsymbol{j}_{ct}}\cdot\left(\frac{w_{i}\boldsymbol{c}_{i}}{c_{s}^{2}}\right)\]
The evolution of \(g_i\) follows (215) with \(\mathcal{G}_i=0\).
Method 2: equilibrium for ADE with source term (\(\mathcal{G}_i\neq0\))
Because the moments of order 0, 1 and 2 are the same than ADE, the second method uses the \(g_i^{eq,ADE}\) and considers the divergence of the counter term as a source term \(\mathcal{G}_i\):
(230)\[\overline{g}_{i}^{eq,ADE}(\boldsymbol{x},t)=w_{i}\phi\left[1+\frac{\boldsymbol{c}_{i}\cdot\boldsymbol{u}}{c_{s}^{2}}\right]-\frac{\delta t}{2}\mathcal{G}_{i}\]
where the source term is defined by
(231)\[\mathcal{G}_{i}=\frac{4}{W}\phi(1-\phi)w_{i}\boldsymbol{c}_{i}\cdot\boldsymbol{n}\]
Let us notice that in the right-hand side, the mobility coefficient does not appear compared to the expression which appears in the equilibrium (229).
Moment
Finally, after collision and streaming, the new \(\phi\) is updated by
(232)\[\phi=\sum_{i}g_{i}+\frac{\delta t}{2}\mathcal{G}_{i}\]
For the first method \(\mathcal{G}_i=0\).
Section author: Alain Cartalade