Equilibrium distribution functions for incompressible Navier-Stokes equations

The equilibrium distribution function (EqDF) \(f_i^{eq}\) is one of the most important functions in LBM because it indicates what macroscopic equation the algorithm simulates. The equilibrium distribution function \(f_i^{eq}\) is involved in the collision operator \(\Omega\) in the discrete LBE:

(183)\[f_{i}(\boldsymbol{x}+\boldsymbol{c}_{i}\delta t,t+\delta t)=f_{i}(\boldsymbol{x},t)+\Omega_{i}(f_{i},f_{i}^{eq})+\mathcal{F}_i\]

One of the most attractive aspect of LBM, is to use the same evolution Eq. (183) and simply change the definition of \(f_i^{eq}\) to simulate either the Navier-Stokes equations or the Conservative Allen-Cahn equation. For low Mach Navier-Stokes (NS) and standard Advection-Diffusion Equation (ADE), the EqDF are presented in Continuous Boltzmann equation. Here we present those for incompressible two-phase NS and several variations of transport equation such as the Cahn-Hilliard equation and the Conservative Allen-Cahn one.

Equilibrium distribution functions in LBM_Saclay

The equilibrium distribution functions are problem dependent. They are defined for each kernel in the LBMScheme_Kernel_Name.h file of folder src/kernels. For example for kernel NSAC_Comp, the equilibrium distribution functions for all PDEs are written in LBMScheme_NS_AC_Comp.h file.

For low Mach Navier-Stokes (NS), the EqDF are presented in Continuous Boltzmann equation. Here are presented some modifications for simulating incompressible single-phase and two-phase flows.

Target macroscopic PDEs

First we remind the target macroscopic PDEs we want to simulate:

(184)\[\boldsymbol{\nabla}\cdot\boldsymbol{u}=0\]
(185)\[\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)\) and the total dynamic viscosity are two functions of position and time with the phase-field \(\phi\). The first one is obtained by a linear interpolation of \(\rho_0\) and \(\rho_1\):

(186)\[\varrho(\phi)=\rho_0\phi(\boldsymbol{x},t)+\rho_1(1-\phi(\boldsymbol{x},t))\]

and the second with a harmonic interpolation of \(\eta_0\) and \(\eta_1\):

(187)\[\frac{1}{\eta(\phi)}=\frac{\phi}{\eta_{1}}+\frac{1-\phi}{\eta_{2}}\]

Standard \(f_i^{eq}\) for single-phase incompressible NS

\(f_i^{eq}\) for single-phase iNS

For simulating incompressible NS equations for single-phase flows, the EqDF writes [1] & [2]:

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

(189)\[p_{h}(\boldsymbol{x},t)=\sum_{i}f_{i}(\boldsymbol{x},t)\]

Next, the moment of order one is the impulsion:

(190)\[\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 (188) yields to the equivalent system of equations:

(191)\[\frac{\partial{\color{red}p_h}}{\partial t}+{\color{red}\rho_{0}c_{s}^{2}}\boldsymbol{\nabla}\cdot\boldsymbol{u}=0\]
(192)\[{\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. (191) 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 (186). The first method is a direct extension of (188). The second one is a slight modification which introduces the dimensionless pressure \(p_h^{\star}\).

Version 1 for variable density

Version 1 for variable density

That first version is a direct adaptation of Eq. (188) where the variable density \(\varrho(\phi)\) simply replaces the constant density \(\rho_0\):

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

(194)\[p_{h}(\boldsymbol{x},t)=\sum_{i}f_{i}(\boldsymbol{x},t)\]

and its moment of order 1:

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

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

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

(198)\[\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. (196) can be rewritten:

(199)\[\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. (196) 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:

(200)\[S_p=\boldsymbol{u}\cdot\boldsymbol{\nabla}[\varrho(\phi)c_s^2]\]

in the LBE.

Version 2 for variable density

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 [3]:

(201)\[p_h^{\star}=\frac{p_{h}}{\varrho(\phi)c_{s}^{2}}\]

The EqDF writes:

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

(203)\[p_h^{\star}(\boldsymbol{x},t)=\sum_{i}f_{i}(\boldsymbol{x},t)\]
(204)\[\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:

(205)\[\frac{\partial p^{\star}}{\partial t}+\boldsymbol{\nabla}\cdot\boldsymbol{u}=0\]
(206)\[\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. (206), the pressure gradient can be expanded:

(207)\[\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. (185). In that equation, The viscous term can be expanded:

(208)\[\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. (206). A supplementary force \(\boldsymbol{F}_v\) has to be added to match Eq. (185).

Version 2: forces to add

To match (206) to Eq. (185), two force terms must be added in the microscopic forcing term of LBE (e.g. see [4]). The first one is the pressure force pressure force \(\boldsymbol{F}_p\) defined by

(209)\[\boldsymbol{F}_{p}=-\frac{p_{h}}{\varrho}\boldsymbol{\nabla}\varrho\]

The second one is the viscosity force \(\boldsymbol{F}_v\) defined by

(210)\[\boldsymbol{F}_{v}=\nu\left[\boldsymbol{\nabla}\boldsymbol{u}+(\boldsymbol{\nabla}\boldsymbol{u})^{T}\right]\cdot\boldsymbol{\nabla}\varrho(\phi)\]

References

Section author: Alain Cartalade