## Abstract

The lattice Boltzmann method (LBM) for computational fluid dynamics benefits from a simple, explicit, completely local computational algorithm making it highly efficient. We extend LBM to recover hydrodynamics of multi-component immiscible fluids, while retaining a completely local, explicit and simple algorithm. Hence, no computationally expensive lattice gradients, interaction potentials or curvatures, that use information from neighbouring lattice sites, need to be calculated, which makes the method highly scalable and suitable for high performance parallel computing. The method is analytical and is shown to recover correct continuum hydrodynamic equations of motion and interfacial boundary conditions. This LBM may be further extended to situations containing a high number (O(100)) of individually immiscible drops. We make comparisons of the emergent non-Newtonian behaviour with a power-law fluid model. We anticipate our method will have a range applications in engineering, industrial and biological sciences.

## 1. Introduction

The lattice Boltzmann method (LBM) [1] solves a minimal form of the differential kinetic Boltzmann equation for a mesoscopic particle distribution function (PDF) whose result matches the continuum Navier–Stokes equations. The LBM numerical scheme is a simple two-step algorithm, consisting of PDF collision and then streaming with several noteworthy differences from traditional schemes [1]. Among these, we highlight that the collision step and all hydrodynamic quantities of interest (pressure, velocity, shear rate, stress) are computed locally, on a grid node, avoiding any need for (i) data communication between nodes and (ii) solving a Poisson equation for pressure. Boundary conditions are readily implemented, from simple, local, microscopic rules. LBM, then, is a high performance, O(2) accurate, simple, explicit, parallel algorithm scalable to greater than 10^{5} processing cores [2], which opens up many new avenues in applied computational science.

LBMs for multi-component flows have popularly been taken up. Benefitting from some of LBM’s intrinsic features, various multi-component lattice Boltzmann (MCLB) models exist. Briefly, we note only the more popular (original) variants. The method of Shan & Chen [3] uses a microscopic interaction potential to model phase changes and interfacial tension behaviour. The method of Swift *et al*. [4] rests upon a total free energy, consisting of Landau phase transition terms and surface tension terms, after Cahn–Hilliard theory. The method of Lee & Fischer [5] uses another variant of the Cahn–Hilliard approach, combining a potential form of the free energy and a high-order finite differencing scheme. The method by Gunstensen *et al*. [6] evolves earlier lattice gas automata models, and has sharp-interfacial immiscibility (a source of lattice pinning error) but an appropriate interfacial tension may be directly imposed. The model of Halliday *et al*. [7] similarly assumes arrested coalescence and immiscibility, but uses a diffuse interface, with directly imposed interfacial tension.

Of these MCLB models, methods of Shan & Chen [3], Swift *et al*. [4] and Lee & Fischer [5] are useful in problems where phase separation kinematics features, whereas the methods of Gunstensen *et al*. [6] and Halliday *et al*. [7] avoid these calculations, by considering completely phase-separated, isothermal systems. Only the method of Lee & Fischer [5] is free from so-called spurious velocities due, in main, to its appropriate, non-local, treatment of both directional and non-directional derivatives separately. This does not exclude the use of the alternative schemes; rather, it limits their use to cases in which the hydrodynamic velocity is larger than any spurious velocities.

The common feature of all MCLB models, and indeed many other multi-component numerical methods, is that they calculate gradient terms, ranging from first to third order, which require information from neighbouring lattice nodes and, sometimes, additional sweeps of the underlying lattice algorithm, which discards some of the key advantages of the LBM concept, in which computations were local.

Here, we aim to demonstrate a new MCLB method that retains this advantage of core LBM, namely, local computations and an efficiently scaling parallel algorithm. Previously, a method by Santos *et al*. [8] had developed a local MCLB method based on the concept of field mediators; however, the method contains a velocity-dependent interfacial tension whereas in the proposed model we avoid this and develop, possibly, a simpler algorithm. Our method is based on the isothermal MCLB method of Halliday *et al*. [7], which has been shown to be analytically tractable, to recover the results of analytical calculations of interfacial deformation [9], to accommodate extra force terms to allow moderate density ratio and viscosity ratio mixtures [10] and to accommodate dynamic wetting conditions [11]. We then apply this new method to study emergent, non-Newtonian properties in suspensions of deformable particles, in a channel flow.

## 2. A local multi-component lattice Boltzmann method

Following Halliday *et al*. [7], the LBM PDF, , is split into separate PDFs, each representing a different fluid type: . Here, the superscript indexes fluid type (species or colour). It should not be confused with orders of a Chapman–Enskog expansion. The subscript refers to the usual lattice links [1], is the discrete lattice coordinate, *t* the discrete time and *T* the number of different fluid species.

For now, let *T*=2. The total *colour-blind* PDF is and the macroscopic observables of species density, total density and momentum are , , *ρ*=*ρ*^{0}+*ρ*^{1}, , where is the lattice link velocity. We also define a local phase field .

The colour-blind PDF undergoes a collision rule with a source or force term, . We present here, for definiteness, the Bhatnagar–Gross–Krook (BGK) collision dynamics (other models, like the multi-relaxation time and incompressible models, equally work)
2.1
where represents the post collision, pre-segregation state, *τ* is the combined fluids’ dimensionless relaxation time, △*t* is the time step and is defined as
2.2
in which *t*_{i} and *c*^{2}_{s} are lattice specific weights (see [1]).

The force term, defined later, is chosen to impose the Laplace law pressure step. The collided distribution function (equation (2.1)) then undergoes segregation and propagation
2.3
in which is the interfacial unit normal vector (UNV) between fluid species 0 and 1, and *β*^{01} is a constant related to the interface. Note that product *ρ*^{0}*ρ*^{1}/(*ρ*^{0}+*ρ*^{1})^{2} is evaluated only in the interfacial region, and all quantities on the right-hand side of equation (2.3) are local to the lattice node. The detailed properties of the analytical fluid segregation rules (equation (2.3)) have been analysed in Halliday *et al*. [7]. Its principal advantages over other segregation methods, such as Tolke *et al*. [12] and Gunstensen *et al*. [6], are that interfaces are isotropic, the dynamics of the phase field are appropriate and lattice pinning effects are eliminated. However, the cost is a more diffuse interface. For a flat interface centred on *x*_{c}, with equation (2.3) gives rise to density profiles
2.4
from which we see *β*^{01} controls interface width, *L*^{01}, and has units of mass per unit length. Importantly, we note that for this form of segregation rule we can write .

Returning to the force term, , we require this to effect a Laplace pressure step (△*P*) across an interface, proportional to the fluids’ interfacial tension (*σ*^{01}) and mean radius of curvature (). Note that to obtain the necessary curvature unfortunately involves non-local gradients. In Halliday *et al*. [7], an interface force is calculated directly, input into the LBM body force term, , and effectively imposed through its first moment. To eliminate the explicit calculation of the curvature, we choose to impose the force in the second moment of the LBM force term only by using a novel form
2.5
Here, *A*^{01} is proportional to the emerging interfacial tension that is also effected by the *β*^{01} parameter and, as with equation (2.3), all quantities are local to the lattice node. We note this form of *F*_{i} does not require redefinition of the velocity moment, which is the case in most other MCLB methods. Also, we retain a simple equation of state for the pressure field, adding to the simplicity of this method. Summarizing, the moments of equation (2.5) are
2.6
Taking the divergence of the second moment, above, results in a Laplace or capillary force term close to that first used and derived by Lafaurie *et al*. [13]. One interpretation of this expression is as a bulk pressure perturbation in the vicinity of the interface and projected onto the interfacial surface.

Executing a Chapman–Enskog analysis for mass, momentum [14] and phase fields [7] identifies the continuum transport equations on standard LBM lattices [1] and up to second order in the velocity as
2.7
2.8
and
2.9
where and the transport coefficients and pressure are
2.10
Here, we note the LBM force term appearing on the right-hand side of equation (2.8) can be viewed in two ways. First, by inserting equation (2.6) and making use of the product rule, it is straightforward to show that it represents a force , identical to the force expression used in Halliday *et al*. [7]. Second, quantity △*σ*_{αβ} may be viewed as a stress perturbation and, as such, should obey interfacial kinematic boundary conditions and of continuity of the fluid velocity. We note the *a priori* interfacial tension expression. A local expression for the viscous stress tensor is now given by
2.11

From equation (2.9), provided the force term generates a correct Laplace law pressure step, the term and, thus, the right-hand side of equation (2.9) is zero indicating the phase field is advected with the underlying fluid velocity. We also note, if using BGK collision dynamics, to reduce compressibility errors the condition *Ma*=*u*/*c*_{s}≪1 should be maintained.

To complete the MCLB method, summarized in equations (2.1), (2.3) and (2.5), we must calculate the interfacial UNV. We show two ways in which to do this.

#### (i) Non-local multi-component lattice Boltzmann method

Here, standard finite differences are used to calculate the lattice gradient of the phase field, which is in turn normalized, . Using this method results in a reliable and robust algorithm but ultimately non-local algorithm, owing to the use of finite differences.

#### (ii) Local multi-component lattice Boltzmann method

To approximate a local interfacial UNV (hence local MCLB behaviour) we define a meso-scale, or *link* phase field, . We then form a crude approximation to the phase field gradient . Although this is a poor approximation of the gradient it remains a good approximation to its direction. A normalization () provides a local estimation of the interfacial UNV. Using this method together with equations (2.1), (2.3) and (2.5) results in a completely local MCLB model with obvious advantages in parallel communication.

## 3. Verification of the local multi-component lattice Boltzmann method

Consider a two-dimensional circular droplet of fluid species *ρ*^{1} embedded in background fluid *ρ*^{0}. The law of similarity is used to match key dimensionless *Reynolds* and *Weber* numbers. In hydrostatic equilibrium, the interfacial UNVs are known analytically. Relative to a polar origin at the centre of mass, , where *θ* is the angular position of the interface. Typical (percentage) error in the local method interfacial UNVs, compared with the analytical values, is less than or equal to 1.5 per cent.

The radial pressure profile shows smooth density variation across the interface, according to the Laplace law, over four decades of interfacial tension value, facilitating a wide range of materials (figure 1). Using the local MCLB method increases the magnitude of the interfacial spurious velocities by a factor of three, owing to the 1.5 per cent inaccuracies in the local UNV calculation. However, we note these values still compare well with all other MCLB and numerical methods (with the exception of [5]). All these methods should always ensure simulation velocities exceed the spurious velocity values, for their effect to be negligible and to underwrite kinematic boundary conditions.

In simulations, one should check that local UNV fields remain continuous within interfacial regions. In tests of droplet break-up and thread formation it can lead to inaccuracies in the local method from which our non-local MCLB counterpart, above, is free. However, continuum methods are known to be inaccurate for droplet break-ups involving threads of a few molecules in width. Moreover, a wide range of applications that do not involve drop scission remain accessible, as we now show.

## 4. Multi-component lattice Boltzmann for suspensions of deformable particles

Particle suspensions give rise to non-Newtonian flow behaviour and flows of deformable particles are less commonly studied than those of solid particles. We follow Zhou & Pozrikidis [15] in studying a pressure-driven, two-dimensional channel suspension flow of identical, deformable particles (drops). Zhou & Pozrikidis [15] studied the effect of relatively few drops *T*≤24, using the computationally expensive boundary integral method. Here, we study a similar system but with a much larger range of particle numbers. We then examine emergent, non-Newtonian behaviour and compare with the well-known non-Newtonian power-law fluid model. The power-law model proposes a shear-dependent shear viscosity where *η*_{0} is the consistency parameter, *n* is the power-law index and is the local shear rate. Despite some of the power-law model shortcomings ( and the units of *η*_{0} are Pa s^{n}) the solution for channel flows is analytically available for shear-thinning (*n*<1) and shear-thickening (*n*>1) fluids: , where *Q* is the discharge. Here, we investigate the relationship of particle concentration and the power-law index.

### (a) Multi-component lattice Boltzmann algorithm extension for more than two fluids

Following Dupin *et al*. [16] we extend the two-component MCLB method of §2 to any number of components *T*≫2. The principle behind this method has a PDF for each immiscible fluid species. However, to avoid excessive computational memory demands when *T*>2, we apply an exclusion principle because, for immiscible fluids with a diffuse interface, any one location in space (lattice node) can only accommodate a critical, finite number, *T*_{C}, of fluid species. In the two-dimensional systems here, this number is empirically *T*_{C}≤6. Consequently, on any link only ≤*T*_{C} PDFs indexed by an associated label identifying the fluid species (which may range ) are tracked. The simplest way to do this is to ensure PDFs stream onto similarly indexed lattices. The resultant algorithm’s memory consumption and computational timings for *T*≫*T*_{C} closely approximate those with *T*≈*T*_{C}.

### (b) High-density deformable particle suspension in a channel flow

To narrow the wide parameter space in suspension flows, we investigate a mono-disperse suspension of *T*−1 identical immiscible deformable particles in a background fluid in a periodic channel flow driven by a fixed pressure gradient. The viscosity and density ratios between fluid species have unit value and drop–drop interfacial tension is 10 times the drop–background fluid value, to mimic surfactant effects and ensure drops remain mono-disperse. The drop radius to channel half height ratio *r*/*H*=0.066 and the channel aspect ratio 2*H*/*L*=0.74 and we use lattice values *H*=148, *τ*=1/2 and *ρ*=1.8. Relevant dimensionless numbers are the no-particle Reynolds=*UH*/*ν*∼35.5 and Weber=2*ρU*^{2}*r*/*σ*∼16.7. The two-dimensional particle concentration is *ϕ*=(*T*−1)*πr*^{2}/2*HL*. The solid walls are modelled with the half link bounce back method and the background fluid preferentially wets the walls.

Figure 2*a* shows a pseudo steady-state snapshot of the flows studied. Particle concentrations *ϕ*,0<*ϕ*<0.5, were studied from an initially random arrangement of particle positions. Single particles experience lift to the channel centre line, where flow is fastest and deformation least. When *ϕ* is increased, all the drops vie for a central channel position and layers of particles form around the channel centre, with more central layers at higher concentration, as evidenced in the short time-averaged cross-sectional averaged concentration profiles shown in figure 2*b*. Increasing concentration decreases the particle free regions near the wall (figure 2*c*).

Figure 2*d* plots short time-averaged velocity profiles at various concentrations. Clearly, higher particle concentrations cause flow blunting. The velocity profiles can be compared with the analytical power-law solution to measure the power-law index. In figure 2*e*, we find *n* decreases with increasing concentration but it tends to plateau, around *ϕ*∼0.32. Defining the effective channel flow viscosity after Zhou & Pozrikidis [15] as *η*_{eff}=*η*_{0}4*UH*/3*Q*, where *U* is the maximum velocity in a system without particles, the effective viscosity increases with the particle concentration as shown in figure 2*f*.

## 5. Conclusion

An isothermal, immiscible MCLB method has been developed. It is simple, explicit, analytical and local. The local computations remove the need for lattice gradients and excessive parallel communication overheads found in previous forms. Tests show this method recovers key features of immiscible multi-component fluids. Density differences between fluids were not presented but an additional force [10] will model this. The MCLB method was extended to include arbitrary numbers of immiscible fluid species and to study high-density suspensions. Channel flows of deformable particle suspensions showed emergent non-Newtonian shear thinning and particle layering at all concentrations. A relationship between a power-law fluid power index and particle concentration was explored.

## Acknowledgements

This work has been supported by a research grant from BBSRC (UK) BB/F013744/1.

## Footnotes

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

- This journal is © 2011 The Royal Society