###### Abstract

The lattice Boltzmann method, now widely used for a variety of applications,
has also been extended to model multi-phase flows through different
formulations. While already applied to many different configurations in the
low Weber and Reynolds number regimes, applications to higher Weber/Reynolds
numbers or larger density/viscosity ratios are still the topic of active
research. In this study, through a combination of the decoupled phase-field
formulation –conservative Allen-Cahn equation– and a cumulants-based collision
operator for a low-Mach pressure-based flow solver, we present an algorithm
that can be used for higher Reynolds/Weber numbers. The algorithm is validated
through a variety of test-cases, starting with the Rayleigh-Taylor instability
both in 2-D and 3-D, followed by the impact of a droplet on a liquid sheet. In
all simulations, the solver is shown to correctly capture the dynamics of the
flow and match reference results very well. As the final test-case, the solver
is used to model droplet splashing on a thin liquid sheet in 3-D with a
density ratio of 1000 and kinematic viscosity ratio of 15 –matching the
water/air system– at We=8000 and Re=1000. The results show that the solver
correctly captures the fingering instabilities at the crown rim and their
subsequent breakup, in agreement with experimental and numerical observations
reported in the literature.

###### keywords:

lattice Boltzmann method; multi-phase flows; Conservative Allen-Cahn; Phase
field.

xx 1 5 Received: date; Accepted: date; Published: date Lattice Boltzmann
solver for multi-phase flows: Application to high Weber and Reynolds numbers
S.A. Hosseini 1,2*, H. Safari 1 and D. Thévenin 1 S.A. Hosseini, H. Safari and
D. Thévenin Correspondence<EMAIL_ADDRESS>
lbLBlattice Boltzmann lbmLBMlattice Boltzmann method acACAllen-Cahn chCHCahn-
Hilliard pdePDEpartial differential equation edfEDFequilibrium distribution
function dvbeDVBEdiscrete velocity Boltzmann equation nsNSNavier-Stokes
acmACMartificial compressibility method eosEoSequation of state mrtMRTmultiple
relaxation time srtSRTsingle relaxation time rhsRHSright hand side
dfDFdistribution function ceCEChapman Enskog

##  1 Introduction

The is a discrete solver for the so-called , initially developed as an
alternative to classical solvers for the incompressible hydrodynamic regime
Krüger et al. (2017); Guo and Shu (2013). Due to the simplicity of the
algorithm, low computational cost of discrete time-evolution equations and
locality of non-linear terms and boundary conditions it has rapidly grown over
the past decades Succi (2002). It is worth noting that while intended for the
incompressible regime, the formally solves the compressible isothermal
equations at a reference temperature. While originally tied to the considered
flow’s temperature, in the context of the solver, the reference temperature is
a numerical parameter allowing to control convergence and consistency of the
results Krüger et al. (2017). The weak compressibility in the formulation
along with the parabolic nature of the governing the evolution of pressure –as
opposed to Chorin’s original – have made the scheme efficient and applicable
to unsteady flows Chorin (1997). Although originally used for single-phase
flows it has since been extended to multi-phase, multi-species and
compressible flows.

While generally based on diffuse-interface formulations, solvers for multi-
phase flows can be categorized as pertaining to one of three major categories:
(a) pseudo-potential Shan and Chen (1993, 1994), (b) free energy Swift et al.
(1996, 1995) and (c) phase-field. Other types of formulations can also be
found in the literature, however they are not as widely spread and/or
developed as these three.
In the context of the free energy formulation, the expression for the non-
local non-ideal pressure tensor is found through the free energy functional.
The appropriate pressure tensor is then introduced into the solver via a
moment-matching approach assigning coefficients to different terms in the
Swift et al. (1995). The interesting point that makes this formulation
consistent and differentiates it from the generic double-well potential-based
Cahn-Hilliard formulation, is that in the minimization process of the free
energy, the is explicitly considered. It is interesting to note that, as is
the case for the pseudo-potential formulation, the explicit intervention of
the within the free functional makes the thickness of the interface tied to
physical parameters, i.e. surface tension, density ratio, etc. As a
consequence, the choice of the and/or tuning of the coefficients in the is a
method of choice to widen the area of accessible density ratios. This approach
was later extended by introducing non-ideal components of the pressure tensor
via an external body force. Introducing these effects with a body force made
the scheme more stable by reducing Galilean invariance issues tied to the
third-order moments of the Wagner and Li (2006).
The pseudo-potential formulation follows more of a bottom-up approach in
introducing non-ideal dynamics into the solver. It follows the general
philosophy of the Boltzmann-Vlasov equation, introducing a non-local
_potential_ to account for non-ideal effects. While the original formulation
relied on what was termed an effective _density_ , actual were introduced into
the pseudo-potential in Kupershtokh et al. (2009); Yuan and Schaefer (2006).
Apart from thermodynamic consistency, the possibility of using different
allowed for higher density ratios to be modeled. As the free energy
formulation, this model is limited to lower Weber number regimes because it
naturally comes with large surface tension values. While more advanced models
allow for independent tuning of the surface tension Sbragaglia et al. (2007),
the spectrum of values covered by the model is rather limited and barely
allows for variations of one order of magnitude Li and Luo (2013).
The last category is based on the free energy functional minimization
approach, just like the free energy approach. However, contrary to the latter,
the surface and bulk energies used in the minimization process are those of a
generic double-well potential Fakhari and Rahimian (2010), allowing to
decouple –among other parameters– the interface thickness from the fluid
physical properties. Another consequence of this choice of functional is a
partial loss of thermodynamic consistency, making the extension of the
formulation to more complex physics such as thermal flows, compressible flows,
or acoustics less straightforward, although a number of attempts have been
documented in the literature Safari et al. (2013, 2014); Yazdi et al. (2018).
Nevertheless, it has been observed to be very effective and robust for multi-
phase flows in the incompressible regime and readily able to deal with larger
Weber numbers. For a more comprehensive overview of the developments of such
models, interested readers are referred to Wang et al. (2019). It is also
worth noting that approaches relying on explicit tracking of the interface
with a consistent energy functional making use of the non-ideal have also been
proposed as ways to improve the stability of the original free energy
formulation He et al. (1999); Inamuro et al. (2004).

Over the past decades a lot of efforts have been put in developing phase
field-based solvers for various applications Safari et al. (2014);
Amirshaghaghi et al. (2016, 2018). Given that in such formulations the local
density is a dependent variable –on the local value of the order parameter,
they have to be coupled to a modified form of the solver for the flow usually
referred to as the incompressible formulation. The so-called low-Mach
formulation is mostly based on the modified distribution function introduced
in He et al. (1999) where the pressure is the zeroth-order moment of the
distribution function. This flow solver has been combined with different forms
of interface tracking formulations, e.g. , conservative or to model multi-
phase flows. The aim of the present study is to introduce a multi-phase solver
relying on the pressure-based formulation of He et al. (1999) and a
realization –for the flow solver– coupled with a solver for the conservative .
The use of the collision operator in cumulants space along with the decoupled
interface tracking allow for simulations in the high Reynolds and Weber
regimes. After a brief introduction of the model, it will be used to simulate
a variety of test-cases proving its ability to reproduce correct physics and
its robustness. It is worth noting that all models were implemented in our in-
house multi-physics solver ALBORZ Hosseini (2020).

##  2 Theoretical background

###  2.1 Target macrosopic system

As briefly stated in the introduction, the aim of the present work is to solve
the multi-phase flow equations within the context of a diffuse interface
formulation in the limit of the incompressible regime, where interface
dynamics are followed and accounted for via an additional indicator field,
$\phi$. As such, at the macroscopic level the low Mach equations are targeted:

${\partial_{t}\rho u_{i}+\partial_{j}\rho
u_{i}u_{j}+\partial_{j}\sigma_{ij}+\mu_{\phi}\partial_{i}\phi+F_{b,i}=0,}$ (1)

where $u_{i}$ is the fluid velocity, $\rho$ the fluid density and $F_{b,i}$
designates external body forces. The stress tensor $\sigma_{ij}$ is defined
as:

${\sigma_{ij}=p_{h}\delta_{ij}-\eta\left(\partial_{i}u_{j}+\partial_{j}u_{i}\right)+\frac{2}{3}\eta\partial_{k}u_{k}\delta_{ij},}$
(2)

where $\eta$ is the fluid dynamic viscosity, tied to the kinematic viscosity
$\nu$ as $\eta=\rho\nu$, and $p_{h}$ the hydrodynamic pressure. The chemical
potential $\mu_{\phi}$ is defined as:

${\mu_{\phi}=2\beta\phi\left(\phi-1\right)\left(2\phi-1\right)\kappa\Delta\phi,}$
(3)

where $\Delta=\nabla^{2}$ is the Laplacian operator and $\beta$ and $\kappa$
are parameters specific to the formulation. It must be noted that the second
term on the of Equation 1 accounts for surface tension effects. For the sake
of clarity the free parameters will be detailed in the next paragraph.
The interface is tracked using the conservative equation, where the order
parameter $\phi$ evolves as Sun and Beckermann (2007); Chiu and Lin (2011):

${\partial_{t}\phi+\partial_{i}u_{i}\phi-\partial_{i}M\left[\partial_{i}\phi-
n_{i}\frac{4\phi(1-\phi)}{W}\right]=0,}$ (4)

where the parameter $\phi$ takes on values between 0 and 1, $M$ is mobility,
$W$ is the interface thickness and $n_{i}$ is the unit normal to the interface
obtained as:

${n_{i}=\frac{\partial_{i}\phi}{\lvert\lvert\nabla\phi\lvert\lvert}.}$ (5)

The interfaces can be found through iso-surfaces of the order parameter, i.e.
$\phi=1/2$. To recover the correct surface tension, the free parameters
appearing in the chemical potential, i.e. $\kappa$ and $\beta$ are tied to the
surface tension $\sigma$ and interface thickness $W$ in the equation via
$\beta=12\sigma/W$ and $\kappa=3\sigma W/2$.

###  2.2  formulation for the conservative phase-field equation

The conservative equation can readily be recovered by appropriately defining
the discrete equilibrium state and relaxation coefficient in the advection-
diffusion model:

$\partial_{t}g_{\alpha}+c_{\alpha,i}\partial_{i}g_{\alpha}+\mathcal{S}_{\alpha}=\Omega^{\phi}_{\alpha},$
(6)

where $g_{\alpha}$ and $c_{\alpha}$ are the populations and velocities in the
discrete velocity kinetic model and the collision operator is defined as:

$\Omega^{\phi}_{\alpha}=\frac{1}{\tau_{\phi}}\left(g_{\alpha}^{(eq)}-g_{\alpha}\right).$
(7)

The is defined as:

${g_{\alpha}^{(eq)}=w_{\alpha}\phi\sum_{n=0}^{2}\frac{1}{n!c_{s}^{2n}}\mathcal{H}_{n}:a^{(eq)}_{n},}$
(8)

where $\mathcal{H}_{n}$ and $a_{n}^{(eq)}$ are the Hermite polynomial and
coefficient of order $n$, $c_{s}$ the lattice sound speed and $w_{\alpha}$
weights tied to each discrete velocity (resulting from the Gauss-Hermite
quadrature). The expressions for these polynomials and corresponding
coefficients are listed in Appendix A. The source term in Equation 6 is
defined as Fakhari et al. (2017):

$\mathcal{S}_{\alpha}=w_{\alpha}\mathcal{H}_{i}n_{i}\frac{4\phi(1-\phi)}{W}.$
(9)

Given that the source term affects the first-order moment –a non-conserved
moment of the distribution function– the distribution function is tied to the
phase parameter as:

$\phi=\sum_{\alpha}g_{\alpha}.$ (10)

The relaxation coefficient is fixed as:

$\tau_{\phi}=\frac{M}{c_{s}^{2}}.$ (11)

After integration in space/time the now-famous collision-streaming form can be
recovered:

${\bar{g}_{\alpha}\left(x+c_{\alpha}\delta_{t},t+\delta_{t}\right)=\left(1-\frac{\delta_{t}}{\bar{\tau}_{\phi}}\right)\bar{g}_{\alpha}\left(x,t\right)+\frac{\delta_{t}}{\bar{\tau}_{\phi}}g^{(eq)}_{\alpha}\left(x,t\right)+\delta_{t}\bar{\mathcal{S}}_{\alpha}\left(x,t\right),}$
(12)

where the source term takes on a new form, i.e.:

$\bar{\mathcal{S}}_{\alpha}=\left(1-\frac{1}{2\tau_{\phi}}\right)w_{\alpha}\mathcal{H}_{i}n_{i}\frac{4\phi(1-\phi)}{W},$
(13)

and:

${\bar{\tau}_{\phi}=\tau_{\phi}+\frac{\delta_{t}}{2}.}$ (14)

It is also worth noting that the derivatives of the order parameter appearing
in the various discrete time-evolution equations are computed using isotropic
finite differences, i.e.:

$\partial_{i}\phi=\frac{1}{c_{s}^{2}}\sum_{\alpha}w_{\alpha}c_{\alpha,i}\phi(x+c_{\alpha}),$
(15)

and:

$\partial^{2}_{i}\phi=\frac{2}{c_{s}^{2}}\sum_{\alpha}w_{\alpha}\left[\phi(x+c_{\alpha})-\phi(x)\right].$
(16)

While the present work makes use of a second-order , one must note that the
same macroscopic , i.e. Equation 4, can also be recovered by using a first-
order and an additional correction term of the following form Wang et al.
(2016):

$C_{\alpha}=\frac{w_{\alpha}}{c_{s}^{2}}\mathcal{H}_{i}\partial_{t}\phi
u_{i},$ (17)

where, as for Equation 13, post-discretization it changes into:

$\bar{C}_{\alpha}=\left(1-\frac{1}{2\tau_{\phi}}\right)\frac{w_{\alpha}}{c_{s}^{2}}\mathcal{H}_{i}\partial_{t}\phi
u_{i}.$ (18)

Such correction terms were first introduced in the context of advection-
diffusion solvers Chopard et al. (2009) and further extended to non-linear
equations in the same context Hosseini et al. (2019). Detailed derivation and
multi-scale analyses are readily available in the literature, e.g. Zu et al.
(2020).

###  2.3  model for flow field

The flow solver kinetic model follows the low-Mach formulation used, among
other sources, in Lee and Lin (2003); Hosseini et al. (2019, 2020) and based
on the original model introduced in He et al. (1999):

$\partial_{t}f^{{}^{\prime}}_{\alpha}+c_{\alpha,i}\partial_{i}f^{{}^{\prime}}_{\alpha}=\Omega_{\alpha}+\Xi_{\alpha},$
(19)

where the collision operator is:

$\Omega_{\alpha}=\frac{1}{\tau}\left({f_{\alpha}^{(eq)}}^{{}^{\prime}}-f^{{}^{\prime}}_{\alpha}\right),$
(20)

and $\Xi_{\alpha}$ is defined as:

$\Xi_{\alpha}=c_{s}^{2}\left(\frac{f_{\alpha}^{(eq)}}{\rho}-w_{\alpha}\right)\left(c_{\alpha,i}-u_{i}\right)\partial_{i}\rho+w_{\alpha}c_{s}^{2}\rho\partial_{i}u_{i}+\left(F_{b,i}+F_{s,i}\right)\left(c_{\alpha,i}-u_{i}\right)\frac{f_{\alpha}^{(eq)}}{\rho},$
(21)

and the relaxation coefficient $\tau$ is tied to the fluid kinematic viscosity
$\nu$ as:

$\tau=\frac{\nu}{c_{s}^{2}}.$ (22)

The forces $F_{b,i}$ and $F_{s,i}$ represent respectively external body forces
and surface tension, i.e.:

${F_{s,i}=\mu_{\phi}\partial_{i}\phi.}$ (23)

The modified distribution function, $f^{{}^{\prime}}_{\alpha}$, is defined as:

$f^{{}^{\prime}}_{\alpha}=w_{\alpha}p_{h}+c_{s}^{2}\left(f_{\alpha}-w_{\alpha}\rho\right),$
(24)

where $f_{\alpha}$ is the classical iso-thermal distribution function. The
modified equilibrium follows the same logic and is defined as:

${{f^{(eq)}_{\alpha}}^{{}^{\prime}}=w_{\alpha}p_{h}+w_{\alpha}\rho
c_{s}^{2}\sum_{n=1}^{2}\frac{1}{n!c_{s}^{2n}}\mathcal{H}_{n}:a^{(eq)}_{n}.}$
(25)

The density is tied to the order parameter as:

$\rho=\rho_{l}+\left(\rho_{h}-\rho_{l}\right)\phi,$ (26)

where $\rho_{h}$ and $\rho_{l}$ are respectively the densities of the heavy
and light fluid. For a detailed analysis of the macroscopic equations
recovered by this model and the derivation of the discrete equations,
interested readers are referred to Hosseini et al. (2019); Hosseini (2020). In
the context of the present study the low-Mach model is wrapped in a moments-
based formulation where the post-collision populations
$f^{{}^{\prime}*}_{\alpha}$, to be streamed as:

${{f^{{}^{\prime}}}_{\alpha}\left(x+c_{\alpha}\delta_{t},t+\delta_{t}\right)={f^{{}^{\prime}*}}_{\alpha}\left(x,t\right),}$
(27)

are computed as:

${f^{{}^{\prime}*}}_{\alpha}=\rho
c_{s}^{2}f_{\alpha}^{p*}+\frac{\delta_{t}}{2}\Xi_{\alpha}.$ (28)

The post-collision pre-conditioned population $f_{\alpha}^{p*}$ is:

${f^{p*}_{\alpha}}=\mathcal{C}^{-1}\left(\mathcal{I}-\mathcal{W}\right)\mathcal{K}^{p}+\mathcal{C}^{-1}\mathcal{W}\mathcal{K}^{p},$
(29)

where $\mathcal{C}$ is the moments transform matrix –from pre-conditioned
populations to the target momentum space, $\mathcal{I}$ the identity matrix
and $\mathcal{W}$ the diagonal relaxation frequency matrix. Following Geier et
al. (2020) prior to transformation to momentum space the populations are pre-
conditioned as:

$f^{p}_{\alpha}=\frac{1}{\rho
c_{s}^{2}}f^{{}^{\prime}}_{\alpha}+\frac{\delta_{t}}{2\rho
c_{s}^{2}}\Xi_{\alpha}.$ (30)

This pre-conditioning accomplishes two tasks, namely normalizing the
populations with the density and thus eliminating the density-dependence of
the moments and introducing the first half of the source term. As such the
moments $\mathcal{K}^{p}$ are computed as:

$\mathcal{K}^{p}_{\beta}=\mathcal{C}_{\alpha\beta}f^{p}_{\alpha}.$ (31)

The transformation from s to cumulants is carried out using the steps
suggested in Geier et al. (2015), which allows for a more efficient algorithm.
The s are first transformed into central moments:

${\widetilde{\Pi}_{\beta}^{p}=\sum_{\alpha}{\left(c_{\alpha,x}-u_{x}\right)}^{n_{x}}{\left(c_{\alpha,y}-u_{y}\right)}^{n_{y}}{\left(c_{\alpha,z}-u_{z}\right)}^{n_{z}}f_{\alpha}^{p}.}$
(32)

where here $\beta=x^{n_{x}}y^{n_{y}}z^{n_{z}}$. The central moments are then
transformed into the corresponding cumulants using the following relations:

$\displaystyle\mathcal{K}^{p}_{x}=$ $\displaystyle\widetilde{\Pi}_{x}^{p},$
(33a) $\displaystyle\mathcal{K}^{p}_{xy}=$
$\displaystyle\widetilde{\Pi}_{xy}^{p},$ (33b)
$\displaystyle\mathcal{K}^{p}_{x^{2}}=$
$\displaystyle\widetilde{\Pi}_{x^{2}}^{p},$ (33c)
$\displaystyle\mathcal{K}^{p}_{xy^{2}}=$
$\displaystyle\widetilde{\Pi}_{xy^{2}}^{p},$ (33d)
$\displaystyle\mathcal{K}^{p}_{xyz}=$
$\displaystyle\widetilde{\Pi}_{xyz}^{p},$ (33e)
$\displaystyle\mathcal{K}^{p}_{x^{2}yz}=$
$\displaystyle\widetilde{\Pi}_{x^{2}yz}^{p}-\left[\widetilde{\Pi}_{x^{2}}^{p}\widetilde{\Pi}_{yz}^{p}+2\widetilde{\Pi}_{xy}^{p}\widetilde{\Pi}_{xz}^{p}\right],$
(33f) $\displaystyle\mathcal{K}^{p}_{x^{2}y^{2}}=$
$\displaystyle\widetilde{\Pi}_{x^{2}y^{2}}^{p}-\left[\widetilde{\Pi}_{x^{2}}^{p}\widetilde{\Pi}_{y^{2}}^{p}+2{(\widetilde{\Pi}_{xy}^{p})}^{2}\right],$
(33g) $\displaystyle\mathcal{K}^{p}_{xy^{2}z^{2}}=$
$\displaystyle\widetilde{\Pi}_{xy^{2}z^{2}}^{p}-\left[\widetilde{\Pi}_{z^{2}}^{p}\widetilde{\Pi}_{xy^{2}}^{p}+\widetilde{\Pi}_{y^{2}}^{p}\widetilde{\Pi}_{xz^{2}}^{p}+4\widetilde{\Pi}_{yz}^{p}\widetilde{\Pi}_{xyz}^{p}+2(\widetilde{\Pi}_{xz}^{p}\widetilde{\Pi}_{y^{2}z}^{p}+\widetilde{\Pi}_{xy}^{p}\widetilde{\Pi}_{yz^{2}}^{p})\right],$
(33h) $\displaystyle\mathcal{K}^{p}_{x^{2}y^{2}z^{2}}=$
$\displaystyle\widetilde{\Pi}_{x^{2}y^{2}z^{2}}^{p}-\left[4{(\widetilde{\Pi}_{xyz}^{p})}^{2}+\widetilde{\Pi}_{x^{2}}^{p}\widetilde{\Pi}_{y^{2}z^{2}}^{p}+\widetilde{\Pi}_{y^{2}}^{p}\widetilde{\Pi}_{x^{2}z^{2}}^{p}+\widetilde{\Pi}_{z^{2}}^{p}\widetilde{\Pi}_{x^{2}y^{2}}^{p}+4(\widetilde{\Pi}_{xy}^{p}\widetilde{\Pi}_{x^{2}yz}^{p}+\right.$
$\displaystyle\left.\widetilde{\Pi}_{xz}^{p}\widetilde{\Pi}_{xy^{2}z}^{p}+\widetilde{\Pi}_{xy}^{p}\widetilde{\Pi}_{xyz^{2}}^{p}+2(\widetilde{\Pi}_{xy^{2}}^{p}\widetilde{\Pi}_{xz^{2}}^{p}+\widetilde{\Pi}_{x^{2}y}^{p}\widetilde{\Pi}_{yz^{2}}^{p}+\widetilde{\Pi}_{x^{2}z}^{p}\widetilde{\Pi}_{y^{2}z}^{p}))+\right.$
$\displaystyle\left.(16\widetilde{\Pi}_{xy}^{p}\widetilde{\Pi}_{xz}^{p}\widetilde{\Pi}_{yz}^{p}+4({(\widetilde{\Pi}_{xz}^{p})}^{2}\widetilde{\Pi}_{y^{2}}^{p}+{(\widetilde{\Pi}_{yz}^{p})}^{2}\widetilde{\Pi}_{x^{2}}^{p}+{(\widetilde{\Pi}_{xy}^{p})}^{2}\widetilde{\Pi}_{z^{2}}^{p})+2\widetilde{\Pi}_{x^{2}}^{p}\widetilde{\Pi}_{y^{2}}^{p}\widetilde{\Pi}_{z^{2}}^{p})\right].$
(33i)

The remainder of the moments can be easily obtained via permutation of the
indices. The collision process is performed in cumulant space according to
Geier et al. (2015). The fluid viscosity is controlled via the collision
factor related to second-order cumulants (e.g. $\mathcal{K}^{p}_{xy}$,
$\mathcal{K}^{p}_{x^{2}}-\mathcal{K}^{p}_{y^{2}}$,
$\mathcal{K}^{p}_{x^{2}}-\mathcal{K}^{p}_{z^{2}}$ etc). The rest of the
collision factors are set to unity for simplicity. Once the collision step has
been applied, cumulants are transformed back into central moments as:

$\displaystyle\widetilde{\Pi}_{x}^{p*}=$ $\displaystyle\mathcal{K}^{p*}_{x},$
(34a) $\displaystyle\widetilde{\Pi}_{xy}^{p*}=$
$\displaystyle\mathcal{K}^{p*}_{xy},$ (34b)
$\displaystyle\widetilde{\Pi}_{x^{2}}^{p*}=$
$\displaystyle\mathcal{K}^{p*}_{x^{2}},$ (34c)
$\displaystyle\widetilde{\Pi}_{xy^{2}}^{p*}=$
$\displaystyle\mathcal{K}^{p*}_{xy^{2}},$ (34d)
$\displaystyle\widetilde{\Pi}_{xyz}^{p*}=$
$\displaystyle\mathcal{K}^{p*}_{xyz},$ (34e)
$\displaystyle\widetilde{\Pi}_{x^{2}yz}^{p*}=$
$\displaystyle\mathcal{K}^{p*}_{x^{2}yz}+\left[\widetilde{\Pi}^{p*}_{x^{2}}\widetilde{\Pi}^{p*}_{yz}+2\widetilde{\Pi}^{p*}_{xy}\widetilde{\Pi}^{p*}_{xz}\right],$
(34f) $\displaystyle\widetilde{\Pi}_{x^{2}y^{2}}^{p*}=$
$\displaystyle\mathcal{K}^{p*}_{x^{2}y^{2}}+\left[\widetilde{\Pi}_{x^{2}}^{p*}\widetilde{\Pi}_{y^{2}}^{p*}+2{(\widetilde{\Pi}_{xy}^{p*})}^{2}\right],$
(34g) $\displaystyle\widetilde{\Pi}_{xy^{2}z^{2}}^{p*}=$
$\displaystyle\mathcal{K}^{p*}_{xy^{2}z^{2}}+\left[\widetilde{\Pi}^{p*}_{z^{2}}\widetilde{\Pi}^{p*}_{xy^{2}}+\widetilde{\Pi}^{p*}_{y^{2}}\widetilde{\Pi}^{p*}_{xz^{2}}+4\widetilde{\Pi}^{p*}_{yz}\widetilde{\Pi}^{p*}_{xyz}+2(\widetilde{\Pi}^{p*}_{xz}\widetilde{\Pi}^{p*}_{y^{2}z}+\widetilde{\Pi}^{p*}_{xy}\widetilde{\Pi}^{p*}_{yz^{2}})\right],$
(34h) $\displaystyle\widetilde{\Pi}^{p*}_{x^{2}y^{2}z^{2}}=$
$\displaystyle\mathcal{K}^{p*}_{x^{2}y^{2}z^{2}}+\left[4{(\widetilde{\Pi}^{p*}_{xyz})}^{2}+\widetilde{\Pi}^{p*}_{x^{2}}\widetilde{\Pi}^{p*}_{y^{2}z^{2}}+\widetilde{\Pi}^{p*}_{y^{2}}\widetilde{\Pi}^{p*}_{x^{2}z^{2}}+\widetilde{\Pi}^{p*}_{z^{2}}\widetilde{\Pi}^{p*}_{x^{2}y^{2}}+4(\widetilde{\Pi}^{p*}_{xy}\widetilde{\Pi}^{p*}_{x^{2}yz}+\right.$
$\displaystyle\left.\widetilde{\Pi}^{p*}_{xz}\widetilde{\Pi}^{p*}_{xy^{2}z}+\widetilde{\Pi}^{p*}_{xy}\widetilde{\Pi}^{p*}_{xyz^{2}}+2(\widetilde{\Pi}^{p*}_{xy^{2}}\widetilde{\Pi}^{p*}_{xz^{2}}+\widetilde{\Pi}^{p*}_{x^{2}y}\widetilde{\Pi}^{p*}_{yz^{2}}+\widetilde{\Pi}^{p*}_{x^{2}z}\widetilde{\Pi}^{p*}_{y^{2}z}))-\right.$
$\displaystyle\left.(16\widetilde{\Pi}^{p*}_{xy}\widetilde{\Pi}^{p*}_{xz}\widetilde{\Pi}^{p*}_{yz}+4({(\widetilde{\Pi}^{p*}_{xz})}^{2}\widetilde{\Pi}^{p*}_{y^{2}}+{(\widetilde{\Pi}^{p*}_{yz})}^{2}\widetilde{\Pi}^{p*}_{x^{2}}+{(\widetilde{\Pi}^{p*}_{xy})}^{2}\widetilde{\Pi}^{p*}_{z^{2}})+2\widetilde{\Pi}^{p*}_{x^{2}}\widetilde{\Pi}^{p*}_{y^{2}}\widetilde{\Pi}^{p*}_{z^{2}})\right].$
(34i)

After this step, the post-collision central moments can be readily transformed
back to populations. All transforms presented here and upcoming simulations
are based on the D3Q27 stencil. It must also be noted that the following set
of 27 moments are used as moments basis:

$\beta\in\\{0,x,y,z,xy,xz,yz,x^{2}-y^{2},x^{2}-z^{2},x^{2}+y^{2}+z^{2},\\\
xy^{2}+xz^{2},xyz,xy^{2}-xz^{2},x^{2}+yz^{2},x^{2}z+y^{2}z,x^{2}y-yz^{2},x^{2}z-y^{2}z,x^{2}y^{2}-2x^{2}z^{2}+y^{2}z^{2},\\\
x^{2}y^{2}+x^{2}z^{2}-2y^{2}z^{2},x^{2}y^{2}+x^{2}z^{2}+y^{2}z^{2},x^{2}yz,xy^{2}z,xyz^{2},x^{2}y^{2}z,x^{2}yz^{2},xy^{2}z^{2},x^{2}y^{2}z^{2}\\},$
(35)

where $\beta=x^{2}-y^{2}$ stands for a central moment of the form
$\widetilde{\Pi}^{p}_{x^{2}}-\widetilde{\Pi}^{p}_{y^{2}}$. Previous systematic
studies of the flow solver have shown second-order convergence under diffusive
scaling Hosseini et al. (2019).

##  3 Numerical applications

The proposed numerical method will be validated through different test-cases
in the present section. All results and simulation parameters are reported in
units, i.e. non-dimensionalized with time-step, grid-size and heavy fluid
density.

###  3.1 Static droplet: Surface tension measurement

As a first test, to validate the hydrodynamics of the model, we consider the
case of a static droplet in a rectangular domain with periodic boundaries all
around. All cases consist of a domain of size $256\times 256$ filled with a
light fluid. A _droplet_ of the heavier fluid is placed at the center of the
domain. Simulations are pursued till the system converges. The pressure
difference between the droplet and surrounding lighter fluid is then
extracted. Using Laplace’s law, i.e.:

$\Delta P=\frac{\sigma}{r},$ (36)

where $\Delta P$ is the pressure difference and $r$ the droplet radius, one
can readily obtain the effective surface tension. Three different surface
tensions, i.e. $\sigma=1\times 10^{-1}$, $1\times 10^{-3}$ and $1\times
10^{-6}$, along with four different droplet radii, i.e. $r=25$, $30$, $35$ and
$45$, were considered here. The obtained results are shown in Figure 1. The
results presented here consider a density ratio of 20 and non-dimensional
viscosity of 0.1.

Figure 1: Changes of pressure difference around droplet for different surface
tensions and droplet radii. Red, blue and black symbols illustrate results
from the present study with respectively $\sigma=10^{-1},10^{-3}$ and
$10^{-6}$.

It is readily observed that the model satisfies Laplace’s law and recovers
correct surface tensions. Furthermore, it is seen that it can span a wide
range of surface tensions, as opposed to other classes of multi-phase solvers
such as the free energy or pseudo-potential formulations Qin et al. (2018);
Mazloomi M et al. (2015) and maintain relatively low spurious currents. For
example, at a density ratio of 1000 and $\sigma=10^{-3}$, the spurious
currents were found to be only of the order of $10^{-6}$, in strong contrast
with the previously cited approaches.

###  3.2 Rayleigh-Taylor instability

The Rayleigh-Taylor instability is a well-known and widely studied gravity-
driven effect occurring when a layer of a heavier fluid lies on top of another
layer of a lighter fluid Yang et al. (2018a, b); Rahmat et al. (2014). A
perturbation at the interface between the two fluids causes the heavier one to
penetrate into the lighter fluid. In general, the dynamics of this system are
governed by two non-dimensional parameters, namely the Atwood and Reynolds
numbers. The former is defined as:

$\hbox{At}=\frac{\rho_{h}-\rho_{l}}{\rho_{h}+\rho_{l}},$ (37)

while the latter is:

$\hbox{Re}=\frac{\rho_{h}U^{*}L}{\mu_{h}},$ (38)

where $\rho_{l}$ and $\rho_{h}$ are the densities of the heavy and light
fluids, $\mu_{h}$ is the dynamic viscosity of the heavy fluid, $L_{x}$ is the
size of the domain in the horizontal direction and $U^{*}$ is the
characteristic velocity defined as:

$U^{*}=\sqrt{gL_{x}},$ (39)

where $g$ is gravity-driven acceleration. The characteristic time for this
case is defined as:

$T=\frac{L_{x}}{U^{*}}.$ (40)

Following the set-up studied in He et al. (1999), we consider a domain of size
$L_{x}\times 4L_{x}$ with $L_{x}=600$. Initially the top half of the domain is
filled with the heavy liquid and the bottom half with the lighter one. The
interface is perturbed via the following profile:

$h_{i}(x)=\frac{L}{10}\cos\left(\frac{2\pi x}{L_{x}}\right)+2L_{x}.$ (41)

While periodic boundaries were applied in the horizontal direction, at the top
and bottom boundaries no-slip boundary conditions were applied using the half-
way bounce-back scheme Krüger et al. (2017). The At number is set to 0.5 while
two different Re numbers are considered, i.e. Re=256 and 2048. In both cases
$g=6\times 10^{-6}$ while the non-dimensional viscosities were respectively
0.1406 and 0.0176. To validate the simulations, the position of the downward-
plunging heavy liquid spike is measured over time and compared to reference
data from He et al. (1999). The results are illustrated in Figure 2.

Figure 2: (Right) Evolution of interface for the Rayleigh-Taylor instability
for (top row) Re=256 and (bottom row) Re=2048 at different times: (from left
to right) $t/T=$1,2,3,4 and 5. (Left) Position of the penetrating spike over
time: (black) Re=256 and (red) Re=2048. (Plain lines) present results and
(symbols) data from He et al. (1999).

It is observed that both simulations agree very well with the reference
solution of He et al. (1999). To showcase the ability of the solver to handles
under-resolved simulations and illustrate the convergence of the obtained
solutions, the simulations were repeated at two additional lower resolutions
with $L_{x}=$300 and 150 with an acoustic scaling of the time-step size. The
results obtained with those lower resolutions are shown in Figures 3 and 4.

Figure 3: (Right) Interface for the Rayleigh-Taylor instability at $t/T=$5 and
Re=256 for three different resolutions (left to right) $L_{x}$=150, 300 and
600. (Left) Position of the penetrating spike over time: (black) $L_{x}$=600,
(red) $L_{x}$=300 and (blue) $L_{x}$=150.

By looking at the position of the plunging spike it can be clearly seen that
while minor differences exist, even the lowest resolution captures the correct
position. Smaller feature however, especially at Re=2048, need higher
resolutions to be correctly captured. At Re=256 for instance, even the
secondary instability is converged as at $L_{x}$=300 no segmentation is
observed. For Re=2056 on the other hand, while larger structure start to
converge, thinner features clearly need more resolutions.

Figure 4: (Right) Interface for the Rayleigh-Taylor instability at $t/T=$5 and
Re=2048 for three different resolutions (left to right) $L_{x}$=150, 300 and
600. (Left) Position of the penetrating spike over time: (black) $L_{x}$=600,
(red) $L_{x}$=300 and (blue) $L_{x}$=150.

###  3.3 Turbulent 3-D Rayleigh-Taylor instability

To further showcase the ability of the solver to deal with complex flows, we
also consider the Rayleigh-Taylor instability in 3-D. The studied
configuration follows those studied in Liang et al. (2016). The definitions of
non-dimensional parameters are similar to those used in the previous section.
The domain is discretized using $100\times 100\times 1200$ grid-points, with
$L=100$. The interface is placed at the center of the domain along the
$z$-axis, and perturbed using:

$h_{i}(x,y)=\frac{L}{10}\left[\cos\left(\frac{2\pi
x}{L}\right)+\cos\left(\frac{2\pi y}{L}\right)\right]+6L,$ (42)

and the Reynolds and Atwood numbers are set to respectively 1000 and 0.15. As
for previous configurations, periodic boundaries were applied in the
horizontal direction and no-slip boundaries at the top and bottom. The body
force was set to $g=3.6\times 10^{-5}$ and viscosity to 0.006. The position of
the downward-plunging spike was measured over time and compared to reference
data from Liang et al. (2016). After the penetration of two liquids into each
other, the Kelvin-Helmholtz instability causes the plunging spike to roll up
and take a mushroom-like shape. As the mushroom-shaped spike further
progresses into the lighter fluid, the cap disintegrates into four finger-like
structures. It is interesting to note that, as will be shown later, these
fingers are reminiscent of the instabilities leading to splashing in the
impact of a droplet on liquid surfaces.

Figure 5: (Right) Evolution of interface for the 3-D Rayleigh-Taylor
instability for Re=1000 at different times: (from left to right) $t/T=$1.9,
3.9, 5.8, 7.8 and 9.7. (Left) Position of the penetrating spike over time:
(Plain lines) present results and (symbols) data from Liang et al. (2016).

Overall, as shown in Figure 5, the results obtained from the present
simulation are in good agreement with reference data.

###  3.4 Droplet splashing on thin liquid film

As the final case, we consider the impact of a droplet on a thin liquid layer.
This configuration is interesting as it involves complex dynamics such as
splashing and is of interest in many areas of science and engineering
Hagemeier et al. (2011, 2012). Immediately after impact the liquid surface is
perturbed. In many instances, at the contact point (line) a thin liquid jet
forms, and then continues to grow and propagate as a _corolla_. As the crown-
like structures propagates radially, a rim starts to form. At high enough
Weber numbers the structure breaks into small droplets via the
Rayleigh–Plateau instability Josserand and Zaleski (2003). A detailed study of
the initial stages of the spreading process have shown that the spreading
radius scales with time regardless of the Weber and Reynolds numbers Josserand
and Zaleski (2003). While widely studied in the literature using different
numerical formulations Hu et al. (2019); Liang et al. (2018); Fakhari et al.
(2017); Sitompul and Aoki (2019), simulations have usually been limited to
lower density and viscosity ratios and/or Weber and Reynolds numbers Hu et al.
(2019); Liang et al. (2018); Fakhari et al. (2017); Qin et al. (2018). As such
we first focus on a 2-D configuration considering three sets of We and Re
numbers, namely: Re=200 and We=220, Re=1000 and We=220 and Re=1000 and
We=2200. In all simulations the density and viscosity ratios are set to
$\rho_{h}/\rho_{l}=1000$ and $\nu_{l}/\nu_{h}=15$ emulating a water/air
system. The geometrical configuration is illustrated in Figure 6.

Figure 6: Geometrical configuration of the droplet impact on liquid sheet case
in 2-D.

The top and bottom boundary conditions are set to walls modeled with the half-
way bounce-back formulation while symmetrical boundaries are applied to the
left and right. The droplet diameter is resolved with 100 grid-points. The
initial velocity in the droplet is set to $U_{0}=$0.05 and $\nu_{L}$ is
determined via the Reynolds number:

$\hbox{Re}=\frac{\rho_{h}U_{0}D}{\mu_{h}}.$ (43)

Furthermore, the We number is defined as:

${\hbox{We}=\frac{\rho_{l}D{U_{0}}^{2}}{\sigma}.}$ (44)

The evolution of the liquid surface as obtained from the simulations is shown
in Figure 7. Following Josserand and Zaleski (2003), breakup of the rims and
splashing occurs for larger impact parameters defined as:

$K=\hbox{We}^{1/2}\hbox{Re}^{1/4}.$ (45)

Accordingly, the impact parameters for the studied 2-D cases are: K=55.7, 83.4
and 263.8. Looking at the evolution of the systems in Figure 7 it can be
clearly observed that in agreement with observations in Josserand and Zaleski
(2003), larger values of the impact parameter lead to droplet detachment from
the rim and splashing.

Figure 7: Impact of circular droplet on liquid sheet at different We and Re
numbers with $\rho_{h}/\rho_{l}=1000$ and $\nu_{l}/\nu_{h}=15$. (black) Re=200
and We=220, (red) Re=1000 and We=220, and (blue) Re=1000 and We=2200.

Furthermore, the evolution of the spreading radii $r_{K}$ over time for
different cases are shown in Figure 8. As shown there the radii scale with
time at the initial stages of the impact, in agreement with results reported
in Josserand and Zaleski (2003).

Figure 8: Evolution of the spreading radius $r_{K}$ as a function of time for
the droplet impact on liquid film case. Circular symbols designate 2-D
simulations: (black) Re=200 and We=220, (red) Re=1000 and We=220 and (blue)
Re=1000 and We=2200. Rectangular symbols belong to the 3-D simulation with
Re=1000 and We=8000. The dashed line is $\frac{r_{K}}{D}=1.1\sqrt{t/T}$.

As a final test-case, to showcase the robustness of the proposed algorithm, a
3-D configuration with Re=1000 and We=8000 was also ran. The evolution of the
liquid surface over time is shown in Figure 9.

Figure 9: Impact of spherical droplet on thin liquid sheet at We=8000 and
Re=1000 at different times with $\rho_{h}/\rho_{l}=1000$ and
$\nu_{l}/\nu_{h}=15$.

After initial impact a thin liquid jet is formed at the contact line between
the droplet and the sheet. Then, the crown evolves and spreads. At later
stages the finger-like structures start to form at the tip of the crown. These
liquid fingers then get detached from the crown and liquid splashing is
observed. The sequence of events is in excellent agreement with those
presented in Josserand and Zaleski (2003). Furthermore, the spreading radius,
as plotted in Figure 8 is in agreement with theoretical predictions.

##  4 Conclusions

A -based solver relying on the conservative equation and a modified
hydrodynamic pressure/velocity-based distribution and collision operator in
cumulants space was presented in the this study with the aim to model multi-
phase flows in the larger Weber/Reynolds regimes. While the stability at high
Weber numbers –i.e. low surface tensions– is achieved through the decoupled
nature of the conservative formulation, the added stability in terms of
kinematic viscosity –i.e. larger Reynolds numbers– is brought about by the
collision operator and the modified pressure-based formulation for the flow.
Compared to other models available in the literature based on the formulation,
the use of cumulants allows for stability at considerably higher Reynolds
numbers –i.e. lower values of the relaxation factor. For instance,
configurations such as the 3-D droplet splashing were not stable with the
formulation for the same choice of non-dimensional parameters, i.e. resolution
and relaxation factor. The algorithm was shown to capture the dynamics of the
flow and be stable in the targeted regimes. The application of the proposed
algorithm to more complex configurations such as liquid jets is currently
being studied and will be reported in future publications.

conceptualization, S.A.H. and H.S.; methodology, S.A.H.; software, S.A.H.;
validation, S.A.H. and H.S.; formal analysis, S.A.H.; investigation, S.A.H.;
data curation, S.A.H.; writing–original draft preparation, S.A.H.;
writing–review and editing, S.A.H., H.S. and D.T.; visualization, S.A.H.;
supervision, D.T.

S.A.H. and H.S. would like to acknowledge the financial support of the
Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) in TRR 287
(Project-ID 422037413).

The authors declare no conflict of interest..

References

yes

## References

  * Krüger et al. (2017) Krüger, T.; Kusumaatmaja, H.; Kuzmin, A.; Shardt, O.; Silva, G.; Viggen, E.M.  The Lattice Boltzmann Method: Principles and Practice; Graduate Texts in Physics, Springer International Publishing: Cham, 2017.  doi:black10.1007/978-3-319-44649-3.
  * Guo and Shu (2013) Guo, Z.; Shu, C.  Lattice Boltzmann Method and Its Applications in Engineering; Vol. 3, Advances in Computational Fluid Dynamics, WORLD SCIENTIFIC, 2013.  doi:black10.1142/8806.
  * Succi (2002) Succi, S.  The Lattice Boltzmann Equation for Fluid Dynamics and Beyond; 2002.
  * Chorin (1997) Chorin, A.J.  A Numerical Method for Solving Incompressible Viscous Flow Problems.  Journal of Computational Physics 1997, 135, 118–125.  doi:black10.1006/jcph.1997.5716.
  * Shan and Chen (1993) Shan, X.; Chen, H.  Lattice Boltzmann model for simulating flows with multiple phases and components.  Phys. Rev. E 1993, 47, 1815–1819.  doi:black10.1103/PhysRevE.47.1815.
  * Shan and Chen (1994) Shan, X.; Chen, H.  Simulation of nonideal gases and liquid-gas phase transitions by the lattice Boltzmann equation.  Phys. Rev. E 1994, 49, 2941–2948.  doi:black10.1103/PhysRevE.49.2941.
  * Swift et al. (1996) Swift, M.R.; Orlandini, E.; Osborn, W.R.; Yeomans, J.M.  Lattice Boltzmann simulations of liquid-gas and binary fluid systems.  Phys. Rev. E 1996, 54, 5041–5052.  doi:black10.1103/PhysRevE.54.5041.
  * Swift et al. (1995) Swift, M.R.; Osborn, W.R.; Yeomans, J.M.  Lattice Boltzmann Simulation of Nonideal Fluids.  Phys. Rev. Lett. 1995, 75, 830–833.  doi:black10.1103/PhysRevLett.75.830.
  * Wagner and Li (2006) Wagner, A.; Li, Q.  Investigation of Galilean invariance of multi-phase lattice Boltzmann methods.  Physica A: Statistical Mechanics and its Applications 2006, 362, 105–110.  doi:black10.1016/j.physa.2005.09.030.
  * Kupershtokh et al. (2009) Kupershtokh, A.; Medvedev, D.; Karpov, D.  On equations of state in a lattice Boltzmann method.  Computers & Mathematics with Applications 2009, 58, 965–974.  doi:black10.1016/j.camwa.2009.02.024.
  * Yuan and Schaefer (2006) Yuan, P.; Schaefer, L.  Equations of state in a lattice Boltzmann model.  Physics of Fluids 2006, 18, 042101.  doi:black10.1063/1.2187070.
  * Sbragaglia et al. (2007) Sbragaglia, M.; Benzi, R.; Biferale, L.; Succi, S.; Sugiyama, K.; Toschi, F.  Generalized lattice Boltzmann method with multirange pseudopotential.  Phys. Rev. E 2007, 75, 026702.  doi:black10.1103/PhysRevE.75.026702.
  * Li and Luo (2013) Li, Q.; Luo, K.H.  Achieving tunable surface tension in the pseudopotential lattice Boltzmann modeling of multiphase flows.  Phys. Rev. E 2013, 88, 053307.  doi:black10.1103/PhysRevE.88.053307.
  * Fakhari and Rahimian (2010) Fakhari, A.; Rahimian, M.H.  Phase-field modeling by the method of lattice Boltzmann equations.  Physical Review E 2010, 81, 036707.
  * Safari et al. (2013) Safari, H.; Rahimian, M.H.; Krafczyk, M.  Extended lattice Boltzmann method for numerical simulation of thermal phase change in two-phase fluid flow.  Physical Review E 2013, 88, 013304.
  * Safari et al. (2014) Safari, H.; Rahimian, M.H.; Krafczyk, M.  Consistent simulation of droplet evaporation based on the phase-field multiphase lattice Boltzmann method.  Physical Review E 2014, 90, 033305.
  * Yazdi et al. (2018) Yazdi, H.; Rahimiani, M.H.; Safari, H.  Numerical simulation of pressure-driven phase-change in two-phase fluid flows using the Lattice Boltzmann Method.  Computers & Fluids 2018, 172, 8–18.
  * Wang et al. (2019) Wang, H.; Yuan, X.; Liang, H.; Chai, Z.; Shi, B.  A brief review of the phase-field-based lattice Boltzmann method for multiphase flows.  Capillarity 2019, 2, 33–52.
  * He et al. (1999) He, X.; Chen, S.; Zhang, R.  A Lattice Boltzmann Scheme for Incompressible Multiphase Flow and Its Application in Simulation of Rayleigh–Taylor Instability.  Journal of Computational Physics 1999, 152, 642–663.  doi:black10.1006/jcph.1999.6257.
  * Inamuro et al. (2004) Inamuro, T.; Ogata, T.; Tajima, S.; Konishi, N.  A lattice Boltzmann method for incompressible two-phase flows with large density differences.  Journal of Computational Physics 2004, 198, 628–644.  doi:black10.1016/j.jcp.2004.01.019.
  * Amirshaghaghi et al. (2016) Amirshaghaghi, H.; Rahimian, M.; Safari, H.  Application of a two phase lattice Boltzmann model in simulation of free surface jet impingement heat transfer.  International Communications in Heat and Mass Transfer 2016, 75, 282–294.
  * Amirshaghaghi et al. (2018) Amirshaghaghi, H.; Rahimian, M.H.; Safari, H.; Krafczyk, M.  Large Eddy Simulation of liquid sheet breakup using a two-phase lattice Boltzmann method.  Computers & Fluids 2018, 160, 93–107.
  * Hosseini (2020) Hosseini, S.A.  Development of a lattice Boltzmann-based numerical method for the simulation of reacting flows.  PhD thesis, Otto-von-Guericke Universität/Universite Paris-Saclay, 2020\.
  * Sun and Beckermann (2007) Sun, Y.; Beckermann, C.  Sharp interface tracking using the phase-field equation.  Journal of Computational Physics 2007, 220, 626–653.
  * Chiu and Lin (2011) Chiu, P.H.; Lin, Y.T.  A conservative phase field method for solving incompressible two-phase flows.  Journal of Computational Physics 2011, 230, 185–204.
  * Fakhari et al. (2017) Fakhari, A.; Bolster, D.; Luo, L.S.  A weighted multiple-relaxation-time lattice Boltzmann method for multiphase flows and its application to partial coalescence cascades.  Journal of Computational Physics 2017, 341, 22–43.  doi:black10.1016/j.jcp.2017.03.062.
  * Wang et al. (2016) Wang, H.; Chai, Z.; Shi, B.; Liang, H.  Comparative study of the lattice Boltzmann models for Allen-Cahn and Cahn-Hilliard equations.  Physical Review E 2016, 94, 033304.
  * Chopard et al. (2009) Chopard, B.; Falcone, J.L.; Latt, J.  The lattice Boltzmann advection-diffusion model revisited.  The European Physical Journal Special Topics 2009, 171, 245–249.
  * Hosseini et al. (2019) Hosseini, S.A.; Darabiha, N.; Thévenin, D.  Lattice Boltzmann advection-diffusion model for conjugate heat transfer in heterogeneous media.  International Journal of Heat and Mass Transfer 2019, 132, 906–919.
  * Zu et al. (2020) Zu, Y.; Li, A.; Wei, H.  Phase-field lattice Boltzmann model for interface tracking of a binary fluid system based on the Allen-Cahn equation.  Physical Review E 2020, 102, 053307.
  * Lee and Lin (2003) Lee, T.; Lin, C.L.  Pressure evolution lattice Boltzmann equation method for two-phase flow with phase change.  Physical Review E 2003, 67, 056703.
  * Hosseini et al. (2019) Hosseini, S.A.; Safari, H.; Darabiha, N.; Thévenin, D.; Krafczyk, M.  Hybrid Lattice Boltzmann-finite difference model for low Mach number combustion simulation.  Combustion and Flame 2019, 209, 394–404.
  * Hosseini et al. (2020) Hosseini, S.A.; Abdelsamie, A.; Darabiha, N.; Thévenin, D.  Low-Mach hybrid lattice Boltzmann-finite difference solver for combustion in complex flows.  Physics of Fluids 2020, 32, 077105.
  * Geier et al. (2020) Geier, M.; Lenz, S.; Schönherr, M.; Krafczyk, M.  Under-resolved and large eddy simulations of a decaying Taylor–Green vortex with the cumulant lattice Boltzmann method.  Theor. Comput. Fluid Dyn. 2020.  doi:black10.1007/s00162-020-00555-7.
  * Geier et al. (2015) Geier, M.; Schönherr, M.; Pasquali, A.; Krafczyk, M.  The cumulant lattice Boltzmann equation in three dimensions: Theory and validation.  Computers & Mathematics with Applications 2015, 70, 507–547.  doi:black10.1016/j.camwa.2015.05.001.
  * Qin et al. (2018) Qin, F.; Mazloomi Moqaddam, A.; Kang, Q.; Derome, D.; Carmeliet, J.  Entropic multiple-relaxation-time multirange pseudopotential lattice Boltzmann model for two-phase flow.  Physics of Fluids 2018, 30, 032104.  doi:black10.1063/1.5016965.
  * Mazloomi M et al. (2015) Mazloomi M, A.; Chikatamarla, S.; Karlin, I.  Entropic Lattice Boltzmann Method for Multiphase Flows.  Phys. Rev. Lett. 2015, 114, 174502.  doi:black10.1103/PhysRevLett.114.174502.
  * Yang et al. (2018a) Yang, X.; He, H.; Xu, J.; Wei, Y.; Zhang, H.  Entropy generation rates in two-dimensional Rayleigh–Taylor turbulence mixing.  Entropy 2018, 20, 738.
  * Yang et al. (2018b) Yang, H.; Wei, Y.; Zhu, Z.; Dou, H.; Qian, Y.  Statistics of heat transfer in two-dimensional turbulent Rayleigh-Bénard convection at various Prandtl Number.  Entropy 2018, 20, 582.
  * Rahmat et al. (2014) Rahmat, A.; Tofighi, N.; Shadloo, M.; Yildiz, M.  Numerical simulation of wall bounded and electrically excited Rayleigh–Taylor instability using incompressible smoothed particle hydrodynamics.  Colloids and Surfaces A: Physicochemical and Engineering Aspects 2014, 460, 60–70.
  * Liang et al. (2016) Liang, H.; Li, Q.X.; Shi, B.C.; Chai, Z.H.  Lattice Boltzmann simulation of three-dimensional Rayleigh-Taylor instability.  Phys. Rev. E 2016, 93, 033113.  doi:black10.1103/PhysRevE.93.033113.
  * Hagemeier et al. (2011) Hagemeier, T.; Hartmann, M.; Thévenin, D.  Practice of vehicle soiling investigations: A review.  International Journal of Multiphase Flow 2011, 37, 860–875.
  * Hagemeier et al. (2012) Hagemeier, T.; Hartmann, M.; Kühle, M.; Thévenin, D.; Zähringer, K.  Experimental characterization of thin films, droplets and rivulets using LED fluorescence.  Experiments in fluids 2012, 52, 361–374.
  * Josserand and Zaleski (2003) Josserand, C.; Zaleski, S.  Droplet splashing on a thin liquid film.  Phys. Fluids 2003, 15, 1650.  doi:black10.1063/1.1572815.
  * Hu et al. (2019) Hu, Y.; Li, D.; Jin, L.; Niu, X.; Shu, S.  Hybrid Allen-Cahn-based lattice Boltzmann model for incompressible two-phase flows: The reduction of numerical dispersion.  Phys. Rev. E 2019, 99, 023302.  doi:black10.1103/PhysRevE.99.023302.
  * Liang et al. (2018) Liang, H.; Xu, J.; Chen, J.; Wang, H.; Chai, Z.; Shi, B.  Phase-field-based lattice Boltzmann modeling of large-density-ratio two-phase flows.  Phys. Rev. E 2018, 97, 033309.  doi:black10.1103/PhysRevE.97.033309.
  * Sitompul and Aoki (2019) Sitompul, Y.P.; Aoki, T.  A filtered cumulant lattice Boltzmann method for violent two-phase flows.  Journal of Computational Physics 2019, 390, 93–120.  doi:black10.1016/j.jcp.2019.04.019.

##  Appendix A Hermite polynomials and coefficients

The Hermite polynomials used in the s of different solvers are defined as:

$\displaystyle\mathcal{H}_{0}$ $\displaystyle=1,$ (46a)
$\displaystyle\mathcal{H}_{i}$ $\displaystyle=c_{\alpha,i},$ (46b)
$\displaystyle\mathcal{H}_{ij}$
$\displaystyle=c_{\alpha,i}c_{\alpha,j}-c_{s}^{2}\delta_{ij},$ (46c)

where $\delta_{ij}$ denotes the Kronecker delta function, while corresponding
equilibrium coefficients are:

$\displaystyle a^{(eq)}_{0}$ $\displaystyle=\rho,$ (47a) $\displaystyle
a^{(eq)}_{i}$ $\displaystyle=\rho u_{i},$ (47b) $\displaystyle a^{(eq)}_{ij}$
$\displaystyle=\rho u_{i}u_{j},$ (47c)