Improved three-dimensional color-gradient lattice Boltzmann model for immiscible multiphase flows

In this paper, an improved three-dimensional color-gradient lattice Boltzmann (LB) model is proposed for simulating immiscible multiphase flows. Compared with the previous three-dimensional color-gradient LB models, which suffer from the lack of Galilean invariance and considerable numerical errors in many cases owing to the error terms in the recovered macroscopic equations, the present model eliminates the error terms and therefore improves the numerical accuracy and enhances the Galilean invariance. To validate the proposed model, numerical simulation are performed. First, the test of a moving droplet in a uniform flow field is employed to verify the Galilean invariance of the improved model. Subsequently, numerical simulations are carried out for the layered two-phase flow and three-dimensional Rayleigh-Taylor instability. It is shown that, using the improved model, the numerical accuracy can be significantly improved in comparison with the color-gradient LB model without the improvements. Finally, the capability of the improved color-gradient LB model for simulating dynamic multiphase flows at a relatively large density ratio is demonstrated via the simulation of droplet impact on a solid surface.


I. Introduction
In the past three decades, the lattice Boltzmann (LB) method [1][2][3][4][5][6][7][8][9], which originates from the lattice gas automaton (LGA) method [10], has been developed into an efficient numerical approach for simulating fluid flow and heat transfer. Different from conventional numerical methods, which are based on the direct discretization of macroscopic governing equations, the LB method is built on the mesoscopic kinetic equation. It tracks the evolution of a particle distribution function and then accumulates the particle distribution function to obtain the macroscopic properties. Owing to its kinetic nature, the LB method exhibits some advantages over conventional numerical methods. For example, in the LB equation the convective operator (the streaming process) is linear, whereas the convective terms of the Navier-Stokes equations are nonlinear [11]. Moreover, in the LB simulations the complex boundary conditions can be formulated with the elementary mechanical rules such as the bounce-back rule according to the interaction of the "LB particles" with the solid walls. Furthermore, the LB method is ideal for parallel computing because of its explicit scheme and the local interactions.
Since the emergence of the LB method, its applications to multiphase flows have always been a very important theme of this method and various multiphase LB models have been developed from different points of view [12]. Generally, most of the existing multiphase LB models can be classified into the following four categories [5][6][7], i.e., the color-gradient LB method, the pseudopotential LB method, the free-energy LB method, and the phase-field LB method. The first color-gradient LB model was proposed by Gunstensen et al. [13], which is also the earliest mulitcomponent extension of the LGA method to the LB method [14]. In the color-gradient LB method, two distribution functions are introduced to represent two different fluids and a color-gradient-based perturbation operator is employed to generate the surface tension as well as a recoloring step for separating different phases or components. The pseudopotential LB method, which is the simplest multiphase LB method, was introduced by Shan and Chen [15,16]. In this method, the fluid interactions are mimicked by an interparticle potential, through which the separation of different phases or components can be achieved naturally. The free-energy LB method was developed by Swift et al. [17,18] based on thermodynamics considerations. They proposed to modify the second-order moment of the particle equilibrium distribution function so as to include a non-ideal thermodynamic pressure tensor. The phase-field LB method is based on the phase-field theory, in which the interface dynamics is described by an order parameter that obeys the Cahn-Hilliard equation or a Cahn-Hilliard-like equation [19].
Each of these multiphase LB methods has its advantages and limitations. A comprehensive review of the pseudopotential LB method and the phase-field LB method can be found in Ref. [7]. In addition, the book by Huang, Sukop and Lu [12] is also dedicated to the multiphase LB methods. In this work, we restrict our study to the color-gradient multiphase LB method, which exhibits very low dissolution for tiny droplets or bubbles [20] in comparison with other multiphase LB methods. In the original color-gradient LB model devised by Gunstensen et al. [13], the work done by the color gradient against the color flux was maximized to force the colored particles to move towards fluids with the same color.
In addition, the model of Gunstensen et al. suffers from the limitation of equal densities for two-phase flows. Some improvements have been conducted to overcome the shortcomings of the original color-gradient model. Grunau et al. [21] modified the form of the particle equilibrium distribution function to allow for variable density and viscosity ratios. Latva-Kokko and Rothman [22] replaced the numerical maximization recoloring step of Gunstensen et al.'s model with a formulaic segregation algorithm, which solves the lattice pinning problem at the interface region and significantly improves the computational efficiency of the color-gradient LB method.
Later, Reis and Phillips [23] proposed a new perturbation operator for generating the surface tension of the color-gradient LB method and derived a theoretical expression for the surface tension through its mechanical definition. Liu et al. [24] extended the model of Reis and Philips to three-dimensional space by deriving a generalized perturbation operator, in which an expression for the surface tension parameter is directly obtained without approximations. However, similar to the free-energy multiphase LB method, the color-gradient multiphase LB method also modifies the equilibrium distribution function [23,24] to incorporate the pressure of fluid. Hence it also suffers from the lack of Galilean invariance [7]. Through the Chapman-Enskog analysis, Huang et al. [25] showed that some error terms exist in the macroscopic momentum equation recovered from the color-gradient multiphase LB method. They demonstrated that for two-phase flows with different densities the error terms significantly affect the numerical accuracy. A scheme has been proposed by Huang et al. [25] to eliminate the error terms, but they emphasized that their scheme just works well for cases of density ratios less than 10.

II. The existing 3D color-gradient LB models
A. The color-gradient LB equation In the color-gradient LB method, the two immiscible fluids are represented by a red fluid and a blue fluid, respectively. The corresponding distribution functions are denoted by k i f , where i is the lattice velocity direction and k R  or B denotes the color ("Red" or "Blue"). The total distribution function is defined as The evolution of the distribution functions is governed by the following LB where C is the color gradient and the free parameter A is proportional to the surface tension.
Another three-dimensional color-gradient LB model can be found in the study of Liu et al. [24], who extended the perturbation operator proposed by Reis and Phillips [23] to three-dimensional space where Moreover, Liu et al. [24] employed the recoloring algorithm proposed by Latva-Kokko and Rothman [22], which can solve the lattice pinning problem and reduce the spurious velocities. According to the recoloring algorithm of Latva-Kokko and Rothman, the recoloring steps for the red and blue fluids can be defined as follows [22,29]: where  is a free parameter controlling the interface thickness, i  is the angle between the color gradient N   and the lattice direction i e , which yields , and * i f is the post-perturbation value of the total distribution function, namely * , , ( With Eqs. (10) and (11), the "streaming" process is implemented as x . In the study of Liu et al. [24], the equilibrium distribution function is also defined by Eq.

C. The error terms
The error terms in the momentum equation recovered from the three-dimensional color-gradient LB models have been identified by Huang et al. [25] through the Chapman-Enskog analysis and are given by where the subscripts  ,  , and  denote the x, y, or z coordinate and   is the Kronecker delta.
For two-phase flows with identical densities, k i i    is usually adopted and then , k eq i f given by Eq.
(6) reduces to the standard equilibrium distribution function in the LB method, which leads to . Accordingly, the error terms disappear for two-phase flows with identical densities.
However, for two-phase flows with different densities, the error terms in Eq. (14) will make the Galilean invariance lost and may affect the numerical accuracy significantly since the density gradient cannot be neglected near the interface.
Recently, Saito et al. [30] constructed a three-dimensional 27-velocity (D3Q27) color-gradient LB model, in which an enhanced equilibrium distribution function devised by Leclaire et al. [31] is adopted where k i  is an additional term, which was originally employed by Che Sidik and Takahiko [32] for a free-energy LB model and was extended to the color-gradient LB method by Leclaire et al. [31].
Nevertheless, it is noticed that both Che Sidik and Takahiko [32] and Leclaire et al. [31] showed that there are still some error terms in the recovered macroscopic momentum equation, which can be found in Eqs. (29)-(33) of Ref. [31]. The main error terms are similar to the aforementioned error terms given by Eq. (14). Ⅲ . Improved 3D color-gradient LB model

A. Theoretical analysis
In this section, the physical origin of the error terms in Eq. (14) is analyzed. Taking the second-order and third-order moments of the equilibrium distribution function given by Eq. (6), we can find that As seen in Eq. (16) However, the symmetry of the standard lattices (such as the D2Q9, D3Q19, and D3Q27 lattices) is insufficient to completely support the replacement of 2 Fortunately, the off-diagonal elements of the third-order moment of the equilibrium distribution function can satisfy the required relationship by adding a high-order term to the equilibrium distribution function, as shown in Ref. [27] for recovering p RT   in a compressible LB model on standard lattices. Following the study of Li et al. [27], the new equilibrium distribution function is defined as which yields For the off-diagonal elements of the third-order moment of To remove the error terms caused by the diagonal elements of the third-order moment of , k eq i f , a correction term can be added to the single-phase collision operator where k i G is the correction term and the coefficient   eliminating the discrete effect of a forcing or source term in the LB equation [33]. The zeroth-and first-order moments of the correction term satisfy the following relationships: The constraints on the second-order moment of the correction term can be derived through the Chapman-Enskog analysis, which can be implemented by introducing the following multi-scale expansions [34]: where  is the expansion parameter.
According to the studies of Reis and Phillips [23] and Liu et al. [24], the recoloring step is not considered in the Chapman-Enskog analysis, and therefore the recoloring operator regarded as a unit operator. Meanwhile, the perturbation operator (2) ( ) k i  only affects the surface tension term and has been well demonstrated in Refs. [23,24]. Therefore, in the present study, the Chapman-Enskog analysis is performed for Eq.
being treated as (1) ( ) Taking the Taylor-series expansion of the left-hand side of Eq. (1) and using the multi-scale expansions given by Eqs. (22) and (23), we can rewrite the LB equation in the consecutive orders of  as follows: With the help of Eq. (25), Eq. (26) can be rewritten as Taking the summations of Eqs. (25) and (27) The continuity equation can be obtained by combining Eq. (28) with Eq. (29). Similarly, the first-order moments of Eqs. (25) and (27) yield, respectively Using the new equilibrium distribution function given by Eq. (18), the off-diagonal elements of the third-order moment can satisfy the required relationship, as shown in Eq. (19) ). Hence we can obtain the following constraints on the correction term k i G : However, the diagonal elements of the third-order moment With these constraints, the error terms caused by the diagonal elements of the third-order moment of the equilibrium distribution function can be removed, and then the following equation can be derived from Eq. (33) by substituting Eqs. (16) and (19) as well as the above constraints into Eq. (33): where the kinematic viscosity k  is given by To sum up, the new equilibrium distribution function given by Eq.

B. Improved model based on the MRT collision operator
Using the MRT collision operator [35][36][37][38], the single-phase collision operator with the correction term can be written as follows: where ij  is the Kronecker delta and   can be implemented in the moment space: where I is the unit matrix, The recoloring steps for the red and blue fluids are still given by Eqs. (10) and (11), respectively.
In the present work, the improved color-gradient MRT-LB model is constructed using the D3Q19 lattice and the following transformation matrix M is employed [39,40] where the relaxation times v  and e  determine the shear and bulk viscosities, respectively, while q  and   are related to non-hydrodynamic moments. The equilibria in the moment space can be obtained where x Q , y Q , and z Q are given by (see also Eqs. (35)-(37)), respectively In numerical implementation, the second-order isotropic difference scheme is applied to the spatial gradients in Eqs. (46)-(48), i.e., where  denotes the x, y, or z coordinate and   The Chapman-Enskog analysis can also be applied to the MRT collision operator, which is similar to that of the BGK collision operator. Readers are referred to Refs. [40][41][42] about the Chapman-Enskog analysis of the three-dimensional MRT-LB method. It can be found that, using the equilibria , k eq m in the Appendix and the correction term k C given by Eq. (45), the following macroscopic momentum equation can be derived in the low Mach number limit: where the dynamic shear viscosity k  and the bulk viscosity b k  are given by, respectively The kinematic viscosity k  is given by k where  is a free parameter related to the interface thickness and is usually set as 0.98   [12], and The surface tension in Eq. (9) depends on the relaxation time. A simple treatment to make the surface tension independent of the relaxation time is to change the perturbation operator from Eq. (8) to (2), new , and then the surface tension is given by the perturbation operator within the framework the MRT-LB method can be redefined as (2), new Similar to the single-phase collision operator in Eq. (40), the perturbation operator given by Eq. (53) can also be executed in the moment space.

IV. Numerical results and discussion
In this section, numerical simulations are carried out to validate the improved three-dimensional color-gradient LB model. First, the test of a moving droplet in a uniform flow field is employed to verify the Galilean invariance of the improved model. Subsequently, the numerical accuracy of the improved model is demonstrated through simulating the layered two-phase flow and three-dimensional Rayleigh-Taylor instability. Finally, the capability of the improved model for simulating dynamic multiphase flows at a large density ratio is validated by the simulation of droplet impact on a solid wall.

A. Moving droplet in a uniform flow field
In LB community [43,44], it has been reported that a circular droplet in a uniform flow field will become an elliptic one when employing a multiphase LB model with broken Galilean invariance. To  being set to 0.9 and 0.7, respectively, which satisfy  Figure 1 shows the simulated snapshots of a moving circular droplet. For comparison, the numerical results of the color-gradient MRT-LB model without the improvements are also presented, which is the MRT version of the three-dimensional color-gradient LB model of Liu et al. [24] and is hereinafter referred to as the original model. When the original model is employed, the shape of the droplet becomes elliptic, as shown in Fig. 1(a), which means that the lack of Galilean invariance leads to deformation of the droplet. On the contrary, from Fig. 1(b) we can see that the improved color-gradient model allows the droplet to retain its circular shape, demonstrating that the Galilean invariance is restored in the improved model.

B. Layered two-phase flow in a channel
In this subsection, the layered two-phase flow between two parallel plates is simulated to validate the numerical accuracy of the improved color-gradient LB model. As shown in Fig. 2, the channel height in the y direction with 0 y  at the center of the channel. The red fluid is initially located in the central region a y a    , whereas the blue fluid is located in the regions a y b   . The layered two-phase flow is driven by a constant body force   , 0, 0 G . By assuming a Poiseuille-type flow in the channel, we can obtain the following the analytical solution for the velocity profile [12]: where the coefficients are defined as   is the dynamic viscosity ratio [12].

FIG. 2.
Schematic of the layered two-phase flow between two parallel plates.
In our simulations, the computational domain is divided into 10 100 4 and 50 b  . The non-slip boundary condition [45] is applied to the two parallel plates, while the periodic boundary condition is employed in the x and z directions. Three cases are investigated:

C. Three-dimensional Rayleigh-Taylor instability
The phenomenon of Rayleigh-Taylor instability is associated with the penetration of a heavy fluid into a light fluid and can be found in a wide range of scientific and environmental fields. This problem involves complex interfacial interactions and has been intensively studied because of its practical and scientific importance [46][47][48]. Rayleigh-Taylor instability develops from the following single mode initial perturbation: where h is the height of the fluid interface. can be seen that the shapes of the spike and girdle predicted by the original model obviously deviate from those obtained by a multiphase flux solver [48], while the shapes simulated by the improved model are in good agreement with those reported in Ref. [48].

D. Droplet impact on a solid surface
Finally, the impingement of a droplet on a flat surface is simulated to validate the capability of the proposed color-gradient MRT-LB model for simulating dynamic multiphase flows at a relatively large density ratio. Impingement of droplets on a solid surface is a very important phenomenon in many engineering applications, such as ink-jet printing and spray cooling. The dynamics of droplet impact on solid surfaces is usually governed by the following two non-dimensional parameters: where Re and We are the Reynolds number and the Weber number, respectively. In Eq. (57), 0 U is the impact speed of the droplet and 0 D is the initial diameter of the droplet.
The computational domain is divided into 300 300 150  In our simulations, the Reynolds number varies from Re 75  to 1000 . Figure 8 displays some snapshots of the droplet impingement process at Re 1000  . As shown in the figure, immediately after the impingement, the shape of the droplet resembles a truncated sphere (Fig. 8a). Later, a lamella is formed as the liquid moves radially outwards (Figs. 8b and 8c). The lamella continues to grow radially ( Fig. 8d) until the maximum spreading diameter is reached and the spreading process ends, during which the kinetic energy is transformed into the surface energy by increasing the area of the droplet [48] Figure 9 shows a comparison of the maximum spreading factor between the experimental correlation formula of Asai et al. [50], the experimental data of Scheller and Bousfield [51], and the numerical results predicted by the proposed color-gradient LB model. From the figure it can be seen that our numerical results are in good agreement with the experimental correlation/data reported in the previous studies, demonstrating that the improved color-gradient LB model is capable of simulating dynamic multiphase flows at a relatively large density ratio. Comparison of the maximum spreading factor between the present numerical results, the experimental correlation in Ref. [50], and the experimental data in Ref. [51].

V. Conclusions
The previous three-dimensional color-gradient LB models usually suffer from the lack of Galilean invariance and considerable numerical errors because of the error terms in the recovered macroscopic equations. In this paper we have theoretically analyzed the physical origin of the error terms in the previous models. Based on the theoretical analysis, we have proposed an improved three-dimensional color-gradient LB model for simulating immiscible multiphase flows. Specifically, a high-order term is added to the equilibrium distribution function, through which the off-diagonal elements of the third-order moment of the equilibrium distribution function can satisfy the required relationship for recovering the correct Navier-Stokes equations. Meanwhile, the deviations of the diagonal elements are corrected via a correction term in the LB equation. Compared with the previous models, the present model eliminates the error terms and therefore improves the numerical accuracy and enhances the Galilean invariance.
To validate the proposed color-gradient LB model, numerical simulation have been performed. The test of a moving droplet in a uniform flow field has been employed to verify the Galilean invariance of the improved model. It has been shown that the shape of the droplet becomes elliptic when the original model is used, while the improved model allows the droplet to retain its circular shape. Numerical simulations have also been carried out for the layered two-phase flow and three-dimensional  (1 ) 2 ( ) 6 (1 ) 2 ( ) 6