Appendix D: Lishchuk multi-component algorithm¶

The second set of LBE exercises make use of a multi-component algorithm devised by Sergey Lishchuk [Phys Rev E 67, 036701 (2003)]. This algorithm assumes all fluid components are immiscible (i.e. do not mix together) and there is continuity of velocity and continuity in tangential stress across any fluid-fluid interface, but a stress is applied along the interfacial normal that is proportional to the interfacial curvature. This algorithm has been observed to give relatively small spurious velocities from fluid-fluid interfaces (microcurrents), enabling larger interfacial tensions and flow speeds than some other multi-component algorithms can manage. It is designed to work purely hydrodynamically and thus it cannot tackle separation kinetics properly.

Without loss of generality, we can look at a two-fluid system with a background fluid (coloured blue) and a fluid making up a suspended drop (red). Each fluid has its own set of distribution functions on a lattice: \(B_i (\vec{x}, t)\) for the blue fluid, \(R_i (\vec{x}, t)\) for red. Summing up these coloured distribution functions over all lattice links gives individual fluid densities, \(\rho_B = \sum_i B_i\) and \(\rho_R = \sum_i R_i\). We can also define achromatic distribution functions, \(f_i = B_i + R_i\), and summing these over lattice links gives the total fluid density \(\rho = \sum_i f_i\).

We define a phase field or phase index for these two fluids:

\[\rho^N = \frac{\rho_R - \rho_B}{\rho},\]

which at each lattice point can range from +1 for purely red fluid to -1 for purely blue fluid. Spatial gradients of this quantity in regions where an interface exists between the two fluids provide normals to that interface:

\[\hat{n} = \frac{\nabla \rho^N}{|\nabla \rho^N|}.\]

We can make use of ‘stencils’ to approximate gradients of a property in an LBE calculation with high accuracy, thanks to the isotropy of the lattices in use. For the phase index gradient, a fourth-order accurate gradient can be obtained using:

\[\nabla \rho^N \approx \frac{1}{c_s^2} \sum_i t_i \rho^N (\vec{x} + \hat{c}_i) \hat{c}_i.\]

where \(c_s = \frac{\Delta x}{\sqrt{3} \Delta t}\) is the lattice speed of sound and \(t_i\) is the weighting factor for lattice link \(i\). The local curvature at an interfacial lattice point can be calculated from the surface derivative of the interfacial normal, which for a two-dimensional system is:

\[K = n_x n_y \left( \frac{\partial n_x}{\partial y} + \frac{\partial n_y}{\partial x} \right) - n_x^2 \frac{\partial n_y}{\partial y} - n_y^2 \frac{\partial n_x}{\partial x}.\]

Finally, the interfacial force between the red and blue fluids can be calculated:

\[\vec{F} = \frac{1}{2} \alpha K \nabla \rho^N,\]

where \(\alpha\) is the interfacial tension between the two fluids. This force can be applied during the collision step using any available forcing scheme: the forcing scheme devised by Guo et al. [Phys Rev E 65, 046308 (2002)] is recommended to apply this force correctly as confirmed by Chapman-Enskog analysis.

To conserve velocity continuity between the fluids, the collisions (including applying the interfacial force) must be carried out achromatically using the combined \(f_i\) distribution functions. To recover the distribution functions for each fluid, a segregation step is required. This step uses the post-collisional achromatic distribution functions, \(f_i (\vec{x}, t^{+})\), the interfacial normals and the previously-calculated fluid densities to find post-collisional fluid distribution functions:

\[R_i (\vec{x}, t^{+}) = \frac{\rho_R}{\rho} f_i (\vec{x}, t^{+}) + \beta t_i \frac{\rho_R \rho_B}{\rho} \hat{c}_i \cdot \hat{n},\]

\[B_i (\vec{x}, t^{+}) = \frac{\rho_R}{\rho} f_i (\vec{x}, t^{+}) - \beta t_i \frac{\rho_R \rho_B}{\rho} \hat{c}_i \cdot \hat{n},\]

where \(\beta\) is a segregation parameter used to control the width of the interface between the two fluids and the numerical stability of calculations. This form of segregation normally produces phase index profiles that resemble continuous hyperbolic tangent functions, leading to phase index gradients that do not change overly rapidly between lattice points and lower spurious velocities (microcurrents) generated at fluid-fluid interfaces.

Variants of this algorithm exist for non-spherical or non-circular drops at rest (e.g. fluid-filled vesicles), to avoid calculating interfacial curvatures between fluids, to directly calculate forcing terms for collsions without calculating interfacial terms (the ‘Spencer tensor’ variant), and to calculate phase index gradients and interfacial normals locally at each lattice site without using stencils.