Lattice Boltzmann simulation of deformable fluid-filled bodies: progress and perspectives

With the rapid development of studies involving droplet microfluidics, drug delivery, cell detection, and microparticle synthesis, among others, many scientists have invested significant efforts to model the flow of these fluid-filled bodies. Motivated by the intricate coupling between hydrodynamics and the interactions of fluid-filled bodies, several methods have been developed. The objective of this review is to present a compact foundation of the methods used in the literature in the context of lattice Boltzmann methods. For hydrodynamics, we focus on the lattice-Boltzmann method due to its specific ability to treat time-and spatial-dependent boundary conditions and to incorporate new physical models in a computationally efficient way. We split the existing methods into two groups with regard to the interfacial boundary: fluid-structure and fluid-fluid methods . The fluid-structure methods are characterised by the coupling between fluid dynamics and mechanics of the flowing body, often used in applications involving membranes and similar flexible solid boundaries. We further divide fluid-structure-based methods into two subcategories, those which treat the fluid-structure boundary as a continuum medium and those that treat it as a discrete collection of individual springs and particles. Next, we discuss the fluid-fluid methods, particularly useful for the simulations of fluid-fluid interfaces. We focus on models for immiscible droplets and their interaction in a suspending fluid and describe benchmark tests to validate the models for fluid-filled bodies.


I. INTRODUCTION
The many-body dynamics of deformable objects, passive and active alike, in confined fluid flows is a central theme of modern soft matter research, with many applications in science, engineering, and medicine.From a fundamental perspective, the main challenge relates to the self-consistent interplay between external degrees of freedom, positions and momenta, and the internal ones describing the shape changes of the deformable body.Such an interplay is expected to spawn new dynamic regimes that are simply inaccessible to rigid particles.Likewise, such new dynamic regimes are expected to lead to new states of soft matter, hence potentially new materials.
In this review, we shall focus on the transport of deformable fluid-filled bodies, an important problem with various applications, from oil to pharmaceutical and food-processing industries [1][2][3][4] .Familiar examples include paints [5][6][7] , milk [8][9][10] , and blood 11,12 , where the internal constituents are ink droplets, fat droplets, and cells, respectively, and all deform when subjected to mechanical stresses, such as the ones resulting from the flow of the surrounding fluid.An analytical approach proves challenging due to the strong coupling between the fluid and the deformable particles.This interaction is reciprocal, for the fluid affects the shape of the particle due to hydrodynamic forces and the deformation in turn changes the boundaries of the fluid flow.While still facing numerous challenges, computational methods have become an indispensable alternative tool to study these systems [13][14][15][16] .
Numerical models for the interplay between the macroscopic fluid properties and its internal deformable constituents usually combine computational fluid dynamics (CFD) with computational mechanics [17][18][19] .This is particularly true for the simulation of fluid-filled bodies with solid boundaries.Deformability alters the flow, hence it affects the way energy is dissipated within the flow, resulting in strong rheological and mechanical nonlinearities, as signalled by the crossover from Newtonian to shear-thinning regimes, such as in most emulsions.A fluid-filled body consists of a body, usually spherical, with a fluid core (see Fig. 1).Common examples include droplets and capsules.The latter is a liquid droplet enclosed by a thin membrane.These bodies are often surrounded by a fluid environment and consequently, they exhibit fluid-fluid boundaries (droplets) or fluid-solid boundaries (capsules).
The modelling approaches for fluid-filled bodies can be divided into two groups, based on the boundary between the fluid-filled body and its environment.These are fluid-structure and fluid-fluid methods.The former consists of modelling a fluid-solid boundary, where the fluid and a deformable solid interact with each other.The second consists of modelling In recent years the LBM has become a useful tool in solving the hydrodynamics of soft matter systems 29 .The fluid flow is described by the Navier-Stokes (NS) equations, where u is the fluid velocity, t is time, p is the pressure field, µ is the dynamic viscosity, ρ is the density and g is the body force.LBM however does not solve the NSE directly.Instead, it is a mesoscopic method which recovers the results of the NS equations.As opposed to molecular dynamics, which solves the dynamics of individual molecules/atoms, in LBM a fluid "particle" consists of a group of molecules and consequently is much larger than the individual scale of molecules.Information about the fluid is contained in the lattice nodes and can only move in specific directions reducing the number of degrees of freedom.
Although it originates from Lattice Gas Automata 30 , it is now considered a separate method.Several works have shown the potential of LBM to solve fluid-related problems 13,[31][32][33][34][35] .It is an efficient algorithm particularly popular for Stokes flow, systems with complex solid boundaries 36 , multiphase flow 37 , hydrodynamic dispersion 38 and convection 39 .
A. Lattice BGK model The Boltzmann equation that contains the fluid particle information reads where f = f (x, v,t) is the distribution function for a fluid particle in phase space, having microscopic velocity, v, at time t, at position x, where a is the acceleration and Ω( f ) is the collision operator which describes the interaction between molecules.Different LBM formulations will have different collision operators.The simplest one is the Bhatnagar-Gross-Krook (BGK) operator given by where τ is the relaxation time which drives the relaxation to the equilibrium distribution f eq (x, v,t).The relaxation time is directly related to the viscosity of the fluid since the fluid kinematic viscosity ν is given by where c s is the speed of sound, τ is the single relaxation time and ∆t is the time step.It should be stated, however, that connecting the kinematic viscosity and the speed of sound is shown to be valid for low Mach number (< 0.3) and nearly incompressible flow.
The equilibrium distribution function is related to the macroscopic fluid velocity u and density ρ where θ is the normalised temperature which is usually assumed to be the unity and D is the spatial dimension.Note that ρ, u and θ depend on time and space in general.The above equation can be expanded using the Hermite series to the second order as f eq α (x i ,t) = w α ρ 1 + where the subscript α represents the discrete lattice speed (see Fig. 2) , ξ α is the microscopic velocity vector, w α is the weight of the lattice and, i is the index of lattice sites.f eq α (x i ,t) is the distribution at equilibrium at position x i and time t.Both ρ and u are a function of position x i and time t.The parameters ξ α and w α depend on the lattice chosen for the discretization.The velocity is discretized into structured lattice velocity (1D, 2D or 3D) vectors.For example, the D1Q3 lattice represents a onedimensional chain with three vectors per node.Usually, more lattice vectors mean increased precision but also increased computational effort.Popular lattices include the D2Q9 (see Fig. 2) and D3Q19.Different lattices have different speeds of sound c s and weights w α for each vector 40 .
The final discrete lattice Boltzmann equation (with the BGK operator) is, where F α is the forcing term which includes both external forces and intra/intermolecular fluid forces.There are different force schemes in LBM used to implement these terms.Several articles and reviews have been written about the advantages/disadvantages of each one (see Ref. [41][42][43] ).A popular one, the Guo forcing scheme 44 , is given by where F is the force.There are two main steps in Eq. 7: collision and streaming (see Fig. 2).The information for both steps comes from the nearest neighbours (usually within the first belt) of the respective fluid node.This makes the method very adaptable to parallel computing such as with graphical processing units (GPUs), when compared to other conventional methods.
The macroscopic density and velocity are calculated through In this review, we focus on the LBM to model the hydrodynamics.The LBM is particularly well suited to multiphase flow, particularly drops, bubbles, and emulsions.Some of the commonly quoted strengths and weaknesses of the LBM method are discussed, for example, in Ref. 40,45 .We also focus on the single relaxation time collision operator due to its simplicity, but we emphasise that more sophisticated ones are used to mitigate spurious effects and may offer enhanced stability, like the MRT collision operator 46,47 and the entropic LBM 48,49 .

III. FLUID-STRUCTURE BASED METHODS
At the microscale, there is a plethora of particles which are bound by membranes.Examples include living cells, artificial capsules (e.g.polymerised membranes) and vesicles (e.g.lipid bilayer membranes).An example of biological membranebound particles is red blood cells (RBCs).We highlight that variations of the models presented here have been used for soft objects like microgels 50 and capsules 51 .While these bodies can suffer strong deformations, for the majority of them volume is relatively unchanged.Consequently, researchers have been interested in the transport of soft matter and other complex fluids flowing at low Reynolds numbers.In order to model biological and synthetic membrane-bound particles, one should consider the membrane elastic, bending and viscous properties along with the inner and outer fluid viscosity.Depending on what is being modelled a combination of the aforementioned properties is necessary to accurately capture the physics.Furthermore, the meshing of the membrane is required, particularly for the methods described in this section.In 2D, the mesh consists of line segments and 3D of surfaces.In general, there are two ways to approach this problem.These approaches depend on whether the boundary (in most cases a membrane) is treated as a continuum medium or as a discrete one.We show in Fig. 3 an illustration of the two different approaches.
The first approach is a continuum one where the membrane is treated as a 2D continuum surface in 3D space.One simulates the solid mechanics of the membrane using appropriate constitutive laws that are descriptions of the energy-strain or stress-strain relations of the membrane.It is then possible to model strain-softening membranes, i.e. the membrane stress increases slower than strain, or strain-hardening membranes, where the membrane stress increases faster than strain 52 .Although constitutive laws are rigorous physical descriptions of the membrane, the numerical implementation is complicated.In addition, due to the continuity of the membrane properties, this approach is limited to modelling at length scales where local differences in the membrane are insignificant 53 .Another consequence of applying a constitutive model in the entire membrane domain is that the shape memory might be neglected in some cases 54 .This approach is discussed in subsection III A.
In the second approach, the membrane is treated as a discrete network of springs and particles.This approach provides some advantages namely, its mathematical description is simpler and local differences or thermal fluctuations can be included in contrast to the continuum approach.However, spring constants can depend strongly on the mesh configuration as shown in Ref. 55 .In fact, Ref. 55 shows that discrete spring models have anisotropic mechanical properties, in general.Furthermore, when studying biological membranes, local area incompressibility has to be enforced which is not possible using simple discrete spring models 55 .Additional terms have to be included.This approach is also known to suffer less numerical instability than the continuum approach 53 .This discrete approach is discussed in subsection III B.
The elastic forces are treated differently in each approach.
In the continuum approach, these are solved using continuumbased solvers such as the finite element method (FEM) while in the discrete approach, the motion and interaction of individual springs are treated in a particle-based manner.Each approach has its own advantages and drawbacks which we highlight in the following sections.The bending along with the volume and area constraints are in common with both the continuum and discrete approaches and are discussed in section III C. Additionally, the fluid can be coupled with membrane dynamics using the immersed boundary method (IBM).In subsection III D, we discuss how to couple the fluid mesh with the membrane mesh.This review focuses on LBM for the fluid which has a Cartesian fixed mesh and for the membrane a (triangular) evolving mesh.The IBM has proven to be useful for partitioned simulations, where the flow and the displacement of the structure are calculated separately 56 and where the elastic bodies have negligible mass 57 .In IBM, the flow is solved on an Eulerian grid while the membrane is treated on a Lagrangian grid 58 .The total force acting on the membrane is treated as a source term in the NS equation.This means that artificial forces are spread as body forces around the membrane boundary in order to mimic the effects of the moving membrane.
In the LBM, this is included in the g term of Eq. ( 1), which ensures no slip and no flow penetration boundaries.Because the IBM is often employed for structures with negligible inertia (overdamped regime), the velocity of a membrane node is interpolated from the surrounding fluid while the forces acting on the membrane node are spread in the surrounding fluid nodes.IBM is discussed in more detail in subsection III D.
Once we define the constitutive models and constraints, it is straightforward to derive the nodal forces as will be shown.However, it is nontrivial to implement the force calculations.The vast majority of meshes used in the literature are triangular meshes in both the continuum and discrete approaches.We aim, in the following subsections, to describe approaches to obtain those forces for elastic, immersed bodies.We illustrate different possible shapes with the fluid-structure-based models in Fig. 6.

A. Continuum approach
The membrane can be treated as a continuum sheet with infinitesimal thickness.The total energy of a membrane is where W s is the elastic strain (in-plane) energy, W b is the bending energy, W A is the area penalty energy, and W V is the volume penalty energy.We highlight that W s is described by a continuum function which describes the in-plane energy.It is possible then to calculate elastic in-plane forces from W s using FEM.Moreover, explicit expressions can be used for the calculation of bending and area/volume penalty forces as this does not necessarily require the FEM.In general, to solve the physics of the continuum membrane we need to discretize it using an appropriate mesh (usually triangular).The elementary surfaces in the mesh are connected by nodes which interact

Discrete mechanics approach Continuum mechanics approach
Boundary treated as a continuous element.
Based on solid mechanics.
Diverse strain-softening or hardening models available.
Rigorous but more numerically unstable.
No local differences or heterogenities in the boundary.
More difficult to implement.
Boundary treated as a discrete network of springs and particles.
Spring constants are mesh dependent.
Most models are strain-softening.
Numerically more stable.
Possible to implement local heterogenities (e.g.thermal noise).
Easier to implement.

Deformation Deformation
FIG. 3. Discrete mechanics approach (left) vs continuum mechanics approach (right).On the discrete mechanics approach the dynamics of the membrane is given by a connected network of springs.On the continuum mechanics approach the membrane is treated as a surface i.e. a continuum medium.The advantages and disadvantages of each approach are highlighted.
with the fluid.The force acting on membrane node i at position x i is given by the principle of virtual work [59][60][61][62] , Details about derivations using the equation for each term can be found in Ref. 63 .For simplicity, we usually consider more straightforward expressions for the force instead of the full expressions derived in the next sections.

Elastic membrane models and finite element method
To determine the elastic forces acting on membrane nodes, it is essential to choose an appropriate constitutive equation that accurately captures the membrane's inherent properties.We refer the reader to Ref. [64][65][66] for commonly used models which we also discuss in this section.The Neo-Hookean (NH) model 64 is one of the most straightforward choices among different constitutive equations, in the form where E s is the surface shear elasticity modulus, λ 1 and λ 2 are the principle stretch ratios (see Fig. 4a) and dA is a surface element area.The NH constitutive equation describes the membrane as a material that doesn't change volume when deformed.If the membrane area increases, it gets thinner to keep the volume constant.An alternative form of Eq. 12 is the Zero-Thickness (ZT) shell representation 67 with the strain energy function given as: A widely used model suggested by Skalak 68 ; see e.g.Ref. 18,69 .The membrane is treated as a thin, 2D elastic surface due to its small thickness compared to its surface area.This allows for an effective modelling of the membrane's behaviour under both low stress and strong deformations.In this model, the strain energy density function is defined as: where E s and C are the membrane elastic shear and area dilation moduli.The same is true for other elastic models.Equations ( 12), (13), and ( 14) can also be written as a function of the strain invariants I 1 and I 2 defined for a 2D membrane (with isotropic and homogeneous properties) as We note that some of the methods, such as that proposed by Ref. 70 , allow one to model a more generic capsule whose interface is also a 2D isotropic sheet with no bending energy.The membrane energy is given by where W 0 is a reference value, I 1 and I 2 are the strain invariants and α 1 , α 2 and α 3 are coefficients which control the elasticity and deformation of the capsule.In fact, one can tune these coefficients to obtain other models such as the Skalak or Neo-Kookean models previously mentioned and even the properties of a droplets 71 .
Once the constitutive law is chosen, we can compute the elastic forces using the linear FEM 59 .The elastic forces are calculated at the Lagrangian nodes of the membrane mesh, using a linear FEM, commonly employed in capsule deformation studies 4,18,72 .These forces result from the planar deformation of discrete surface elements deviating from their original configurations.To enable a comparison between the current and reference configurations, each deformed surface element is projected onto a common plane, aligning with its initial configuration, as shown in Fig. 4b.Triangular surface elements are often used to discretize the membrane thus we also use triangular surface elements in this review.Next, the in-plane displacements, denoted as u x and u y , can be calculated at every vertex.The procedure, as described below, is motivated by Ref. 18,59,60,63,[73][74][75] Let us denote a undeformed membrane marker described by the reference coordinates x and y and after deformation by X and Y .The relation between the undeformed and deformed state is given by: As the membrane is modelled using an 2D energy law W , the nodal forces, F s,i x and F s,i y , are derived using the principle of virtual work for a discrete element.This is expressed in terms of virtual displacements, δ u i x and δ u i y , given by: where, as in Fig. 4b, we use the superscripts i = 1, 2, 3 to identify the variables vertices of the triangle.δW s is the first order variation in the strain energy of the discrete surface element.
There are two assumptions one can make for simplicity.Firstly, a discrete element can be assumed to undergo uniform deformation which leads to the stretch ratios being constant on that discrete element.Consequently, the virtual work for an discrete element is given by where A 0 is the original area of the discrete element and W strain energy density.Secondly, one can assume that the strain energy is related to the principal stretch ratios, λ 1 and λ 2 due to the membrane being isotropic and incompressible.This leads to the first variation of the strain energy function, δW , computed using the chain rule resulting in the following equation for virtual work of a discrete element due to displacements, δ u i x and δ u i y : To obtain expressions for the nodal forces, F s,i x and F s,i y , Eqs. ( 18) and ( 20) are substituted in Eq.( 19) and using the fact that displacements, δ u i x and δ u i y , are arbitrary gives: The forces F s,i x and F s,i y are in the local coordinate system, which after computing are then transformed back to the global Cartesian coordinates and added to the node's Lagrangian force term as shown in Ref. 59 .As an example, in a triangular mesh, every node serves as the vertex for 5 or 6 triangles as shown in Fig. 5 where we highlight a few Lagrangian nodes in red.Consequently, the total elastic force at each node is the combined force from each of these 5 or 6 triangles.
The derivatives in Eq. ( 21) are computed using a linear FEM, with linear shape functions, N i , i = 1, 2, 3, where N 1 , is defined as 59,74,75 and N 2 and N 3 can be obtained by cycling the indices from 1 → 2 → 3 → 1. 2A 0 is twice the area of the undeformed element and is computed as follows: In turn, the spatial dependence of the displacement on the element can be expressed as x + N 2 (x, y)u x + N 3 (x, y)u The relation between (x, y) and (X,Y ) (see Eq. ( 17)) is given in matrix notation by dX = F • dx where the matrix of the Cartesian components of the deformation gradient tensor is 59,60 .Lastly, one just needs to compute the principle stretch ratios, λ 1 and λ 2 , from the eigenvalues of the two-dimensional right Cauchy-Green tensor, G = F T F, whose components are given by FIG. 4. a) Example of a deformation of a 2D sheet due to forces F 1 and F 2 .Deformation in the x direction is characterised by the stretch ratio λ 1 or by the displacement u 1 (the same in the y direction).λ 1 and λ 2 are the stretch ratios in two principle directions, defined as the ratio of the final to the initial length.b) Computation of the deformation is done by mapping both the original and the deformed triangles to a common plane.
The derivatives of the in-plane displacements, u and v, in Eq. ( 25) for G 21 = G 12 are computed by differentiating their linear representations given by Eq. ( 24).Finally, the principle stretch ratios are expressed in terms of the components of G as follows:

B. Discrete approach
A different approach to model membrane elasticity is based on molecular dynamics, as initially suggested in Ref. 76 , and subsequently widely used 77,78  flows, the computational cost of molecular dynamics is prohibitively expensive to reach the length and time scales relevant to the motion and deformation of flexible objects.Addressing this regime requires a mesoscale method.One strategy involves using a discrete arrangement of springs and (coarsegrained) particles.As one represents multiple atoms as a single particle, the computational cost is reduced.In addition, this simplification helps in studying larger spatial scales or longer timeframes while sacrificing some fine-scale details.Furthermore, this allows for defining an effective interaction potential between each particle to simulate the mechanical behaviour of the system.By adopting this approach, the number of degrees of freedom is minimised, and the time and length scales are increased.The force acting in each particle can be obtained by taking the derivative of the potential energy.Other forces such as fluid forces and particle-particle repulsion are also taken into account.
While membranes have different constituents, most share common physical properties.The in-plane stretching energy can be modelled through several different potentials.These can range from the simple linear such as the harmonic potential to a complex nonlinear ones such as worm-like-chain (WLC) potential or the finite extensible nonlinear elastic (FENE) potential.This approach has proven effective in making predictions that consistently match the energy density functions used in the continuum model 55,76 .The total energy of a membrane is where W volume is included only for 3D simulations.Once again, the nodal force at node i at position x i derived from the potential energy is given by Full derivations are not shown but can be found in Ref. 79 .

Elastic membrane models
A harmonic stretching potential can be used to model the stretching energy as where k s is the linear spring stretching constant, N S is the number of springs, l j and l j0 are length and the equilibrium length of the j-th spring, respectively.Eq. ( 29) has been used to study diseased RBCs 80 , blood flow in capillaries 81 and has been validated against modern experimental techniques such as atomic force microscopy 82 .Eq. ( 29) is one of the simplest form of in-plane to model in-plane streching energy, however it is restricted in its ability to model nonlinear characteristics.Nonetheless, some attempts have been made to account for an increase in stretching resistance with elongation of the spring by using a nonlinear spring constant k s as k s = k s0 e 2(ι−1) where the bond stretch ratio (ι) is included 83 .The stretch ratio ι is defined as ι = l/l 0 where l and l 0 are length and the equilibrium length of the spring, respectively.It has been shown to achieve good results between experiments and simulations 81 .Similarly, other nonlinear potentials may be used to model the stretching energy.The two most popular nonlinear potentials 77 are the WLC potential and FENE potential.Explicitly, they are defined as where k B is the Boltzmann constant, p is the persistence length, T is the temperature, x = l/l m ∈ (0, 1), l m is the maximum spring extension and k s is the FENE spring constant.Notice that in both models the maximum extension is limited to l m as the force will approach infinity as the spring length approaches l m .
Note that these springs (FENE and WLC) represent attractive potentials, thus they tend to reduce the area (area compression).A repulsive force field should be combined to restrict this reduction.Ref. 77 proposed two potentials, one of which is an inverse power repulsive potential (W POW ) that is based on the spring length i.e it will restrict the length of the spring.This means where the inverse power repulsive potential is given by and k p is the repulsive stiffness.Finally, the forces for Eqs. ( 30) and ( 32) are where x = l/l m and li j = ⃗ l i j /l is the vector of unit length between nodes i and j.

Bending
The resistance to bending and the prevention of buckling is relevant in many membrane-bounded bodies, particularly biological ones.The implementation of bending is not trivial and it is computationally intricate.In fact, it has been the subject of several reviews [84][85][86] .Since there are many ways to implement bending, for simplicity, we focus on the two most common approaches, one for the continuum approach and another for the discrete approach.We refer the reader to Ref. 84,85 for bending models and numerical implementations.In line with the continuum approach, one of the most popular bending formulations is that of Helfrich 87 where κ is the mean curvature and κ G is the Gaussian curvature, c 0 is the spontaneous curvature (curvature for which the bending energy is minimal).It has been shown that the bending force acting on node i can be written as 88,89 where n is the outwards normal vector to the surface and ∆ LB is the Laplace-Beltrami operator.Calculation of Eq. ( 35) and ∆ LB κ can be found in Ref. 84,85,90 .Another popular and simplified expression for the bending energy 86 often used with the discrete approach is, where k b is the spring constant for bending, N b is the number of bending springs, θ k is the angle between two outward surface normals of two neighbouring triangular elements that share the edge k and θ k0 is the equilibrium angle.The corresponding force Eq. ( 36) is given by Details on the calculation of ∂ θ k ∂ x i can be found in Ref. 63,86 .

Area and volume constraints
In order to impose area and volume conservation it may be necessary to include appropriate constraints.This will depend on the problem at hand.In some cases like a capsule in shear flow, the volume change is negligible and it is not necessary to impose volume conservation.The area and volume energy constraints can be enforced by adding the energy terms 63 where k a and k v are the area and volume coefficients that control the strength of the force which regulates the change in total area A and the volume V from the initial area A 0 and volume V 0 .While A 0 and V 0 are input parameters, A and V need to be computed at every time step.From Eqs. (38) we find that the corresponding forces are More details and computation of ∂ A ∂ x i and ∂V ∂ x i may be found in Ref. 63 .While A can be calculated by summing the elemental areas (e.g.triangles), the calculation of volume requires a more sophisticated approach.An efficient way to calculate V is described in Ref. 91 .We can also express Eqs.(38) as 92 , where n is the outwards normal to the surface.These constraints for area and volume are also used in models which employ a discrete approach.Finally, since a node can be the vertex of multiple surface elements, the total force where j runs over each surrounding element.

D. Linking fluid to structure
The field of fluid-structure interaction (FSI) revolves around the intricate dynamics and interdependence between two fundamental components: a fluid and a solid.This dynamic interplay gives rise to a range of challenges and opportunities, offering insights into complex phenomena across various scientific and engineering domains, encompassing areas like paint pigments, polymers, gels, proteins, red blood cells flow in the heart, and more [93][94][95] .The essence of FSI lies in understanding how fluids and solids interact dynamically.When a fluid flows around or interacts with a solid structure, it induces forces and deformations in the solid, while the solid, in turn, influences the behaviour of the fluid.This mutual influence creates a complex coupling that is prevalent in numerous natural and engineered systems.However, they pose challenges due to their nonlinear and multiscale nature.The nonlinear interactions between fluid and solid components, coupled with the need to consider multiple scales of phenomena, make analytical solutions often impossible.leading to the development of numerical methods to simultaneously address fluid and structure dynamics.In these methods, interface boundary conditions are of crucial as they ensure consistency and continuity between the fluid and solid domains.They dictate how forces, velocities, and deformations are transmitted across the interface.When considering interfacial conditions, there are two main approaches.The first is the Arbitrary Lagrangian-Eulerian (ALE) method 93 .In the ALE method, the mesh used to discretize the fluid domain is allowed to move.Unlike the traditional Eulerian approach where the mesh is fixed in space or the Lagrangian approach where the mesh moves with the fluid, the ALE method allows the motion of the mesh.On the interface between fluid and solid, the velocity and stress should be continuous.The fluid mesh follows the shape of the solid mesh, making it easy to set accurate and efficient boundary conditions at the fluid-solid interface.Yet, when solids undergo significant deformations, re-meshing is often needed to prevent mesh tangling which can affect the mesh motion.This process of mesh regeneration can be challenging during simulations.As an alternative the (IBM) was proposed 58 using a non-boundary fitting method.The IBM uses two independent meshes for fluid and solid components.The solid boundaries are immersed in the fluid mesh, and forces from the solid are spread into the fluid.This eliminates the need for the fluid mesh to conform to the solid boundaries.
The mechanics of the solid and fluid flow need to be solved.This can be done in two ways: a partitioned way where the structural mechanics and flow equations are solved separately or a monolithic way where they are solved simultaneously.FSI problems using the LBM often use the IBM as a coupling algorithm.This allows software modularity and different, possibly more efficient algorithms can be used to solve the structural mechanics.This section will discuss the technical details of IBM and its coupling with the LBM and is primarily driven by Ref. 40,58,96,97 .Other coupling techniques such as the stress integration approach 98 and friction coupling approach 29,[99][100][101] are briefly mentioned.The IBM was selected due to its efficiency and popularity.It was first proposed by Peskin to study blood flow in the heart 58,102 using vortex methods.The fluid is solved on a stationary mesh, while the immersed solids are modelled using a flexible moving mesh.This mesh is not tied to the fixed structure of the fluid mesh.Information is exchanged between the fluid and solid mesh through nodal interpolation.Two-way coupling implies that both the fluid and the immersed structures influence each other.The motion and forces of the fluid impact the immersed structures, and, conversely, the motion and forces of the immersed structures affect the fluid through an effect force density and can prevent fluid penetration in the solid.IBM has been used in the simulation of jellyfish 103 , blood flow 18,104,105 and platelet migration 106 .We refer the reader to Ref. 58,107 for a more detail description of IBM and applications.
Firstly, one can start a description of IBM through a Lagrangian perspective.The immersed structure is characterised as a parametric surface, denoted as X(p, q, r,t).This means that the shape and position of the structure are defined by curvilinear coordinates (p, q, r) varying in a 3D space, with the additional temporal parameter t representing time.Curvilinear coordinates are particularly useful for describing complex shapes or deformations that may not align with a Cartesian coordinate system (e.g.fluid domain).In contrast, The fluid domain is described by coordinates which are fixed in space and provide a frame of reference independent of the motion of the fluid i.e.Eulerian coordinates denoted by x.Ref. 58,97 defines the force density g(x,t) exerted by the structure on the fluid as g(x,t) = F (p, q, r,t)δ (x − X(p, q, r,t))dp dq dr.(42)   This force is distributed then used by the momentum equation of the surrounding fluid nodes.δ (x) is the delta function δ (x) δ (y) δ (z) where x, y, z are the Cartesian components of the position vector x.Likewise, the structural velocity v(X(p, q, r,t),t) is updated by interpolating the neighbouring fluid velocities u(x,t) as v(X(p, q, r,t),t) = u(x,t)δ (x − X(p, q, r,t))dx.(43)   Note that Eq. ( 43) represents the continuity of velocity on the fluid-solid boundary.
The idea of IBM coupling is shown in Fig. 7 where the black dots represent fluid nodes, and magenta dots represent the solid structure.How many nodes are used for the interpolation and their weights depend on δ (x).For instance, Fig. 7(a) illustrates that the velocity of the magenta node X will be interpolated from the surrounding four black fluid nodes within the shaded square frame.Likewise, in Fig. 7(b) the structural force will spread to the four surrounding nodes.
x FIG. 7. IBM two-way fluid-solid coupling.a) The solid velocity v(X,t) is interpolated from neighbouring fluid nodes within the shaded square box.The contribution from a fluid node is weighted by the δ (x) function.b) The solid force g(x,t) will spread to the local fluid nodes as a force density.We see that the fluid node x will receive a force contribution from the magenta nodes within the shaded area as the magenta nodes will spread the force to the nearest four black nodes.As x is the nearest fluid node to the magenta nodes in the shaded area, it will receive the force from them.The spread of the force among fluid nodes is given by the δ (x) function.The number of neighbouring nodes to spread the force can vary.

Spatial and temporal discretization
Implementing the IBM requires both spatial and temporal discretization of Eqs. ( 42) and (43), is needed for the numerical implementation of the IBM.Following Ref. 58,96 , we address the discretization in space first, then proceed to discuss the selection of δ (x), and, finally the discretization in time.
The discretization in space of the Eulerian grid, denoted by g x is a uniform and regular grid in x, y, z coordinates, e.g., x = (x j , y j , z j ) ∆x, where ∆x is the spacing, (x j , y j , z j ) are the components of the position in each direction.This discretization also matches the discretization used in the LBM, which leads to ∆x = δ x and ∆t = δt.Equally, the discretization of the Lagrangian grid, denoted by G s , is the set of (p, q, r) of the form (p k δ p, q k δ q, r k δ r), where (p k , q k , r k ) are integers.Ref. 58 recommended that δ s < ∆x 2 , s ∈ {p, q, r} to prevent fluid penetration/leaking.However, Ref. 63 systematically studied mesh ratio between solid and fluid for small deformations, e.g., δ s ∆x , and concluded that a range of approximately (0.5, 1.5) is enough to prevent fluid penetration/leading without notably affecting the results.Nonetheless, if large deformations occurs that the mesh should be refined accordingly to capture these deformations.Resuming, the force spreading equation (42) becomes g(x,t) = ∑ (p,q,r)∈G s F (p, q, r,t)δ ∆ (x − X(p, q, r,t))∆p ∆q ∆r, (44)  where F (p, q, r,t) is the force density from the solid structures.As in Ref. 96 , let us define F = F (p, q, r,t)∆p ∆q ∆r.F can be viewed as the integration of the force density F over the element of volume dv = ∆p ∆q ∆r, which is the force term applied to each node.The Lagrangian nodes of the solid can be indexed as i without losing generality.Equation ( 44) is then simplified to where X i is the position of the i-th Lagrangian node.Eq. ( 43) is simplified to In the LBM, the time and spatial step in LB units are usually assumed unity, e.g., ∆t = 1, ∆x = 1, respectively.Thus, Eq. ( 46) is simplified to Ref. 58 demonstrated that the δ ∆ (x) function must adhere to specific constraints and properties to ensure consistent calculation of mass, force, and torque from both Eulerian and Lagrangian point of view.Here, we assumed that δ ∆ (x) expressed through the scalar function φ (x) where (x, y, z) are the three components of the position vector x.
We do not delve into the specifics of the φ (x) function.Rather, our intention is to highlight the frequently utilised four-point interpolation function, which is, ) Ref. 58 showed that φ (x) can be accurately estimated using a straightforward formula However, Eq. ( 49) is much faster as it involves evaluating polynomial functions rather than the cos(x) function.Given that the delta function is used to transfer quantities between the solid and the fluid, the choice of its regularisation determines the accuracy of the IBM.Certain properties have to be satisfied in a discrete way by the discrete version of the delta function, to ensure the conservation of force, mass and torque, or to prevent void jumps of variables on the grid nodes 108 .Often three or four-grid approximations are sufficient to satisfy the necessary discrete conditions.Linear interpolation is too simplistic of an approximation and does not lead to the conservation of the physical quantities.Studies of the stability of different regularised delta functions and smoothing techniques can be found in Ref. [108][109][110] .
Because the fluid and solid are solved in an alternating fashion, the IBM follows a partitioned approach.To achieve a temporal discretization scheme with second-order accuracy, Ref. 58,97 suggested an integration scheme based midpoint rule 58 .For the sake of simplicity, the solution at time step n will be denoted by a superscript on the variable.We can then consider the solid nodal position X n i at time step n, the intermediate position at time step n + 1 2 is calculated using Using the updated solid nodal position at the n + 1 2 time step, the force acting on the structure can be evaluated as where we take the gradient of the energy function W . Next, the structure force F n+ 1 2 i will be spread out into the fluid through For specific 2D simulations (e.g.sedimentation) an area scaling factor A b can be used in Eq. ( 53) 111 .Using the force density g n+ 1 2 in the LBM forcing scheme (e.g.Eq. ( 8)), we can then stream and collide using Eq. ( 7) and then calculate velocity u n+ 1 2 using Eq. ( 9).Finally, the solid position at time n + 1 is updated as We recall that we use ∆t = 1 and ∆x = 1 in the above equations.

Lubrication corrections
When the gap between particles or between particles and the wall is equal to or less than the resolution of the discretization, the fluid flow can no longer be resolved in this region.In most cases this happens when the gap is less than one lattice spacing.To solve this, there are two alternative solutions.The first one is to use a refined mesh for the entire fluid system (or just near the region in question if possible).This would lead to a more expensive simulation.Depending on the problem a slightly finer mesh might prove to be enough.The second approach is to introduce a lubrication force that mimics the repulsive force of the fluid being squeezed out of the gap.This second approach is computationally more manageable.This idea has been incorporated into the LBM [112][113][114] , using lubrication theory.The lubrication force between two rigid spheres is µ is the dynamic viscosity of the fluid, where r is the radius of the spheres, s is the dimensionless gap s = R/r − 2 where R is the central distance between two spheres.x i j is the position vector difference between sphere i and j, defined as x i j = x i − x j , xi j is the unit vector.u is the sphere's velocity.Equation ( 55) can also be extended to the case where a sphere approaches a stationary wall or object by setting u j = 0.For deformable bodies, a repulsive force between any two nodes of two meshes in close proximity is sufficient 115 .This force is zero for node-to-node distances larger than one lattice constant and increases as 1/r 2 at smaller distances where d i j is the distance between the two nodes i and j, κ int is a constant regulating the strength and g int ji = −g int i j .

Other coupling schemes
Here we discuss only the no-slip boundary conditions between fluid and structure, motivated by Ref. 96 , which consists of imposing velocity and force continuity along the interface where u f and u s are the fluid and solid velocity, ϒ is the fluidsolid boundary and σ i j is the stress tensor with superscript f and s meaning fluid and solid, respectively.n j is the surface normal.In the IBM, the solid velocity is interpolated from the fluid, as seen in Eq. ( 43), while the force is spread into the fluid, as shown in Eq. (42).Alternatively, the inverse method can be used 98,116 .In this approach, we initially apply the fluid stress to the structure, subsequently solve for the structural response, and ultimately enforce the structural velocity as a boundary condition on the fluid.This technique is referred to as the stress integration approach 98,117 .The total force acting on the structure is where dA is the differential area over the interface.The stress tensor σ i j is given by The pressure term p is usually assumed to follow an ideal gas law p = ρc 2 s .Following Ref. 118 , the deviatoric shear stress in the LBM τ i j := ρν (u i, j + u j,i ) can be evaluated as where ω = ∆t τ , τ is the relaxation time for the LBM.α represents the discrete lattice speed, f neq α is the non-equilibrium part of the density distribution defined as f neq α = f α − f eq α with f eq α calculated from Eq. ( 6).ξ α is the discretized velocity vector.
Using the stress computed from the LBM through Eqs.(58)   and ( 59), the force exerted on the structure can be calculated by taking the dot product between the surface normal and the stress tensor, using the second equation in (56).
In soft matter systems, the flows are often at low Reynolds numbers.Therefore, it is appropriate to Stokes friction force to couple the immersed solids and fluids [99][100][101]119 . In Rf. 100 , polymer monomers are treated as points with a friction force where ζ is the friction coefficient, which might differ from the one provided by Einstein's diffusion relation.u s , u f are solid point and local fluid velocities, respectively.Ref. 99 notes that u f can be estimated through interpolation from the nearest neighbouring fluid nodes.The solid will experience a friction force F ζ , among other forces such as elastic, to determine its motion.Notice that in order to conserve momentum, a force of equal magnitude but opposite direction −F ζ is applied to the neighbouring fluid nodes.For example, Ref. 99,100 uses a force density −F ζ /δ x 3 .Tuning might be necessary for an appropriate choice of the friction coefficient to reduce the effect of slip.Lastly, while in this section we focus on the simulation of membrane-bound particles such as capsules using primarily the IBM, recent works have also used the IBM to simulate droplets 71,120,121 .

IV. FLUID-FLUID BASED METHODS
This section gives an overview of widely multicomponent (different fluids) models for studying droplets and emulsions using the LBM as seen in Ref. 35,40,122,123 .These have applications in several fields, such as for example in injet printing 124 .The LBM facilitates the relatively seamless integration of interparticle forces, forming interfaces between components.However, these forces are often simplified to save time in numerical simulations.Various models have been modified and improved to explore different parameters and reduce issues.Despite these efforts, no single model is perfect for all situations.First, we present the multicomponent formulation of each method and then we explain how frustrated coalescence can be achieved.

A. Colour Modelling
Colour modelling was among the early multicomponent models developed as seen in Ref. 125 .This model is an extension of a basic lattice-gas cellular automaton (LGCA) that handles two immiscible components with surface tension.Instead of the usual indistinguishable particles, it uses two distinct particles labelled as red and blue.The model follows rules similar to the standard LGCA, ensuring that collisions maintain the count of red and blue particles.This modification enables the simulation of immiscible fluids and incorporates surface tension effects within the LGCA framework.Yet, additional rules are introduced to promote the clustering of similar colours and generate surface tension effects.Following this, the model was integrated into the LBM in Ref. 126 and later, in Ref. 127 the model underwent changes to accommodate variations in viscosity and density ratios.Another key development was proposed in Ref. 128 , namely, to add a recolouring step which significantly reduces the computational requirements while minimising the spurious currents and removing the lattice pinning effect 129,130 .The model works by introducing the distribution function for each of the fluid components f α,k (x,t), where k is the component index.The discrete LBM equation is given by where Ω α,k is the multistage collision operator expressed as Each has a distinct aim: • α,k is the normal single-phase collision operator for each component; • Ω (2) α,k is the perturbation operator, generating the interfacial tension; • Ω (3) α,k is the recolouring operator responsible for phase separation.
Macroscopic quantities are obtained by calculating the standard moments of f α,k (x,t) : where ρ = Σ k ρ k .Consequently, the velocity u is referred to as the colour-blind velocity.The perturbation step, responsible for creating surface tension at the interface between two fluids i.e. components, is expressed as: where ∇ρ is the colour gradient, A k regulates the intensity of the surface tension and B α are constants unique to the DNQn set.For the D2Q9 velocity set, B α are given by Likewise, for the D3Q19 velocity set, B α are given by where χ is a free parameter: Values for B α and χ are derived in Ref. 131 .Calculation of partial derivatives is necessary for ∇ρ and it has been shown 132 that the following isotropic central difference contributes to enhancing numerical stability and lowering the discretization error, where j runs through the spatial coordinates.Lastly, one just needs to do the recolouring.In Ref. 125 this involves the minimisation of the work, W , which is given by where q is the local or colour flux Nevertheless, this recolouring method may lead to pinning.This was corrected in Ref. 128 by enabling the blending of red and blue components through symmetric colour distribution concerning the colour gradient.The particle distributions after recolouring are given as where B ∈ (0, 1) is used to manage the thickness of the interface, f ′ α = Σ k f ′ α,k where f ′ α,k is the post-collision state after Ω The total zero-velocity equilibrium distribution function f eq α (ρ, u = 0) = ∑ k f α,k (ρ, u = 0) eq is given by where the weights φ α (along with w α ) are lattice dependent 133,134 .Further modification of the recolouring step in Eq. ( 70) to improve spurious velocities, lattice pinning and convergence was done first in Ref. 133 .Furthermore, θ α is the angle between ∇ρ and ξ α .Finally, the streaming step is carried for each colour.
The advantages of the colour model include its ease of extension beyond two components.Moreover, it provides free control over both the surface tension and thus, the interface thickness.Recent advancements have enabled the simulations of fluids with high density 135 and viscosity ratios 136 .One possible limitation, contingent on the modelling needs, is the absence of a direct connection to thermodynamics.Another downside involves the creation of undesired velocities at curved interfaces, a phenomenon commonly observed in various other fluid-fluid models.

Frustrated coalescence
In colour modelling, a repulsive force has been used to prevent the coalescence of immiscible droplets [137][138][139] .This repulsive force term (operating exclusively on the fluids interfaces) can be introduced to account for any near-contact forces acting between droplets that are closer than the lattice spacing where A h [h(x)] is the parameter controlling the strength of this repulsive force, n is a unit normal to the interface and δ I = which confines the repulsive force to the interface.A graphical representation of this force is shown in Fig. 8.
For example in Ref. 139 , h min = 2 and h max = 4 lattice spacing.Despite the possibility of other functional forms, this one is generally adequate to prevent droplet coalescence and, more significantly, to accurately describe the physics at various length scales.The repulsive force leads to which is the NS momentum equation for a multicomponent system, with a surface-localised repulsive term, expressed through the potential function ∇π.Including extra forces to frustrate coalescence can in principle be used for any multicomponent method, as will be discussed in the following sections.

B. Free-Energy Modelling
We introduce a straightforward free-energy LBM for binary liquid mixtures 140,141 .In this method, discrete density distributions f k,α for each k-component model the hydrodynamics and the evolution of a phase field, respectively.This approach utilises a Cahn-Hilliard phase field method to represent a scalar determining the fluid's composition.We will describe a model where both the fluid flow and evolution of the scalar phase field are solved with the LBM.However, it is also possible to use a hybrid scheme such that the fluid flow is solved with LBM while the phase field is solved with finite differences which leads to reduced memory requirement and usually a marginal loss of numerical stability [142][143][144][145] .The scalar not only moves with the flow but also undergoes diffusion to minimise a freeenergy functional.The functional includes two components: a double-well potential leading to phase separation and an energy penalty discouraging composition gradients, giving tension to interfaces.The thickness of the interface and its interfacial tension can be specified independently (in terms of α and κ φ as discussed below Eq. ( 83)) which is an important advantage over the pseudopotential models described in the next section.
The distribution function f k,α evolves according to where τ k specify the relaxation rates.By relating macroscopic values to the density distributions and appropriately choosing the equilibrium distributions f eq k,α , the simulations model the required continuum equations.For the f k,α field, the macroscopic density and momentum are and the equilibrium distribution is where index notation has been used for the last term and summation over repeated xy indices is implied.The coefficients A α , B α ,C α , D α and G α needs to adhere to conservation constraints, yet these constraints do not uniquely define the coefficients.It is possible to use coefficients that reduce spurious currents 146 .For a D2Q9 velocity set, the full list of coefficients is given in Ref. 45 and details for construction on a D3Q19 velocity set are given in Ref. 146 .
For the phase field, the scalar φ specifies the composition of the fluid, and it varies between -1 (continuum phase) and 1 (droplet phase).It is determined from the density distribution by The continuum equation for φ is the Cahn-Hilliard one where M is the mobility of the chemical potential µ.This diffusivity is determined by the relaxation time τ φ and a free parameter P according to M = P τ g − 1 2 , while the chemical potential is determined by the free-energy of the system.The free-energy functional F[φ (x)] is 141 The first term gives an ideal gas equation of state, the second term is a double-well potential that causes phase separation with minima at compositions φ 0 = ± A B , and the third term sets the interfacial tension by associating an energy penalty for changes in φ with an elastic constant κ φ .The parameters A and B specify the shape of the double-well potential.Two symmetrical phases with φ 0 = ±1, are obtained when A = B.The strength of the free energy resulting from concentration gradients is governed by the parameter κ φ .The chemical potential is then given by The one-dimensional steady-state solution for φ between two infinite domains offers crucial insights into the interface, particularly regarding its thickness and excess free energy.The solution is as follows 141 where ℓ φ = 2κ φ /B is the interfacial thickness of the droplets.The free energy model is also used to study active droplets 143,147,148 .

Frustrated coalescence
Literature on frustrated coalescence using the free energy LBM is scarce.We present a method used by Ref. 149 to study the shear thinning of immiscible droplets in a 2D channel driven by flow.It tracks two fields: the phase-field variables representing each droplet's density φ i where i = 1, . . ., N is the number of droplets and the fluid velocity field u.The free energy is given by: where the first term represents the double well potential which ensures droplet stability.It yields two minima for φ i = φ 0 and φ i = 0, which represent the inside and outside region of the i-th droplet, respectively.The droplet deformability properties are determined by their surface tension Γ = 8κ φ α/9, and can FIG.9. Droplets moving in a channel.This was obtained using a free energy LBM and solving Eq. (84).Colour depicts ∑ i φ i which is ≈ 2 for droplets (yellow) and ≈ 0 for the background fluid (black).therefore be tuned by changing the value of κ φ in Eq. ( 83).The parameters κ φ and α also determine the interfacial thickness of the droplets ℓ φ = 2κ φ /α.The third, final term describes a soft repulsion pushing the droplets apart when they overlap, therefore preventing coalescence.The strength of this repulsion is regulated by the value of the positive constant ε.The dynamics of φ i is given by where M is the mobility and is the chemical potential of the i-th droplet.Equations ( 83) and ( 84) are solved using a combination of the LBM (for the fluid flow) and finite differences (for the phase field) similar to Ref. 140,150 .This may be used for a small number of droplets.
As the number of times that Eqs. ( 83) and ( 84) are solved increases the computational cost rapidly as the number of droplets increases.This method is also useful to simulate droplets inside droplets 151 .

C. Pseudopotential Method
The pseudopotential (or Shan-Chen models) first introduced in Ref. 152,153 as a way to simulate multicomponent and multiphase flows.It works by adding inter-particle interactions resulting in the separation of components.Some instances of applications of pseudopotential models include gravity driven droplets in confined environments 154 and transport of H2O in air 155,156 .Different fluids in a multicomponent mixture have different molecular forces.The force acts on pairs of molecules located at x and x ̸ = x.It is also assumed that a higher density ρ(x) of molecules leads to larger forces.The total force is given as an integral over all possible effective interactions as where G(x, x) strength of force between two fluid elements at x and x ′ and ψ(x) is the pseudopotential, which is a function of the fluid density.Assuming nearest-neighbour interactions at the mesoscopic scale, the effective force on the kth component is given by where G kk is a k × k matrix for the molecular interaction strength between fluid components and its sign regulates whether the forces are attractive (negative) or repulsive (positive).Subscript r stands for repulsive.To create immiscible droplets F k r needs to be repulsive.k and k are just two different components of the fluid system.For a multicomponent mixture (no multiphase), the diagonal components of the matrix are zero and G AB = G BA for two components A and B. ψ k is the pseudopotential (or effective density).The pseudopotential can take different forms.A commonly used form is where ρ 0 is a freely chosen reference density, usually unity.ψ k (x,t) = ρ k can also be used and it is an approximation of Eq. ( 87) at low-density.However, at higher densities, this assumption will lead to numerical instabilities 157 .In the absence of attraction between components, the can often be reduced to the physical density for simplicity.Nevertheless, using Eq. ( 87) instead of ψ k (x,t) = ρ k leads to improved numerical stability.The pseudopotential model allows for multicomponent multiphase (MCMP) simulations, but, for simplicity, we do not present the variants of pseudopotential MCMP LBM.In what follows, we focus only on the multicomponent variant as it is the most widely used method for immiscible droplets.In Ref. 152,153 , they proposed a simple forcing scheme to include these interactions.Similar to the body force, it is also possible to incorporate the effective force of the interaction potential in the fluid velocity field.Nevertheless, for a multicomponent system, we introduce u eq which is a combined velocity (different from the physical velocity) of different components given by: where ρ k u k is the kth component of momentum and τ k is the relaxation time of the BGK collision operator for the component k.The barycentric velocity u of the fluid mixture, i.e., the physical velocity and the total density are given by Also, the density and momentum of each component can be calculated using, For a multi-component system, the time evolution of f k,α (x,t) is given by the discrete Boltzmann equation for each component k (similar to Eq. ( 7)), The equilibrium distribution is given by The Guo forcing scheme 44 can be used to implement these forces acting in a fluid as it yields a viscosity-independent surface tension 158 , particularly with the frustrated coalescence mechanism described in the next section.Thus similar to Eq. ( 8), where F k is the total force acting on component k.In this case F k = F r k .Different forcing schemes can be used.Additionally, as shown in the next section, the frustrated coalescence mechanism can be introduced as an extra force term.Although we described a model based on the BGK operator, extensions for MRT are possible 159 .
The pseudopotential model provides isotropy of the surface tension and spontaneous component separation which improves the numerical efficiency of the simulation.The method is successful because of its simplicity and flexibility in largescale simulations of complex fluids.The implementation is relatively simple when compared to the other models reported so far.However, it also comes with some nonphysical features such as a diffusing interface and spurious currents 40 .Some attempts have also been made to derive a thermodynamically consistent pseudopotential model [160][161][162] .Moreover, it is possible to obtain an effective free energy in the pseudopotential models 163,164

Frustrated coalescence
To prevent coalescence a multi-range repulsive force is used which acts between the nearest neighbours and the next nearest neighbours 158,165 .This force is used in addition to the repulsive component separation force F r k .This means that in Eq. ( 93), k is a competing force given by where G k,1 and G k,2 regulate the strength of force.It is a competing mechanism (superscript c) where the first term is an attractive force and the second term is a repulsive force between same components.This multi-range potential uses 2 different lattices L1 and L2.In 2D these are the D2Q9 and D2Q25 165 , respectively while in 3D they are D3Q41 for L1 and D3Q39 for L2 166 .In 2D, it is possible to completely separate the range of the forces i.e. the L1 lattice runs over only on the first belt of neighbours while L2 runs over on the second belt of neighbours.However, in 3D both lattices (D3Q41 and D3Q39) run over to the third belt of neighbours but on different lattice points which is sufficient to prevent droplet coalescence 167 .In the competing force F c k , the attractive force must overcome the repulsive force to stabilise droplets, so we set G k,1 > G k,2 .To obtain the desired repulsive and attractive forces: G k,1 < 0, G k,2 > 0, and G k,k > 0. For simplicity, we assume: The computational cost does not increase with the number or size of the droplets by contrast to the colour gradient.Compared to the free energy model where each droplet is a different component, the pseudopotential model performs computationally better.This provides an efficient algorithm for simulating a large number of droplets 168 .This is because the repulsion between droplets is based on the density of each LB lattice node.The extra cost comes from using a second lattice to calculate the force.However, this model suffers from spurious velocities caused by an imbalance between the forces.Using a higher order lattice for the streaming can help to reduce these velocities 167 .These velocities increase as the viscosity ratio deviates from unity or the surface tension increase.Proper and careful calibration of the G parameters is necessary otherwise the droplet may increase or decrease in size.With this method, it is not possible to vary the surface tension and the interface thickness independently.

V. BENCHMARKS
We now summarise a few simple benchmarks that can be used to test the model implementation in hydrostatic and hydrodynamic conditions.The benchmarks are for one capsule/droplet in shear flow, capsule sedimentation and a test of the Laplace law.We refer the reader to Ref. 158 for benchmarks on multiple droplets.

A. Capsule and droplet in shear flow
Shear flow is a common validation flow to test for both droplets 169,170 and capsules 117,171 (fluid-structure interactions).Here we use capsules as a generic term for membrane-bound particles.This benchmark consists of placing a spherical (circular in 2D) droplet/capsule in shear flow (between two parallel walls with opposite velocities) and allowing it to deform until a steady state regime is reached 170 (see Fig. 11).In what follows, we have separated into two subsections the shear flow of a droplet and a capsule for clarity.FIG.11.The final steady-state shape of droplet/capsule with the major axis, a, the minor axis of deformation, b, and the third axis, c.An additional reference frame (x ′ , y ′ , z ′ ) is defined, with the axes corresponding to the deformed droplet principal axes.The different regions of the droplet have been also highlighted for ease of reference: tips (green), belly (red) and sides (blue).

Droplet
As a droplet deforms in shear flow, the final shape is a result of viscous forces competing with surface tension.The ratio of these forces is the Capillary number, where u w is the wall velocity, H is the distance between walls, R is the radius of the droplet, µ the dynamic viscosity and Γ is the surface tension.Ca is the main parameter that controls the deformation and final shape of the droplet.The final shape is evaluated using the Taylor deformation parameter given by where a and b are the lengths of the major and minor axis, respectively.In the limit of small Capillary numbers (thus small deformations) a steady state is always achieved.In the limit of nearly spherical droplets an analytical solution was first reported in Ref. 172 by solving the momentum and continuity equation using the Lorentz reflection method.Ref. 172 gives the first order solution D SH (SH stands for the initials of the original authors in Ref. 172 ) for small confinements ratio R/H between two parallel walls 172,173 where p is the viscosity ratio between the droplet and surrounding fluid, H is the gap between the walls, R is the droplet radius and parameter C SH which depends on the position of the droplet relative to the walls and this value grows when closer to the walls.Ref. 172 performed analytical calculations for which they found that a droplet placed in the middle of the walls (R/H = 0.5) has C SH = 5.699 while a droplet at R/H = 0.125 has C SH = 193.3.Other values at different droplet positions are also reported.Some values have been tested against experiments as shown in Ref. 173 .The theoretical model given Eq. ( 97) and Ref. 172 has an additional term D Taylor obtained from Taylor small-deformation bulk theory 174,175 and is given by where Ca is the Capillary number and θ is the angle between reference axis x and major axis a (see Fig. 11) and usually assumed to be π/4 radians as in Ref. 173,174 .Taylor bulk theory 174,175 is valid for small deformations and viscosity ratios p close to 1. Furthermore, while this theoretical model predicts that the magnitude of the deformation will scale with the confinement ratio, the shape should remain ellipsoidal.We highlight that this test can be performed in both 2D and 3D for droplets.For small Ca the deformation parameter is the same for a 2D and 3D droplet (same parameters), while for higher Ca results for 2D and 3D droplets start to deviate.An extensive study of 2D vs 3D droplets can be found in Ref. 176 .The methodology presented here is valid for all the fluid-fluid models mentioned in this review.

Capsule
A capsule in shear will also reach a steady state and the main parameter defining the deformation is also the Capillary number now defined as Ca = u w Rµ k s H where k s is the stretching coefficient and has the same physical units as the surface tension Γ.The Taylor deformation parameter is also used to characterise capsules in shear flow.We highlight that the correct constitutive model needs to be used to compare with other literature or analytical results.Results and validation for Hookean constitutive laws can be found in Ref. [177][178][179] .For Skalak models we refer the reader to Ref. 18,180,181 .It is not trivial to obtain analytical results from constitutive models.In most cases, models are validated against other numerical implementations.An analytical result for the Taylor deformation parameter has been found for the neo-Hookean and Skalak laws, Eqs ( 12) and ( 14) and is given by 18,70,182 , where α 1 = 0, α 2 = 2/3, and α 3 = 1/3 are coefficients extracted from the constitutive model.Equation ( 99) is only valid for small deformation in Stokes flow (Re < 1).More details can be found in Ref. 18

B. Capsule sedimentation
The process of capsule sedimentation involves placing a capsule in a static fluid and allowing it to accelerate downward under gravity until it reaches a terminal velocity, where the resultant drag force balances the gravitational load (see Fig. 12).We present a simple validation for a confined 2D capsule using the harmonic potential given in Eq. (29).The capsule should be rigid.The elastic and bending coefficients can be adjusted to achieve a quasi-rigid state.It is clear that simulating sedimentation of a 2D capsule i.e. a cylinder, will encounter Stokes' paradox, which states that there is no drag force on the cylinder in the Stokes regime.However, experiments and approximations have been derived for cylinder sedimentation in the presence of walls.Therefore, this benchmark consists of placing the 2D capsule between 2 vertical parallel walls.The top and bottom boundaries can be walls or periodic but the vertical walls need to be long enough to minimize boundary effects.
According to Ref. 111 , the analytical wall-corrected settling velocity of a 2D particle can be calculated from where D is the diameter of the particle, g is gravitational acceleration, ρ is the density where s and f stand for solid and fluid, respectively, µ is the dynamic viscosity and Λ is the wall correction factor.Subscript c stands for corrected.Notice that Eq. ( 100) is obtained by equating the gravitational force per unit length F g = π 4 gD 2 ρ s − ρ f with the drag force per unit length F d = V c µΛ.Moreover, Eq. ( 100) is valid for low Reynolds numbers.There are several approximations for Λ which can be found in Ref. 111,183 .A commonly cited approximation can be found in Ref. 184,185 where Stokes equation for Newtonian fluids using asymptotic expansion with no-slip boundary conditions on the walls given by Λ = −4π 0.9157 + ln(k) − 1.724k 2 + 1.730k 4 − 2.406k 6 + 4.591k 8 (101) where k = D/H is the dimensionless particle size.D is the diameter of the particle and H is the distance between the two vertical walls i.e., the width of the channel.Ref. 184,185 applied an integral transform which produces a convergent series and allows the boundary conditions to be applied simultaneously.It is generally known that Eq. ( 101) is valid for k ∈ [0, 0.5].Usually, Λ is only a function of k for low or very high Reynolds numbers Re.For intermediate values, it is also a function of the Reynolds number.The results for a 2D sedimenting capsule are shown in Fig. 12 and are in agreement up to k ≈ 0.4 in the range k ∈ [0, 0.5].Note that for a 3D spherical capsule, Stokes' law can be used for validation.

C. Droplet Laplace test
The Laplace test consists of placing a spherical droplet in the centre of a fluid domain and verifying the Laplace law.Laplace law dictates that when the system reaches equilibrium, the pressure difference across the droplet interface ∆P is related to the interfacial tension where R is the radius of the droplet.In 2D the Laplace pressure is ∆P = Γ R .One should measure the pressure at the droplet centre and then as far away from the centre as possible to prevent undesired effects from affecting the measurements such as spurious velocities.The domain needs to be large enough to prevent the effect of periodic copies (periodic boundaries) or walls (solid boundaries).A linear relationship should be obtained and the slope is the surface tension.As an example, we plot results for a droplet using the pseudopotential method depicted in Fig. 13.Additionally, it is possible to compare the surface tension value to theoretical and numerical values 186 .Moreover, other hydrostatic tests such as contact angle test 186 (for wetting boundary conditions) can be useful to further check the accuracy of the model.

VI. SUMMARY AND OUTLOOK
We reviewed different methods to model fluid-filled objects in flow and addressed some of the advantages, disadvantages, and trade-offs of the different methodologies, as perceived based on our own (necessarily biased) experience and the literature survey.Below, we summarise the main items.
Fluid-fluid models are usually computationally more economical than fluid-structure ones because the interfaces emerge spontaneously from the underlying mesoscale physics and do not require any explicit tracking.The price for this flexibility is a lack of specific control of the mechanical properties of the interfaces, as often required in biological applications.They are also generally affected by spurious currents resulting from lack of sufficient isotropy of the underlying lattice; even though several efficient remedies have been developed over the years, spurious currents remain an issue to be critically monitored.
Both pseudo-potential and free energy must be supplemented with a model for near-contact interactions in order to prevent droplet coalescence, especially under conditions typical of microfluidic experiments.Such near-contact interactions can be implemented very effectively within high-performance LB codes, yielding excellent performances on parallel GPU clusters.For instance, optimized codes running on GPU accelerators can reach GLUPS performance (one billion lattice updates per second) on single GPU machines with a few thousands CUDA cores 187,188 .To convey a concrete idea means that one can simulate one millimeter cube material at micron resolution in space and nanosecond resolution in time, over one millisecond timespan in about 10 6 seconds, namely about two weeks wall clock time.On massively parallel GPU clusters, these figures can be boosted by two orders of magnitude, meaning that the above simulation can be performed in just a few hours of wall-clock time.These numbers open up exciting perspectives for the computational design of new soft materials 24,189 .Fluid-fluid models readily extend to three dimensions by simply enlarging the set of discrete velocities, although at a correspondingly increased computational cost (especially on account of memory access) 190 .
As mentioned above, fluid-structure models are computationally intensive but offer a more accurate control of the mechanical properties of the interfaces, hence of the shape of the resulting droplets and capsules.They also provide a significant latitude in shape-space and do not suffer from significant spurious current effects.Since the structural dynamics is explicitly taken into account, the computational cost raises with the number of droplets and the number of degrees of freedom per droplet/capsule, setting the case for a judicious trade-off between physical fidelity and computational viability.Significant progress over the recent years has led to near-GLUPS performance on multi-million grid sizes with embedded propellers 191 .
Finally, three-dimensional extensions are demanding not only in terms of computational resources, but also in terms of model and programming complexity, such as the switch to finite-element representations of the internal degrees of freedom.
By and large, it appears fair to say that the major progress of LB models over the last decade, both in its fluid-fluid and fluid-structure versions has come to a point of enabling the computational study of new smart materials, both passive and active.Furthermore, the capability of tracking droplets and capsules in realistically complex geometries, such as the human body, may disclose new exciting perspectives in computational microphysiology and medicine 192 .

FIG. 1 .
FIG. 1.Comparison of flow around a fluid-solid (left) and a fluidfluid (right) interface.On the left, the red line represents a fluid-solid interface such as a membrane on a rigid capsule.There is no internal circulation because the velocity at the interface is zero (considering no-slip boundary conditions).On the right, we see internal circulation due to the continuity of the velocity across the fluid-fluid interface.Points A and B indicate stagnation points and the colour gradient of the external fluid represents a pressure gradient which drives the flow.Both images are in the capsule reference frame (fluid velocity relative to the capsule/droplet).The internal fluid in both images is blue.Black arrows represent the fluid flow.

FIG. 2 .
FIG. 2. Two main steps of the LBM: (a) collision and (b) streaming using a D2Q9 lattice.The collision is described by the distribution functions going from one node after the streaming step propagates them to the neighbouring nodes according to the direction.The arrows represent the nine distribution functions.

1 FIG. 5 .FIG. 6 .
FIG. 5. Loop subdivision mesh.The mesh is started from a solid such as an icosahedron (or other initial shapes such as a pyramid) which has 20 sides and 12 nodes.A triangular surface element is divided into 4.There are 3 levels of refinement shown: a) to b), b) to c) and c) to d).For clarification, we highlight an initial triangle 1 which is divided into triangles 1a, 1b, 1c, and 1d.The red points just indicate the nodes of the initial face element.
have been applied.The term cos (θ α ) is conveniently expressed as cos

FIG. 8 .
FIG. 8. Near contact forces representation.Mesoscale modelling of the interactions at the point of contact of two immiscible fluid droplets.The repulsive force is denoted by F rep , and the unit vector n = −∇ρ/|∇ρ| is perpendicular to the fluid interface.The dotted line follows the fluid interface outermost boundary, and the vectors x and y represent the coordinates taken there.The repulsive parameter maximum and minimum values are A hM and A hm , respectively.The colours represent the density field i.e. yellow is the continuum fluid and blue is the dispersed (droplets) fluid.

FIG. 10
FIG. 10.(a) The D3Q41 lattice and (b) the D3Q39 lattice.The D3Q41 is a zero-one-three lattice meaning that it is analogous to the Brillouin zones (also referred to as lattice-belts, or lattice-stencils) covered by ξ α .The D3Q39 is a zero-one-two-three lattice.(c) The D2Q9 lattice for short-range nodes (red).The D2Q25 lattice for short and mid-range nodes (blue).

FIG. 12 .
FIG. 12. (a) Schematic of an initially circular particle of diameter D moving in a channel of length L and width H due to the force of gravity g.The particle is placed at x 0 = 0.5H.(b) Terminal velocity V t against dimensionless size k.The particle has a diameter of 12.5 l.u and the fluid has a relaxation parameter τ = 1.The analytical approximation corresponds to Eq. (100).

FIG. 13 .
FIG.13.Laplace test using the pseudopotential method at equilibrium for G k,1 = −7.9,G k,2 = 4.9 (see Eq. (94)) and τ A = τ B = 1.From the linear fit we obtain the surface tension Γ = 0.00367.The images show the steady state of droplets as we vary the radius.