Royal Society Publishing

A local lattice Boltzmann method for multiple immiscible fluids and dense suspensions of drops

Timothy J. Spencer, Ian Halliday, Chris M. Care


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 105 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, Embedded Image, is split into separate PDFs, each representing a different fluid type: Embedded Image. 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], Embedded Image 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 Embedded Image and the macroscopic observables of species density, total density and momentum are Embedded Image, Embedded Image, ρ=ρ0+ρ1, Embedded Image, where Embedded Image is the lattice link velocity. We also define a local phase field Embedded Image.

The colour-blind PDF undergoes a collision rule with a source or force term, Embedded Image. We present here, for definiteness, the Bhatnagar–Gross–Krook (BGK) collision dynamics (other models, like the multi-relaxation time and incompressible models, equally work) Embedded Image 2.1 where Embedded Image represents the post collision, pre-segregation state, τ is the combined fluids’ dimensionless relaxation time, △t is the time step and Embedded Image is defined as Embedded Image 2.2 in which ti and c2s 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 Embedded Image 2.3 in which Embedded Image 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, Embedded Image 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 xc, with Embedded Image equation (2.3) gives rise to density profiles Embedded Image 2.4 from which we see β01 controls interface width, L01, and has units of mass per unit length. Importantly, we note that for this form of segregation rule we can write Embedded Image.

Returning to the force term, Embedded Image, 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 (Embedded Image). Note that to obtain the necessary curvature unfortunately involves non-local gradients. In Halliday et al. [7], an interface force Embedded Image is calculated directly, input into the LBM body force term, Embedded Image, 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 Embedded Image 2.5 Here, A01 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 Fi 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 Embedded Image 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 Embedded Image 2.7 Embedded Image 2.8 and Embedded Image 2.9 where Embedded Image and the transport coefficients and pressure are Embedded Image 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 Embedded Image, 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 Embedded Image and Embedded Image 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 Embedded Image 2.11

From equation (2.9), provided the force term generates a correct Laplace law pressure step, the term Embedded Image 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/cs≪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, Embedded Image which is in turn normalized, Embedded Image. 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, Embedded Image. We then form a crude approximation to the phase field gradient Embedded Image. Although this is a poor approximation of the gradient it remains a good approximation to its direction. A normalization (Embedded Image) 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, Embedded Image, 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.

Figure 1.

Range of interfacial tensions recovered with the local MCLB method. Inset (a) interfacial pressure profile, (b) analytical versus numerical comparison. Filled circles, data; solid line, linear slope 1/r. (Online version in colour.)

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 Embedded Image where η0 is the consistency parameter, n is the power-law index and Embedded Image is the local shear rate. Despite some of the power-law model shortcomings (Embedded Image and the units of η0 are Pa sn) the solution for channel flows is analytically available for shear-thinning (n<1) and shear-thickening (n>1) fluids: Embedded Image, 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, TC, of fluid species. In the two-dimensional systems here, this number is empirically TC≤6. Consequently, on any link only ≤TC PDFs indexed by an associated label identifying the fluid species (which may range Embedded Image) 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 TTC closely approximate those with TTC.

(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 2H/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ρU2r/σ∼16.7. The two-dimensional particle concentration is ϕ=(T−1)πr2/2HL. The solid walls are modelled with the half link bounce back method and the background fluid preferentially wets the walls.

Figure 2a 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 2b. Increasing concentration decreases the particle free regions near the wall (figure 2c).

Figure 2.

(a) Channel flow of 180 identical deformable particles, (b) time-averaged concentration profiles, (c) decreasing wall particle free region, (d) velocity profiles and power-law fitted data (plus symbols, ϕ=0; cross symbols, ϕ=0.053; asterisks, ϕ=0.106; open squares, ϕ=0.159; filled squares, ϕ=0.212; open circles, ϕ=0.265; filled circles, ϕ=0.318; open triangles, ϕ=0.371; filled triangles, ϕ=0.424; open inverted triangles, ϕ=0.475), (e) concentration versus power index relationship, (f) effective viscosity increase (filled circles, simulation; solid line, 1+0.158ϕ−1.799ϕ2+ 9.931ϕ3). (Online version in colour.)

Figure 2d 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 2e, 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=η04UH/3Q, where U is the maximum velocity in a system without particles, the effective viscosity increases with the particle concentration as shown in figure 2f.

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.


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



View Abstract