Phase-field-based lattice Boltzmann modeling of large-density-ratio two-phase flows

# Phase-field-based lattice Boltzmann modeling of large-density-ratio two-phase flows

Hong Liang Department of Physics, Hangzhou Dianzi University - Hangzhou 310018, China    Jiangrong Xu Department of Physics, Hangzhou Dianzi University - Hangzhou 310018, China    Jiangxing Chen Department of Physics, Hangzhou Dianzi University - Hangzhou 310018, China    Huili Wang School of Mathematics and Statistics, Huazhong University of Science and Technology, Wuhan, 430074, China    Zhenhua Chai School of Mathematics and Statistics, Huazhong University of Science and Technology, Wuhan, 430074, China Hubei Key Laboratory of Engineering Modeling and Scientific Computing, Huazhong University of Science and Technology, Wuhan 430074, China    Baochang Shi School of Mathematics and Statistics, Huazhong University of Science and Technology, Wuhan, 430074, China Hubei Key Laboratory of Engineering Modeling and Scientific Computing, Huazhong University of Science and Technology, Wuhan 430074, China
July 10, 2019
###### Abstract

In this paper, we present a simple and accurate lattice Boltzmann (LB) model for immiscible two-phase flows, which is able to deal with large density contrasts. This model utilizes two LB equations, one of which is used to solve the conservative Allen-Cahn equation, and the other is adopted to solve the incompressible Navier-Stokes equations. A novel forcing distribution function is elaborately designed in the LB equation for the Navier-Stokes equations, which make it much simpler than the existing LB models. In addition, the proposed model can achieve superior numerical accuracy compared with previous Allen-Cahn type of LB models. Several benchmark two-phase problems, including static droplet, layered Poiseuille flow, and Spinodal decomposition are simulated to validate the present LB model. It is found that the present model can achieve relatively small spurious velocity in the LB community, and the obtained numerical results also show good agreement with the analytical solutions or some available results. At last, we use the present model to investigate the droplet impact on a thin liquid film with a large density ratio of 1000 and the Reynolds number ranging from 20 to 500. The fascinating phenomena of droplet splashing is successfully reproduced by the present model and the numerically predicted spreading radius exhibits to obey the power law reported in the literature.

###### pacs:
47.11.-j 47.55.-t 68.03.-g
preprint: APS/123-QED

r

## I Introduction

Two-phase flows are ubiquitous in nature and engineering applications. Numerical modeling of such flows becomes an important complement to experimental studies with the rapid development of computational science, while it may face some certain challenges owing to complex interfacial dynamics involving multiple space and time scales. Physically, the interfacial phenomenon can be recognized as a natural consequence of intermolecular interactions. In this regard, the lattice Boltzmann (LB) method Succi (); Guo (), based on the mesoscopic kinetic theory, becomes a suitable candidate to model and simulate two-phase flows.

Over the past three decades, the LB method has received great success in modeling multi-phase fluid systems Succi (); Guo () and some nonlinear equation systems Shi (); Chai (). The reasons behind its success lie in the algorithmic simplicity, nature parallelization and easy implementation of complex boundary. Additionally, thanks to the kinetic nature, the LB method can deal with the intermolecular interactions in a straightforward manner, which is also regarded as its unique advantage that distinguishes from the traditional CFD methods. Up to now, a variety of LB models for multiphase flows have been proposed, which mainly fall into four categories, including color-gradient model Gunstensen (), pseudo-potential model Shan (), free energy model Swift () and phase-field based model He (); Lee1 (); Lee2 (); Liang1 (); Liang2 (). For the detailed expositions, the readers can refer to the recent reviews Liu (); Li1 () on LB approaches for multiphase flows and the references therein.

Most of the previously proposed LB models are only able to handle two-phase flows with small or moderate density contrasts. Generally, the density ratio of liquid and vapor phases is larger than 100, and it even could approach 1000 for a realistic water-air two-phase system. Within this context, to develop a multiphase model that can simulate large-density-ratio flows is an attractive topic in LB community. Inamuro et al. Inamuro () proposed a first LB model based on the free-energy method that can tolerate large density differences. While they needs to solve an additional Poisson equation for the pressure to enforce the incompressible condition, which seems to be complex, and undermines the simplicity of LB method. Besides, an empirical cutoff value is used to determine fluid density, which could lead to the violation of the mass conservation, like the level set method Sethian (). The extension of the original pseudo-potential model Shan () to large density ratio cases was attributed to Yuan and Schaefer Yuan (). They evaluated the performances of different equations of state in the pseudo-potential model and found that large density ratio can be reached with a suitable choice of equation of state. While it is noticed that their studies only focus on the stationary two-phase problems, and the model will suffer from some limitations more or less when it is readily applied to dynamic problems. To remove this limitation, Li et al. Li2 () presented an improved pseudo-potential model that can satisfy thermodynamic consistency. Meanwhile it can improve numerical stability of the pseudo-potential method at a large density ratio for dynamic flows, which was demonstrated by the simulation of droplet splashing on a liquid film with the largest density ratio of 700. But they also declared in the literature that it would induce numerical instability when the density ratio was increased to 1000 in this case. Ba et al. Ba () also developed a color-gradient-based LB model for simulating two-phase flows with high density ratio, in which a modified equilibrium distribution function and a simple source term are introduced. They significantly improved the performance of the color-gradient method and achieved satisfactory results in the simulation of droplet splashing problem with the largest density ratio of 100. On the other hand, several researchers have also attempted to develop a large-density-ratio LB model based on the phase-field theory, which has become increased popular in modeling multiphase flows Jacqmin (). Zheng et al. Zheng () proposed a LB two-phase model based on the Cahn-Hilliard phase field equation and claimed that their model can simulate large-density-ratio flows. Actually, it is noted that they only consider the Navier-Stokes equations on the average density of binary fluids instead of the real fluid density, and therefore their model in theory is only able to deal with density-matched binary fluids, which is also numerically proved by Fakhari and Rahimian Fakhari1 (). Lee et al. Lee1 (); Lee2 () also presented another LB model for large-density-ratio two-phase flows from the phase field viewpoint. The key point of their model in achieving large density ratio is the use of a stable mixing difference scheme for computing gradient terms, while it will induce the violation of mass and momentum conservation Lou (). Besides, the inconsistency between the recovered interfacial equation and the target equation in their models was also found Zu (); Liang1 (). Wang et al. Wang () proposed an interesting LB flux model for two-phase flows with large density ratios, in which a stable high-order WENO difference scheme is used to solve the Cahn-Hilliard equation, and a like finite volume method for particle distribution function is utilized to solve the incompressible Navier-Stokes equations. Recently, Ren et al. Ren () proposed a LB model from the perspective of the Allen-Cahn phase field equation, while they only concentrated on two-phase flows limited to small or moderate density ratio, and whether it can be applicable for large-density-ratio flows has not been discussed. Besides, their model contains many complex gradient terms, and seems to be implemented difficultly. More recently, Fakhari and Bolster Fakhari2 () developed a simple LB model based on the Allen-Cahn phase field equation that can simulate large-density-ratio two-phase flows. However, it is found that the model contains some artificial terms in the recovered interfacial equation, which may affect the numerical accuracy HWang (); Ren ().

In this paper, we intend to present a simple, accurate and also robust two-phase model for large-density-ratio flows in the framework of LB method. The proposed LB model is based on the Allen-Cahn phase field theory, which only contains a most second-order gradient term. Therefore the present model can achieve a lower numerical dispersion in interface tracking, compared with the previous LB models He (); Lee1 (); Lee2 (); Fakhari1 (); Zu (); Liang1 () based on the four-order Cahn-Hilliard equation. This is the main contribution to obtain a better numerical stability at a large density ratio. In addition, a novel force distribution function is introduced in this model, which can be much simpler than those of the existing LB models Lee1 (); Lee2 (); Zu (); Liang1 (); Ren (); Fakhari2 (). The inconsistency of the recovered interfacial equation with the target equation in Fakhari’s model Fakhari2 () is also remedied in this model by the incorporation of a proper source term Ren (); HWang (). Through the Chapman-Enskog analyis, our model can recover both the conservative Allen-Cahn and incompressible Navier-Stokes equations correctly, which can be more accurate than the previous Allen-Cahn based LB models Ren (); Fakhari2 (). The rest of the paper is arranged as follows. In Sec. II, the macroscopic governing equations are first given, and a novel LB model for two-phase flows based on the Allen-Cahn phase field theory is then presented. Numerical experiments to validate the present model and a detailed comparison with some previous LB models can be found in Sec. III, and finally we made a brief summary in Sec. IV.

## Ii LB model for two-phase flows

In this section, we first give a brief introduction on the governing equations in the framework of the Allen-Cahn phase-field theory Sun (); Chiu (), and then present a LB model for two-phase incompressible flows. Based on the collision operator used, the LB method can be roughly divided into three categories: the single-relaxation-time or so-called BGK method Qian (), the two-relaxation-time method Ginzburg (), and the multiple-relaxation-time (MRT) method Lallemand (). Considering its simplicity and high computational efficiency, the present model is constructed based on the BGK collision operator and its extension to the advanced MRT version can be conducted directly, which constitutes one of our future research branches.

### ii.1 Governing equations

The conservative Allen-Cahn equation can be expressed by Sun (); Chiu ()

 ∂ϕ∂t+∇⋅(ϕu)=∇⋅[M(∇ϕ−λn)], (1)

where is the mobility, is the unit vector normal to the interface,

 n=∇ϕ|∇ϕ|, (2)

and is a function of defined by

 λ=4ϕ(1−ϕ)W, (3)

where is the interface thickness, taking and represents the liquid and gas phase fluids, respectively, and the interface is marked by the contour level of . Here we consider the incompressible two-phase flows, and the fluid velocity in Eq. (1) is governed by the following Navier-Stokes equations with the force Unverdi (),

 ∇⋅u=0, (4a) ∂(ρu)∂t+∇⋅(ρuu)=−∇p+∇⋅[μ(∇u+∇uT)]+Fs+G, (4b)

where is the fluid density, is the hydrodynamic pressure, is the dynamic viscosity by , is the kinematic viscosity, is the surface tension force, and is the possible body force. In the literature Kim (), there exists several different forms of the surface tension force, here we choose the widely used one of the potential form in the phase field methods Jacqmin (); Zu (); Liang1 (); Liang3 (),

 Fs=μϕ∇ϕ, (5)

where is the chemical potential defined by

 μϕ=4βϕ(ϕ−1)(ϕ−0.5)−k∇2ϕ, (6)

where and are physical parameters that depend on the interface thickness and the surface tension (),

 k=32σW,  β=12σW. (7)

### ii.2 LB model for the conservative Allen-Cahn equation

The LB evolution equation with the BGK collision operator for the conservative Allen-Cahn equation can be written as Shi (); HWang ()

 fi(x+ciδt,t+δt)−fi(x,t)=−1τf[fi(x,t)−feqi(x,t)]+δtFi(x,t), (8)

where is the particle distribution function, is the non-dimensional relaxation time related to the mobility, is the source term, and a simple form of the equilibrium distribution function is given by

 feqi=ωiϕ(1+ci⋅uc2s) (9)

where is the sound speed, are the discrete velocities, and are the weighting coefficients. and depend on the choice of the lattice model. For the two-dimensional flows considered here, the D2Q5 or D2Q9 lattice model can be applied in the LB algorithm for Allen-Cahn equation. Considering the consistency with the LB algorithm for Navier-Stokes equations, in this work we adopt the popular D2Q9 lattice model Qian (); Liang1 (); Liang2 (); Wei1 (). Then the weighting coefficients can be given by , , , and the discrete velocities are

 ci=⎧⎪ ⎪⎨⎪ ⎪⎩(0,0)c, i=0,(cos[(i−1)π/2],sin[(i−1)π/2])c, i=1−4,√2(cos[(i−5)π/2+π/4],sin[(i−5)π/2+π/4])c, i=5−8, (10)

where / is the lattice speed with , representing the grid spacing and the time increment respectively, and . By convention, and are set as the length and time units, i.e., .

To recover the Allen-Cahn equation exactly with the multi-scale analysis, the source term in Eq. (8) should be defined as Ren (); HWang ()

 Fi=(1−12τf)ωici⋅[∂t(ϕu)+c2sλn]c2s, (11)

where the time derivative term is introduced to eliminate the artificial term in the recovered equation, which is similar to the technique used in LB models Liang1 (); Liang2 (); Liang3 () for the Cahn-Hilliard Equation. One notice that in the existing LB model Fakhari2 () based on the Allen-Cahn theory, the term is not included, which results in the deviation between the recovered equation and the target equation Ren (); HWang ().

The order parameter in the present model can be computed by

 ϕ=∑ifi. (12)

The distribution of fluid density in a two-phase system physically is consistent with that of the order parameter. To satisfy this physical property, the fluid density should take the linear interpolation,

 ρ=ϕ(ρl−ρg)+ρg, (13)

where and represent the densities of the liquid and gas phases. Following the Chapman-Enskog analysis in Ref. HWang (), it is found that the conservative Allen-Cahn equation can be recovered correctly from the LB equation (8) and the mobility can be determined by

 M=c2s(τf−0.5)δt. (14)

### ii.3 LB model for the Navier-Stokes Equations

The LB equation with the BGK collision operator for the Navier-Stokes Equations can be expressed as Guo1 (); Wei2 ()

 gi(x+ciδt,t+δt)−gi(x,t)=−1τg[gi(x,t)−geqi(x,t)]+δtGi(x,t), (15)

where is the distribution function for solving the flow field, is its corresponding equilibrium distribution function, is the dimensionless relaxation time related to the viscosity, is the force distribution function. To satisfy the divergence-free condition of velocity, should be elaborately designed as Liang1 (); Liang3 ()

 geqi=⎧⎪⎨⎪⎩pc2s(ωi−1)+ρsi(u), i=0,pc2sωi+ρsi(u), i≠0 (16)

with

 si(u)=ωi[ci⋅uc2s+(ci⋅u)22c4s−u⋅u2c2s]. (17)

For the two-dimensional flows, the D2Q9 lattice model is also adopted for flow field and the related physical coefficients , are also chosen as those of the previous section.

Different from the previous LB models Lee1 (); Lee2 (); Zu (); Liang1 (); Ren (); Fakhari2 (), a novel force distribution function is given by

 Gi=(1−12τg)ωi[u⋅∇ρ+ci⋅Fc2s+u∇ρ:(cici−c2sI)c2s], (18)

where is the total force,

 F=Fs+G. (19)

We would like to point out that the force distribution function given in Eq. (18) can be much simpler than those of the previous models Lee1 (); Lee2 (); Zu (); Liang1 (); Ren (); Fakhari2 (). In addition, the present LB model with the force term (18) can recover the Navier-Stokes equations correctly using the Chapman-Enskog analysis (see Appendix A for the detais). Substituting Eqs. (5) and (13) into Eq. (18), one can further simplify Eq. (18) as

 Gi=(1−12τg)ωi[ci⋅(μϕ∇ϕ+G)c2s+(ρl−ρg)u∇ϕ:cicic2s]. (20)

Taking the zeroth- and the first order moments of the distribution function , the macroscopic quantities and can be evaluated as Liang1 (); Liang3 ()

 ρu=∑icigi+0.5δtF, (21a) p=c2s(1−ω0)⎡⎣∑i≠0gi+δt2u⋅∇ρ+ρs0(u)⎤⎦, (21b)

which can be further recast as

 ρu=∑icigi+0.5δt(μϕ∇ϕ+G), (22a) p=c2s(1−ω0)⎡⎣∑i≠0gi+δt2(ρl−ρg)u⋅∇ϕ+ρs0(u)⎤⎦, (22b)

with the substitution of Eqs. (5) and (13). Based on the Chapman-Enskog analysis, the fluid kinematic viscosity can be determined by

 ν=c2s(τg−0.5)δt. (23)

In a two-phase system, the viscosity is no longer a uniform value due to its jump at the liquid-gas interface. There are several manners to treat the viscosity across the interface. To be smooth across the interface, the viscosity in the diffusion-interface methods is usually supposed to be a linear function of the order parameter He (),

 ν=ϕ(νl−νg)+νg, (24)

where and are the kinematic viscosities of the liquid and gas phases. In addition to Eq. (24), another popular treatment to determine the viscosity is the inverse linear form Lee2 (); Fakhari2 (),

 1ν=ϕ(1νl−1νg)+1νg. (25)

Oftentimes, to avoid the sharp-interface limit of the phase-field methods, a step function is also applied for the dynamic viscosity Ren (),

 μ={μl,   ϕ≥0.5,μg,   ϕ<0.5, (26)

where and are the dynamic viscosities of two different phases. The scheme (26) can achieve a considerable accuracy in tracking the interface, while similar to the sharp-interface methods, it could be unstable when it is applied to interfacial dynamic problems with large topology change Anderson (). In this work, a simple linear form as used for density is adopted, if not specified.

For numerical iterations, the derivative terms in the model should be discretized with suitable difference schemes. As commonly used in LB literatures Zu (); Liang1 (); Liang2 (); Liang3 (), the gradient term is computed by the second-order isotropic central scheme,

 ∇ϕ(x)=∑i≠0ωiciϕ(x+ckδt)c2sδt (27)

and the laplace operator is calculated by

 ∇2ϕ(x)=∑i≠02ωi[ϕ(x+ciδt)−ϕ(x)]c2sδ2t. (28)

Occasionally, the gradient term can be computed with the nonequilibrium part in some certain LB approaches Chai (); HWang () for the convection-diffusion equations. Through the multiscale analysis, it is shown that the gradient of and its gradient norm in this model can be computed by the local nonequilibrium schemes HWang ()

 |∇ϕ|=−|C|−BA, (29a) ∇ϕ=CA+B|∇ϕ|, (29b)

where , , and . When the and are computed by Eqs. (29a) and (29b), the unit normal vector n then can be obtained. This treatment on the gradient term enables the collision process being implemented locally, in addition to the computation of , which is one of the striking features in LB approaches. Therefore the scheme will be adopted in our simulations, unless otherwise stated. However, we find that the velocity actually satisfies an implicit equation, when Eq. (29a) and (29b) are applied for the statistics of the velocity, and for simplicity we just applied Eqs. (27) and (28) in this step.

At the end of this section, we intend to give some remarks on the present model for two-phase flows. Firstly, the present model is developed based on the conservative Allen-Cahn Equation, which contains a lower-order diffusion term compared with the fourth-order Cahn-Hilliard equation. Therefore the Allen-Cahn based model in theory can achieve a lower numerical dispersion in solving the index function and also the density field via Eq. (13) than the Cahn-Hilliard type of LB models. The lower-dispersion solution of plays a significant role in simulating large density ratio flows since a small deviation could be more likely to lead a unphysically negative value of fluid density, causing numerical instabilities. It is worth noting that a type of large-density-ratio LB models Lee1 (); Lee2 () have been proposed based on the Cahn-Hilliard equation, which is attributed to the use of a mixed scheme that combines the central and biased differences. However, it can lead to the violations of mass and momentum conservation Lou (). On the contrary, in our model the isotropic central scheme and the local nonequilibrium scheme are applied, which not only preserve a secondary-order accuracy in space, but also can ensure the global mass conservation of a two-phase system. Secondly, a novel force distribution function for flow field is proposed in the present model, which can be much simpler than those of the existing Allen-Cahn based LB models Ren (); Fakhari2 (). It is noted that our model only contains one type of non-local gradient term for the order parameter, which is much smaller than those of the previous model Ren (). In addition, the gradient term and its Laplacian in our model can be computed with local nonequilibrium schemes, which makes the collision process to be conducted locally if has been given. While in the previous models Ren (); Fakhari2 (), only the central difference schemes are applied, and thus the collision process cannot be conducted locally. Thirdly, both the conservative Allen-Cahn Equation and the incompressible Navier-Stokes equations can be recovered exactly from the present model with the multi-scale analysis. While the model of Fakhari et al. Fakhari2 () contains some artificial terms in the recovered interfacial equation. The numerical experiments conducted below will demonstrate that the present model is more accurate than the previous Allen-Cahn based LB models Ren (); Fakhari2 (). At last, we would like to stress that the present model is a standard LB scheme for simulating large-density-ratio two-phase flows without the use of an advanced finite difference or finite volume method Wang (), therefore it can naturally inherit the advantages of LB method in dealing with complex physical boundary and parallel computing.

## Iii Numerical Results and discussions

In this section, several typical benchmark problems, including static droplet, layered Poiseuille flows, and Spinodal decomposition are used to validate the present LB model for large density ratio flows. We attempt to conduct a detailed comparison between the present results and the analytical solutions or some available results. At last, we also investigated droplet impact on a thin liquid film, where the effect of the Reynolds number are discussed in detail.

### iii.1 Static droplet

The static droplet is a basic two-phase problem, which has been widely used to verify the developed numerical methods Lee1 (); Fakhari1 (); Zu (); Liang1 (); Ba (); Li2 (). In this subsection, we will simulate this problem with large density ratio to validate the present LB model. Initially, a liquid droplet with the radius of surrounded by the gas phase is located at the center of the square domain and the periodic boundary conditions are applied at all boundaries. The distribution profile of the order parameter is initialized by

 ϕ(x,y)=0.5+0.5tanh2[R−(x−100)2−(y−100)2]W, (30)

which makes its value to be smooth across the interface. In the simulation, we set the density ratio to be , and some other physical parameters are given as , , . Fig. 1(a) depicts the interface pattern of the droplet at the equilibrium state, together with the initial one given by Eq. (30). It can be found that they line up over each other exactly, which indicates that the present model has a high accuracy in track the interface. Furthermore, we quantitatively plotted in Fig. 1(b) the density distributions along the horizontal center line with different values of the mobility . It is shown that numerical predictions of the density field are all in good agreement with the analytical solution.

The spurious velocity around the interface is a commonly concerned problem in LB approaches for two-phase flows, and cannot be completely eliminated in the framework of the LB method Guo2 (). In Fig. 1(a), we also display the velocity distribution in the whole computational domain obtained by the present model. It can be found that the spurious velocities indeed exist at the vicinity of the interface, and their maximum magnitude computed by is about . We also compared the present model with some previously improved LB models in terms of the spurious velocity. It is shown that the maximum amplitude of spurious velocities in an improved Shan-Chen model Yu () has the order of . Recently, Ba et al. Ba () developed an improved color-gradient based model for high density ratio, which produced spurious velocities with the order of . As for the Cahn-Hilliard type of LB model Liang1 (), they can obtain spurious velocities at the level of . Therefore we can conclude that the present LB model is able to produce relatively small spurious velocities.

### iii.2 Layered Poiseuille flow

The layered Poiseuille flow is a classical two-phase problem, which can provide a good benchmark for validating the developed LB approaches Zu (); Ba (); Wang (); Ren (); Huang (); Liang4 (). While to our knowledge, most of the previous studies are limited to the small density ratio less than 10, due to the numerical instability problem. In this subsection, we will simulate the layered two-phase Poiseuille flows with the largest density ratio of 1000, and also conduct comparisons of the present model with the existing Allen-Cahn based LB models Ren (); Fakhari2 (). Consider a channel flow of two immiscible fluids driven by a constant body force . Initially, the gas phase fluid is placed in the upper region of and the region of is fill with the liquid phase fluid. Periodic boundary conditions are applied in the -direction, and the bottom and top boundaries are the solid walls, which are treated by the halfway bounce-back boundary conditions. Based on these boundary conditions, one can derive the analytical solution for the horizontal velocity profile (),

 ux(y)=⎧⎪ ⎪⎨⎪ ⎪⎩Gxh22μg[−(yh)2−yh(μg−μlμg+μl)+2μgμg+μl], 0

where , which provides a steady horizontal velocity of at the center. To quantitatively describe the accuracy of the present model and also conveniently compare with the existing LB models, the following relative error is used,

 Eu=∑y|unx(y,t)−uax(y)|∑y|uax(y)|, (32)

where the subscripts and denote the numerical and analytical solutions.

In the simulation, the computational grid is chosen to be , and the initial distribution of the order parameter is set as

 ϕ(x,y)=0.5+0.5tanh2(0.5Ny−y)W, (33)

which gives the profile of the planar interface. is fixed as a small value of , which ensures that the incompressible limit can be satisfied, and some other related parameters are given as , , and . Four different cases of the density ratios are considered, where the corresponding kinematic viscosity ratios are 1 for the former three cases, and for the latter case. Here the dynamic viscosity is given by Eq. (26) in the present test of our model. Note that the two-phase system with and considered here is very close to the realistic water-air system at the room temperature and normal atmospheric pressure. Figure 2 shows the profiles of the horizontal velocity () with various density ratios obtained by the present model, together with the corresponding analytical solutions. For comparisons, we also simulated the above cases with the previous Allen-Cahn based LB models Ren (); Fakhari2 () under identical computational conditions, and the obtained numerical results are also presented in Fig. 2. It can be observed from Fig. 2 that numerical results of the present model agree well with the analytical solutions for all density ratios, while some obvious discrepancies with the analytical solutions are found in the results of the existing Allen-Cahn based LB models Ren (); Fakhari2 (), especially at high density ratios. We also conducted a quantitative comparison between the present model and the previous LB models Ren (); Fakhari2 (). The relative errors of the velocity with these LB models were measured and the results are summarized in Table 1. It is found that the previous models produce large relative errors, and they all increase significantly with the density ratio. In contrast, a much smaller relative error can be derived by the the present model, which also seems to be independent of the density ratio. Based on the above discussion, we can see that the present model is more accurate than the previous Allen-Cahn based LB models Ren (); Fakhari2 ().

### iii.3 Spinodal decomposition

Spinodal decomposition Cahn () is a fundamental property of a fluid mixture with two different species. For suitable compositions and quenches, the initial homogeneous mixture is unstable in the presence of small fluctuations, and then the spinodal decomposition phenomenon will take place. This phenomenon, ubiquitous in physics and chemistry, has been studied extensively. Several researchers have also investigated the spinodal decomposition problem using the LB approaches Shan (); Swift (); Zu (); Gan (), while they mainly focus on the process of phase separation with small or moderate density ratios. In this subsection, we intend to simulate this problem with the large density ratio of 1000 by the present LB model, where the gradient terms are computed by Eqs. (27) and (28). This exercise is devoted to the demonstration of the capability of our method in studying complex high-density-ratio two-phase flows. The computational mesh used here is chosen to be . The periodic boundary conditions are applied at all boundaries. In the simulation, the initial distribution of the order parameter with small fluctuations can be given by

 ϕ(x,y)=13+rand(x,y), (34)

where is a random function with the maximum amplitude of 0.01. Then a small perturbation can be imposed on a homogeneous density field via Eq. (13), where and are set to be 1000 and 1. We only consider binary fluids with the viscosity ratio of , which approaches that of a water-air system. The remaining parameters in the simulation are fixed as , and . Figure 3 depicts the time evolution of the density distribution during the phase separating process, where the time () has been nondimensionalized by the viscous time of the liquid phase . It can be found that the early stage of phase separation induces small fluctuations of the density into large scale inhomogeneties. Then some tiny droplets with random shapes are formed in the system. The droplet sizes increase with time, and some of them also coalesce into the larger ones, which leads to the eventual separation of binary fluid components. The above phase separating processes are results of the hydrodynamics and surface tension action, which conform to the expectation.

### iii.4 Droplet impact on a thin liquid film

At last, to show the capacity of the present model, we consider a complex problem of droplet impact dynamics with large density ratio. Droplet impact on liquid surfaces Yarin () is a familiar spectacle in natural event of falling raindrop on the wet ground or puddle. Further, it plays a prominent role in many technical applications, such as ink jet printing, spray cooling and and coating. In spite of its ubiquity and extensive researches Yarin (); Coppola (); JLee (); Thoraval (), numerical simulation of such flows still poses some challenges due to complex interfacial changes in topology, and yet there exists a large density difference for a water-air system. In addition, a numerical singularity may be produced at the impact point. In this section, we will simulate a two-dimensional droplet impact on a preexisting thin liquid film with a large density ratio of 1000 by the present LB model, in the absence of the gravitational field.

The simulations are performed on a uniform computational mesh with the size of , as illustrated in Fig. 4. A wetting liquid film with the height of is initially located at the bottom wall, and a circular droplet with the radius () of lattice units is just placed on the upper of the liquid film. In the simulation, the distribution of the order parameter can be initialized by

 ϕ(x,y)=0.5+0.5tanh2(Hw−y)W, (35)

and also

 ϕ(x,y)=0.5+0.5tanh2[R−(x−0.5L)2−(y−R−Hw)2]W, y>Hw, (36)

where is the interface thickness and is set to be 5. The velocity field at initial time can be assigned by

 (u,v)={(0,−ϕU), y>Hw,(0,0),     y≤Hw, (37)

where is the impact velocity with a fixed value of . The periodic boundary conditions are applied at the left and right boundaries, while the no-slip bounce back boundary condition is imposed at the bottom wall and the open boundary condition is utilized at the top boundary. Two major dimensionless parameters governing droplet impact are the Reynolds number and the Weber number, which are respectively defined by Yarin ()

 Re=ρlDUμl, (38)

and

 We=ρlDU2σ, (39)

where is the droplet diameter. The Weber number has been taken fixed and equals to be , as commonly used in other studies  Josseranda (); Lee1 (); Coppola (). The density ratio of the liquid and gas phases is set to . Three typical Reynolds numbers , 100, and 20 are considered in this work, which are derived by tuning the kinematic viscosity of the liquid phase while keeping the gas kinematic viscosity as a constant. With this strategy, it is found that the lowest achievable liquid kinematic viscosity is 0.02 at the largest of 500, and the viscosity ratio at this situation is 10, which is very close to that of a realistic water-air two-phase system. With the driving of the impact velocity, the system is released and the droplet will instantly impact onto the underneath film. Here we mainly concentrate on the interfacial dynamics and the variation law of the spreading radius versus time. Figs. 5-7 depict typical scenic representations of droplet impact process at three different Reynolds numbers of 20, 100 and 500, where the time instants is the normalized time defined by , is the iteration step. For high of 500, the droplet moves downward instantly with slight deformation at initial stage, and some tiny bubble-ring entrapments are visible in the neck connecting the droplet and film. The bubble entrapment phenomenon has also been reported in the recent studies on the droplet impact  JLee (); Thoraval (). Then the droplet continues to spread, followed by the formation of the ejecta sheet at the intersection region between the droplet and liquid layer. The ejecta sheet grows into a splashing lamella (known as crown in axisymmetric or 3D geometry) propagating radially with increasing time and tends to bend at its end rim. The splashing phenomenon is also observed in a moderate of 100, as shown in Fig. 6, while the extent is significantly reduced. This is ascribed to the increasing frictional force between the liquid phase and its ambient vapor phase at a larger viscosity, and then the interface layer can be more stable. As the Reynolds number is lowered to a small value of 20, we do not observe the splashing behavior of the droplet in Fig. 7, as expected. The droplet only merge with the thin liquid film, which evolves in a manner of the outward moving surface wave. This process of droplet impact is oftentimes named as deposition, which is in line with the results of the previous studies Lee1 (); Li2 ().

We also conducted a quantitative study on the spreading radius, which is a concerning physical quantity in droplet impact dynamics  Yarin (). The previous researches Josseranda (); Coppola (); Ba (); Li2 (); Lee1 () have indicated that the growth of the spreading radius generally can be described by the power law , where is a coefficient that depends on the flow geometry. For the axisymmetric or 3D modeling of the droplet impact, the coefficient is found to be 1.1 by Josseranda and Zaleskib Josseranda (). While for the plane two-dimensional situation, the scaled prefactor is found to be larger than 1.1, as reported in several literatures Lee1 (); Ba (); Li1 (); Coppola (). Figure 8 shows the time variation of the numerically predicted spreading radius by the present model. For a comparison, the theoretical results of the fitting power formula is also presented. The comparison between them shows a good agreement in general, except for slight deviation at initial instants. In addition, the spreading radius in the present simulation exhibits to obey the power law .

## Iv Summary

Numerical modeling of two-phase flows with large density ratios is still a challenging task in the framework of LB approach. In this paper, we propose a simple and accurate LB model for two-phase systems, which is capable of simulating large-density-ratio flows. The proposed LB model is based on the conservative phase-field equation, which involves a lower-order diffusion term compared with the commonly used Cahn-Hilliard equation in interface capturing. Therefore, the present model is expected to achieve a better numerical accuracy and stability. In addition, a novel force distribution function is also is elaborately designed in this model such that it contains only one nonlocal macroscopic quantity, which is much simpler than the previous phase-field-based LB models Lee1 (); Lee2 (); Zu (); Liang1 (); Ren (); Fakhari2 (). The multi-scale analysis also demonstrates that both the conservative Allen-Cahn equation and the incompressible Navier-Stokes equations can be derived correctly from the present model. To validate the present model, we firstly simulated two basic steady problems of static droplet and layered Poiseuille flows, which have their own analytical solutions. In the former test, it is found that the present can accurately capture the density field distributions in the bulk regions and also across the interface at the density ratio of 1000. In addition, it is also shown that the present model can obtain relatively small spurious velocities in the LB community, with the maximum magnitude of the order of . In the latter test, we simulated the channel flow with a density ratio ranging from 10 to 1000, and also conducted detailed comparisons with the previous Allen-Cahn based LB models Ren (); Fakhari2 (). It is found that the present model can obtain satisfactory results in the velocity predictions, and is also more accurate than the previous LB models Ren (); Fakhari2 (). Next, we consider two dynamic problems of Spinodal decomposition and droplet impact on a thin liquid film with a large density ratio of 1000. The phase separation process can be clearly observed in the system, which is in line with the expectation. The present model also successfully reproduces the classical splashing phenomenon, and the predicted spreading radius is found to exhibit the power law reported in the literature, which provides a good validation of the present LB model in dealing with complex high-density-ratio two-phase flows. Finally, we anticipate that our numerical method will be useful to scientific applications, such as micro-fluidics, material science, and oil recovery industry.

## Acknowledgments

One of the authors (Hong Liang) gratefully thanks insightful discussions with Prof. Qing Li in the study of droplet impact problem. This work is financially supported by the National Natural Science Foundation of China(Grant Nos. 11602075 and 51576079).

## Appendix A Chapman-Enskog analysis of the present model

The Chapman-Enskog analysis is now performed to demonstrate the consistency of the LB evolution equation (15) with the incompressible Navier-Stokes equations. The moment conditions is first given based on the expressions of the equilibrium and force distribution functions:

 ∑igeqi=0, ∑iciαgeqi=ρuα,∑iciαciβgeqi=ρuαuβ+pδαβ, ∑iciαciβciγgeqi=ρc2sΔαβγθuθ, (40)
 ∑iGi=(1−12τg)uα∂αρ, ∑iciαGi=(1−12τg)Fα, (41) Λ=:∑iciαciβGi=(1−12τg)[uα∂β(ρc2s)+uβ∂α(ρc2s)+(uγ∂γρc2s)δαβ],

where is the Kronecker delta function, . To derive the macroscopic equations, we expand the particle distribution function, the time and space derivatives, and the force in consecutive scales of ,

 gi=g(0)i+ϵg(1)i+ϵ2g(2)i+⋅⋅⋅, (42a) ∂t=ϵ∂t1+ϵ2∂t2, ∂α=ϵ∂1α, (42b) Fα=ϵF(1)α, (42c)

where is a small expansion parameter. Applying the Taylor expansion to Eq. (15), and substituting Eq. (A3) into the expanded result, we can obtain the following multi-scale equations,

 ϵ0: g(0)i=g(eq)i, (43a) ϵ1: D1ig(0)i=−1τgδtg(1)i+G(1)i, (43b) ϵ2: ∂t2g(0)i+D1ig(1)i+δt2D21ig(0)i=−1τgδtg(2)i, (43c)

where . The substitution of Eq. (A4b) into Eq. (A4c) yields

 ∂t2g(0)i+D1i(1−12τg)g(1)i=−1τgδtg(2)i−δt2D1iG(1)i. (44)

Following Refs. Liang1 (); Liang3 (), the zero-order moment of can be defined as,

 ∑kgk=−δt2uα∂αρ. (45)

Applying the expansion formula (A3) to Eqs. (21a) and (A6), one can easily derive

 ∑ig(1)i=−δt2uα∂1αρ, ∑ig(n)i=0, (n≥2), (46)
 ∑iciαg(1)i=−δt2F(1)α, ∑iciαg(n)i=0, (n≥2). (47)

The recovered equations at scale can be obtained by summing Eq. (A4b) and Eq. (A4b) over , respectively,

 ∂1αuα=0, (48)
 ∂t1(ρuβ)+∂1α(ρuαuβ+pδαβ)=F(1)β. (49)

Similarly, the recovered equations at scale can be derived from Eq. (A5),

 ∂t1(−δt2uα∂1αρ)+∂1α(−δt2F(1)α)=−δt2[∂t1(uα∂1αρ)+∂1αF(1)α] (50)
 ∂t2(ρuβ)+(1−12τg)∂1αΠ(1)=−δt2∂1αΛ(1), (51)

where is the first-order momentum flux tensor determined below, and . From Eq. (A4b), one can get

 Π(1)=∑iciαciβg(1)i=−τgδt∑iciαciβ[D1ig(0)i−G(1)i]=−τgδtc2s[∂1α(ρuβ)+∂1β(ρuα)+(∂1γρuγ)δαβ]+τgδtΛ(1), (52)

where the terms of have been neglected under the incompressible limit. Substituting Eq. (A13) into Eq. (A12), one can simplify Eq. (A12) as

 ∂t2(ρuβ)−∂1α[νρ(∂1αuβ+∂1βuα)]=0, (53)

where is the kinematic viscosity. Combining Eqs. (A9) and (A11) at and scales, together with Eqs. (A10) and (A14), we have

 ∂αuα=0, (54)
 ∂t(ρuβ)+∂α(ρuαuβ+pδαβ)=∂1α[νρ(∂1αuβ+∂1βuα)]+Fβ, (55)

which clearly shows that the incompressible Navier-Stokes equations can be exactly recovered from the present LB model.

## References

• (1) S. Succi, The Lattice Boltzmann Equation for Fluid Dynamics and Beyond, Oxford University Press, Oxford, 2001.
• (2) Z. L. Guo, C. Shu, Lattice Boltzmann method and its applications in engineering, World Scientific Singapore, 2013.
• (3) B. C. Shi and Z. L. Guo, Lattice Boltzmann model for nonlinear convection-diffusion equations, Rev. E 79 (2009) 016701.
• (4) Z. H. Chai, B. C. Shi, Z. L. Guo, A multiple-relaxation-time lattice Boltzmann model for general nonlinear anisotropic convection-diffusion equations, J. Sci. Comput. 69 (2016) 355-390.
• (5) A. K. Gunstensen, D. H. Rothman, S. Zaleski, and G. Zanetti, Lattice Boltzmann model of immiscible fluids, Phys. Rev. A 43 (1991) 4320.
• (6) X. Shan, H. Chen, Lattice Boltzmann model for simulating flows with multiple phases and components, Phys. Rev. E 47 (1993) 1815.
• (7) M. Swift, W. Osborn, and J. Yeomans, Lattice Boltzmann simulation of nonideal fluids, Phys. Rev. Lett. 75 (1995) 830.
• (8) X. He, S. Chen, and R. Zhang, A lattice Boltzmann scheme for incom- pressible multiphase flow and its application in simulation of Rayleigh- Taylor instability, J. Comput. Phys. 152 (1999) 642-663.
• (9) T. Lee and L. Liu, Lattice Boltzmann simulations of micron-scale drop impact on dry surfaces, J. Comput. Phys. 229 (2010) 8045-8063.
• (10) H. Liang, B. C. Shi, Z. L. Guo, Z. H. Chai, Phase-field-based multiple- relaxation-time lattice Boltzmann model for incompressible multiphase flows, Phys. Rev. E 89 (2014) 053320.
• (11) H. Liang, Z. H. Chai, B. C. Shi, Z. L. Guo, and T. Zhang, Phase-field- based lattice Boltzmann model for axisymmetric multiphase flows, Phys. Rev. E 90 (2014) 063311.
• (12) H. Liu, Q. J. Kang, C. R. Leonardi et. al., Multiphase lattice Boltzmann simulations for porous media applications, Comput. Geosci. 20 (2016) 777- 805
• (13) Q. Li, H. K. Luo, Q. J. Kang, Y. L. He, Q. Chen, Q. Liu, Lattice Boltzmann methods for multiphase flow and phase-change heat transfer, Prog. Energy Combust. Sci. 52 (2016) 62-105.
• (14) T. Inamuro, T. Ogata, S. Tajima, N. Konishi, A lattice Boltzmann method for incompressible two-phase flows with large density differences, J. Comput. Phys. 198 (2004) 628.
• (15) J. A. Sethian, Level set methods and fast marching methods: evolving interfaces in computational geometry, fluid mechanics, computer vision, and materials science, Cambridge university press, 1999.
• (16) P. Yuan, L. Schaefer, Equations of state in a lattice Boltzmann model, Phys. Fluids, 18 (2006) 042101.
• (17) Q. Li, K. H. Luo, X. J. Li, Lattice Boltzmannmodeling of multiphase flows at large density ratio with an improved pseudopotential model, Phys. Rev. E 87 (2013) 053301.
• (18) Y. Ba, H. Liu, Q. Li, Q. Kang, J. Sun, Multiple-relaxation-time color-gradient lattice Boltzmann model for simulating two-phase flows with high density ratio, Phys. Rev. E 94 (2016) 023310.
• (19) D. Jacqmin, Calculation of two-phase Navier-Stokes flows using phase-field modeling, J. Comput. Phys. 155 (1999) 96-127.
• (20) H. W. Zheng, C. Shu, Y. T. Chew, A lattice Boltzmann model for multiphase flows with large density ratio, J. Comput. Phys. 218 (2006) 353.
• (21) A. Fakhari, M.H. Rahimian, Phase-field modeling by the method of lattice Boltzmann equations, Phys. Rev. E 81 (2010) 036707.
• (22) T. Lee, C L. Lin, A stable discretization of the lattice Boltzmann equation for simulation of incompressible two-phase flows at high density ratio, J. Comput. Phys. 206 (2005) 16-47.
• (23) T. Lee, L. Liu, Lattice Boltzmann simulations of micron-scale drop impact on dry surfaces, J. Comput. Phys. 229 (2010) 8045.
• (24) Q. Lou, Z. L. Guo, B. C. Shi, Effects of force discretization on mass conservation in lattice Boltzmann equation for two-phase flows, Europhys. Lett. 99 (2012) 64005.
• (25) Y. Q. Zu, S. He, Phase-field-based lattice Boltzmann model for incompressible binary fluid systems with density and viscosity contrasts, Phys. Rev. E 87 (2013) 043301.
• (26) Y. Wang, C. Shu, H. B. Huang, C. T. Teo, Multiphase lattice Boltzmann flux solver for incompressible multiphase flows with large density ratio, J. Comput. Phys. 280 (2015) 404-423.
• (27) F. Ren, B. W. Song, M. C. Sukop, and H. B. Hu, Improved lattice Boltzmann modeling of binary flow based on the conservative Allen-Cahn equation, Phys. Rev. E 94 (2016) 023311.
• (28) A. Fakhari, D. Bolster, Diffuse interface modeling of three-phase contact line dynamics on curved boundaries: A lattice Boltzmann model for large density and viscosity ratios, J. Comput. Phys. 334 (2017) 620-638.
• (29) H. L. Wang, Z. H. Chai, B. C. Shi, and H. Liang, Comparative study of the lattice Boltzmann models for Allen-Cahn and Cahn-Hilliard equations, Phys. Rev. E 94 (2016) 033304.
• (30) Y. Sun, C. Beckermann, Sharp interface tracking using the phase-field equation, J. Comput. Phys. 220 (2007) 626-653.
• (31) P. H. Chiu, Y. T. Lin, A conservative phase field method for solving incompressible two-phase flows, J. Comput. Phys. 230 (2011) 185-204.
• (32) Y. H. Qian, D. d¡¯Humires, and P. Lallemand, Lattice BGK models for Navier-Stokes equation, Europhys. Lett. 17 (1992) 479-484.
• (33) I. Ginzburg, F. Verhaeghe, D. d¡¯Humieres, Two-relaxation-time lattice Boltzmann scheme: About parametrization, velocity, pressure and mixed boundary conditions, Commun. Comput. Phys. 3 (2008) 427-478.
• (34) P. Lallemand, L. S. Luo, Theory of the lattice Boltzmann method: dispersion, dissipation, isotropy, Galilean invariance, and stability, Phys. Rev. E 61 (2000) 6546.
• (35) S. O. Unverdi, G. Tryggvason, A front-tracking method for viscous, incompressible, multi-fluid flows, J. Comput. Phys. 100 (1992) 25-37.
• (36) J. Kim, A continuous surface tension force formulation for diffuse-interface models, J. Comput. Phys. 204 (2005) 784-804.
• (37) H. Liang, B. C. Shi, and Z. H. Chai, Lattice Boltzmann modeling of three-phase incompressible flows, Phys. Rev. E 93 (2016) 013308.
• (38) Y. K. Wei, Z. D. Wang, H. S. Dou, Y. H. Qian, A novel two-dimensional coupled lattice Boltzmann model for incompressible flow in application of turbulence Rayleigh-Taylor instability, Comput. Fluids 156 (2017) 97-102.
• (39) Z. L. Guo, C. G. Zheng, and B. C. Shi, Discrete lattice effects on the forcing term in the lattice Boltzmann method, Phys. Rev. E 65 (2002) 046308.
• (40) Y. K. Wei, Z. D. Wang, J. F. Yang, H. S. Dou, Y. H Qian, A simple lattice Boltzmann model for turbulence Rayleigh-B¨¦nard thermal convection, Comput. Fluids 118 (2015) 167-171.
• (41) D. M. Anderson, G. B. McFadden, and A. A. Wheeler, Diffuse-interface methods in fluid mechanics, Annu. Rev. Fluid Mech. 30 (1998) 139-165.
• (42) Z. L. Guo, C. G. Zheng, and B. C. Shi, Force imbalance in lattice Boltzmann equation for two-phase flows, Phys. Rev. E 83 (2011) 036707.
• (43) Z. Yu, L. S. Fan, Multirelaxation-time interaction-potential-based lattice Boltzmann model for two-phase flow, Phys. Rev. E 82 (2010) 046708.
• (44) H. B. Huang, X. Y. Lu, Relative permeabilities and coupling effects in steady-state gas-liquid flow in porous media: A lattice Boltzmann study, Phys. Fluids 21 (2009) 092104.
• (45) H. Liang, B. C. Shi, Z. H. Chai, An efficient phase-field-based multiple-relaxation-time lattice Boltzmann model for three-dimensional multiphase flows, Comput. Math. Appl. 73 (2017) 1524-1538.
• (46) J. W. Cahn, Phase separation by spinodal decomposition in isotropic systems, J. Chem. Phys. 42 (1965) 93-99.
• (47) Y. B. Gan, A. G. Xu, G. C. Zhang, and S. Succi, Discrete Boltzmann modeling of multiphase flows: hydrodynamic and thermodynamic non-equilibrium effects, Soft Matter 11 (2015) 5336-5345.
• (48) A. L. Yarin, Drop impact dynamics: splashing, spreading, receding, bouncing…, Annu. Rev. Fluid Mech. 38 (2006) 159-192.
• (49) G. Coppola, G. Rocco, and L. Luca, Insights on the impact of a plane drop on a thin liquid film, Phys. Fluids 23 (2011) 022105.
• (50) J. S. Lee, B. M. Weon, J. H. Je, and K. Fezzaa, How does an air film evolve into a bubble during drop impact?, Phys. Rev. Lett. 109 (2012) 204501.
• (51) M. J. Thoraval, K. Takehara, T. G. Etoh and S. T. Thoroddsen, Drop impact entrapment of bubble rings, J. Fluid Mech. 724 (2013) 234-258.
• (52) C. Josseranda and S. Zaleskib, Droplet splashing on a thin liquid film, Phys. Fluids 15 (2003) 1650.
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters