On the lattice Boltzmann method for multiphase flows with large density ratios

Accepted Manuscript On the lattice Boltzmann method for multiphase flows with large density ratios

Seung Hyun Kim, Heinz Pitsch

PII: DOI: Reference:

S0021-9991(15)00622-1 http://dx.doi.org/10.1016/j.jcp.2015.09.029 YJCPH 6132

To appear in:

Journal of Computational Physics

Received date: Revised date: Accepted date:

29 January 2015 9 July 2015 15 September 2015

Please cite this article in press as: S.H. Kim, H. Pitsch, On the lattice Boltzmann method for multiphase flows with large density ratios, J. Comput. Phys. (2015), http://dx.doi.org/10.1016/j.jcp.2015.09.029

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

On the lattice Boltzmann method for multiphase flows with large density ratios Seung Hyun Kima , Heinz Pitschb a

Department of Mechanical and Aerospace Engineering, The Ohio State University, Columbus, OH 43210, USA b Institut f¨ ur Technische Verbrennung, RWTH Aachen University, Aachen, Germany

Abstract An analysis of the lattice Boltzmann (LB) method for multiphase flows with large density ratios is presented. It is shown that for incompressible, multiphase LB methods, the divergence-free condition is not satisfied within the formal accuracy of the LB method, when the density ratio between the two phases is large enough. The discrete differentiation-by-parts rule is responsible for this error. A new multiphase LB method to resolve this issue is proposed. Keywords: Lattice Boltzmann method, Kinetic theory, Multiphase flow

1. Introduction The lattice Boltzmann (LB) method [1, 2, 3, 4, 5] is a reduced-order kinetic model to reproduce the Navier-Stokes hydrodynamics (and beyond) at a macroscopic level. Since the pioneering work of McNamura and Zanetti [1], the LB method has been extended to various complex flows involving, for example, multicomponent and interfacial phenomena, and has been particularly successful for multiphase flows in complex geometry [3, 5]. While there are exceptions [6], most multiphase LB methods that have been developed so far [3, 7, 8, 9] belong to a class of the diffuse interface model [10], where an interface is numerically resolved with a few grid points across it. A specific multiphase formulation of the discrete kinetic model can be derived from a continuum kinetic equation such as the Enskog equation or it can be constructed a posteriori in such a way that the corresponding macroscopic equation follows a desired hydrodynamic equation for multiphase flows. The former approach includes the model Email address: [email protected] (Seung Hyun Kim) Preprint submitted to J. Comput. Phys.

September 21, 2015

of Luo and Girimaji [11] and of He et al. [9], while the latter includes the Shan-Chen model [7] and the free energy model [8]. Typically, the density ratio in simulations using these methods is restricted to O(10). Several approaches have been developed to allow for application to flows with large density ratios [12, 13, 14]. He et al. [9] presented an incompressible multiphase LB method where the evolution of the phase interface is separated from that of mass and momentum. Lee and Lin [13] showed that with a stable discretization scheme for the forcing terms, the mean-field LB method of He et al. [9] can be applied to flows with density ratios of up to 1000, while numerical stability also depends on the shear rate [15]. Inamuro et al. [12] employed the projection method, for which the pressure Poisson equation is solved at each time step to enforce the divergence-free condition for the velocity field, and showed that their method can consider density ratios up to 1000. Wang et al. [16] presented a multiphase LB flux solver for incompressible multiphase flows with large density ratios. Here, an analysis of the multiphase LB method for large density ratios is presented. The focus is on the satisfaction of the divergence-free velocity field in incompressible multiphase flows. An origin of the incompressibility error in the multiphase LB method is identified and discussed, and a new LB formulation is presented to overcome this issue. The term, “multiphase flows,” here refers to two-phase interfacial flows.

2. Formulation and analysis 2.1. Macroscopic equations The macroscopic governing equations recovered by the incompressible and multiphase LB method [9, 13] can be written as, ∂uβ 1 ∂p +ρ = 0, c2s ∂t ∂xβ      ∂uβ ∂uα ∂uα ∂p ∂ ∂uα + uβ =− ρν + ρFαs , + + ρ ∂t ∂xβ ∂xα ∂xβ ∂xα ∂xβ   ∂C ∂C ∂φ ∂ + uβ M , = ∂t ∂xβ ∂xβ ∂xβ

(1) (2) (3)

where uα is the velocity, p is the pressure, ρ is the density, and ν is the kinematic viscosity. In the phase-field method, the dynamics of interfaces between different phases are described using the 2

order parameter C, which is an indicator of different phases. In the equation for the order parameter, M and φ are the mobility and the chemical potential, respectively. The surface tension force Fαs can be implemented in different ways, and the difference between the different implementations of the surface tension can be absorbed in the “pseudo pressure,” p, in the momentum equation. Here, the stress form of the surface tension is used, which is given by   ∂2φ ∂φ 1 ∂ Fαs = − δαβ . ρ ∂xβ ∂xα ∂xβ ∂xγ ∂xγ

(4)

2.2. Lattice Boltzmann method: Momentum-based formulation The LB method that reproduces the above macroscopic equations has been first derived in He et al. [9] by adopting a transformation of the distribution function that separates the density evolution from the pressure and the momentum evolution in the LB equations. The discrete velocity equations (DVEs) for the multiphase LB method can be written as ciα − uβ ∂ρc2s 1 ∂fi ∂fi + ciα = − (fi − fie ) + Γiαβ + Ψiα ρFαs , ∂t ∂xα τ c2s ∂xα

(5)

∂gi ∂gi 1 + ciα = − (gi − gie ), ∂t ∂xα τg

(6)

where

 Ψiα = wi Γiαβ =

 (ciα − uα ) uγ ciα ciγ , − c2s c4s

ciα uα (ciα ciβ − c2s δαβ )uα uβ + . c2s 2c4s

(7) (8)

ciα is the discrete velocity for the distribution functions fi and gi . The DVE for the distribution function fi describes the mass and momentum conservation, while that for gi describes the evolution of the order parameter. The DVE for gi is different from that in He et al. [9] and Lee and Lin [13]. τ and τg are the relaxation time for fi and for gi , respectively. The equilibrium functions are given by

  p ciα uα (ciα ciβ − c2s δαβ )uα uβ , + ρ + c2s c2s 2c4s   Cuα ciα + G (mφδ + Cu u ) , gie = wi0 C + wi α β iαβ αβ c2s fie = wi



where Giαβ = wi

ciα ciβ − c2s δαβ , 2c4s

w10 = 1, wi0 = 0(i = 1). 3

(9) (10)

(11) (12)

wi is the weight for the i-th distribution function. The discrete velocities and weights are from the standard scheme, D2Q9 for 2-D problems. The lower moments of fi and gi are 

fi =

i

 i

p , c2s

ciα fi = ρuα ,



gi = C.

(13) (14) (15)

i

The LB equation can be obtained by the spatial and time integration of the DVB equation. The second-order integration along a direction of a discrete velocity gives fˆi (x + ci δt, t + δt) − fˆi (x, t) = −

1 τ (fˆi − fie ) + δtFi , τ /δt + 0.5 τ + 0.5δt

(16)

where δt fˆi = fi − (Ωi + Fi ), 2 ciα − uβ ∂ρc2s Fi = Γiαβ + Ψiα ρFαs . c2s ∂xα

(17) (18)

Ωi is the collision operator for the i-th distribution function. The equation for the order parameter distribution function can be obtained in the same way. 2.3. Incompressibility error To illustrate the error in the continuity equation for the above LB formulation, a special case of τ /δt = 0.5 is considered. When τ /δt = 0.5, the equation for the distribution function fˆi can be written as 1 fˆi (x, t) = fie (x − ci δt, t − δt) + δtFi (x − ci δt, t − δt), 2

(19)

which can be transformed to fi (x, t) − fie (x, t − δt) = fie (x − ci δt, t − δt) − fie (x, t − δt) − [fi (x, t) − fie (x, t)] 1 + δtFi (x, t − δt) + δt [Fi (x − ci δt, t − δt) − Fi (x, t − δt)] . 2

(20)

Summing over the discrete directions leads to p(x, t) − p(x, t − δt) =



[fie (x − ci δt, t − δt) − fie (x, t − δt)]

i

+ δtuα

∂ρ 1  (x, t − δt) + δt [Fi (x − ci δt, t − δt) − Fi (x, t − δt)] . ∂xα 2 i

4

(21)

The first and the third term on the r. h. s. correspond to the approximation of the spatial derivatives. The explicit expressions for these terms can be obtained by using the Taylor expansion of a scalar φ: φ(x + δtci ) = φ(x) + δtciα

∂φ ∂2φ 1 + δt2 ciα ciβ + O(δt3 ). ∂xα 2 ∂xα xβ

(22)

From the isotropy of the quadrature, the first-order derivatives and the Laplacian can be evaluated by

δφ 1  1  ∂φ ≈ = 2 wi ciα [φ(x + ci δt) − φ(x)] = 2 wi ciα φ(x + ci δt), ∂xα δxα cs δt cs δt

(23)

∂2φ δ2φ 2  ≈ = 2 2 wi [φ(x + ci δt) − φ(x)] . 2 2 cs δt ∂xβ δxβ i

(24)

i

i

Note that the discretization scheme in [17], which is shown to reduce the spurious currents, is identical to the above. Using the above equations, Eq. (21) can be written as δρuβ ∂ρ 1 p(x, t) − p(x, t − δt) = − uβ + O(δt2 ). 2 cs δt δxβ ∂xβ

(25)

The last term on the r. h. s. of Eq. (20) is of O(δt3 ), when the use is made of the discretization for the Laplacian, Eq. (24). The terms on the r.h.s. are evaluated at x and t − δt. Equation (25) is the leading-order, discrete continuity equation reproduced from the LB equation with τ /δt = 0.5. To impose the divergence-free velocity field within the error of O(Ma2 ), the following relation is to be satisfied within the formal accuracy of the LB method: ρ

δρuβ ∂uβ ∂ρ ≈ − uβ . ∂xβ δxβ ∂xβ

(26)

Note that the first term on the r. h. s. is originated from the propagation step of the LB method. This term, which is the approximation of the divergence of the momentum vector, is fixed by the LB propagation scheme, and given by Eq. (23). When applying the same discretization scheme as Eq. (23) to the density derivative term, the following differentiation-by-parts has to be satisfied: δρuβ δuβ δρ ≈ − uβ . δxβ δxβ δxβ

(27)

δuβ δρ 1 δρuβ 1 ≈ − uβ . δxβ ρ δxβ ρ δxβ

(28)

ρ Dividing by the density gives

In Eq. (28), the l. h. s. should be O(Ma2 ) in the LB method. Near the phase interface, each term on the r.h.s. is usually much larger than the l.h.s., when the density ratio is large, especially on 5

the lower-density side of the interface. For incompressibility, these two large terms have to cancel within the error of O(Ma2 ). This condition is not generally satisfied. Indeed, since the evolution of the pressure field to satisfy the incompressibility condition is driven by the divergence multiplied by the density, the differentiation-by-part error is amplified in the momentum-based formulation. As a result, the truncation error is proportional to the density ratio near the interface and can be much larger than what is acceptable for the pseudo incompressibility. Even when other discretization schemes are applied to the density gradient in Eq. (26), the analysis still remains valid. While the above analysis is restricted to a special case of τ /δt = 0.5, a general case can be considered using the Chapman-Enskog expansion or the method of the successive approximation [8]. Using the method of the successive approximation [8], we obtain  (0)  (0) (∂t + ciα ∂α ) fi − δtτ ∂t2 + 2∂t ciα ∂α + ciα ciβ ∂αβ fi + δtτ (∂t Fi + ∂α ciα Fi ) + O(δt2 ) =−

(0)

where fi

is the equilibrium distribution function for the discrete velocity ciα . By summing

Eq. (29) over the discrete velocities, we obtain ∂t pc−2 s

(29)

1 τ (0) (fi − fi ) + Fi , δt(τ + 1/2) τ + 1/2

+ ∂α ρuα − δtτ

∂t2 pc−2 s

+ 2∂t ∂α ρuα + ∂αβ



(0) fi ciα ciβ

+ δtτ ∂t

i

 τ 1 uα ∂α ρ + Fi = uα ∂α ρ. = 2(τ + 1/2) τ + 1/2



Fi + ∂α ciα Fi + O(δt2 )

i

i

(30) The first-order continuity equation can then be written as ∂t pc−2 s + ∂α ρuα − uα ∂α ρ = O(δt).

(31)

The enforcement of the divergence-free condition is made through the use of differentiation-byparts. The present analysis is related to the specific LB formulation that employes the transformation proposed by He et al. [9]. However, the discrete form of differentiation-by-parts results in similar errors for most multiphase LB models. Also, such a rule is used when deriving the viscous term, which can result in the error in momentum conservation. The origin of this error is that the leadingorder contribution in the macroscopic hydrodynamic equations is obtained by the subtraction of two large terms that are proportional to the density ratio. The error is significant when the density ratio between two phases is of O(100) or larger. 6

2.4. Lattice Boltzmann method: Velocity-based formulation As a remedy to the aforementioned problem, we propose to adopt the equilibrium function that is based on the velocity uα , not the momentum ρuα . The DVB equation can be written as ∂fi ∂fi 1 + ciα = − (fi − fie ) + Ψiα (Fαs + Fαp + Fαv ) + wi0 F d , ∂t ∂xα τ where fie

=



p wi0 2 cs

+ wi

ciα uα (ciα ciβ − c2s δαβ )uα uβ + c2s 2c4s Fαp = −

Fαv =



ν ∂ρ ρ ∂xβ

1 ∂p , ρ ∂xα

∂uβ ∂uα + ∂xα ∂xβ

F d = −(cd ρ − 1)

(32)

 ,

(33) (34)

 ,

∂uγ . ∂xγ

(35) (36)

The second-order integration along a direction of a discrete velocity gives 1 (fˆi − fie ) τ /δt + 0.5  δt Ψiα (Fαs + Fαp + Fαv ) + wi0 F d .

fˆi (x + ci δt, t + δt) − fˆi (x, t) = − +

τ τ + 0.5δt

(37)

The velocity and pressure are obtained by uα =



ciα fˆi +

i

p = c2s

δt 2

 i



 Fαs + Fαp + Fαu , ρ

δt fˆi + F d . 2

(38)

(39)

In order to solve the coupled equations for the pressure and velocity in Eqs. (38) and (39), the predictor-corrector method is employed. At the predictor step, p∗ = c2s

 i

u∗α =



(n)

ciα fˆi

i

+

 δt s(n) F + F p(∗) + F v(n−1) . 2

At the corrector step,



δt (n) fˆi + F d(∗) , 2 i 

 δt (n) F s(n) + F p(n) + F v(∗) . = ciα fˆi + 2 p(n) = c2s

u(n) α

δt (n) fˆi + F d(n−1) , 2

i

7

(40)

(41)

(42) (43)

Before evaluating the moments, a filtering operation is applied to the pressure field to enhance the numerical stability. The second-order filtering operation used here is ˜ φ(x) =



wi φ(x + ci δt).

(44)

i

Higher-order filters [18] are given by ˜ φ(x) = φ(x) − a

D M  

dm φ(x + ξ j ),

(45)

j=1 m=−M

where ξj is the unit vector in the j-direction, and D is the dimensionality of the problem. The coefficients for the fifth-order filter are d0 = 6/16, d1 = −4/16, and d2 = 1/16, while those for the seventh-order filter are d0 = 5/16, d1 = −15/64, d2 = 2/32, and d3 = −1/64. Because the filtering operation is applied only to the pressure field, it smoothes out high-frequency oscillations only for the pressure field and does not directly affect the velocity field and viscous dissipation. In the velocity-based formulation, the divergence term on the r.h.s. of the pressure evolution equation is evaluated directly, not by the subtraction of the two large terms that are proportional to the density ratio as in the momentum-based formulation. As a result, the enforcement of the divergence-free condition is not subject to the differentiation-by-part error. The present formulation resembles that of Inamuro [12]. The primary difference is that in the present method, the pressure evolution is determined through the LB equation, while the pressure Poisson equation is solved in Inamuro’s method. To stabilize computation, a filtering operation is applied to the pressure field only.

3. Results and discussion As a test case to demonstrate the incompressibility error in the LB formulation for which the equilibrium function is based on the momentum ρuα , rather than the velocity uα , a single bubble rising in a viscous fluid is considered. A stationary-state bubble rise case in Unverdi and Tryggvason [19], where a reference solution obtained by the front tracking method is available, is simulated. The density ratio ρh /ρl and the viscosity ratio μh /μl are 40 and 493, respectively. The subscripts, h and l, represent a high density fluid and a low density fluid, respectively. The E¨ otv¨os number and the Morton number are 10 and 0.1, respectively. The simulation is performed with 150 × 75 lattices. The relaxation time τ is given as τ = [ρh τh C + ρl τl (1 − C)]/ρ. To match 8

(a)

(b)

(c)

(d)

Figure 1: Axial velocity and divergence fields for a bubble rising in a viscous fluid. (a) axial velocity predicted by vLBM; (b) divergence predicted by vLBM; (c) axial velocity predicted by mLBM; (d) divergence predicted by mLBM. The thick solid line represents a bubble shape predicted by the LB method, while the circles represent a bubble shape predicted by the front tracking method [19].

9

No filter

F1

F2(a = 0.25)

F2(a = 0.5)

F3 (a = 0.25)

F3 (a = 0.5)

Case 1

Unstable

1.8539×10−2

1.8550×10−2

1.8561×10−2

1.8561×10−2

1.8565×10−2

Case 2

Unstable

1.1188×10−2

1.1195×10−2

1.1193×10−2

1.1197×10−2

1.1195×10−2

Table 1: Rising velocity of a single bubble (F1: Eq. (44), F2: fifth-order filter, F3: seventh-order filter).

the non-dimensional numbers, the relaxation time for a high density fluid is set as τh /δt = 0.26. The gravitational force is added as an external forcing term. The relaxation time for the order parameter, τg /δt, and the mobility coefficient, m, are set to be 0.2. Below, all quantities are presented in lattice units. Figure 1 shows the velocity and divergence fields predicted by the velocity-based and the momentum-based formulations. Also compared with the reference solution is the predicted bubble shape. The bubble shape predicted by the velocity-based formulation, vLBM, is in very good agreement with the reference solution by the front tracking method [19]. The momentum-based formulation, mLBM, leads to a larger bubble due to the incompressibility error. Note that for mLBM, the error in the divergence field is much larger than the truncation error in the asymptotic expansion of the LB equation, O(Ma2 ), in Fig. 1. The divergence error in mLBM increases as the density ratio increases. The divergence error for the velocity-based formulation, vLBM, is within the truncation error of the LB method regardless of the density ratio. For incompressible flow simulations, the LB method can be regarded as an artificial compressibility method, where √ the divergence field depends on the sound speed. In vLBM, the sound speed is cd ρcs , and the divergence error depends on cd . With increasing cd , the sound speed increases and the error in the divergence field decreases. However, a large value of cd can make computation unstable. The value of cd used for the results in Fig. 1 is 6. Figure 2 shows the divergence of the velocity along the vertical line that passes through the center of the bubble. The underlying velocity field is predicted by the momentum based formulation. The error in the divergence-free condition is most significant near the interface where the density changes sharply. It can also be noticed that the difference between the divergence directly evaluated from the velocity, the l. h. s. of Eq. (28), and that evaluated using the differentiationby-parts, the r. h. s. of Eq. (28), is significant near the interface, for this moderately large density ratio. 10

0.005

velocity divergence

0.004 0.003 0.002 0.001 0

-0.001 -0.002

-20

-10

0

10

20

x

Figure 2: Divergence of the velocity field along the vertical line that passes through the center of a bubble, predicted by mLBM (solid line: δuβ /δxβ , dashed line: (δρuβ /δxβ − uβ δρ/δxβ ) /ρ).

11

The velocity-based formulation is not numerically stable without filtering, and the effects of filtering are investigated for rising bubble cases. Table 1 shows the rising velocity of a single bubble in a viscous fluid, predicted using different filters. Two cases are considered. The simulation parameters for Case 1 are those used for Fig. 1. For Case 2, a higher density ratio of 1000 is considered. The viscosity ratio μh /μl is 100. The E¨otv¨os number and the Morton number are 16 and 7.7×10−4 , respectively. The computational domain is discretized into 256×128 lattices. Without filtering, the computation is unstable for both cases. The rising velocity of the bubble predicted using different filters is almost identical. As in Fig. 3, the velocity field and the bubble shape are hardly affected by the use of different filters. There is, however, a noticeable effect on the divergence field in Fig. 3. When the second-order filter in Eq. (44) is used, the incompressibility error is more pronounced than using the fifth- and seventh-order filters. However, the difference between the fifth- and seventh-order filters is minimal. For the results in Fig. 1, the fifth-order filter with a = 0.25 is used. In the present method, the velocity magnitude can reach a typical upper bound for the LB method, O(0.1) in the lattice unit, while the stability domain depends on non-dimensional parameters. 4. Conclusions An analysis of the incompressible multiphase LB method is presented. The analysis shows that the incompressible multiphase LB method can generate substantial divergence in the velocity field when density ratios between two phases are large. This error in the incompressibility condition arises from the fact that the derivative operator used in the LB method does not satisfy the differentiation-by-parts condition. To remedy this problem, a new formulation where the equilibrium function is expanded using the velocity rather than the momentum is proposed. The continuity error in the new formulation is shown to be within the truncation error of the LB method. Acknowledgement Initial work at Stanford University was financially supported by Honda R&D Co., Ltd. Fundamental Technology Research Center, Wako, Japan. [1] G. R. McNamara, G. Zanetti, Use of the Boltzmann equation to simulate lattice-gas automata, Phys. Rev. Lett. 61 (1988) 2332.

12

(a)

(b)

(c)

(d)

(e)

(f)

Figure 3: Effects of filtering on the vLBM prediction of a bubble rising in a viscous fluid (Case 2). Axial velocity fields predicted using (a) the second-order filter, Eq. (44); (b) the fifth-order filter with a = 0.25; (c) the seventh-order filter with a = 0.25 are shown (red: 0.012, blue: -0.006). Also shown are divergence fields predicted using (d) the second-order filter, Eq. (44); (e) the fifth-order filter with a = 0.25; (f) the seventh-order filter with a = 0.25 (red: 2×10−5 , blue: -2×10−5 ). The black line represents the shape of a bubble.

13

[2] Y. H. Qian, D. d’Humieres, P. Lallemand, Lattice BGK models for Navier–Stokes equation, Europhys. Lett. 17 (1992) 479–484. [3] S. Chen, G. D. Doolen, Lattice Boltzmann method for fluid flows, Annu. Rev. Fluid Mech. 30 (1998) 329–364. [4] R. Benzi, S. Succi, M. Vergassola, The lattice Boltzmann-equation - Theory and applications, Phys. Reports 222 (1992) 145–197. [5] S. Succi, Lattice Boltzmann Equation for Fluid Dynamics and Beyond, Clarendon Press, 2001. [6] G. Th¨ ommes, J. Becker, M. Junk, A. K. Vaikuntam, D. Kehrwald, A. Klar, K. Steiner, A. Wiegmann, A lattice Boltzmann method for immiscible multiphase flow simulations using the level set method, J. Comput. Phys. 228 (2009) 1139–1156. [7] X. Shan, H. Chen, Lattice Boltzmann method for simulating flows with multiple phases and components, Phys. Rev. E 47 (1993) 1815–1819. [8] M. R. Swift, E. Orlandini, W. R. Osborn, J. M. Yeomans, Lattice Boltzmann simulation of liquid-gas and binary fluid systems, Phys. Rev. E 54 (1996) 5041. [9] X. He, S. Chen, R. Zhang, A lattice Boltzmann scheme for incompressible multiphase flow and its application in simulation of Rayleigh-Taylor instability, J. Comput. Phys. 152 (1999) 642–663. [10] D. Jacqmin, Calculation of two-phase Navier-Stokes flows using phase-field modeling, J. Comput. Phys. 155 (1999) 96–127. [11] L. S. Luo, S. S. Girimaji, Theory of the lattice Boltzmann method: Two-fluid model for binary mixtures, Phys. Rev. E 67 (2003) 036302. [12] 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–644. [13] T. Lee, C.-L. Lin, A stable discretization of the lattice Boltzmann equation for simulation of incompressible two-phase flows at large density ratio, J. Comput. Phys. 206 (2005) 16–47. [14] 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–371. [15] V. L. Chenadec, H. Pitsch, A monotonicity preserving conservative sharp interface flow solver for high density ratio two-phase flows, J. Comput. Phys. 249 (2013) 185–203. [16] Y. Wang, C. Shu, H. B. Huang, C. J. Teo, Multiphase lattice boltzmann flux solver for incompressible multiphase flows with large density ratio, J. Comput. Phys. 280 (2015) 404–423. [17] C. M. Pooley, K. Furtado, Eliminating spurious velocities in the free-energy lattice Boltzmann method, Phys. Rev. E 77 (2008) 046702. [18] D. Ricot, S. Marie, P. Sagaut, C. Bailly, Lattice Boltzmann method with selective viscosity filter, J. Comput. Phys. 228 (2009) 4478–4490. [19] S. O. Unverdi, G. Tryggvason, A front-tracking method for viscous, incompressible, mulit-fulid flows, J. Comput. Phys. 100 (1992) 25–37.

14