{'fraction_non_alphanumeric': 0.023848895748355567, 'fraction_numerical': 0.00598978429427112, 'mean_word_length': 4.673129038982697, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Lattice Boltzmann for Binary Fluids with Suspended Colloids

A new description of the binary fluid problem via the lattice Boltzmann method is presented which highlights the use of the moments in constructing two equilibrium distribution functions. This offers a number of benefits, including better isotropy, and a more natural route to the inclusion of multiple relaxation times for the binary fluid problem. In addition, the implementation of solid colloidal particles suspended in the binary mixture is addressed, which extends the solid-fluid boundary conditions for mass and momentum to include a single conserved compositional order parameter. A number of simple benchmark problems involving a single particle at or near a fluid-fluid interface are undertaken and show good agreement with available theoretical or numerical results.


Introduction
The lattice Boltzmann equation (LBE) offers a very attractive way to study complex fluid flows governed by the Navier-Stokes equations (for a recent review see, e.g., [1]). This is particularly the case for flows involving irregular and/or changing geometry, e.g., flows in porous media, multiphase systems, or multicomponent fluid mixtures. In the case of fluid mixtures where flows can have highly convoluted evolving structure, the LBE approach avoids the need to track the time evolution of the interfaces between two or more components -pinch-offs and topological reorganisations occur naturally on the lattice. The situation becomes more complex still if solid particles, or colloids, are suspended in the flow. Such colloidal suspensions exhibit a range of interesting properties, and are of considerable interest in many technological and industrial applications (e.g., the food, cosmetic, and pharmaceutical industries).
A number of LBE approaches to multiphase and multicomponent fluids have been advanced [2,3], which adopt different strategies to capture the important thermodynamic interactions which give rise to surface tension between the phases. While numerical LBE results for single phase flow are very robust, multiphase flow is more problematic: if the thermodynamics is not consistent, spurious unphysical currents can be generated in interfacial regions owing to lack of detailed balance. While a complete, thermodynamically consistent, LBE approach to the binary fluid problem is still an open issue, the current methods can be used to provide useful results for many problems of interest, provided the errors incurred are monitored carefully, e.g., for phase separation dynamics [4,5] and droplet motion and break-up in shear flow [6]. An alternative approach [7] is to use the LBE in the momentum sector, but employ a standard finite difference technique for the evolution of a scalar order parameter describing the composition.
Solid-fluid boundary conditions for stationary objects of complex geometry are readily implemented within the LBE. In addition, a general method for the inclusion of moving solid objects in a single phase is available [8] which allows the representation of, e.g., moving spheres [9] and ellipsoids [10]. These objects are generally resolved on the fluid lattice (they are at least several lattice spacings in size) without the need for complicated curved boundary treatments. Long range hydrodynamic interactions between particles are captured by the LBE and, even when particles are separated by distances small compared to the lattice size, potentially important short-range lubrication interactions can be included for simple objects such as spheres by adding back any unresolved contribution using an analytical expression [11,12].
In this work, a number of these different strands are combined and extended to address a new problem: that of colloids in a binary solvent. First, a binary fluid LBE following [2] is adopted in which the equilibrium distributions are cast in a different, and more intuitive form. This leads to a number of improvements in the behaviour of fluid-only problems. Second, solid-fluid boundary conditions are extended to include colloidal particles. As one important advantage of the LB approach is the relative ease with which solid-fluid boundary conditions are included, it is useful to retain LB in both the momentum and thermodynamic sectors. The boundary conditions for the two sectors can then be implemented in a consistent fashion.
The paper is organised as follows. The following section gives an overview of the basic LBE algorithm used and the extension of solid-fluid boundary conditions for moving particles to the order parameter sector. Section 3 presents a number of results for simple test problems for the fluid only, and including a single particle.

The Lattice Boltzmann Approach
In this section, the LB approach to the binary fluid problem is described. As this differs from previous approaches [2,5], so the formalism is first set out for the lattice-Boltzmann equation (LBE) for the single phase in a way which makes clear the extension to the lattice kinetic equation (LKE) describing the second phase.

LBE for Navier-Stokes
The LBE is commonly used to solve the conservation laws for mass and momentum describing the hydrodynamics of an isothermal fluid and These equations describe the dynamics of the mass density ρ and the momentum density g = ρv for fluid velocity v. The momentum density can be written as the mass flux g α and the momentum flux as Π αβ so that for a Newtonian fluid where Greek subscripts denote Cartesian components. Together, these equations lead to the Navier-Stokes equations for isothermal flow, namely where p is the pressure, while η is the shear viscosity and ζ is the bulk viscosity.
The Lattice Boltzmann equation is derived from the (continuous time) discrete velocity equation where the right-hand side corresponds to a linearised collision operator. The discrete velocity set {c i } are the nodes of a Gauss-Hermite quadrature [3,14] which ensure that the conserved moments of the distribution function, f i , have the same values as in the continuum. For the isothermal LBE, these are the mass and momentum densities The Boltzmann kinetic description is restricted to a dilute gas with an ideal equation of state p = nk B T = ρc 2 s , where n is the number density and c s is the isothermal speed of sound. It is convenient to subtract the trivial kinetic contribution to the pressure from the momentum flux tensor to define a deviatoric momentum flux where Π αβ = i f i c iα c iβ is the momentum flux and Q iαβ = c iα c iβ − c 2 s δ αβ is referred to as the kinetic projector [1].
In the most commonly used DdQn models, with n velocities (or quadrature nodes) in d dimensions, the equilibrium distribution functions are given by where the w i are weights defining the quadrature and repeated Greek indices are understood to be summed over. This form is obtained by a truncation of the Hermite polynomial expansion of the Maxwell-Boltzmann distribution at second order [14]. Finally, the BGK collision matrix L ij must satisfy the constraints imposed by the conservation laws and rotational symmetry. In the single relaxation time approximation L ij = δ ij /τ , and the shear and bulk viscosities are related to the relaxation time τ by η = ρc 2 s τ and ζ = (2/d)ρc 2 s τ , respectively. For multiple relaxation time models, the above relations remain valid if τ is replaced by separate relaxation time for the shear and bulk viscous stress.
The final LBGK equation is obtained by a further discretisation of the continuous time equation. A second-order characteristic-based method [15] can be followed to obtain This does not introduce lattice artefacts at second order so that the definition of the viscosity remains as in the continuum, unlike that in a first order scheme, where the lattice error is absorbed into the viscosity to give η = ρc 2 s (τ − 1/2).

LKE for Cahn-Hilliard
The starting point for the binary fluid mixture in this work is the (Landau-Ginzburg-Wilson) free energy functional which describes the total energy of a system of fixed volume as a functional F [φ] of a single compositional order parameter φ(r, t). The order parameter measures the ratio of the number density of particles of the two components of the mixture n 1 and n 2 so that φ = (n 1 − n 2 )/(n 1 + n 2 ). The choice of the free energy uniquely defines the physical quantities of interest such as the fluid-fluid interfacial tension σ and the interfacial width ξ 0 . The chemical potential is obtained from the free energy functional via so that the thermodynamic force density acting on the fluid is then −φ∇µ.
The equation of motion for the order parameter is the Cahn-Hilliard equation This is a conservation law involving the single conserved quantity φ, the flux of which φv − M∇µ is made up of an advective component related to the fluid velocity v and a diffusive component related to the gradient of the chemical potential by an order parameter mobility M. The mobility is assumed to be constant and independent of φ. A kinetic relaxation scheme for the Cahn-Hilliard equation is obtained by introducing distributions g i obeying the kinetic relaxation equation The discrete velocities c i are the same as for the LBE, and physical quantities are again related to moments of the distribution via φ = i g i and j α = i g i c iα . In choosing the form of the equilibrium distribution, constraints similar to those of [2] are employed, namely While the physical interpretation of the constraints on the zeroth and first moments is clear, that for the second moment is somewhat less so. There are two contributions, the first of which ensures a diffusive contribution to the evolution of the order parameter related to the chemical potential, while the second represents an advective flux related to the velocity field. It is then possible to write down an equilibrium distribution by analogy with Eq. 9, i.e., While this choice satisfies the above constraints, it is not unique. As it turns out, a slightly modified version is required in practice (see following section).
Further discretisation again provides an LBGK equation with a single relaxation time for the order parameter The order parameter mobility is related to the relaxation time via M = c 2 s τ φ , which again can be adjusted to absorb the discrete lattice correction. Note that the relaxation time is not fixed in this approach, in contrast to previous work.

Implementation
The choice of free energy functional follows that of [5], where This describes a symmetric binary mixture in which the bulk free energy is related to two parameters A and B, while a term in a third parameter κ penalises gradients in the order parameter, i.e., it tends to minimise curvature of the interface. The fluid-fluid interfacial tension is then If the further constraint that that A = −B < 0 is added, then the order parameter lies predominantly on the interval [−1, 1]. The equilibrium interfacial profile is a tanh with characteristic width For the results presented in this work, the thermodynamic force arising from the order parameter sector is introduced via a correction to the equilibrium stress from a chemical pressure tensor Having made these choices, the lattice Boltzmann algorithm of collision followed by propagation can be implemented in the normal way. The postcollision distributions, denoted f ⋆ i , are based on Eq. 9 so that where the density and momentum are unchanged by the collision process so that ρ ⋆ = ρ and g ⋆ α = g α . The deviatoric stress is relaxed with a single relaxation time satisfying η = ρc 2 s τ where S eq αβ = ρu α u β +P αβ . Note that this reprojection of the physical properties to the distribution readily admits the use of multiple relaxation times, e.g., if separate values for the shear and bulk viscosities are required.
For the order parameter, there is only one conserved quantity in the collision φ ⋆ = φ, where again the star refers to post-collision quantities. The order parameter flux is relaxed with single relaxation time satisfying M = c 2 s τ φ so that with j eq α = φv α . Finally, the values of µ and v are known, allowing a reprojection of the physical quantities to the g ⋆ i distributions. If the equilibrium distributions Eq. 17 are used the model is, in general, numerically unstable. A simple finite-difference expansion [13] of the evolution equation (not shown) shows a spurious up-gradient diffusion arises which is the cause of the instability. Instead, it is possible to use where the δ i0 has the effect of moving most of the order parameter into the non-propagating rest distribution g 0 . This choice is found to have very good stability properties, while also satisfying the constraints on the moments of g i .

Colloidal particles
A very general method for the representation of solid particles within the LB approach was first put forward by Ladd [8]. Solid particles (of any shape) can be defined by a boundary surface which intersects a set of vectors {c b ∆t} which connect lattice nodes inside and outside the surface. These are referred to as boundary links. In the original approach, this set of links defined a (spherical) shell with fluid occupying all lattice nodes both inside and outside the "solid" object. While this so-called internal fluid can be criticised as unphysical, it actually exerts significant effect on the dynamics of the particle in a single phase flow only on the time scale that it takes sound waves to cross the particle [17]. However, in a binary fluid mixture, it becomes essential to have truly solid particles to prevent possible unphysical transfer of fluid across particle surfaces. Such unphysical transfer would impact on the net composition of the real fluid (i.e., that outside the particles). This is particularly true for particles which wet one species of fluid preferentially (non-neutral wetting), where a particle might capture one species of fluid and hence unbalance the net composition and/or generate spurious thermodynamic forces in the region of the surface. A number of implementations of truly solid particles have been developed [16,17]; this work extends that of Nguyen and Ladd [11] to the binary fluid problem.
A schematic diagram showing the distribution of links for a section of a spherical particle is shown in Figure 1. The centre of a sphere of radius a at r 0 defines the position of the links, and the particle moves smoothly across the lattice with linear velocity U and angular velocity Ω. Boundary nodes are defined to be half way along the links, which the set of vectors joining the centre of the sphere to the boundary nodes denoted {r b }. For a single integral r b Fig. 1. A schematic showing the distribution of boundary links for a section of a spherical particle. The nominal radius of the particle is a, while the vector joining the centre to the boundary nodes is r b . If the particle moves relative to the lattice, fluid nodes may be exposed (right), in which case fluid with appropriate properties must be added back. The actual lattice used in calculations is D3Q15.

particle (having a full complement of links)
where w c b are again the quadrature weights appropriate for the boundary links. These results will be useful in the next section.

Boundary conditions and particle dynamics
The boundary condition developed by Ladd is usually referred to as bounceback on links (BBL). For a moving particle, this takes the incoming distribution at the solid-fluid boundary, f b , and reflects it back along the incoming direction where c b ′ = −c b . The correction is related to the local solid body velocity at the boundary node and has the effect of transferring mass from the leading edge of the particle to the trailing edge. If the approximation is made that ρ ≈ ρ 0 , i.e., that the density is approximately constant around the particle, then Equations 27-30 can be combined to show that the BBL applied around the whole particle conserves mass.
Conservation of linear and angular momentum entails computing the momentum transfer from the fluid to the particle as a result of BBL, and updating the particle linear and angular velocity appropriately. For this, the fully implicit method described by Nguyen and Ladd [11] is used.
The BBL condition for density is extended to the order parameter [18] so that, at the surface, the distributions g b are again reflected with a correction However, the order parameter φ around the particle cannot, in general, be approximated by a constant as is the case for density. This means that the order parameter is not conserved at BBL by an amount for a given particle. To ensure the composition of the fluid does not drift over time, the deficit or excess δφ is added back as a correction at the following time step: The size of the additional term is generally small compared with the term in u b .c b , so this represents a negligible perturbation to the motion.
Perhaps more serious are the corrections that arise when the particle changes shape as it moves across the lattice. When the particle motion exposes a fluid lattice node or recovers a solid node, fluid with the appropriate properties must be added or removed. For density, this gives rise to a correction to the BBL Eq. 29 that ensures the mean fluid density ρ 0 does not drift [11]. A similar correction can be made in the order parameter sector which ensures that the mean fluid composition φ 0 does not change. For example, if fluid is added at a newly exposed lattice node, then new properties are determined by a linear interpolation of the distribution functions at adjacent fluid sites.
For order parameter, an extra correction to the BBL at the following step of δφ = (φ − φ 0 ) is required to maintain mean fluid composition, φ being the order parameter added or removed. These corrections lead to small unphysical fluctuations in the order parameter near a moving colloid surface. At present, these are accepted as the cost of preventing drift in the fluid composition.
Finally, the gradient in the order parameter field is required to compute P αβ from Eq. 22. Near the solid particles, these gradients are computed following [18] where the order parameter is extrapolated along links to the boundary nodes. In this work only neutral wetting is considered where the contact angle between the solid and the fluid-fluid interface is 90 degrees. It should be noted that in imparting the thermodynamic force to the fluid via the stress P αβ , there is no direct thermodynamic force on the colloid. The order parameter only then affects the particle motion indirectly via the fluid velocity fluid.

Particles close to contact
The BBL for the density distribution functions f i allows the net hydrodynamic force and torque on a particle to be computed, from which the particle velocity and position can be updated in turn. For an isolated particle this is straightforward. However, if two or more moving particles are close enough on the lattice scale that there are no fluid sites in the interstice, the full lubrication force between the particles will not be resolved. In this case, it is possible to add back the unresolved part from the analytical expressions available for pairwise forces between spheres [21] and, after appropriate calibration, recover the correct lubricating behaviour at close approach [11].
The same effect can occur in the order parameter sector when, for example, two particles are close together at a fluid-fluid interface. The particles should experience capillary forces owing to the curvature of the interface in the gap between them, but, again, if no fluid is present this force will not be resolved. Unfortunately, there is no simple way to add back the unresolved force as is done for the hydrodynamic lubrication as the capillary force depends in a complex way on the curvature of the interface in the vicinity of the particles. While these capillary forces may be important for some problems, they do not have the same impact as the lubrication forces in, for example, preventing the particles overlapping. However, the results presenting in the following section only consider a single particle. Calibration of the hydrodynamic radius of different particles for the single fluid is carried out following [11]; these are assumed to remain unchanged for the binary fluid.

Fluid only
As a demonstration that the equilibrium distributions presented in the previous section lead to better behaviour in the fluid sector than those used  Fig. 2. An initially steady spherical interface between two fluids is allowed to relax for some 100s of time steps. The LBE method described in [2] ("SOY") is compared with the new method. While the two methods give similar results at surface tension 0.0042(a and b), anisotropies are exposed using the SOY distributions at higher surface tension 0.055 (c and d).
previously, Fig. 2 shows a comparison of results for a simple test problem. An initially steady spherical fluid droplet of one phase is initialised surrounded by the second phase. The system, here 64 3 lattice sites with periodic boundary conditions, is then allowed to relax for a few hundred time steps. This problem is repeated for two parameter sets taken from Table 3 of [5] corresponding to surface tensions of σ = 0.0042 and σ = 0.055. While there is little visible difference between the results for the lower surface tension, anisotropy in the old distributions is exposed at higher surface tensions which do not infect the new distributions. This improved behaviour is also manifest in improved numerical stability. While by no means unconditionally stable, use of the new distributions allows a wider range of parameter space to be investigated safely.

Single particle at fluid-fluid interface
A spherical particle subject to gravity may be suspended at a fluid-fluid interface by interfacial tension. For a particle which has a density difference ∆ρ with the surrounding fluid, the key dimensionless number in this situation is the Bond number where g is the gravitational acceleration. The Bond number expresses the balance between the downward (or upward for buoyant particles) force and the opposing interfacial tension. If Bo << O(1) the interfacial tension can support the particle with little deformation. As Bo increases, the interface bows downward until the particle can no longer be supported. For a neutrally wetting particle with a contact angle of 90 o between solid and fluid-fluid interface, the critical Bond number is 3/4; particles break away from the interface for higher values.
In the LB model, a single particle of radius a is placed at rest across an initially flat interface and subject to a downward force. The particle is allowed to fall The curve is the theoretical result [19] is valid for Bo << 1. The particle detaches from the interface for Bond numbers higher than the critical value B = 3/4.
until it comes to rest, in which case the equilibrium displacement below the level interface is measured, or it detaches from the interface. The width of the system in the horizontal direction is at least 10a so that there is no significant curvature of the interface at the periodic boundary conditions. In addition, the second interface required in the system is placed at sufficient distance that it does not affect the motion of the particle. Figure 3 shows the normalised displacement h/a as a function of Bond number for a number of different particle sizes. For Bo < 0.01 the interface is essentially rigid with no significant displacement, while the displacement increases to a significant fraction of the particle size for higher Bond numbers. The agreement between the numerical results and the theoretical result is generally excellent. As Bo approaches the critical value of 3/4 the smaller particles drop somewhat below the theoretical curve and can detach from the interface at Bond numbers as low as 0.55. This is likely to be caused by a poor representation of the contact line for smaller discrete particles.

Particle approaching a fluid-fluid interface
As a second test of particle motion in a binary mixture, the drag on a single sphere sedimenting toward a stationary fluid-fluid interface is computed (Fig. 4). First, the mean drag on the sphere is obtained in a single phase LB calculation to provide the calibration 6πηa as for the hydrodynamic radius (the nominal radius used here is a = 2.3). The drag is then recomputed for the same particle in the binary fluid as it sediments vertically toward a hori- zontal fluid-fluid interface. The Bond number is kept small enough that there is little deformation in the interface as the particle approaches. The drag on the sphere is measured as a function of the normal separation h/a; as the exact value depends upon the position of the particle relative to the lattice, the symbols in Fig. 4 represent an average over 20 different particle configurations and binned at intervals in h/a.
The results show that, far from the interface, the drag is very similar to that seen in the single phase. This suggests the redistribution of order parameter as the particle moves does not strongly influence the dynamics. As the particle approaches the interface, the drag increases but then decreases sharply as the leading edge of the particle touches the interface. The effect of differing interfacial width is shown: the sharper interface induces a steeper rise in the drag as the particle approaches. However, in both cases the particle is finally captured by the interface. For comparison, the exact results for a similar problem in which a sphere approaches a flat interface of zero thickness is shown [20]. In this case, the particle never actually reaches the interface.

Closing Remarks
This work has demonstrated the use of a new binary fluid lattice Boltzmann approach applied to the problem of colloidal particles The representation of colloidal particles within an LBE binary fluid allows a large number of interesting physical problems to be investigated. The results presented shown an excellent agreement with exact results for a number of simple test problems, and provides confidence that the approximations made in the approach do not result in undue errors.
A number of potential improvements can be identified. First, the thermodynamic force follows [2] in entering via the equilibrium stress. This can be added directly as an additional force on the fluid in the LBE. The implementation of this additional force in the presence of solid particles requires some care and is addressed elsewhere [22]. Second, the problem of resolving nearcontact capillary forces would be most elegantly addressed via the use of a finite-volume like approach where there is always some fluid retained in the interstice between particles. This approach would also have the benefit of preventing abrupt changes in the discrete shape of a moving particle. However, this would represent a considerable increase in complexity over the current approach.
In the meantime, the work demonstrates the flexibility of the lattice Boltzmann method in addressing problems that are very difficult to address at all, and currently intractable using other methods.