Open AccessArticle

Lattice Boltzmann Solver for Multiphase Flows: Application to High Weber and Reynolds Numbers

Seyed Ali Hosseini

^1,2,*,

Hesameddin Safari

¹ and

Dominique Thevenin

Laboratory of Fluid Dynamics and Technical Flows, University of Magdeburg “Otto von Guericke”, D-39106 Magdeburg, Germany

Department of Mechanical and Process Engineering, ETH Zürich, 8092 Zürich, Switzerland

Author to whom correspondence should be addressed.

Entropy 2021, 23(2), 166; https://doi.org/10.3390/e23020166

Submission received: 26 December 2020 / Revised: 14 January 2021 / Accepted: 21 January 2021 / Published: 29 January 2021

(This article belongs to the Special Issue Computational Fluid Dynamics and Conjugate Heat Transfer)

Download

Browse Figures

Versions Notes

Abstract

The lattice Boltzmann method, now widely used for a variety of applications, has also been extended to model multiphase flows through different formulations. While already applied to many different configurations in 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 a decoupled phase-field formulation—the conservative Allen–Cahn equation—and a cumulant-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 was validated through a variety of test cases, starting with the Rayleigh–Taylor instability in both 2D and 3D, followed by the impact of a droplet on a liquid sheet. In all simulations, the solver correctly captured the flow dynamics andmatched reference results very well. As the final test case, the solver was used to model droplet splashing on a thin liquid sheet in 3D with a density ratio of 1000 and kinematic viscosity ratio of 15, matching the water/air system at We = 8000 and Re = 1000. Results showed that the solver correctly captured 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; multiphase flows; conservative Allen–Cahn; phase field

1. Introduction

The lattice Boltzmann method (LBM) is a discrete solver for the so-called discrete velocity Boltzmann equation (DVBE), initially developed as an alternative to classical solvers for the incompressible hydrodynamic regime [1,2]. Due to the simplicity of the algorithm, low computational cost of discrete time-evolution equations, and locality of nonlinear terms and boundary conditions, it has rapidly grown over the past few decades [3]. While intended for the incompressible regime, the LBM formally solves the compressible isothermal Navier-Stokes (NS) equations at a reference temperature. While originally tied to the considered flow’s temperature, in the context of the lattice Boltzmann (LB) solver, the reference temperature is a numerical parameter allowing for control over convergence and consistency [1]. Weak compressibility in the formulation along with the parabolic nature of the partial differential equation (PDE) governing the evolution of pressure, as opposed to Chorin’s original artificial compressibility method (ACM), made the scheme efficient and applicable to unsteady flows [4]. Although originally used for single-phase flows, it has since been extended to multiphase, multispecies, and compressible flows.

While generally based on diffuse-interface formulations, LB solvers for multiphase flows can be categorized as pertaining to one of three major categories: (a) pseudopotential [5,6], (b) free energy [7,8], and (c) phase field. Other types of formulations can also be found in the literature, but they are not as widely spread and/or developed as these three.

In the context of the free-energy formulation, the expression for the nonlocal nonideal pressure tensor is found through the free-energy functional. The appropriate pressure tensor is then introduced into the LB solver via a moment-matching approach assigning coefficients to different terms in the equilibrium distribution function (EDF) [8]. This formulation is consistent and differentiated from the generic double-well potential-based Cahn–Hilliard formulation because, in the minimization process of free energy, the equation of state (EoS) is explicitly considered. As is the case for the pseudopotential formulation, the explicit intervention of the EoS within the free functional ties the thickness of the interface to physical parameters, e.g., surface tension, density ratio, and EoS. As a consequence, the choice of the EoS and/or tuning of the coefficients in the EoS is a method of choice to widen the area of accessible density ratios. This approach was later extended by introducing nonideal components of the pressure tensor via external body forces. 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 EDF [9].

The pseudopotential formulation follows more of a bottom–up approach in introducing nonideal dynamics into the solver. It follows the general philosophy of the Boltzmann–Vlasov equation, introducing a nonlocal potential to account for nonideal effects. While the original formulation relied on what was termed effective density, actual EoS were introduced into the pseudopotential in [10,11]. Apart from thermodynamic consistency, the possibility of using different EoS allowed for higher density ratios to be modelled. 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 the independent tuning of surface tension [12], the spectrum of values covered by the model is rather limited and barely allows for variations of one order of magnitude [13].

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 [14], allowing for decoupling, among other parameters, the interface thickness from the fluid physical properties. Another consequence of this choice of functional is the partial loss of thermodynamic consistency, making the extension of the formulation to more complex physics such as thermal flows, compressible flows and acoustics less straightforward, although a number of attempts were documented in the literature [15,16,17]. Nevertheless, it was observed to be very effective and robust for multiphase 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 [18]. Approaches relying on the explicit tracking of the interface with a consistent energy functional making use of the nonideal EoS were also proposed as ways to improve the stability of the original free-energy formulation [19,20].

Over the past few decades, much effort has been put into developing phase-field-based LB solvers for various applications [16,21,22]. Given that in such formulations 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 LB solver for the flow usually referred to as incompressible formulation. The so-called low-Mach formulation is mostly based on the modified distribution function introduced in [19], where pressure is the zeroth-order moment of the distribution function. This flow solver was combined with different forms of interface-tracking formulations, e.g., Allen-Cahn (AC), conservative AC, or Cahn-Hilliard (CH) to model multiphase flows. The aim of the present study is to introduce a multiphase solver relying on the pressure-based formulation of [19] and a multiple relaxation time (MRT) realization for the flow solver coupled with a LB solver for the conservative AC. The use of the MRT collision operator in cumulant space with the decoupled interface tracking allows for simulations in high Reynolds and Weber regimes. After a brief introduction of the model, it is used to simulate a variety of test cases, proving its ability to reproduce correct physics and its robustness. All models were implemented in our in-house multiphysics solver, ALBORZ [23].

2. Theoretical Background

2.1. Target Macrosopic System

As briefly stated in the introduction, the aim of the present work is to solve multiphase flow equations within the context of the diffuse interface formulation in the limit of an incompressible regime, where interface dynamics are followed and accounted for via an additional indicator field,

ϕ

. As such, at the macroscopic level, low Mach NS equations are targeted:

\partial_{t} ρ u_{i} + \partial_{j} ρ u_{i} u_{j} + \partial_{j} σ_{i j} + μ_{ϕ} \partial_{i} ϕ + F_{b, i} = 0,

(1)

where

u_{i}

is fluid velocity,

ρ

the fluid density, and

F_{b, i}

designates external body forces. The stress tensor

σ_{i j}

is defined as:

σ_{i j} = p_{h} δ_{i j} - η (\partial_{i} u_{j} + \partial_{j} u_{i}) + (\frac{2}{3} η - ξ) \partial_{k} u_{k} δ_{i j},

(2)

where

η

is the fluid dynamic viscosity tied to kinematic viscosity

ν

η = ρ ν

ξ

the bulk viscosity and

p_{h}

the hydrodynamic pressure. The chemical potential

μ_{ϕ}

is defined as

μ_{ϕ} = 2 β ϕ (ϕ - 1) (2 ϕ - 1) κ Δ ϕ,

(3)

where

Δ = \nabla^{2}

is the Laplacian operator, and

β

and

κ

are parameters specific to the AC formulation. The second term on the right hand side (RHS) of Equation (1) accounts for surface-tension effects. For the sake of clarity, free parameters are detailed in the next paragraph.

The interface was tracked using the conservative AC equation, where order parameter

ϕ

evolved as [24,25]:

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

(4)

where the order parameter

ϕ

takes on values between 0 and 1, M is mobility, W is interface thickness, and

n_{i}

is the unit normal to the interface, obtained as

n_{i} = \frac{\partial_{i} ϕ}{| | \nabla ϕ | |} .

(5)

Interfaces can be found through isosurfaces of the order parameter, i.e.,

ϕ = 1 / 2

. To recover the correct surface tension, free parameters appearing in the chemical potential, i.e.,

κ

and

β

, are tied to surface tension

σ

and interface thickness W in the AC equation via

β = 12 σ / W

and

κ = 3 σ W / 2

2.2. LB Formulation for Conservative Phase-Field Equation

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

\partial_{t} g_{α} + c_{α, i} \partial_{i} g_{α} + S_{α} = Ω_{α}^{ϕ},

(6)

where

g_{α}

and

c_{α}

are populations and velocities in the discrete velocity kinetic model, and the collision operator is defined as

Ω_{α}^{ϕ} = \frac{1}{τ_{ϕ}} (g_{α}^{(e q)} - g_{α}) .

(7)

The EDF is defined as

g_{α}^{(e q)} = w_{α} ϕ \sum_{n = 0}^{2} \frac{1}{n! c_{s}^{2 n}} H_{n} : a_{n}^{(e q)},

(8)

where

H_{n}

and

a_{n}^{(e q)}

are the Hermite polynomial and coefficient of order n,

c_{s}

is lattice sound speed, and

w_{α}

are 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 [26]

S_{α} = w_{α} H_{i} n_{i} \frac{4 ϕ (1 - ϕ)}{W} .

(9)

Given that the source term affects the first-order moment, a nonconserved moment of the distribution function, the distribution function is tied to the phase parameter as

ϕ = \sum_{α} g_{α} .

(10)

The relaxation coefficient is fixed as

τ_{ϕ} = \frac{M}{c_{s}^{2}} .

(11)

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

{\bar{g}}_{α} (x + c_{α} δ_{t}, t + δ_{t}) = (1 - \frac{δ_{t}}{{\bar{τ}}_{ϕ}}) {\bar{g}}_{α} (x, t) + \frac{δ_{t}}{{\bar{τ}}_{ϕ}} g_{α}^{(e q)} (x, t) + δ_{t} {\bar{S}}_{α} (x, t),

(12)

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

{\bar{S}}_{α} = (1 - \frac{1}{2 τ_{ϕ}}) w_{α} H_{i} n_{i} \frac{4 ϕ (1 - ϕ)}{W},

(13)

and:

{\bar{τ}}_{ϕ} = τ_{ϕ} + \frac{δ_{t}}{2} .

(14)

The derivatives of the order parameter appearing in the various discrete time-evolution equations are computed using isotropic finite differences, i.e.,

\partial_{i} ϕ = \frac{1}{c_{s}^{2}} \sum_{α} w_{α} c_{α, i} ϕ (x + c_{α}),

(15)

and

\partial_{i}^{2} ϕ = \frac{2}{c_{s}^{2}} \sum_{α} w_{α} [ϕ (x + c_{α}) - ϕ (x)] .

(16)

While the present work makes use of a second-order EDF, the same macroscopic PDE, i.e., Equation (4), can also be recovered by using a first-order EDF and an additional correction term of the following form [27]:

C_{α} = \frac{w_{α}}{c_{s}^{2}} H_{i} \partial_{t} ϕ u_{i},

(17)

which as for Equation (13), postdiscretization changes into

{\bar{C}}_{α} = (1 - \frac{1}{2 τ_{ϕ}}) \frac{w_{α}}{c_{s}^{2}} H_{i} \partial_{t} ϕ u_{i} .

(18)

Such correction terms were first introduced in the context of advection–diffusion LB solvers [28], and further extended to nonlinear equations in the same context [29]. Detailed derivation and multiscale analyses are readily available in the literature, e.g., [30].

2.3. LB Model for Flow Field

The flow solver kinetic model follows the low-Mach formulation used, among other sources, in [31,32,33], and is based on the original model introduced in [19]

\partial_{t} f_{α}^{^{'}} + c_{α, i} \partial_{i} f_{α}^{^{'}} = Ω_{α} + Ξ_{α},

(19)

where the collision operator is

Ω_{α} = \frac{1}{τ} ({f_{α}^{(e q)}}^{^{'}} - f_{α}^{^{'}}),

(20)

Ξ_{α}

is defined as

Ξ_{α} = c_{s}^{2} (\frac{f_{α}^{(e q)}}{ρ} - w_{α}) (c_{α, i} - u_{i}) \partial_{i} ρ + w_{α} c_{s}^{2} ρ \partial_{i} u_{i} + (F_{b, i} + F_{s, i}) (c_{α, i} - u_{i}) \frac{f_{α}^{(e q)}}{ρ},

(21)

and the relaxation coefficient

τ

is tied to fluid kinematic viscosity

ν

τ = \frac{ν}{c_{s}^{2}} .

(22)

Forces

F_{b, i}

and

F_{s, i}

represent external body forces and surface tension, respectively, i.e.,

F_{s, i} = μ_{ϕ} \partial_{i} ϕ .

(23)

The modified distribution function

f_{α}^{^{'}}

is defined as

f_{α}^{^{'}} = w_{α} p_{h} + c_{s}^{2} (f_{α} - w_{α} ρ),

(24)

where

f_{α}

is the classical isothermal distribution function. The modified equilibrium follows the same logic and is defined as

{f_{α}^{(e q)}}^{^{'}} = w_{α} p_{h} + w_{α} ρ c_{s}^{2} \sum_{n = 1}^{2} \frac{1}{n! c_{s}^{2 n}} H_{n} : a_{n}^{(e q)} .

(25)

Density is tied to the order parameter as

ρ = ρ_{l} + (ρ_{h} - ρ_{l}) ϕ,

(26)

where

ρ_{h}

and

ρ_{l}

are the densities of the heavy and light fluid, respectively. For detailed analysis of the macroscopic equations recovered by this model and the derivation of the discrete equations, interested readers are referred to [23,32]. In the context of the present study, the low-Mach model was wrapped in a moment-based formulation where postcollision populations

f_{α}^{^{'} *}

to be streamed as

{f^{^{'}}}_{α} (x + c_{α} δ_{t}, t + δ_{t}) = {f^{^{'} *}}_{α} (x, t),

(27)

are computed as

{f^{^{'} *}}_{α} = ρ c_{s}^{2} f_{α}^{p *} + \frac{δ_{t}}{2} Ξ_{α} .

(28)

The postcollision preconditioned population

f_{α}^{p *}

f_{α}^{p *} = C^{- 1} (I - W) K^{p} + C^{- 1} W K^{p},

(29)

where

C

is the moment transform matrix from preconditioned populations to the target momentum space,

I

is the identity matrix, and

W

is the diagonal relaxation frequency matrix. Following [34], prior to transformation to momentum space, populations are preconditioned as

f_{α}^{p} = \frac{1}{ρ c_{s}^{2}} f_{α}^{^{'}} + \frac{δ_{t}}{2 ρ c_{s}^{2}} Ξ_{α} .

(30)

This preconditioning accomplishes two tasks, namely, normalizing the populations with density and thus eliminating the density dependence of the moments, and introducing the first half of the source term. As such, moments

K^{p}

are computed as

K_{β}^{p} = C_{α β} f_{α}^{p} .

(31)

The transformation from distribution function (DF)s to cumulants is carried out using the steps suggested in [35], which allows for a more efficient algorithm. The DFs are first transformed into central moments:

{\tilde{Π}}_{β}^{p} = \sum_{α} {(c_{α, x} - u_{x})}^{n_{x}} {(c_{α, y} - u_{y})}^{n_{y}} {(c_{α, z} - u_{z})}^{n_{z}} f_{α}^{p} .

(32)

Here,

β = x^{n_{x}} y^{n_{y}} z^{n_{z}}

. The central moments are then transformed into the corresponding cumulants using the following relations:

\begin{matrix} K_{x}^{p} = & {\tilde{Π}}_{x}^{p}, \end{matrix}

(33a)

\begin{matrix} K_{x y}^{p} = & {\tilde{Π}}_{x y}^{p}, \end{matrix}

(33b)

\begin{matrix} K_{x^{2}}^{p} = & {\tilde{Π}}_{x^{2}}^{p}, \end{matrix}

(33c)

\begin{matrix} K_{x y^{2}}^{p} = & {\tilde{Π}}_{x y^{2}}^{p}, \end{matrix}

(33d)

\begin{matrix} K_{x y z}^{p} = & {\tilde{Π}}_{x y z}^{p}, \end{matrix}

(33e)

\begin{matrix} K_{x^{2} y z}^{p} = & {\tilde{Π}}_{x^{2} y z}^{p} - [{\tilde{Π}}_{x^{2}}^{p} {\tilde{Π}}_{y z}^{p} + 2 {\tilde{Π}}_{x y}^{p} {\tilde{Π}}_{x z}^{p}], \end{matrix}

(33f)

\begin{matrix} K_{x^{2} y^{2}}^{p} = & {\tilde{Π}}_{x^{2} y^{2}}^{p} - [{\tilde{Π}}_{x^{2}}^{p} {\tilde{Π}}_{y^{2}}^{p} + 2 {({\tilde{Π}}_{x y}^{p})}^{2}], \end{matrix}

(33g)

\begin{matrix} K_{x y^{2} z^{2}}^{p} = & {\tilde{Π}}_{x y^{2} z^{2}}^{p} - [{\tilde{Π}}_{z^{2}}^{p} {\tilde{Π}}_{x y^{2}}^{p} + {\tilde{Π}}_{y^{2}}^{p} {\tilde{Π}}_{x z^{2}}^{p} + 4 {\tilde{Π}}_{y z}^{p} {\tilde{Π}}_{x y z}^{p} + 2 ({\tilde{Π}}_{x z}^{p} {\tilde{Π}}_{y^{2} z}^{p} + {\tilde{Π}}_{x y}^{p} {\tilde{Π}}_{y z^{2}}^{p})], \\ K_{x^{2} y^{2} z^{2}}^{p} = & {\tilde{Π}}_{x^{2} y^{2} z^{2}}^{p} - [4 {({\tilde{Π}}_{x y z}^{p})}^{2} + {\tilde{Π}}_{x^{2}}^{p} {\tilde{Π}}_{y^{2} z^{2}}^{p} + {\tilde{Π}}_{y^{2}}^{p} {\tilde{Π}}_{x^{2} z^{2}}^{p} + {\tilde{Π}}_{z^{2}}^{p} {\tilde{Π}}_{x^{2} y^{2}}^{p} + 4 ({\tilde{Π}}_{x y}^{p} {\tilde{Π}}_{x^{2} y z}^{p} + \\ {\tilde{Π}}_{x z}^{p} {\tilde{Π}}_{x y^{2} z}^{p} + {\tilde{Π}}_{x y}^{p} {\tilde{Π}}_{x y z^{2}}^{p} + 2 ({\tilde{Π}}_{x y^{2}}^{p} {\tilde{Π}}_{x z^{2}}^{p} + {\tilde{Π}}_{x^{2} y}^{p} {\tilde{Π}}_{y z^{2}}^{p} + {\tilde{Π}}_{x^{2} z}^{p} {\tilde{Π}}_{y^{2} z}^{p})) + \end{matrix}

(33h)

\begin{matrix} (16 {\tilde{Π}}_{x y}^{p} {\tilde{Π}}_{x z}^{p} {\tilde{Π}}_{y z}^{p} + 4 ({({\tilde{Π}}_{x z}^{p})}^{2} {\tilde{Π}}_{y^{2}}^{p} + {({\tilde{Π}}_{y z}^{p})}^{2} {\tilde{Π}}_{x^{2}}^{p} + {({\tilde{Π}}_{x y}^{p})}^{2} {\tilde{Π}}_{z^{2}}^{p}) + 2 {\tilde{Π}}_{x^{2}}^{p} {\tilde{Π}}_{y^{2}}^{p} {\tilde{Π}}_{z^{2}}^{p})] . \end{matrix}

(33i)

The remainder of the moments can be easily obtained via permutation of the indices. The collision process was performed in cumulant space according to [35]. The fluid viscosity is controlled via the collision factor related to second-order cumulants (e.g.,

K_{x y}^{p}

K_{x^{2}}^{p} - K_{y^{2}}^{p}

K_{x^{2}}^{p} - K_{z^{2}}^{p}

etc.). The rest of the collision factors were set to unity for simplicity. Once the collision step had been applied, cumulants were transformed back into central moments as

\begin{matrix} {\tilde{Π}}_{x}^{p *} = & K_{x}^{p *}, \end{matrix}

(34a)

\begin{matrix} {\tilde{Π}}_{x y}^{p *} = & K_{x y}^{p *}, \end{matrix}

(34b)

\begin{matrix} {\tilde{Π}}_{x^{2}}^{p *} = & K_{x^{2}}^{p *}, \end{matrix}

(34c)

\begin{matrix} {\tilde{Π}}_{x y^{2}}^{p *} = & K_{x y^{2}}^{p *}, \end{matrix}

(34d)

\begin{matrix} {\tilde{Π}}_{x y z}^{p *} = & K_{x y z}^{p *}, \end{matrix}

(34e)

\begin{matrix} {\tilde{Π}}_{x^{2} y z}^{p *} = & K_{x^{2} y z}^{p *} + [{\tilde{Π}}_{x^{2}}^{p *} {\tilde{Π}}_{y z}^{p *} + 2 {\tilde{Π}}_{x y}^{p *} {\tilde{Π}}_{x z}^{p *}], \end{matrix}

(34f)

\begin{matrix} {\tilde{Π}}_{x^{2} y^{2}}^{p *} = & K_{x^{2} y^{2}}^{p *} + [{\tilde{Π}}_{x^{2}}^{p *} {\tilde{Π}}_{y^{2}}^{p *} + 2 {({\tilde{Π}}_{x y}^{p *})}^{2}], \end{matrix}

(34g)

\begin{matrix} {\tilde{Π}}_{x y^{2} z^{2}}^{p *} = & K_{x y^{2} z^{2}}^{p *} + [{\tilde{Π}}_{z^{2}}^{p *} {\tilde{Π}}_{x y^{2}}^{p *} + {\tilde{Π}}_{y^{2}}^{p *} {\tilde{Π}}_{x z^{2}}^{p *} + 4 {\tilde{Π}}_{y z}^{p *} {\tilde{Π}}_{x y z}^{p *} + 2 ({\tilde{Π}}_{x z}^{p *} {\tilde{Π}}_{y^{2} z}^{p *} + {\tilde{Π}}_{x y}^{p *} {\tilde{Π}}_{y z^{2}}^{p *})], \\ {\tilde{Π}}_{x^{2} y^{2} z^{2}}^{p *} = & K_{x^{2} y^{2} z^{2}}^{p *} + [4 {({\tilde{Π}}_{x y z}^{p *})}^{2} + {\tilde{Π}}_{x^{2}}^{p *} {\tilde{Π}}_{y^{2} z^{2}}^{p *} + {\tilde{Π}}_{y^{2}}^{p *} {\tilde{Π}}_{x^{2} z^{2}}^{p *} + {\tilde{Π}}_{z^{2}}^{p *} {\tilde{Π}}_{x^{2} y^{2}}^{p *} + 4 ({\tilde{Π}}_{x y}^{p *} {\tilde{Π}}_{x^{2} y z}^{p *} + \\ {\tilde{Π}}_{x z}^{p *} {\tilde{Π}}_{x y^{2} z}^{p *} + {\tilde{Π}}_{x y}^{p *} {\tilde{Π}}_{x y z^{2}}^{p *} + 2 ({\tilde{Π}}_{x y^{2}}^{p *} {\tilde{Π}}_{x z^{2}}^{p *} + {\tilde{Π}}_{x^{2} y}^{p *} {\tilde{Π}}_{y z^{2}}^{p *} + {\tilde{Π}}_{x^{2} z}^{p *} {\tilde{Π}}_{y^{2} z}^{p *})) - \end{matrix}

(34h)

\begin{matrix} (16 {\tilde{Π}}_{x y}^{p *} {\tilde{Π}}_{x z}^{p *} {\tilde{Π}}_{y z}^{p *} + 4 ({({\tilde{Π}}_{x z}^{p *})}^{2} {\tilde{Π}}_{y^{2}}^{p *} + {({\tilde{Π}}_{y z}^{p *})}^{2} {\tilde{Π}}_{x^{2}}^{p *} + {({\tilde{Π}}_{x y}^{p *})}^{2} {\tilde{Π}}_{z^{2}}^{p *}) + 2 {\tilde{Π}}_{x^{2}}^{p *} {\tilde{Π}}_{y^{2}}^{p *} {\tilde{Π}}_{z^{2}}^{p *})] . \end{matrix}

(34i)

After this step, postcollision central moments could be readily transformed back into populations. All transforms presented here and upcoming simulations are based on the D3Q27 stencil. The following set of 27 moments were used as the basis for the moments:

\begin{matrix} β \in {0, x, y, z, x y, x z, y z, x^{2} - y^{2}, x^{2} - z^{2}, x^{2} + y^{2} + z^{2}, \\ x y^{2} + x z^{2}, x y z, x y^{2} - x z^{2}, x^{2} + y z^{2}, x^{2} z + y^{2} z, x^{2} y - y z^{2}, x^{2} z - y^{2} z, x^{2} y^{2} - 2 x^{2} z^{2} + y^{2} z^{2}, \\ x^{2} y^{2} + x^{2} z^{2} - 2 y^{2} z^{2}, x^{2} y^{2} + x^{2} z^{2} + y^{2} z^{2}, x^{2} y z, x y^{2} z, x y z^{2}, x^{2} y^{2} z, x^{2} y z^{2}, x y^{2} z^{2}, x^{2} y^{2} z^{2}}, \end{matrix}

(35)

where

β = x^{2} - y^{2}

stands for a central moment of form

{\tilde{Π}}_{x^{2}}^{p} - {\tilde{Π}}_{y^{2}}^{p}

. Previous systematic studies of the flow solver showed second-order convergence under diffusive scaling [32].

3. Numerical Applications

In this section, the proposed numerical method is validated through different test cases. All results and simulation parameters are reported in LB units, i.e., nondimensionalized 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 considered the case of a static droplet in a rectangular domain with periodic boundaries all around. All cases consisted of a domain of

256 \times 256

size filled with a light fluid. A droplet of the heavier fluid was placed at the center of the domain. Simulations were pursued till the system had converged. The pressure difference between the droplet and surrounding lighter fluid was then extracted. Using Laplace’s law, i.e.,

Δ P = \frac{σ}{r},

(36)

where

Δ P

is the pressure difference, and r the droplet radius, one can readily obtain the effective surface tension. Three different surface tensions, i.e.,

σ = 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. Obtained results are shown in Figure 1. Results presented here consider a density ratio of 20 and nondimensional viscosity of 0.1.

The model satisfied Laplace’s law and recovered the correct surface tensions. Furthermore, it could span a wide range of surface tensions, as opposed to other classes of multiphase solvers, such as free energy or pseudopotential formulations [36,37], and maintain relatively low spurious currents. For example, at a density ratio of 1000 and

σ = 10^{- 3}

, spurious currents were found to be only of the order of

10^{- 6}

, in strong contrast with 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 [38,39,40]. Perturbation at the interface between the two fluids causes the heavier one to penetrate the lighter fluid. In general, the dynamics of this system are governed by two nondimensional parameters, namely, the Atwood and Reynolds numbers. The former is defined as

At = \frac{ρ_{h} - ρ_{l}}{ρ_{h} + ρ_{l}},

(37)

while the latter is:

Re = \frac{ρ_{h} U^{*} L}{μ_{h}},

(38)

where

ρ_{l}

and

ρ_{h}

are densities of the heavy and light fluids, respectively,

μ_{h}

is the dynamic viscosity of the heavy fluid,

L_{x}

the size of the domain in the horizontal direction and

U^{*}

the characteristic velocity, defined as

U^{*} = \sqrt{g L_{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 setup studied in [19], we considered a domain sized

L_{x} \times 4 L_{x}

with

L_{x} = 600

. Initially, the top half of the domain was filled with the heavy liquid, and the bottom half with the lighter one. The interface was perturbed via the following profile:

h_{i} (x) = \frac{L}{10} cos (\frac{2 π x}{L_{x}}) + 2 L_{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 [1]. The At number was set to 0.5, while two different Re numbers were considered, i.e., Re = 256 and 2048. In both cases,

g = 6 \times 10^{- 6}

, while the nondimensional viscosities were 0.1406 and 0.0176, respectively. To validate the simulations, the position of the downward-plunging heavy liquid spike was measured over time and compared to the reference data from [19]. Results are illustrated in Figure 2.

Both simulations agreed very well with the reference solution of [19]. To showcase the ability of the solver to handle under-resolved simulations, and illustrate the convergence of the obtained solutions, simulations were repeated at two additional lower resolutions with

L_{x}

= 300 and 150, with an acoustic scaling of the time-step size. Results obtained with those lower resolutions are shown in Figure 3 and Figure 4.

The position of the plunging spike clearly shows that, while minor differences exist, even the lowest resolution captures the correct position. Smaller features, however, especially at Re = 2048, need higher resolutions to be correctly captured. At Re = 256 for instance, even the secondary instability was converged as, at

L_{x}

= 300, no segmentation was observed. For Re = 2056, on the other hand, while a larger structure started to converge, thinner features clearly needed more resolutions.

3.3. Turbulent 3D Rayleigh–Taylor Instability

To further showcase the ability of the solver to deal with complex flows, we also considered the Rayleigh–Taylor instability in 3D. The studied configuration followed those studied in [41]. The definitions of nondimensional parameters were similar to those used in the previous section. The domain was discretized using

100 \times 100 \times 1200

grid points, with

L = 100

. The interface was placed at the center of the domain along the z axis and perturbed using

h_{i} (x, y) = \frac{L}{10} [cos (\frac{2 π x}{L}) + cos (\frac{2 π y}{L})] + 6 L,

(42)

and Reynolds and Atwood numbers were set to 1000 and 0.15, respectively. 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 [41]. After the penetration of the two liquids into each other, the Kelvin–Helmholtz instability caused the plunging spike to roll up and take a mushroomlike shape. As the mushroom-shaped spike further progressed into the lighter fluid, the cap disintegrated into four fingerlike structures. As is shown later, these fingers were reminiscent of instabilities leading to splashing in the impact of a droplet on liquid surfaces.

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

3.4. Droplet Splashing on Thin Liquid Film

As the final case, we considered the impact of a droplet on a thin liquid layer. This configuration is interesting, as it involves complex dynamics such as splashing, and it is of interest in many areas of science and engineering [42,43]. Immediately after impact, the liquid surface is perturbed. In many instances, at the contact point (line), a thin liquid jet forms, and it continues to grow and propagate as a corolla. As the crownlike structure radially propagates, a rim starts to form. At high-enough Weber numbers, the structure breaks into small droplets via the Rayleigh–Plateau instability [44]. A detailed study of the initial stages of the spreading process showed that the spreading radius scales with time regardless of Weber and Reynolds numbers [44]. While widely studied in the literature using different numerical formulations [26,45,46,47], simulations are usually limited to lower density and viscosity ratios, and/or Weber and Reynolds numbers [26,36,45,46]. As such, we first focused on a 2D 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, density and viscosity ratios were set to

ρ_{h} / ρ_{l} = 1000

and

ν_{l} / ν_{h}

= 15, emulating a water/air system. The geometrical configuration is illustrated in Figure 6.

The top- and bottom-boundary conditions were set to walls modelled with the half-way bounce-back formulation, while symmetrical boundaries were applied to the left and right. The droplet diameter was resolved with 100 grid points. Initial velocity in the droplet was set to

U_{0} =

0.05, and

ν_{L}

was determined via the Reynolds number:

Re = \frac{ρ_{h} U_{0} D}{μ_{h}} .

(43)

Furthermore, the We number is defined as

We = \frac{ρ_{l} D {U_{0}}^{2}}{σ} .

(44)

The evolution of the liquid surface, as obtained from the simulations, is shown in Figure 7. Following [44], rim breakup and splashing occurred for larger impact parameters, defined as

K = {We}^{1 / 2} {Re}^{1 / 4} .

(45)

Accordingly, impact parameters for the studied 2D cases were K = 55.7, 83.4, and 263.8. The evolution of the systems in Figure 7 clearly shows that, in agreement with observations in [44], larger values of the impact parameter led to droplet detachment from the rim and splashing.

Furthermore, the evolution of spreading radii

r_{K}

over time for different cases is shown in Figure 8. The radii scaled with time at the initial stages of the impact, in agreement with results reported in [44].

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

After the initial impact, a thin liquid jet was formed at the contact line between the droplet and sheet. Then, the crown evolved and spread. At later stages, the fingerlike structures started to form at the tip of the crown. These liquid fingers then became detached from the crown, and liquid splashing was observed. This sequence of events was in excellent agreement with those presented in [44]. Furthermore, the spreading radius, as plotted in Figure 8, agreed with the theoretical predictions.

4. Conclusions

An LB-based solver relying on the conservative AC equation, and a modified hydrodynamic pressure/velocity-based distribution and MRT collision operator in cumulant space was presented in this study with the aim to model multiphase flows in larger Weber/Reynolds regimes. While stability at high Weber numbers, i.e., low surface tensions, is achieved through the decoupled nature of conservative AC formulation, the added stability in terms of kinematic viscosity, i.e., larger Reynolds numbers, is brought about by the collision operator and modified pressure-based LB formulation for the flow. Compared to other models available in the literature based on AC 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 3D droplet splashing were not stable with single relaxation time (SRT) formulation for the same choice of nondimensional parameters, i.e., resolution and relaxation factor. The algorithm was shown to capture flow dynamics 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.

Author Contributions

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. All authors have read and agreed to the published version of the manuscript.

Funding

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).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this study are available on request from the corresponding author.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. Hermite Polynomials and Coefficients

Hermite polynomials used in EDFs of different solvers, defined as

\begin{matrix} H_{0} & = 1, \end{matrix}

(A1a)

\begin{matrix} H_{i} & = c_{α, i}, \end{matrix}

(A1b)

\begin{matrix} H_{i j} & = c_{α, i} c_{α, j} - c_{s}^{2} δ_{i j}, \end{matrix}

(A1c)

where

δ_{i j}

denotes Kronecker delta function, while corresponding equilibrium coefficients are

\begin{matrix} a_{0}^{(e q)} & = ρ, \end{matrix}

(A2a)

\begin{matrix} a_{i}^{(e q)} & = ρ u_{i}, \end{matrix}

(A2b)

\begin{matrix} a_{i j}^{(e q)} & = ρ u_{i} u_{j}, \end{matrix}

(A2c)

References

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, Switzerland, 2017. [Google Scholar] [CrossRef]
Guo, Z.; Shu, C. Lattice Boltzmann Method and Its Applications in Engineering; Advances in Computational Fluid Dynamics; World Scientific: Singapore, 2013; Volume 3. [Google Scholar] [CrossRef] [Green Version]
Succi, S. The Lattice Boltzmann Equation for Fluid Dynamics and Beyond; Oxford University Press: Oxford, UK, 2002. [Google Scholar]
Chorin, A.J. A Numerical Method for Solving Incompressible Viscous Flow Problems. J. Comput. Phys. 1997, 135, 118–125. [Google Scholar] [CrossRef]
Shan, X.; Chen, H. Lattice Boltzmann model for simulating flows with multiple phases and components. Phys. Rev. E 1993, 47, 1815–1819. [Google Scholar] [CrossRef] [PubMed] [Green Version]
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. [Google Scholar] [CrossRef] [PubMed] [Green Version]
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. [Google Scholar] [CrossRef]
Swift, M.R.; Osborn, W.R.; Yeomans, J.M. Lattice Boltzmann Simulation of Nonideal Fluids. Phys. Rev. Lett. 1995, 75, 830–833. [Google Scholar] [CrossRef] [Green Version]
Wagner, A.; Li, Q. Investigation of Galilean invariance of multi-phase lattice Boltzmann methods. Phys. A Stat. Mech. Its Appl. 2006, 362, 105–110. [Google Scholar] [CrossRef] [Green Version]
Kupershtokh, A.; Medvedev, D.; Karpov, D. On equations of state in a lattice Boltzmann method. Comp. Math. Appl. 2009, 58, 965–974. [Google Scholar] [CrossRef] [Green Version]
Yuan, P.; Schaefer, L. Equations of state in a lattice Boltzmann model. Phys. Fluids 2006, 18, 042101. [Google Scholar] [CrossRef]
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. [Google Scholar] [CrossRef] [Green Version]
Li, Q.; Luo, K.H. Achieving tunable surface tension in the pseudopotential lattice Boltzmann modeling of multiphase flows. Phys. Rev. E 2013, 88, 053307. [Google Scholar] [CrossRef] [Green Version]
Fakhari, A.; Rahimian, M.H. Phase-field modeling by the method of lattice Boltzmann equations. Phys. Rev. E 2010, 81, 036707. [Google Scholar] [CrossRef]
Safari, H.; Rahimian, M.H.; Krafczyk, M. Extended lattice Boltzmann method for numerical simulation of thermal phase change in two-phase fluid flow. Phys. Rev. E 2013, 88, 013304. [Google Scholar] [CrossRef]
Safari, H.; Rahimian, M.H.; Krafczyk, M. Consistent simulation of droplet evaporation based on the phase-field multiphase lattice Boltzmann method. Phys. Rev. E 2014, 90, 033305. [Google Scholar] [CrossRef]
Yazdi, H.; Rahimiani, M.H.; Safari, H. Numerical simulation of pressure-driven phase-change in two-phase fluid flows using the Lattice Boltzmann Method. Comput. Fluids 2018, 172, 8–18. [Google Scholar] [CrossRef]
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. [Google Scholar] [CrossRef] [Green Version]
He, X.; Chen, S.; Zhang, R. A Lattice Boltzmann Scheme for Incompressible Multiphase Flow and Its Application in Simulation of Rayleigh–Taylor Instability. J. Comput. Phys. 1999, 152, 642–663. [Google Scholar] [CrossRef]
Inamuro, T.; Ogata, T.; Tajima, S.; Konishi, N. A lattice Boltzmann method for incompressible two-phase flows with large density differences. J. Comput. Phys. 2004, 98, 628–644. [Google Scholar] [CrossRef]
Amirshaghaghi, H.; Rahimian, M.; Safari, H. Application of a two phase lattice Boltzmann model in simulation of free surface jet impingement heat transfer. Int. Commun. Heat Mass Transf. 2016, 75, 282–294. [Google Scholar] [CrossRef]
Amirshaghaghi, H.; Rahimian, M.H.; Safari, H.; Krafczyk, M. Large Eddy Simulation of liquid sheet breakup using a two-phase lattice Boltzmann method. Comput. Fluids 2018, 160, 93–107. [Google Scholar] [CrossRef]
Hosseini, S.A. Development of a Lattice Boltzmann-Based Numerical Method for the Simulation of Reacting Flows. Ph.D. Thesis, Otto-von-Guericke Universität/Universite Paris-Saclay, Gif-sur-Yvette, France, 2020. [Google Scholar]
Sun, Y.; Beckermann, C. Sharp interface tracking using the phase-field equation. J. Comput. Phys. 2007, 220, 626–653. [Google Scholar] [CrossRef]
Chiu, P.H.; Lin, Y.T. A conservative phase field method for solving incompressible two-phase flows. J. Comput. Phys. 2011, 230, 185–204. [Google Scholar] [CrossRef]
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. J. Comput. Phys. 2017, 341, 22–43. [Google Scholar] [CrossRef] [Green Version]
Wang, H.; Chai, Z.; Shi, B.; Liang, H. Comparative study of the lattice Boltzmann models for Allen-Cahn and Cahn-Hilliard equations. Phys. Rev. E 2016, 94, 033304. [Google Scholar] [CrossRef] [Green Version]
Chopard, B.; Falcone, J.L.; Latt, J. The lattice Boltzmann advection-diffusion model revisited. Eur. Phys. J. Spec. Top. 2009, 171, 245–249. [Google Scholar] [CrossRef]
Hosseini, S.A.; Darabiha, N.; Thévenin, D. Lattice Boltzmann advection-diffusion model for conjugate heat transfer in heterogeneous media. Int. J. Heat Mass Transf. 2019, 132, 906–919. [Google Scholar] [CrossRef] [Green Version]
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. Phys. Rev. E 2020, 102, 053307. [Google Scholar] [CrossRef]
Lee, T.; Lin, C.L. Pressure evolution lattice Boltzmann equation method for two-phase flow with phase change. Phys. Rev. E 2003, 67, 056703. [Google Scholar] [CrossRef]
Hosseini, S.A.; Safari, H.; Darabiha, N.; Thévenin, D.; Krafczyk, M. Hybrid Lattice Boltzmann-finite difference model for low Mach number combustion simulation. Combust. Flame 2019, 209, 394–404. [Google Scholar] [CrossRef]
Hosseini, S.A.; Abdelsamie, A.; Darabiha, N.; Thévenin, D. Low-Mach hybrid lattice Boltzmann-finite difference solver for combustion in complex flows. Phys. Fluids 2020, 32, 077105. [Google Scholar] [CrossRef]
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. [Google Scholar] [CrossRef]
Geier, M.; Schönherr, M.; Pasquali, A.; Krafczyk, M. The cumulant lattice Boltzmann equation in three dimensions: Theory and validation. Comput. Math. Appl. 2015, 70, 507–547. [Google Scholar] [CrossRef] [Green Version]
Qin, F.; Mazloomi Moqaddam, A.; Kang, Q.; Derome, D.; Carmeliet, J. Entropic multiple-relaxation-time multirange pseudopotential lattice Boltzmann model for two-phase flow. Phys. Fluids 2018, 30, 032104. [Google Scholar] [CrossRef]
Mazloomi M, A.; Chikatamarla, S.; Karlin, I. Entropic Lattice Boltzmann Method for Multiphase Flows. Phys. Rev. Lett. 2015, 114, 174502. [Google Scholar] [CrossRef]
Yang, X.; He, H.; Xu, J.; Wei, Y.; Zhang, H. Entropy generation rates in two-dimensional Rayleigh–Taylor turbulence mixing. Entropy 2018, 20, 738. [Google Scholar] [CrossRef] [Green Version]
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. [Google Scholar] [CrossRef] [Green Version]
Rahmat, A.; Tofighi, N.; Shadloo, M.; Yildiz, M. Numerical simulation of wall bounded and electrically excited Rayleigh–Taylor instability using incompressible smoothed particle hydrodynamics. Coll. Surf. A Physicochem. Eng. Asp. 2014, 460, 60–70. [Google Scholar] [CrossRef]
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. [Google Scholar] [CrossRef]
Hagemeier, T.; Hartmann, M.; Thévenin, D. Practice of vehicle soiling investigations: A review. Int. J. Multiph. Flow 2011, 37, 860–875. [Google Scholar] [CrossRef]
Hagemeier, T.; Hartmann, M.; Kühle, M.; Thévenin, D.; Zähringer, K. Experimental characterization of thin films, droplets and rivulets using LED fluorescence. Exp. Fluids 2012, 52, 361–374. [Google Scholar] [CrossRef]
Josserand, C.; Zaleski, S. Droplet splashing on a thin liquid film. Phys. Fluids 2003, 15, 1650. [Google Scholar] [CrossRef]
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. [Google Scholar] [CrossRef]
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. [Google Scholar] [CrossRef] [Green Version]
Sitompul, Y.P.; Aoki, T. A filtered cumulant lattice Boltzmann method for violent two-phase flows. J. Comput. Phys. 2019, 390, 93–120. [Google Scholar] [CrossRef]

Figure 1. Changes in pressure difference around droplet for different surface tensions and droplet radii. Red, blue, and black symbols illustrate results from present study with

σ = 10^{- 1}, 10^{- 3}

, and

10^{- 6}

, respectively.

Figure 1. Changes in pressure difference around droplet for different surface tensions and droplet radii. Red, blue, and black symbols illustrate results from present study with

σ = 10^{- 1}, 10^{- 3}

, and

10^{- 6}

, respectively.

Figure 2. (Left) Evolution of interface for 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. (Right) Position of penetrating spike over time: (black) Re = 256 and (red) Re = 2048. (plain lines) Results and (symbols) data from [19].

Figure 2. (Left) Evolution of interface for 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. (Right) Position of penetrating spike over time: (black) Re = 256 and (red) Re = 2048. (plain lines) Results and (symbols) data from [19].

Figure 3. (Left) Interface for Rayleigh–Taylor instability at

t / T =

5 and Re = 256 for three different resolutions (left to right)

L_{x}

= 150, 300, and 600. (Right) Position of penetrating spike over time: (black)

L_{x}

= 600, (red)

L_{x}

= 300, and (blue)

L_{x}

= 150.

Figure 3. (Left) Interface for Rayleigh–Taylor instability at

t / T =

5 and Re = 256 for three different resolutions (left to right)

L_{x}

= 150, 300, and 600. (Right) Position of penetrating spike over time: (black)

L_{x}

= 600, (red)

L_{x}

= 300, and (blue)

L_{x}

= 150.

Figure 4. (Left) Interface for Rayleigh–Taylor instability at

t / T

= 5 and Re = 2048 for three different resolutions (left to right)

L_{x}

= 150, 300, and 600. (Right) Position of penetrating spike over time: (black)

L_{x}

= 600, (red)

L_{x}

= 300, and (blue)

L_{x}

= 150.

Figure 4. (Left) Interface for Rayleigh–Taylor instability at

t / T

= 5 and Re = 2048 for three different resolutions (left to right)

L_{x}

= 150, 300, and 600. (Right) Position of penetrating spike over time: (black)

L_{x}

= 600, (red)

L_{x}

= 300, and (blue)

L_{x}

= 150.

Figure 5. (Left) Evolution of interface for 3D 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. (Right) Position of penetrating spike over time: (plain lines) Results and (symbols) data from [41].

Figure 5. (Left) Evolution of interface for 3D 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. (Right) Position of penetrating spike over time: (plain lines) Results and (symbols) data from [41].

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

Figure 7. Impact of circular droplet on liquid sheet at different We and Re numbers with

ρ_{h} / ρ_{l} = 1000

and

ν_{l} / ν_{h} = 15

. (black) Re = 200 and We = 220, (red) Re = 1000 and We = 220, and (blue) Re = 1000 and We = 2200.

Figure 7. Impact of circular droplet on liquid sheet at different We and Re numbers with

ρ_{h} / ρ_{l} = 1000

and

ν_{l} / ν_{h} = 15

. (black) Re = 200 and We = 220, (red) Re = 1000 and We = 220, and (blue) Re = 1000 and We = 2200.

Figure 8. Evolution of spreading radius

r_{K}

as function of time for droplet impact on liquid film case. Circular symbols designate 2D simulations: (black) Re = 200 and We = 220, (red) Re = 1000 and We = 220, and (blue) Re = 1000 and We = 2200. Rectangular symbols belong to 3D simulation with Re = 1000 and We = 8000. Dashed line is

\frac{r_{K}}{D} = 1.1 \sqrt{t / T}

Figure 8. Evolution of spreading radius

r_{K}

\frac{r_{K}}{D} = 1.1 \sqrt{t / T}

Figure 9. Impact of spherical droplet on thin liquid sheet at We = 8000 and Re = 1000 at different times with

ρ_{h} / ρ_{l} = 1000

and

ν_{l} / ν_{h} = 15

Figure 9. Impact of spherical droplet on thin liquid sheet at We = 8000 and Re = 1000 at different times with

ρ_{h} / ρ_{l} = 1000

and

ν_{l} / ν_{h} = 15

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Hosseini, S.A.; Safari, H.; Thevenin, D. Lattice Boltzmann Solver for Multiphase Flows: Application to High Weber and Reynolds Numbers. Entropy 2021, 23, 166. https://doi.org/10.3390/e23020166

AMA Style

Hosseini SA, Safari H, Thevenin D. Lattice Boltzmann Solver for Multiphase Flows: Application to High Weber and Reynolds Numbers. Entropy. 2021; 23(2):166. https://doi.org/10.3390/e23020166

Chicago/Turabian Style

Hosseini, Seyed Ali, Hesameddin Safari, and Dominique Thevenin. 2021. "Lattice Boltzmann Solver for Multiphase Flows: Application to High Weber and Reynolds Numbers" Entropy 23, no. 2: 166. https://doi.org/10.3390/e23020166

APA Style

Hosseini, S. A., Safari, H., & Thevenin, D. (2021). Lattice Boltzmann Solver for Multiphase Flows: Application to High Weber and Reynolds Numbers. Entropy, 23(2), 166. https://doi.org/10.3390/e23020166

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu

Lattice Boltzmann Solver for Multiphase Flows: Application to High Weber and Reynolds Numbers

Abstract

1. Introduction

2. Theoretical Background

2.1. Target Macrosopic System

2.2. LB Formulation for Conservative Phase-Field Equation

2.3. LB Model for Flow Field

3. Numerical Applications

3.1. Static Droplet: Surface-Tension Measurement

3.2. Rayleigh–Taylor Instability

3.3. Turbulent 3D Rayleigh–Taylor Instability

3.4. Droplet Splashing on Thin Liquid Film

4. Conclusions

Author Contributions

Funding

Institutional Review Board Statement

Informed Consent Statement

Data Availability Statement

Conflicts of Interest

Appendix A. Hermite Polynomials and Coefficients

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI