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

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.


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

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: where u i is fluid velocity, ρ the fluid density, and F b,i designates external body forces. The stress tensor σ ij is defined as: where η is the fluid dynamic viscosity tied to kinematic viscosity ν as η = ρν, ξ the bulk viscosity and p h the hydrodynamic pressure. The chemical potential µ φ is defined as where ∆ = ∇ 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]: 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 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.

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: where g α and c α are populations and velocities in the discrete velocity kinetic model, and the collision operator is defined as The EDF is defined as where H n and a (eq) n 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] 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 The relaxation coefficient is fixed as After integration in space/time, the now-famous collision-streaming form can be recovered: where the source term takes on a new form, i.e., and:τ The derivatives of the order parameter appearing in the various discrete time-evolution equations are computed using isotropic finite differences, i.e., and 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]: which as for Equation (13), postdiscretization changes intō 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].

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] where the collision operator is Ξ α is defined as and the relaxation coefficient τ is tied to fluid kinematic viscosity ν as Forces F b,i and F s,i represent external body forces and surface tension, respectively, i.e., The modified distribution function f α is defined as where f α is the classical isothermal distribution function. The modified equilibrium follows the same logic and is defined as Density is tied to the order parameter as 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 are computed as The postcollision preconditioned population f p * α is 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 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 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: Here, β = x n x y n y z n z . The central moments are then transformed into the corresponding cumulants using the following relations: 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 p xy , K p ). 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 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: where β = x 2 − y 2 stands for a central moment of form Π p x 2 − Π p y 2 . Previous systematic studies of the flow solver showed second-order convergence under diffusive scaling [32].

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.

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 × 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., 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 × 10 −1 , 1 × 10 −3 , and 1 × 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 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.

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 while the latter is: 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 where g is gravity-driven acceleration. The characteristic time for this case is defined as Following the setup studied in [19], we considered a domain sized L x × 4L 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: 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 × 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 Figures 3 and 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.

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 × 100 × 1200 grid points, with L = 100. The interface was placed at the center of the domain along the z axis and perturbed using 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 × 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.

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 halfway 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: Furthermore, the We number is defined as The evolution of the liquid surface, as obtained from the simulations, is shown in Figure 7.
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.

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.