## Abstract

Spurious currents near an interface between different phases are a common undesirable feature of the lattice Boltzmann equation (LBE) method for two-phase systems. In this paper, we show that the spurious currents of a kinetic theory-based LBE have a significant dependence on the parity of the grid number of the underlying lattice, which can be attributed to the chequerboard effect. A technique that uses a Lax–Wendroff streaming is proposed to overcome this anomaly, and its performance is verified numerically.

## 1. Introduction

The capability of the lattice Boltzmann equation (LBE) in simulating multi-phase flows has been recognized as one of its main advantages that distinguishes it from other numerical methods. The kinetic nature enables it to handle fluid–fluid interactions more easily than traditional methods based on hydrodynamic equations. Several types of LBE models have been proposed based on different physical pictures, such as colour-fluid models, pseudo-potential models, free-energy models and kinetic theory-based models, all of which have found their successful applications in multi-phase systems [1–3].

Although the LBE method has gained much success in the simulation of multi-phase flows, some limitations were also reported [2]. For instance, almost all of the existing multi-phase LBE models suffer from a common undesirable feature, i.e. the appearance of ‘spurious currents’ near an interface at equilibrium, which should vanish totally. Gunstensen [4] was the first to report this phenomenon, and later some other researchers have attempted to suppress the spurious currents from different viewpoints [5–8]. These studies revealed that spurious currents can be influenced by the discrete streaming operator, the force formulation and the order of isotropy of the discrete gradient operator in the evaluation of the interaction force. These studies are helpful in understanding the origin of the spurious current better.

The aim of this work is not to pursue a particular technique for reducing spurious currents, but to report an interesting and fundamental phenomenon related to the spurious current in LBE that has never been noticed in previous studies: it is found that the spurious currents in a kinetic theory-based LBE depend severely on the parity of the lattice employed, and the final equilibrium state of the two-phase system could be qualitatively different. After a simple analysis, it is found that this phenomenon is related to the chequerboard effects of the LBE. A simple technique is then proposed that can effectively remove such effects.

## 2. Lattice Boltzmann equation for two-phase flows

Although the two-phase LBE models mentioned in §1 were developed based on different physical pictures, it can be shown that they can all be reformulated in the following form [9]:
2.1
for *i*=0,…,*b*−1. Here *f*_{i}(** x**,

*t*) is the distribution function at position

**and time**

*x**t*associated with the discrete velocity

*c*_{i},

*δ*

_{t}is the time step and is the equilibrium distribution function, given by 2.2 where

*w*

_{i}is the weight,

*c*

_{s}is a model constant, and

*ρ*and

**are the fluid density and velocity defined by and , with**

*u***being the body force due to molecular interactions. The forcing term**

*F**F*

_{i}represents the effects of the interaction force

**and takes different formulations in different two-phase models. Here, we choose the multi-phase LBE based on the extended Boltzmann equation [9], in which**

*F**F*

_{i}is expressed as 2.3 It is noted that part of

*F*

_{i}can be regrouped into the relaxation term in equation (2.1), which was adopted in some previous works [7,9]. Using the Chapman–Enskog expansion, we can derive the hydrodynamic equations from the LBE (equation (2.1)) with a shear viscosity

*ν*=

*c*

^{2}

_{s}

*ρ*(

*τ*−0.5)

*δ*

_{t}.

For a van der Waals fluid, the force due to intermolecular interactions can be expressed as [9]
2.4
where *p*_{0} is the thermodynamic pressure dependent on the equation of state, *κ* is a parameter relating to the surface tension, and *μ*=*μ*_{0}−*κ*∇^{2}*ρ* is the chemical potential with *μ*_{0} being the chemical potential in the bulk phase. Both the thermodynamic pressure *p*_{0} and the bulk chemical potential *μ*_{0} can be found from the bulk free-energy density *ψ*_{0}, *p*_{0}=*ρ*∂_{ρ}*ψ*_{0}−*ψ*_{0} and *μ*_{0}=∂_{ρ}*ψ*_{0}. In some previous studies, the first part of equation (2.4) is termed the ‘pressure form’ while the second part is termed the ‘potential form’ [6,7]. In the pressure form, the quantity *c*^{2}_{s}*ρ*−*p*_{0} is treated as a scalar and discretized using some numerical schemes, while in the potential form the scalar *μ* is adopted for discretization. Although the two formulations are totally identical mathematically, their discrete versions may differ owing to some discretization errors, and these numerical differences may have significant influences on the spurious currents [6,7]. Both the pressure and potential forms of the interaction force have been adopted in some kinetic theory-based LBE models [7,10,11]. Furthermore, besides the kinetic theory-based models, the potential form of interaction force is also employed in some other types of LBE such as those constructed from phase-field theory [6,12]. It should be noted that the terms *κρ*∇∇^{2}*ρ* and *ρ*∇*μ* in equation (2.4) are not expressed in conservative formulations. Their discrete versions may lead to a slight violation of total momentum conservation [6], which has little influence on the overall interfacial dynamics.

In applications, the spatial gradients in the interaction force (2.4) should be discretized with suitable numerical schemes. A widely used method is the second-order isotropic central scheme (ICS),
2.5
where *ϕ* is an arbitrary function. Some other ICSs involving more neighbouring grids were also proposed in order to achieve higher isotropy [8]. Alternatively, Lee & Fischer [7] proposed an alternative second-order mixed difference scheme. In the present work, we will adopt the above ICS to discretize the interaction force in potential form.

## 3. Chequerboard effects on spurious currents

Chequerboard effects were first found in the lattice gas automata (LGA) method [13], the predecessor of the LBE. Owing to the finite symmetry of the discrete velocity set and the simple collision and streaming dynamics on the lattice, it was found that some spurious invariants, other than the mass, momentum and energy, exist in such lattice-based fluid models. Although the LBE method has overcome some shortcomings of the LGA method, it cannot completely eliminate the spurious invariants. Actually, chequerboard effects have been found in different LBE models (e.g. [14–17]).

Previous studies have shown that the chequerboard invariants can be excited by a solid boundary (e.g. [14]), and these unphysical spurious invariants can cause spatial oscillations and induce numerical instability. Although these studies were made for single-phase LBE models, the findings are also instructive for understanding some special behaviours of two-phase LBE models. For instance, the phase interface may act like a boundary and excite the chequerboard invariants, which further affect the spurious currents.

Now, we show that the chequerboard effects do influence the spurious currents. First, we consider a flat liquid–vapour interface. The bulk free energy of the fluid is assumed to be
3.1
where *β* is a constant, and *ρ*_{l} and *ρ*_{v} are the densities of liquid and vapour phases at saturation, respectively. We will use the two-dimensional nine-velocity square (D2Q9) LBE model to simulate this system. The discrete velocities of the D2Q9 model are given by *c*_{i}=*c**e*_{i} with *e*_{0}=(0,0), *e*_{1}=−*e*_{3}=(1,0), *e*_{2}=−*e*_{4}=(0,1), *e*_{5}=−*e*_{7}=(1,1) and *e*_{6}=−*e*_{8}=(−1,1); *c*=*δ*_{x}/*δ*_{t} is the lattice speed with *δ*_{x} the lattice spacing. Other parameters of this model are given as follows: and *w*_{0}=4/9, *w*_{1∼4}=1/9 and *w*_{5∼8}=1/36. We will take *δ*_{x}, *δ*_{t} and *c* as the units of length, time and velocity, respectively.

Initially, a liquid slab with density *ρ*_{l}=1.0 and thickness 50 is put at the centre of the lattice with *N*_{x}×*N*_{y} nodes, where *x* and *y* are set to be the directions parallel and normal to the interface, respectively. Other parts of the lattice are filled with the gas with *ρ*_{v}=0.2. In our simulations, periodic boundary conditions are applied in all directions. It is known that chequerboard invariants depend on the parity of the number of grids with periodic boundary conditions [16]. When both *N*_{x} and *N*_{y} are even numbers (even lattice), a chequerboard invariant may be decoupled into two sub-populations; on the other hand, the sub-populations will mix on the boundaries if *N*_{x} and *N*_{y} are odd numbers (odd lattice).

In figure 1, the profiles of the velocity component in the *y*-direction, *v*, are shown when two lattices are used. It is seen from figure 1*a* that the velocity oscillates between two constant values in each phase on the even lattice, which is a clear chequerboard pattern. It is also observed that the spurious currents in the liquid phase are weaker that those in the gas phase, and the overall magnitude of the spurious currents in the whole domain is of the order of 10^{−4}. Figure 1*b* shows the velocity profile on the odd lattice. Although the velocity also oscillates in the whole domain, the magnitude is much weaker than that on the even lattice. These results confirm the above speculation that the chequerboard invariants do have significant influences on the spurious currents, and therefore it can be expected that the spurious currents can be effectively reduced if the chequerboard effects are controlled. We also noted that the magnitude of the velocity depends on the initial condition of the density field, but the spatial oscillations always exist on the even lattice.

Some efforts have been proposed to suppress spurious invariants in previous studies. For instance, a time-average method has been proposed to remove the staggered invariant [14]. However, it was found that spatial oscillations can still exist even after time averaging [16]. On the other hand, Qian [15] proposed a fractional propagation scheme that can cancel the spurious invariants in both space and time. An LBE with more general propagation was developed later [18]. The principle of these methods is to introduce some spatial averaging during the time evolution so that the chequerboard invariants can be cancelled out. In this work, we will adopt a Lax–Wendroff (LW) scheme to control the chequerboard effects on the spurious currents, which can be expressed as two sub-steps [18]:
3.2
and
3.3
where *α*_{0}=1−*A*^{2}, *α*_{−1}=*A*(*A*+1)/2 and *α*_{1}=*A*(*A*−1)/2, with 0<*A*≤1 being a parameter. Clearly, as *A*=1 the LW–LBE reduces to the standard LBE. In LW–LBE, the lattice speed is defined as *c*′=*Aδ*_{x}/*δ*_{t}′, so if we take *δ*_{x} and *c*′ as the length and velocity units, respectively, the time step should be given by *δ*_{t}′=*Aδ*_{x}=*Aδ*_{t}. The Chapman–Enskog analysis shows that the viscosity of the D2Q9 LW–LBE is *ν*=*c*^{2}_{s}(*τ*′−0.5)*δ*_{t}′. Therefore, in order to keep the same viscosity as the standard LBE, the relaxation time *τ*′ should be taken as *τ*′=0.5+(*τ*−0.5)/*A*, where *τ* is the relaxation time in the standard LBE.

The LW–LBE with *A*=0.999 is applied to the same planar interface mentioned above. It is noted that in this case the LW–LBE is nearly identical to the standard LBE. But it will be shown later that the small but non-zero values of *α*_{0} and *α*_{1} in equation (3.3) play an important role in controlling the chequerboard effects. The velocity profiles at steady state are shown in figure 2. Clearly, the velocity profiles are nearly identical in both cases, and the chequerboard patterns disappear. Furthermore, unlike the standard LBE, the spurious velocity predicted by the LW–LBE only appears near the interface. These observations indicate that the proposed LW–LBE can suppress the chequerboard effects on the spurious currents effectively.

The standard LBE and the LW–LBE are also applied to a two-dimensional liquid–vapour system where a liquid drop with radius 25 is put at the centre of a square periodic lattice of size *N*×*N* with *N*=100 or 101. In this case, the chequerboard effects are more complicated than in the planar interface. The bulk free-energy density *ψ*_{0} is still given by equation (3.1), and the interaction force is discretized by the ICS with potential form. In figure 3, the time histories of the average kinetic energy, , and the velocity profiles across the domain centre are shown for the cases of the even and odd lattices. In the simulations, the iteration is stopped as the maximum relative error in the density field between two successive 10 steps is less than 10^{−10}. The results show that the kinetic energies in the two cases change similarly with time, but the velocity profiles show some differences; some oscillations again appear in the velocity profile when the even lattice is used. However, the overall velocity distributions are still consistent with each other, which can also be observed from the flow patterns shown in figure 4. It is also observed from figure 4 that the radius of the drop at the final steady state is slightly larger than the initial size, and some spurious currents appear near the interface.

The LW–LBE, however, gives parity-independent results for the drop system, as shown in figure 5. The time histories of the kinetic energy are still similar, like that shown in figure 3, and the oscillations in the velocity are completely removed in the case of the even lattice. Actually, the velocity profiles are indistinguishable on the two lattices.

We have also tested the standard and LW–LBE models with other formulations of bulk free-energy density and discretization methods (e.g. the ICS with pressure form and the Lee–Fischer scheme). Similar results are obtained, and the LW versions of these models are found again to be able to suppress the chequerboard effects significantly and yield parity-independent results.

## 4. Summary

Spurious currents are numerical artefacts of LBE for two-phase flows, while chequerboard invariants always exist in LBE as a result of the finite symmetry of the discrete velocity set and the collision/streaming dynamics. The spurious invariants may be excited by the liquid–vapour interface and then affect the spurious currents. In this work, we have investigated the chequerboard effects on the spurious currents of a kinetic theory-based LBE, and found that chequerboard patterns do appear when an even lattice is used with periodic boundary conditions, and the chequerboard effects have a significant influence on the spurious currents. An LW–LBE was then proposed, which has been shown to be able to suppress the chequerboard effects effectively.

## Acknowledgements

This work is supported by the National Natural Science Foundation of China (10972087, 50936001) and the National Basic Research Programme of China (2011CB707305). Z.L.G. is supported by the Programme for NCET in the University (NCET-07-0323).

## Footnotes

One contribution of 26 to a Theme Issue ‘Discrete simulation of fluid dynamics: methods’.

- This journal is © 2011 The Royal Society