## Abstract

We describe some scaling issues that arise when using lattice Boltzmann (LB) methods to simulate binary fluid mixtures—both in the presence and absence of colloidal particles. Two types of scaling problem arise: physical and computational. Physical scaling concerns how to relate simulation parameters to those of the real world. To do this effectively requires careful physics, because (in common with other methods) LB cannot fully resolve the hierarchy of length, energy and time-scales that arise in typical flows of complex fluids. Care is needed in deciding what physics to resolve and what to leave unresolved, particularly when colloidal particles are present in one or both of two fluid phases. This influences steering of simulation parameters such as fluid viscosity and interfacial tension. When the physics is anisotropic (for example, in systems under shear) careful adaptation of the geometry of the simulation box may be needed; an example of this, relating to our study of the effect of colloidal particles on the Rayleigh–Plateau instability of a fluid cylinder, is described. The second and closely related set of scaling issues are computational in nature: how do you scale-up simulations to very large lattice sizes? The problem is acute for systems undergoing shear flow. Here one requires a set of blockwise co-moving frames to the fluid, each connected to the next by a Lees–Edwards like boundary condition. These matching planes lead to small numerical errors whose cumulative effects can become severe; strategies for minimizing such effects are discussed.

## 1. Introduction

The lattice Boltzmann equation (LBE) is a widely used lattice formulation of fluid mechanics (Succi 2001). It offers a faithful discretization of the Navier–Stokes equation of isothermal, near-incompressible fluid flow, and is very well adapted to parallel computation (Amati *et al*. 1997). The LBE approach is particularly adapted to simulating mesoscopic problems (Coveney & Succi 2002). These include, for example, porous medium flows, and flows of complex and multi-component fluids with microstructure (Ladd 1994*a*,*b*; Swift *et al*. 1995; Kendon *et al*. 1999; 2001; Manz *et al*. 1999; Nekovee *et al*. 2000; Pagonabarraga *et al*. 2001; Love *et al*. 2003).

For binary fluids, the LBEs read (with time-step Δ*t*=1)(1.1)(1.2)Here * x* is a lattice position;

**c**_{i}is the local velocity set corresponding to a set of near-neighbour vectors of the lattice (selected to recover rotational and Galilean invariance properties);

*f*

_{i}is the local distribution function for molecular velocities at the given site and its equilibrium value (dependent on local fluid density and flow velocity). The distribution function

*g*

_{i}controls the evolution of a binary fluid order parameter

*ϕ*; its equilibrium value encodes both a well-chosen free energy functional, , and an order parameter mobility (Swift

*et al*. 1995; Kendon

*et al*. 2001). The matrices

*L*

_{ij}and

*M*

_{ij}are collision operators;

*L*

_{ij}is often chosen as a single-relaxation time (lattice BGK) matrix −

*δ*

_{ij}/

*τ*, in which case the fluid viscosity is with

*c*

_{s}the sound speed (controlled by the choice of

**c**_{i}). In our work

*M*

_{ij}is chosen of BGK form with

*τ*=1, so that

*g*

_{i}is reset to every time-step. The density fields and are the zeroth moments of the distribution functions; the fluid mass fluxes obeys and . Finally, the second moment of determines the kinetic stress tensor (or momentum flux tensor).

These standard LBEs neglect thermal noise and lead to deterministic (zero temperature) evolution of the fluid velocity field * v* and order parameter

*ϕ*; the continuum level equations to which they correspond are the isothermal Navier–Stokes equation for

*, properly coupled to an advection–diffusion equation for*

**v***ϕ*. LB must only be used at low Mach number when the fluid density

*ρ*remains nearly constant. It is possible to add thermal noise to the fluid either at the level of low wavevector hydrodynamics via the stress tensor (Ladd 1994

*a*,

*b*), or consistently across all wavevectors by developing discrete Langevin equations of which the above are the deterministic limit (Adhikari

*et al*., in press). When noise is referred to below, it has been included by the latter route, but in the fluid momentum sector only (

*f*

_{i}but not

*g*

_{i}), corresponding to a thermal fluid in which the order parameter dynamics are ‘cold’, so that thermal capillary waves at the fluid–fluid interface are neglected. Work on the full incorporation of Langevin noise into a binary fluid model is underway. Meanwhile, this method allows incorporation of Brownian colloidal particles within a binary fluid system.

## 2. Physical scaling issues

We first address issues of physical scaling, that is, how to choose parameter values in LB that map onto those of the real world. At one level this is trivial, at least in our own implementation of LB (Desplat *et al*. 2001) in which thermodynamic quantities such as interfacial tension *σ*=(8*κ*|*A*^{3}|/9*B*^{2})^{1/2} between two fluids follow from a well-defined free energy functional. The thermodynamic and kinetic properties of an LB fluid (or fluid mixture) are then defined, through the choices made above for *L*_{ij}, *M*_{ij} and , in terms of lattice units for mass, length and time. The mass unit contains an arbitrary scale factor (we set *ρ*=1) whereas lattice units for length and time are directly linked to the lattice constant and time-step in the simulation. In principle there is no problem, therefore, choosing all the thermodynamic and kinetic parameters in LB to correspond to those of a real physical system, by an appropriate mapping of lattice units onto real ones.

In practice, of course, this is rarely possible. Setting aside various stability issues (Kendon *et al*. 2001), the run-times and system sizes required would be simply too big, as illustrated by the examples that follow. Accordingly, LB is generally interpreted as a ‘mesoscale’ simulation method which throws away local information so as to achieve much larger length- and time-scales than could be achieved with, say, molecular dynamics. On the other hand, LB does not throw away as much information as does computational fluid dynamics (CFD), which discretizes directly the Navier–Stokes equation and treats any interfaces as structureless.

### (a) What is the difficulty?

Consider the problem of disconnection, in which a droplet of one fluid surrounded by another (for simplicity, of equal viscosity) breaks into two droplets via a pinchoff event. (This was studied carefully in two dimensions by Wagner *et al*. 2003*b*.) At the CFD level this can be modelled using a Navier–Stokes continuum, with a structureless interface between the fluids, characterized only by an interfacial tension *σ*. But this problem is singular (Eggers 1997) so that direct CFD methods fail, unless explicit measures are taken to identify incipient reconnection of the interface and deal with it ‘by hand’. Such failure of CFD is inevitable, because the ultimate pinchoff of a thin neck is controlled by molecular, not continuum, physics. On the other hand, LB maintains a coarse-grained representation of the molecular physics through its treatment of distribution functions on the lattice. This leads to the smooth pinchoff of fluid necks at the point where their diameter approaches the interfacial width (generally a couple of lattice spacings). When LB is used for binary fluid simulations, this removes the singularities that would obstruct a direct approach, and does so in a physically ‘realistic’ way (effectively by molecular dissolution of the minority fluid under the influence of Laplace pressure).

On the other hand, the LB treatment of this and similar problems are rarely *fully* realistic (Cates *et al*. 2004): the mechanism is valid, but the scale at which it sets in is not. To see this, note that in the real world the molecular cut-off point is not reached until the diameter of a fluid neck is at sub-nanometre level. In LB the molecular physics takes over at the scale of the lattice constant, and if this was always linked to a physical scale of (say) 0.1 nm, there would be no possibility of simulating systems larger than 100 nm across, even with the largest computers available. (These can run 1024^{3} lattices (Harting *et al*. 2004) although 128^{3} is routinely used in our own work.) Fortunately this limitation, which would restrict LB to broadly the same parameter domain as molecular dynamics, is more illusory than real. For, so long as the pinchoff process itself does not dominate the physics, it does not matter whether or not this is resolved in a *fully* realistic manner. The advantage of having a realistic mechanism remains: because the pinchoff process is physically admissible (albeit for some set of fluid parameters that may remain strictly hypothetical in the laboratory), it is less likely to lead to gross artefact than an arbitrary but well-intentioned numerical discretization. Having said this, gross artefacts are sometimes easier to detect and control than subtle ones; LB is quite prone to the latter, and careful testing for systematic errors is always needed (see e.g. Kendon *et al*. 2001 for a fuller discussion).

### (b) Binary fluid demixing

A good example of the mesoscale use of LB simulation is in the coarsening of fluid domains after spinodal decomposition in a pair of fully symmetric binary fluids (Kendon *et al*. 2001). The domain length-scale *L*(*t*) (appropriately defined; see below) at time *t* after quench in such a system is widely held to obey a scaling of the form (Bray 1994; 2000; Onuki 2002)(2.1)where *L*_{0}=*η*^{2}/(*ρσ*) is a characteristic length, *T*_{0}=*η*^{3}/(*ρσ*^{2}) is a characteristic time and *η* and *ρ* are the (equal) viscosities and interfacial tension of the fluids. Since *L*_{0} and *T*_{0} involve only macroscopic quantities, the observation of such scaling in experiments (Onuki 2002) is evidence that the molecular physics of pinchoff does not matter. Moreover, such scaling implies that results for LB runs with widely different *η* and *σ* should collapse onto each other when scaled onto an appropriate plot. This is indeed seen (Pagonabarraga *et al*. 2002), showing that in LB, as in the real world, the details of pinchoff are irrelevant. However, careful measures must first be taken to eliminate artefacts such as, for example, excessive interdiffusion of the two fluids on the domain scale *L*(*t*). These measures amount to demanding that a specific hierarchy of physical effects is obeyed: on the time-scale of coarsening, diffusion must be easily accomplished across the width of a fluid–fluid interface (a few lattice sites), so as to maintain an equilibrium surface tension, but simultaneously diffusion by the same mechanism must be negligible at the domain scale.

### (c) Dimensionless control parameters

In applying LB method to a wider class of problems in complex fluids modelling, it is wise to be explicitly conscious of the hierarchy of scales that the real-world problem presents. Only by respecting this hierarchy within the simulation can one expect realistic results. In fact, reckoning with the hierarchy in physical terms is needed, even to decide whether the problem lies within the scope of LB for any given computational resource. The best way to develop these ideas is using the dimensionless numbers beloved of fluid-mechanicists, who have long since covered a similar territory in developing scaling relations between one set of experimental parameters and another. When done well, these analyses pay as much attention to discussing what can be neglected as to what must be included (Batchelor 1967; Faber 1995).

In the binary fluid coarsening problem, the two relevant numbers are the capillary number and the Reynolds number . The first measures the ratio of viscous to interfacial forces; the second the ratio of inertial to viscous forces. (A third combination *Re* *Ca* therefore measures the ratio of inertial to interfacial forces.) At early times *t*≪*t**, viscous and interfacial forces are in balance so that *Ca*∼1 (giving *L*(*t*)/*L*_{0}∼*t*/*T*_{0}) and *Re*≪1. At late times *t*≫*t**, inertial and interfacial forces balance so that *Re* *Ca*∼1 (giving *L*/*L*_{0}∼(*T*/*T*_{0})^{2/3}) and *Ca*≪1. Here *t** is a crossover time ‘of order’ *T*_{0}; our LB work fixes the prefactors in the two coarsening laws and shows that in fact *t**≃10^{4}*T*_{0} (Kendon *et al*. 1999), perhaps stretching the interpretation of the term ‘of order’.

Intriguingly, careful investigation of one relatively simple-looking extension to this problem—binary fluid demixing under gravity in the presence of a continuous temperature ramp—throws up at least three new dimensionless numbers (Cates *et al*. 2003) to add to the many that already arise in the scaling of complex fluids. Indeed, not just Reynolds, but Peclet, Schmidt, Rayleigh, Bond, Prandtl, Benard, Weissenberg, and more, can feature under different experimental conditions. (Such surnames give no direct clue as to the actual physics of these numbers; names like ‘inertial/viscous ratio’ would be a lot more informative than ‘Reynolds number’, but convention precludes them.) The relatively subtle hierarchies that arise in problems such as the temperature-ramped binary fluid pose a severe test on the LB methodology. In that particular case, LB work remains in progress, with preliminary results described by Wagner *et al*. (2003*a*).

### (d) Colloids in binary fluids

The addition of colloidal particles to a binary fluid system raises further scaling considerations. Our LB algorithm for simulating colloids is based on that of Nguyen & Ladd (2002); it is outlined briefly by Cates *et al*. (2004) and full details will be published elsewhere (Stratford *et al*. in preparation).

Consider first a simple geometry in which a suspension of colloids, at low volume fraction, sediments under gravity within a stratified pair of layers of two immiscible fluids of equal viscosity. Figure 1 shows this process. The colloidal particles are initially placed at random; they then fall under gravity and, for the parameters selected here, end up attached to the interface. The interfacial tension between fluids is *σ*; the solid–fluid tensions are equal, so that particles at the interface between fluids have a 90 degree wetting angle; both fluids have viscosity *η*; the colloid radius is *a*. The relevant dimensionless numbers are those of Reynolds, Peclet and Bond(2.2)(2.3)(2.4)Here *v*_{sed}=Δ*mg*/6*πηa* (with Δ*m*=4*πa*^{2}Δ*ρ*/3) is the sedimentation speed of a colloid, and *D*=*k*_{B}*T*/6*πηa* its diffusion constant; *v*_{0}=*σ*/*η* is a characteristic velocity associated with coarsening (it obeys *v*_{0}=*L*_{0}/*T*_{0} as defined above). The sedimentation velocity is fixed by balancing gravity against viscosity; accordingly the Bond number (ratio of gravity to interfacial force) is essentially the same as the capillary number, which we defined above as the ratio of viscous to interfacial force. (Note that Cates *et al*. (2004) used an inverted definition of capillary number, so that *Ca* in that paper means 9*Bo*^{−1}/2; to avoid further confusion, we use the Bond number, rather than capillary, in what follows.)

For typical colloids in binary fluid mixtures one has *Re*≃10^{−7}(*a*/*a*_{μ})^{3}; *Pe*≃(*a*/*a*_{μ})^{4} and *Bo*≃10^{−8}(*a*/*a*_{μ})^{2}, where *a*_{μ}≡1 μm. Thus a colloid of, say, *a*=10 μm has negligibly small inertial effects (small *Re*); negligibly small diffusion (large *Pe*); and a very large capillary force attaching it to the fluid–fluid interface, easily capable of holding it there against gravity (very small *Bo*). For *a*≤100 nm, *Re* and *Bo* are still negligible, although *Pe* is now also small: diffusion is important. But even for 5 nm colloids, the combination *Pe*/*Bo* remains large; this represents the ratio of interfacial to thermal forces. Thus, for all reasonable interfacial tensions and colloid sizes, neutrally wetting colloidal particles adsorb quasi-irreversibly to the fluid–fluid interface, and *diffusive* detachment does not take place (Binks 2002; Aveyard *et al*. 2003).

Suppose, therefore, that we wish to simulate *a*=0.1 μm radius neutrally wetting colloids in an oil water mixture. Thermal noise can be neglected; and as outlined above in the case of gravity (but generally true) for typical colloidal problems *Re* is extremely small. However, LB works by solving dynamically a discretization of the full hydrodynamic equations. Thus, as explained in more detail by Cates *et al*. (2004), small Reynolds number can only be achieved by having a colloidal velocity that is, at most, of order *Re* in lattice units. For a run such as that of figure 1 in which colloids move a distance of order the simulation box (say 100 lattice sites), the number of time-steps required for *a*/*a*_{μ}=0.1 then obeys *N*∼100/*Re*≃10^{12}. This run-time is roughly a million times longer than anything practically achievable on the world's fastest computers.

Fortunately, as previously explained, ‘fully’ realistic simulations (in which lattice parameter values map directly onto those of the real world) are not the goal of mesoscale LB. Indeed, for the present problem, fluid mechanical lore (Batchelor 1967) asserts that for an isolated sphere, all Reynolds numbers less than unity are practically equivalent. This is not quite enough to show that LB is accurate, since the time discretization error in LB enters formally at the same order as the inertial terms that *Re* describes, but with a different prefactor. However, Cates *et al*. (2004) confirm that any *Re*≤0.02 is clearly negligible and *Re*≤0.1 probably adequate, bearing in mind typical systematic errors of a few percent (Kendon *et al*. 2001) from each of several other sources. In summary, however small *Re* is in a specified real-world problem, we should not waste time making it smaller than we need to, to gain accurate simulation results for that problem.

By the same token, for many mesoscale simulation purposes one expects that any Bond number less than some small value is equivalent to any other. Of course, this is not true if one looks at the detailed deformation of the fluid interface close to a particle suspended there against gravity, but for multi-colloid simulations that is unlikely to be well resolved in LB anyway. Indeed, in such simulations it is normal to use colloidal radii of about three lattice units or even less (Ladd 1994*a*,*b*; Nguyen & Ladd 2002), whereas the full width of the fluid–fluid interface is usually one or two lattice units (Kendon *et al*. 2001). Although the interfacial tension in LB cannot be raised to very large values without risk of rapid instability, it is not difficult to achieve Bond numbers of order 3×10^{−4} (as used in figure 1). The stability issue is controlled mainly by fluid parameters, not colloid size; also the choice of collision operator makes a difference, with better stability from multi-relaxation methods rather than a single relaxation time. Note that, as discussed by Kendon *et al*. (2001), any LB run using our code is liable to become unstable eventually; what matters in practice is that time-scales of interest can be reached well before this happens.

Figure 2 shows the dependence on *Bo* of the mean particle displacement from the zero-gravity equilibrium position, for particles suspended from an interface (as in figure 1, but now with a single particle only) on varying the strength of gravity and/or interfacial tension. (Stability is not a problem in this range of *Bo* for any of the particle sizes studied.) This is compared with a theoretical prediction based on the result of Derjaguin (1946) that should be accurate at low *Bo*. Agreement with the latter is excellent; moreover, bearing in mind other sources of systematic error, so long as *Bo*≤10^{−2} or certainly *Bo*≤10^{−3}, such displacements are close enough to zero for these *Bo* values to serve as an adequate numerical substitute for fully realistic values of order 10^{−8}. (Interestingly, one such source of error involves spurious momentum currents in the plane of the interface; these are small—≤10^{−5} in lattice units—but mean that particles do not come to a complete halt on the interface even, as here, in the absence of thermal noise.) Figure 3 shows a particle at the interface for *Bo*=0.006 and 0.9. In the first case, the interface is effectively flat; in the second, the particle is en-route to detachment. Our preliminary simulations find a critical Bond number for detachment between 0.6 and 0.9, with a small dependence on particle radius *a*; the theoretical value is 3/4. Much higher accuracy would require larger particles, and certainly one would expect to use these in any problem where the local physics of detachment was dominant.

## 3. Computational scaling issues: sheared binary fluids

Computational scaling issues are in one sense purely algorithmic: the concern is how to scale-up the simulation to deal with larger systems sizes and thus reduce discretization errors and finite-size effects. We deal with these questions here, using as an example the important problem of symmetric binary fluids under shear. This problem has been widely studied by simulation (Corberi *et al*. 1998; Cates *et al*. 1999; Wagner & Yeomans 1999; Lamura *et al*. 2002; Harting *et al*. 2004) and also by theory (Doi & Ohta 1991; Cates *et al*. 1999; Bray 2003). But so far definitive LB results have proved elusive, in large part because these require very large system sizes before finite-size effects can be brought under control.

The scientific issue is relatively simple: does steady shear prevent the demixing of binary fluids, and if so, how do the resulting steady-state length-scales for the fluid domains depend on shear rate? ‘Length-scales’ is plural, because theory shows that in principle the characteristic domain sizes measured in the three directions of velocity (*x*), velocity gradient (*y*) and vorticity (*z*) can differ by arbitrarily large factors. Indeed, in one type of approximate theory, in which the shear flow is assumed perfectly laminar at all times (Rapapa & Bray 1999; Bray 2003) all three of these length-scales diverge in time, but with different power laws. Simple scaling arguments, which ignore anisotropy but allow for a fluctuating velocity profile, do not rule out a divergence with time of the (now single) domain scale *L*(*t*) under shear. However, they also allow at least two possible scenarios for its saturation. The argument of Doi & Ohta (1991) proposes that saturation occurs when the capillary number is of order unity; this gives so that the domain scale in steady-state obeys . In contrast, Cates *et al*. (1999) argue that coarsening might cease only when the Reynolds number is of order unity; this gives or .

There are two separate reasons why these arguments still have not been definitively tested by LB (in either two or three dimensions). The first is purely to do with numerical scale-up and involves the need to subdivide the LB lattice to achieve large shear rates. The second concerns finite-size effects: specifically, what shape of simulation box should be used, and how can one know this *in advance* of doing the simulation? We address these issues in turn.

### (a) Multiple Lees–Edwards planes

Until recently the best way to apply steady shear to a binary LB fluid was to create a simulation box with solid walls on two parallel faces (and periodic boundary conditions on the rest). At these solid walls, a specific velocity can be imposed. Unfortunately, this method of applying shear poses a severe obstacle to computational scale-up, as follows. In LB, the local fluid velocity relative to the lattice cannot be allowed to get too large: *v*_{max}=*c*_{s}/10≃0.06 (in lattice units) is a good rule of thumb. (Beyond this, instability and/or severe breakdown of Galilean invariance sets in.) Applying this maximum velocity with opposite signs at the two walls gives a maximum shear rate , with *Λ*_{y} the separation of the walls. This makes it impossible to increase the system size at fixed shear rate; indeed the achievable *Re* increases only linearly with *Λ*, while *Ca* does not increase at all!

A breakthrough for this problem was made by Wagner & Pagonabarraga (2002) who developed Lees–Edwards boundary conditions (LEBC) for LB. In molecular dynamics, LEBC (Lees & Edwards 1972) impose a mean shear rate by applying shifted periodic boundary conditions across planes normal to the shear gradient direction, *y*; particles crossing such a plane are incremented in velocity by an amount . (They are also incremented in *x*-position by a corresponding, time-dependent amount.) The boundary conditions across sample boundaries normal to *x* and *z* are the standard periodic ones. LEBC have the effect of setting up a mean shear rate : but note that the system itself, not the boundary conditions, decide whether this shear is uniformly distributed within the sample. Specifically, LEBC in molecular dynamics does not prevent ‘shear banding’ (i.e. formation of layerwise regions of different ) so long as these respect the imposed *mean* shear rate .

Wagner & Pagonabarraga (2002) not only showed how to extend the LEBC scheme to LB—a highly non-trivial task in view of the underlying lattice—but also showed how to include additional Lees–Edwards planes. These again lie normal to *y*; if there are *n* such planes in total, the velocity increment across each is now . Within molecular dynamics this would be a pointless exercise in coding, but the effect in LB is to refer the physical fluid velocity within each of the *n* slabs to a lattice that is, for a linear flow profile, moving with the velocity at the centre of that slab. This increases by a factor *n* the maximum velocity gradient that can be modelled in LB, circumventing the scaling problem raised above. Crucially, although slabs are co-moving with the *linear* flow profile ) that would arise for a Newtonian fluid, the method makes no direct assumption about the *actual* flow profile . However, if this deviates strongly from a uniform shear rate, the fluid velocity in a certain slab could violate the requirement of *v*<*v*_{max} in the lattice frame. In systems undergoing shear banding, it might therefore be necessary to have an adaptive algorithm controlling the velocities of different slabs. But for the problem of sheared binary fluids of equal viscosity, shear banding is neither expected nor observed in simulation, and the problem does not occur.

The main limitation of this approach then arises from the lattice implementation of the LEBC themselves. The velocity jump across an LE plane is small, and thus does not give an integer displacement of the lattice in each time-step. Hence there is an interpolation required to allow information from the distribution functions to stream across a given LE plane onto the mismatched lattice plane on the far side. Additionally one has to somehow increment the fluid velocity within the information that gets passed. Since the velocity set is discrete, this cannot easily be done at the level of the individual populations *f*_{i} and velocities **c**_{i}. Wagner & Pagonabarraga (2002) simplify matters by applying the required shift only to the equilibrium distribution, that is, they write(3.1)where is a distribution function after passing across the plane, *f*_{i} its original value, and its equilibrium value. (Note that there is no shift for velocities **c**_{i} which do not cut the plane.) A similar shift is applied to the distribution function *g*_{i} that governs the binary fluid order parameter. The errors that arise from using equation (3.1) are, at any given time-step, demonstrably small (Wagner & Pagonabarraga 2002). However, they are not zero and, in combination with the interpolation effects, can build up into artefacts when integrated over long time periods. The latter problem seems to be more severe for the *g*_{i} than for the *f*_{i}: they mainly affect the order parameter sector. Note that this sort of problem is not completely absent in systems where shear is imposed by moving solid walls at the boundaries, but clearly it is of greater concern in the case of scale-up when one wants to add additional boundary planes within the bulk of the sample.

A first-principle approach to this scale-up issue involves redesigning the implementation of the LEBC, so that the streaming of information across LE planes more accurately reflects the Lees–Edwards prescription at the molecular level. This work is ongoing. Meanwhile, we have tried various mitigating measures; below we describe both the artefacts found, and some of those measures.

We noticed two inter-related artefacts associated with the LE planes, one seeming to restrict, and the other to promote, the coarsening of fluid domains. The first occurred when using relatively high values of the velocity jump Δ*v* across the planes (say Δ*v*>0.1) which gave, nevertheless, fluid velocities on either side of the plane that were less than the commonly accepted maximum (*v*<*v*_{max}=*c*_{s}/10≃0.058). In such cases the locations of the LE planes were clearly visible in the order parameter field and it appeared that the fluid was being mixed more vigorously along these planes than elsewhere in the simulation. This effect is illustrated in figure 4*a*. The growth of the domain sizes in such two dimensional simulations reached a saturating steady-state in each direction, but this cannot be considered physically reliable due to the visible artefact of increased fluid–fluid mixing occurring along the LE planes. Broadly similar effects (not shown) were found in three dimensions.

The second, closely related, artefact was more evident when Δ*v* was reduced, so that say Δ*v*<*v*_{max}/10≃0.06. In these cases the growing domains had a tendency to align themselves along the LE planes, not just locally (as already visible in figure 4*a*) but across the full extent of the simulation box. This additional alignment was enough to promote domain wrap-around in the direction of the shear flow. The latter is a finite-size artefact (discussed below) but the interface lock-in onto LE planes greatly increases the damage done by it. Indeed, in some extreme cases, with an even number of LE planes, the fluid reached a steady-state in which each domain wrapped once around the lattice in the direction of the flow to connect with its own tail, with the interfaces aligned closely to the positions of the LE planes. (Using an odd number of LE planes helps to reduce, but does not eliminate, this effect.) The resulting multi-layer steady-state, such as that of figure 4*b*, could be a valid result for a finite sized system under sheared periodic boundary conditions—but it must be treated with great suspicion when, as here, the number and location of the interfaces closely matches those of the LE planes used.

The LB parameters in figure 4*a*,*b* are chosen to emphasize, rather than minimize, the visibility of the lock-in problem. For other parameter values the effects can be much smaller (figure 4*c*) and might easily be overlooked, but their cumulative effect still remains a matter of concern for quantitative work. We noticed that these artefacts became more pronounced as simulations progressed, but always took quite some time to develop.

A first strategy for minimizing the problem was simply to increase the number of planes while simultaneously reducing the values of Δ*v* such that remains constant. This was not satisfactory since, as foreseen by Wagner & Pagonabarraga (2002), the numerical error created by their scheme grows rather than vanishes in this limit. (Each plane slightly degrades the streaming of *f*_{i} information even if Δ*v*=0.) Thus, although the strength of the lock-in on any particular plane was reduced by this method, we could not be confident of the overall accuracy of the simulation.

Secondly, we improved the accuracy of the interpolation of the populations *f*_{i} across the LE plane from linear to cubic spline, but unfortunately this had little or no beneficial effect on the interface lock-in.

Thirdly, in order to reduce the time available for the artefacts to develop, we periodically translated the positions of the LE planes normal to their own direction through a specified number of lattice units (leaving the fluid untranslated). This was done by re-mapping the *f*_{i}, *g*_{i} data, which itself introduces some interpolation errors. (But since the jumps of the planes are infrequent such additional errors are probably minor.) This plane-jumping technique reduced the time available for alignment of fluid interfaces with LE planes, since no such plane remains in the same place for very long. We found that even on very long and thin lattices, for example with (*Λ*_{x}, *Λ*_{y})=(1000, 160), the simulations still suffered from the finite-size effect of fluid domain wrap-around. Of course, this could be the correct physics for the chosen periodic simulation box and may not be an artefact from lock-in; indeed, for lattices as large as (*Λ*_{x}, *Λ*_{y})=(2000, 160) we observed the longest length-scale still to be growing, albeit slowly, right up to our maximum simulation time of *t*=300 000. The plane-jumping technique therefore remains a possibly useful component in overcoming the scale-up problem.

Fourthly, and in a similar vein, tests were run in which a small, constant overall velocity was imposed in the direction normal to the LE planes by initializing the fluid with a non-zero *y* momentum. (This induces a small steady *x* acceleration, which could be removed by a suitable body force.) The hope was that, by drifting the entire fluid configuration sideways across the LE planes at constant velocity, there would again not be long enough for any lock-in between fluid–fluid interfaces and the LE planes to develop. Instead though, these simulations exposed a more general drawback with the LEBC scheme of Wagner & Pagonabarraga (2002): we found that this was not Galilean invariant. Indeed, the observed velocity profile rapidly changed from the constant shear rate expected (and seen, in the absence of transverse momentum) into something more akin to a step function. This has motivated further work on the LEBC issue, which is still ongoing. Preliminary indications are that the Galilean invariance problem can be solved within the broad framework of the Wagner–Pagonabarraga scheme, so long as some specific improvements are made to the streaming across planes of *f*_{i}, *g*_{i} information. In addition, these improvements decrease (but do not wholly eliminate) the errors arising at the LE planes themselves; in combination with the drifting or the plane-hopping method, we hope that this development may soon allow production run work to begin.

As a preliminary indication of what might be achieved once our scale-up problem with LE planes is surmounted, we show in figure 5 results selected from that part of our two-dimensional data that seemed least affected by the above-mentioned lock-in artefacts and/or the resulting finite-size limitations. In the best two-dimensional runs, one can discern near-steady-state length-scales *L*_{x} and *L*_{y}, defined through the order parameter field as follows:(3.2)where the area integral is over the simulation domain, and *Λ*_{y}*m* is the value of the same integral for an infinitely long domain of width *L*_{y}, containing a single interface perpendicular to the *x*-direction. (The expression for *L*_{y} is found by interchanging *x* and *y*; we restrict attention to the symmetric case with equal amounts of the two fluids. The generalization to three dimensions is straightforward but messy.) With these definitions, the right-hand side of equation (3.2) is a measure of the number of interfaces traversed on passing through the domain along the *x*-direction, and *L*_{x} accordingly gives a characteristic distance between these. (Geometrical factors involving the orientation of the interface complicate this interpretation slightly but this need not concern us here.)

The resulting length-scales are only ‘near’ steady-state because there are large temporal fluctuations; and in many cases nucleation into a ‘wrapped’ state, clearly limited by finite-size effects, occurs subsequent to the period of near-steady flow during which the *L*s are defined. We have plotted in figure 5 the resulting dimensionless lengths against the inverse of the dimensionless shear rate , where the scaling parameters *T*_{0} and *L*_{0} are those of the coarsening problem as defined previously. We make no attempt at power-law fits, but draw attention to the extremely wide range of scaled lengths and times spanned by this plot. Our data reduction method follows that of Kendon *et al*. (2001) although we have only two data points per run rather than a segment of the coarsening curve *l*(*t*/*T*_{0}) as considered there for the unsheared case.

### (b) Aspect ratio scaling

As mentioned previously, there are theories that predict three growing length-scales with divergent ratios for binary fluids under shear. That is, they predict (modulo some logarithmic corrections) *L*_{x}∼*t*^{a}, *L*_{y}∼*t*^{b}, *L*_{z}∼*t*^{c} with *a*, *b*, *c* different exponents; see Rapapa & Bray (1999). Whether these theories are correct, or whether coarsening instead ceases at some finite capillary or Reynolds number (Doi & Ohta 1991; Cates *et al*. 1999), such predictions point to a second scale-up problem distinct from the limit on lattice-frame fluid velocities addressed above by the LEBC approach.

Specifically, the shape of the simulation box required to efficiently minimize finite-size effects in a system with three independently growing length-scales is highly anisotropic. Figure 6 shows a typical example of a three-dimensional simulation on an 800×80×80 lattice. (See also the work of Harting *et al*. 2004.) The domains wrap-around the periodic boundary condition (so that a single domain connects with its own periodic image, creating in effect an infinite correlation length in that direction) first in the vorticity direction; this occurs about 1/3 of the way through the run, at around 8000 time-steps. At around 20 000 time-steps wrap-around occurs along the velocity direction; this run was then terminated at 25 000 time-steps.

Clearly a non-square cross section is required, but the optimal anisotropy for minimizing finite-size effects depends on both shear-rate and the simulation size. Thus, doubling *Λ*_{y} at fixed requires for best efficiency that *Λ*_{x} and *Λ*_{z} are increased by some factor other than two. Moreover, unless a theory like that of Rapapa & Bray (1999) can be relied upon, it is impossible to know in advance what these scale-up factors should be. Since our ability to get meaningful results is critically dependent on controlling finite-size effects, these problems (combined with the LEBC scale-up issue described above) have so far prevented us from gaining a clear quantitative insight into binary demixing under shear to match that of Kendon *et al*. (2001) for the unsheared case in three dimensions.

Meanwhile, figure 7 shows typical *L*(*t*) curves for preproduction runs in two-dimensional and three-dimensional. The curves are shown both as raw data and after smoothing. In the smoothed curves, a plateau can be seen at intermediate times in two-dimensional, followed by a finite-size driven increase of the lengths as domain wrap-around occurs. (Note that, with our definitions, a wrapped configuration has *L*≫*Λ* in the direction of the wrap.) Hence the identifications of near steady-state length-scales, as presented in figure 5, are relatively unambiguous in two dimensions. In three-dimensional there is a hint (only) of a similar plateau in the smoothed curve. We expect that, if saturation of *L* is to be seen in three dimensions, it will be by growth of this plateau as one moves to larger system sizes. But already, the development of disparate length-scales in the three directions is clear from this figure.

### (c) Resource usage; other scalability issues

The sheared binary fluid simulation of figure 6 (80×80×400, run for 25 000 time-steps) took roughly 24 h using MPI parallelization on a SUN E15000 machine using 32 processors. This represents a typical production run on University-owned shared-memory hardware (whether parallelized with MPI or OpenMP). Run-times are increased, though not drastically, in problems where dilute colloidal particles are introduced; so far, colloids are implemented only in the OpenMP version of our codes. If colloids become concentrated enough to form clusters in close hydrodynamic (lubrication) contact, scaling becomes poor because the velocity-dependent lubrication force (whose normal component we add explicitly at short distances, following Nguyen & Ladd 2002) requires a fully implicit update of colloid velocities. Good scaling can be restored by adding a short-range thermodynamic repulsion which prevents the surface-to-surface separation between colloids getting too small. We do this in the run of figure 1 where particles end up closely packed on the interface; the difference between the effective thermodynamic and hydrodynamic radii is about 0.1 lattice unit (see Cates *et al*. 2004). Note that the occurrence or otherwise of hydrodynamic clusters depends on all aspects of the physics, not just the global volume fraction. For example, in figure 1 the initial volume fraction is less than 1% but it is obviously much higher, locally, by the end. We have not attempted to simulate dense colloidal suspensions with our codes as yet, but regularly do run at 20% volume fraction using the additional repulsion just detailed. Further information on scaling with shared memory architectures is given by Stratford *et al*. (2004); for distributed memory, the relevant discussion is in Desplat *et al*. (2001). Imposition of shear has only limited effects on scalability in either case.

## 4. Mid-run alteration of sample geometry

One tool that would make life significantly easier for such simulations would be a way of altering the size and shape of a simulation mid-run, for example if incipient wrap-around of fluid domains across (say) the vorticity direction were detected, as in figure 6. Ultimately one would like to do this as a real-time steering operation using tools developed by the RealityGrid project; but first one has to work out how to accomplish it at all. This is possible in principle, by replicating the data of a simulation on a *Λ*_{x}×*Λ*_{y}×*Λ*_{z} lattice to create one of size *Λ*_{x}×*Λ*_{y}×2*Λ*_{z}. However, since the LB algorithm without noise is deterministic, nothing would be gained: the new simulation would simply calculate everything twice and deliver results with the same periodicity as before. On the other hand, the same is not true in problems where thermal or other noise (perhaps amplified by the chaotic nature of the demixing process) is dominant. In such cases the replicated system would quickly develop different behaviour in its two halves, and start to genuinely behave like a system of twice the original size.

### (a) Example: arrest of Rayleigh–Plateau instability by adsorbed particles

To exemplify the principle of mid-run alteration of sample geometry, we apply the idea here to a specific problem. The question is as follows: if a cylindrical interface between two fluids is coated with a random arrangement of colloidal particles with neutral wetting characteristics (see Cates *et al*. 2004 and figure 1), does this arrest the interfacial dynamics?

A good test of this is to see whether the Rayleigh–Plateau instability, whereby a long cylinder of fluid breaks into droplets under the influence of interfacial tension, is suppressed by the presence of the adsorbed particles (figure 8). A technical problem then lies in creating a realistic initial state. It would be possible to simply place the particles on the interface by hand, but there is a danger of creating a state that is not the outcome of any physically realisable process (for example, it might have an unrealistically high particle density). Thus it seems preferable to mimic the dynamics of figure 1, and allow the colloidal particles to coat the surface under the influence of a body force acting on them. A reasonable choice, for a cylinder, is a body force that is radially directed towards the cylinder surface and linearly proportional to the distance from this surface.

However, to allow this to happen, the cylinder of fluid has first to be stabilized *against* the Rayleigh–Plateau effect for long enough that the particles can arrive at its surface and coat this uniformly, before departures arise from a cylindrical shape. This can be done by mid-run intervention as follows. It is known that the instability has a minimum wavelength *Λ*_{c}=2*R* and a most unstable wavelength *Λ*_{m}=9.02*R*, where *R* is the cylinder radius (Faber 1995). Therefore, by imposing periodic boundary conditions with period *Λ*<*Λ*_{c} along the cylinder axis, the instability is switched off. (In practice, *Λ*<2*Λ*_{c} is sufficient since growth rates are very slow away from *Λ*_{m}.) In such a simulation, it is then a simple matter to deliver a coating of particles to the interface in the manner just outlined. This coated cylinder can then be replicated as described above, by gluing two (or more) identical copies side-by-side so that the length of the simulation box is now much greater than *Λ*_{c}. As mentioned previously, noiseless evolution of this state would retain periodicity; but the instability, if present, ensures that any small-amplitude extrinsic perturbations with wavelength above *Λ*_{c} gets amplified. Our simulations are done with thermal noise, so in principle one could wait for the required perturbation to arise spontaneously but in practice this is slow. (This fact suggests that we might get away with using a long cylinder initially rather than our replication method; but the perturbation caused by arriving particles is significantly larger than that of thermal noise so that this is not automatic.) Instead we can add a one-off perturbation to the replicated system by hand, and watch the result. We carried this out for an initial lattice of 32^{3} with cylinder radius *R*=8 in lattice units; the initial particle volume fraction was 20%. About two-thirds of these particles reached the surface of the cylinder within 50 000 time-steps. The system was then replicated (doubled in length) and initialized with a small amplitude perturbation (*ca*. 1 lattice unit in the radius of the cylinder) at wavelength 64=0.87*Λ*_{c}. This was then run on for 60 000 time-steps, during which there was little growth of the instability. Initializing to the same state without particles led instead to pinchoff after 56 000 time-steps. Figure 8 shows equal-time snapshots for the system with and without particles. On running the particle-laden simulation further pinchoff finally does occur; we are currently investigating the effect of particle density and size on this effect to see if arrest is complete for a high density of well-resolved (i.e. large) colloidal particles.

## 5. Conclusions

In this paper, we have considered various physical and computational scaling issues that arise in the LB simulation of binary fluid mixtures, and illustrated these with examples from our recent work on sheared binary fluids and colloidal particles in binary solvents. To make good use of mesoscale simulation methods one needs to appreciate the hierarchy of length- and time-scales that are present in the physical system of interest, and try to achieve the same hierarchy within the simulation. It is, however, neither necessary nor generally practical to fully resolve all the time- and length-scale separations present in this hierarchy, which in the real world can span many decades. Instead, one should pay attention to the major physical scales and at least ensure that these occur in the right order in the simulation.

Computational scale-up of LB work is hampered by the various special measures that are needed to reconcile large-scale shear (which requires local Galilean invariance in the fluid frame) with the existence of an underlying computational lattice. While progress has been made, this issue is not yet fully understood and we hope to report further on it in future work. For systems under shear, strongly anisotropic simulation domains are also needed; in some cases it would be very helpful if the size and aspect ratio of these could be changed interactively. It is not clear yet whether mid-run intervention to create larger simulation boxes by replicating the data will ever be useful in that context; this seems more promising for noise-dominated and/or chaotic systems than any exhibiting simpler dynamics. Such interventions can however already be useful, for example in the creation of a complicated initial condition that would otherwise lie on the wrong side of an instability boundary.

## Acknowledgments

This work was funded in part by EPSRC grants GR/R67699 (Reality Grid) and GR/S10377. We thank Patrick Warren for useful discussions.

## Footnotes

One contribution of 27 to a Theme ‘Scientific Grid computing’.

- © 2005 The Royal Society