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

Abstract

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–4 Familiar examples include paints,^5–7 milk,^8–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–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–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).

Fig. 1 Comparison of flow around a fluid–solid (left) and a fluid–fluid (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.

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 the interactions between two fluids, usually separated by a fluid–fluid interface. The nature of the boundaries requires computational methods to describe both the fluid inside and outside of the particles. Furthermore, tracking the interface is often required, which can make the overall process computationally expensive. In this review, we focus on cases where the internal and external fluids are both Newtonian. We point out that for industrial applications, capsules can have solid cores (even multiple small solid cores) enclosed by a membrane; however this will not be addressed in this review.²⁰

Computational methods can often illuminate regions of parameter space not easily accessible by experiments or analytical methods. Several methods can be used to model the fluid flow. Here, we focus on the lattice Boltzmann method (LBM), as it provides a versatile way to introduce fluid–fluid interactions and to couple fluid–solid boundaries as well. Various methods (continuum, discrete, mesoscopic) have been developed to simulate droplets, cells, capsules and so on, and we refer to existing literature on this.^21,22 We restrict ourselves to those methods that employ or pertain to LBM. Previous reviews and books^23–25 have covered some fluid–solid and fluid–fluid methods in a different context. For instance, Karimnejad et al.²³ have reviewed fluid–structure (FS) methods for rigid particles. Here our focus is on methods for deformable particles in fluids. We present a comprehensive review of those methods and compare them.

Finally, we wish to underline that the methods described in this review straddle across the passive versus active soft matter divide, just because, once the coupling between internal and external degrees of freedom is accounted for, such divide largely blurs out.²⁶ The methods described here may be used to investigate problems in the active matter systems such as the emergent nematic order in epithelial cells, which share many similarities with the emulsions of Section IV, and the shape changes in bacteria, which become more elongated when swarming.^27,28

The review is organised as follows. The LBM is summarised in Section II. In Section III, we review fluid–structure methods which often use the immersed boundary method. We highlight the differences and similarities between the continuum and the discrete approaches. In Section IV, we review the various multicomponent methods for immiscible droplets. In Section V, we provide benchmarks in 2D and 3D for hydrostatic and hydrodynamic conditions. In Section VI, we conclude by summarising our observations and make remarks on open questions and future challenges.

II. Lattice Boltzmann method

In recent years the LBM has become a useful tool in solving the hydrodynamics of soft matter systems.²⁹ 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,³⁰ it is now considered a separate method. Several works have shown the potential of LBM to solve fluid-related problems.^13,31–35 It is an efficient algorithm particularly popular for Stokes flow, systems with complex solid boundaries,³⁶ multiphase flow,³⁷ hydrodynamic dispersion³⁸ and convection.³⁹

A. Lattice BGK model

The Boltzmann equation that contains the fluid particle information readswhere 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 bywhere τ 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 bywhere 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 aswhere 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 one-dimensional 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.²⁵

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.

The final discrete lattice Boltzmann equation (with the BGK operator) is,

where 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. 40–42). A popular one, the Guo forcing scheme,⁴³ is given bywhere F is the force. There are two main steps in eqn (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. 25 and 44. 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^45,46 and the entropic LBM.^47,48

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 membrane-bound particles is red blood cells (RBCs). We highlight that variations of the models presented here have been used for soft objects like microgels⁴⁹ and capsules.⁵⁰ 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.

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.

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.⁵¹ 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.⁵² Another consequence of applying a constitutive model in the entire membrane domain is that the shape memory might be neglected in some cases.⁵³ This approach is discussed in Section IIIA.

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. 54. In fact, ref. 54 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.⁵⁴ Additional terms have to be included. This approach is also known to suffer less numerical instability than the continuum approach.⁵² This discrete approach is discussed in Section IIIB.

The elastic forces are treated differently in each approach. In the continuum approach, these are solved using continuum-based 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 IIIC. Additionally, the fluid can be coupled with membrane dynamics using the immersed boundary method (IBM). In Section IIID, 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⁵⁵ and where the elastic bodies have negligible mass.⁵⁶ In IBM, the flow is solved on an Eulerian grid while the membrane is treated on a Lagrangian grid.⁵⁷ 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 eqn (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 Section IIID.

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 sections, 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

W = W_s + W_b + W_A + W_V, (10)

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 with the fluid. The force acting on membrane node i at position x_i is given by the principle of virtual work,^58–61Details about derivations using the equation for each term can be found in ref. 62. For simplicity, we usually consider more straightforward expressions for the force instead of the full expressions derived in the next sections.

1. 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. 63–65 for commonly used models which we also discuss in this section. The Neo-Hookean (NH) model⁶³ is one of the most straightforward choices among different constitutive equations, in the formwhere E_s is the surface shear elasticity modulus, λ₁ and λ₂ 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 eqn (12) is the zero-thickness (ZT) shell representation⁶⁶ with the strain energy function given as:

Fig. 4 (a) Example of a deformation of a 2D sheet due to forces F₁ and F₂. Deformation in the x direction is characterised by the stretch ratio λ₁ or by the displacement u₁ (the same in the y direction). λ₁ and λ₂ 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. Inspired by ref. 67.

A widely used model suggested by Skalak;⁶⁸ see e.g. ref. 18 and 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. Eqn (12)–(14) can also be written as a function of the strain invariants I₁ and I₂ defined for a 2D membrane (with isotropic and homogeneous properties) asWe note that some of the methods, such as that proposed by ref. 43, 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 W = W(I₁, I₂, α₁, α₂, α₃) withwhere W₀ is a reference value, I₁ and I₂ are the strain invariants and α₁, α₂ and α₃ 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.⁷⁰

Once the constitutive law is chosen, we can compute the elastic forces using the linear FEM.⁵⁸ 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,71 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, 58, 59, 62 and 67,72,73.

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ⁱ_x and δuⁱ_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

δW_s = A₀δW, (19)

where A₀ 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, λ₁ and λ₂ 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ⁱ_x and δuⁱ_y:To obtain expressions for the nodal forces, F^s,i_x and F^s,i_y, eqn (18) and (20) are substituted in eqn (19) and using the fact that displacements, δuⁱ_x and δuⁱ_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. 58. 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.

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)–(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.

The derivatives in eqn (21) are computed using a linear FEM, with linear shape functions, Nⁱ, i = 1, 2, 3, where N¹, is defined as^58,67,73

and N² and N³ can be obtained by cycling the indices from 1 → 2 → 3 → 1. 2A₀ is twice the area of the undeformed element and is computed as follows:

2A₀ = |(y⁽²⁾ − y⁽³⁾)x⁽¹⁾ + (x⁽³⁾ − x⁽²⁾)y⁽¹⁾ + x⁽²⁾y⁽³⁾ − x⁽³⁾y⁽²⁾|. (23)

In turn, the spatial dependence of the displacement on the element can be expressed asThe relation between (x, y) and (X, Y) (see eqn (17)) is given in matrix notation by dX = F·dx where the matrix of the Cartesian components of the deformation gradient tensor is as shown in ref. 58 and 59. Lastly, one just needs to compute the principle stretch ratios, λ₁ and λ₂, from the eigenvalues of the two-dimensional right Cauchy–Green tensor, G = F^TF, whose components are given byThe derivatives of the in-plane displacements, u and v, in eqn (25) for G₂₁ = G₁₂ are computed by differentiating their linear representations given by eqn (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. 74, and subsequently widely used.^75,76 Although theoretically possible to use molecular dynamics to model fluid–structure interaction 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 (coarse-grained) 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.

Fig. 6 Examples of different shapes that can be used in the lattice-Boltzmann method with finite elements and immerse boundary. (a) is a sphere, (b) is a biconcave shell and (c) is an ellipsoid. An icosahedron was used for the subdivision scheme and there were 5 subdivisions, so the total number of elements is N_e⁵ = 20 × 4⁵ = 20 480 faces.

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.^54,74 The total energy of a membrane is

W = W_stretch + W_bending + W_area + W_volume, (27)

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 byFull derivations are not shown but can be found in ref. 77.

1. Elastic membrane models. A harmonic stretching potential can be used to model the stretching energy aswhere 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. Eqn (29) has been used to study diseased RBCs,⁷⁸ blood flow in capillaries⁷⁹ and has been validated against modern experimental techniques such as atomic force microscopy.⁸⁰Eqn (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_s0e^2(ι−1) where the bond stretch ratio (ι) is included.⁸¹ The stretch ratio ι is defined as ι = l/l₀ where l and l₀ are length and the equilibrium length of the spring, respectively. It has been shown to achieve good results between experiments and simulations.⁷⁹ Similarly, other nonlinear potentials may be used to model the stretching energy. The two most popular nonlinear potentials⁷⁵ are the WLC potential and FENE potential. Explicitly, they are defined aswhere 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. 75 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

W_stretch = W_WLC + W_POW or W_stretch = W_FENE + W_POW, (31)

where the inverse power repulsive potential is given byand k_p is the repulsive stiffness. Finally, the forces for eqn (30) and (32) arewhere x = l/l_m and l̂_ij = l⃑_ij/l is the vector of unit length between nodes i and j.

C. Bending and area & volume constraints

1. 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.^82–84 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. 82 and 83 for bending models and numerical implementations. In line with the continuum approach, one of the most popular bending formulations is that of Helfrich⁸⁵where κ is the mean curvature and κ_G is the Gaussian curvature, c₀ 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^86,87

F^b_i = E_B[(2κ + c₀)(2κ² − 2κ_g − c₀κ) + 2Δ_LBκ]n (35)

where n is the outwards normal vector to the surface and Δ_LB is the Laplace–Beltrami operator. Calculation of eqn (35) and Δ_LBκ can be found in ref. 82, 83 and 88.

Another popular and simplified expression for the bending energy⁸⁴ 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 eqn (36) is given byDetails on the calculation of can be found in ref. 62 and 84.

2. 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⁶²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₀ and volume V₀. While A₀ and V₀ are input parameters, A and V need to be computed at every time step. From eqn (38) we find that the corresponding forces areMore details and computation of and may be found in ref. 62. 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. 89. We can also express eqn (38) as,⁹⁰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 F_i = F^s_i + F^b_i + F^A_i + F^V_i acting on node i is,

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.^91–93 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.⁹¹ 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⁵⁷ 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. 25, 57, 94 and 95. Other coupling techniques such as the stress integration approach⁹⁶ and friction coupling approach^29,97–99 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^57,100 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,¹⁰¹ blood flow^18,102,103 and platelet migration.¹⁰⁴ We refer the reader to ref. 57 and 105 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. 57,95 defines the force density g(x, t) exerted by the structure on the fluid as

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) asNote that eqn (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.

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. Inspired by ref. 62.

1. Spatial and temporal discretization. Implementing the IBM requires both spatial and temporal discretization of eqn (42) and (43), is needed for the numerical implementation of the IBM. Following ref. 57 and 94, 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. 57 recommended that to prevent fluid penetration/leaking. However, ref. 62 systematically studied mesh ratio between solid and fluid for small deformations, e.g., 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 eqn (42) becomes

where F(p, q, r, t) is the force density from the solid structures. As in ref. 94, let us define . 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. Eqn (44) is then simplified towhere X_i is the position of the i-th Lagrangian node. Eqn (43) is simplified toIn the LBM, the time and spatial step in LB units are usually assumed unity, e.g., Δt = 1, Δx = 1, respectively. Thus, eqn (46) is simplified toref. 57 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)

δ_Δ(x) = ϕ(x)ϕ(y)ϕ(z), (48)

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. 57 showed that ϕ(x) can be accurately estimated using a straightforward formulaHowever, eqn (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.¹⁰⁶ 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. 106–108.

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. 57,95 suggested an integration scheme based midpoint rule.⁵⁷ 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ⁿ_i at time step n, the intermediate position at time step is calculated using

Using the updated solid nodal position at the time step, the force acting on the structure can be evaluated aswhere we take the gradient of the energy function W. Next, the structure force will be spread out into the fluid throughFor specific 2D simulations (e.g. sedimentation) an area scaling factor A_b can be used in eqn (53).¹⁰⁹ Using the force density in the LBM forcing scheme (e.g.eqn (8)), we can then stream and collide using eqn (7) and then calculate velocity using eqn (9). Finally, the solid position at time n + 1 is updated asWe recall that we use Δt = 1 and Δx = 1 in the above equations.

2. 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,^110–112 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_ij is the position vector difference between sphere i and j, defined as x_ij = x_i − x_j, x̂_ij is the unit vector. u is the sphere's velocity. Eqn (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.¹¹³ This force is zero for node-to-node distances larger than one lattice constant and increases as 1/r² at smaller distances

where d_ij is the distance between the two nodes i and j, κ_int is a constant regulating the strength and g^int_ji = −g^int_ij.

3. Other coupling schemes. Here we discuss only the no-slip boundary conditions between fluid and structure, motivated by ref. 94, which consists of imposing velocity and force continuity along the interfacewhere u_f and u_s are the fluid and solid velocity, is the fluid–solid boundary and σ_ij 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 eqn (43), while the force is spread into the fluid, as shown in eqn (42). Alternatively, the inverse method can be used.^96,114 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.^96,115 The total force acting on the structure iswhere dA is the differential area over the interface. The stress tensor σ_ij is given by

σ_ij = −pδ_ij + ρν(u_i,j + u_j,i). (58)

The pressure term p is usually assumed to follow an ideal gas law p = ρc_s². Following ref. 116, the deviatoric shear stress in the LBM τ_ij := ρν(u_i,j + u_j,i) can be evaluated aswhere 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 eqn (6). ξ_α is the discretized velocity vector. Using the stress computed from the LBM through eqn (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.^97–99,117 In ref. 98, polymer monomers are treated as points with a friction force

F_ζ = −ζ(u_s − u_f), (60)

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. 97 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. 97 and 98 uses a force density −F_ζ/δx³. 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.^70,118,119

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. 25, 35, 120 and 121. These have applications in several fields, such as for example in injet printing.¹²² The LBM facilitates the relatively seamless integration of inter-particle 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. 123. 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. 124 and later, in ref. 125 the model underwent changes to accommodate variations in viscosity and density ratios. Another key development was proposed in ref. 126, namely, to add a recolouring step which significantly reduces the computational requirements while minimising the spurious currents and removing the lattice pinning effect.^127,128 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

f_α,k(x + ξ_αΔt,t + Δt) = f_α,k(x,t) + Ω_α,k(x,t), (61)

where Ω_α,k is the multistage collision operator expressed as

Ω_α,k = [Ω⁽¹⁾_α,k + Ω⁽²⁾_α,k]Ω⁽³⁾_α,k. (62)

Each has a distinct aim:

• Ω⁽¹⁾_α,k is the normal single-phase collision operator for each component;

• Ω⁽²⁾_α,k is the perturbation operator, generating the interfacial tension;

• Ω⁽³⁾_α,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 byLikewise, for the D3Q19 velocity set, B_α are given bywhere χ is a free parameter: values for B_α and χ are derived in ref. 129. Calculation of partial derivatives is necessary for ∇ρ and it has been shown¹³⁰ 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. 123 this involves the minimisation of the work, W, which is given by

W(f^r,f^b) = −∇ρ·q(f^r,f^b), (68)

where q is the local or colour fluxNevertheless, this recolouring method may lead to pinning. This was corrected in ref. 126 by enabling the blending of red and blue components through symmetric colour distribution concerning the colour gradient. The particle distributions after recolouring are given aswhere is used to manage the thickness of the interface, where is the post-collision state after Ω⁽¹⁾_α,k and Ω⁽²⁾_α,k have been applied. The term cos(θ_α) is conveniently expressed asThe total zero-velocity equilibrium distribution function is given bywhere the weights ϕ_α (along with w_α) are lattice dependent.^131,132 Further modification of the recolouring step in eqn (70) to improve spurious velocities, lattice pinning and convergence was done first in ref. 131. 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¹³³ and viscosity ratios.¹³⁴ 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.^135–137 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

F_rep = −A_h[h(x)]nδ_I, (73)

where A_h[h(x)] is the parameter controlling the strength of this repulsive force, n is a unit normal to the interface and is a function proportional to the phase field which confines the repulsive force to the interface. A graphical representation of this force is shown in Fig. 8. A_h[h(x)] can be set equal to a constant A if h ≤ h_min, it decays as f⁻³ if h_min < h ≤ h_max and it is equal to zero if h > h_max. For example in ref. 137, 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 towhich 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.

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.

B. Free-energy modelling

We introduce a straightforward free-energy LBM for binary liquid mixtures.^138,139 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.^140–143 The scalar not only moves with the flow but also undergoes diffusion to minimise a free-energy 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 eqn (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 areand the equilibrium distribution is

f^eq_k,α = A_α + B_αu·ξ_α + C_αu·u + D_α(u·ξ_α)² + G_α,xyξ_α,xξ_α,y, (77)

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.¹⁴⁴ For a D2Q9 velocity set, the full list of coefficients is given in ref. 44 and details for construction on a D3Q19 velocity set are given in ref. 144.

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 , while the chemical potential is determined by the free-energy of the system. The free-energy functional [ϕ(x)] is¹³⁹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 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 ϕ₀ = ±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 byThe 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¹³⁹where is the interfacial thickness of the droplets. The free energy model is also used to study active droplets.^141,145,146

Frustrated coalescence. Literature on frustrated coalescence using the free energy LBM is scarce. We present a method used by ref. 147 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 = ϕ₀ 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 and can therefore be tuned by changing the value of κ_ϕ in eqn (83). The parameters κ_ϕ and α also determine the interfacial thickness of the droplets . 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 bywhere M is the mobility and is the chemical potential of the i-th droplet. Eqn (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. 138 and 148. This may be used for a small number of droplets. As the number of times that eqn (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.¹⁴⁹

C. Pseudopotential method

The pseudopotential (or Shan–Chen models) first introduced in ref. 150 and 151 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¹⁵² and transport of H₂O in air.^153,154 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 aswhere 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 bywhere G_k̃k 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 iswhere ρ₀ is a freely chosen reference density, usually unity. ψ_k(x,t) = ρ_k can also be used and it is an approximation of eqn (87) at low-density. However, at higher densities, this assumption will lead to numerical instabilities.¹⁵⁵ In the absence of attraction between components, the pseudopotential can often be reduced to the physical density for simplicity. Nevertheless, using eqn (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. 150 and 151, 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 ρ_ku_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 byAlso, 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 eqn (7)),

The equilibrium distribution is given byThe Guo forcing scheme¹⁵⁶ can be used to implement these forces acting in a fluid as it yields a viscosity-independent surface tension,¹⁵⁷ particularly with the frustrated coalescence mechanism described in the next section. Thus similar to eqn (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.¹⁵⁸

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 large-scale 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.²⁵ Some attempts have also been made to derive a thermodynamically consistent pseudopotential model.^159–161 Moreover, it is possible to obtain an effective free energy in the pseudopotential models.^162,163

Frustrated coalescence. To prevent coalescence a multi-range repulsive force is used which acts between the nearest neighbours and the next nearest neighbours.^157,164 This force is used in addition to the repulsive component separation force F^r_k. This means that in eqn (93), F_k = F^r_k + F^c_k where F^c_k is a competing force given bywhere 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,¹⁶⁴ respectively while in 3D they are D3Q41 for L1 and D3Q39 for L2.¹⁶⁵ 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.¹⁶⁶ 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: G_A,1 = G_B,1, G_A,2 = G_B,2, G_A,B = G_B,A for the two components A and B.

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.¹⁶⁷ 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.¹⁶⁶ 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 (Fig. 9 and 10).

Fig. 9 Droplets moving in a channel. This was obtained using a free energy LBM and solving eqn (84). Colour depicts which is ≈2 for droplets (yellow) and ≈0 for the background fluid (black).

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).

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. 157 for benchmarks on multiple droplets.

A. Capsule and droplet in shear flow

Shear flow is a common validation flow to test for both droplets^168,169 and capsules^115,170 (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¹⁶⁹ (see Fig. 11). In what follows, we have separated into two sections 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). Inspired by ref. 169.

1. 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 bywhere 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. 171 by solving the momentum and continuity equation using the Lorentz reflection method. ref. 171 gives the first order solution D_SH (SH stands for the initials of the original authors in ref. 171) for small confinements ratio R/H between two parallel walls^171,172where 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. 171 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. 172. The theoretical model given eqn (97) and ref. 171 has an additional term D_Taylor obtained from Taylor small-deformation bulk theory^173,174 and is given bywhere 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. 172 and 173. Taylor bulk theory^173,174 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. 175. The methodology presented here is valid for all the fluid–fluid models mentioned in this review.

2. 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 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. 176–178. For Skalak models we refer the reader to ref. 18, 179 and 180. 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, eqn (12) and (14) and is given by,^18,43,181where α₁ = 0, α₂ = 2/3, and α₃ = 1/3 are coefficients extracted from the constitutive model. Eqn (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 eqn (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.

Fig. 12 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.5H. 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 eqn (100).

According to ref. 109, 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 eqn (100) is obtained by equating the gravitational force per unit length with the drag force per unit length F_d = V_cμΛ. Moreover, eqn (100) is valid for low Reynolds numbers. There are several approximations for Λ which can be found in ref. 109 and 182. A commonly cited approximation can be found in ref. 183 and 184 where Stokes equation for Newtonian fluids using asymptotic expansion with no-slip boundary conditions on the walls given bywhere 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. 183,184 applied an integral transform which produces a convergent series and allows the boundary conditions to be applied simultaneously. It is generally known that eqn (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 . 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.¹⁸⁵ Moreover, other hydrostatic tests such as contact angle test¹⁸⁵ (for wetting boundary conditions) can be useful to further check the accuracy of the model.

Fig. 13 Laplace test using the pseudopotential method at equilibrium for G_k,1 = −7.9, G_k,2 = 4.9 (see eqn (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.

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.^186,187 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⁶ 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,188 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).¹⁸⁹

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.¹⁹⁰

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.¹⁹¹

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

We acknowledge financial support from the Portuguese Foundation for Science and Technology (FCT) under the contracts: PTDC/FIS-MAC/28146/2017 (LISBOA-01-0145-FEDER-028146), EXPL/FIS-MAC/0406/2021 (DOI: 10.54499/EXPL/FIS-MAC/0406/2021), UIDB/00618/2020 (DOI: 10.54499/UIDB/00618/2020), UIDP/00618/2020 (DOI: 10.54499/UIDP/00618/2020), 2020.08525.BD and DL57/2016/CP1479/CT0057 (DOI: 10.54499/DL57/2016/CP1479/CT0057). SS wishes to acknowledge financial support from the ERC-PoC grants DROPTRACK (contract number 101081171).

References

1. S. L. Anna and H. C. Mayer, Phys. Fluids, 2006, 18, 121512.
2. H. Liu and Y. Zhang, Phys. Fluids, 2011, 23, 082101.
3. J. R. Clausen and C. K. Aidun, Phys. Fluids, 2010, 22, 123302.
4. Y. Sui, Y. T. Chew and H. T. Low, Int. J. Mod. Phys. C, 2007, 18, 993–1011.
5. P.-G. Gennes, F. Brochard-Wyart, D. Quéré, et al., Capillarity and wetting phenomena: drops, bubbles, pearls, waves, Springer, 2004.
6. J. N. Israelachvili, The Handbook of Surface Imaging and Visualization, CRC Press, 2022, pp. 793–816.
7. W. Qi and P. B. Weisensee, Phys. Fluids, 2020, 32, 067110.
8. T. Truong, N. Bansal and B. Bhandari, Food Bioprocess Technol., 2014, 7, 3416–3428.
9. E. Dickinson, Food Hydrocolloids, 2012, 28, 224–241.
10. D. J. McClements, Sustainable Food Proteins, 2023, 1, 101–132.
11. D. A. Fedosov, M. Peltomäki and G. Gompper, Soft Matter, 2014, 10, 4258–4267.
12. J. Dupire, M. Socol and A. Viallat, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 20808–20813.
13. C. K. Aidun and J. R. Clausen, Annu. Rev. Fluid Mech., 2010, 42, 439–472.
14. R. Benzi, S. Succi and M. Vergassola, Phys. Rep., 1992, 222, 145–197.
15. T. Ye, D. Pan, C. Huang and M. Liu, Phys. Fluids, 2019, 31, 011301.
16. H. Zhu, Z. Zhou, R. Yang and A. Yu, Chem. Eng. Sci., 2008, 63, 5728–5770.
17. T. Gao and H. H. Hu, J. Comput. Phys., 2009, 228, 2132–2151.
18. T. Krüger, F. Varnik and D. Raabe, Comput. Math. Appl., 2011, 61, 3485–3505.
19. K. Sugiyama, S. Ii, S. Takeuchi, S. Takagi and Y. Matsumoto, Comput. Mech., 2010, 46, 147–157.
20. W. Gao and Y. Chen, Int. J. Heat Mass Transfer, 2019, 135, 158–163.
21. G. Hou, J. Wang and A. Layton, Commun. Comput. Phys., 2012, 12, 337–377.
22. U. D. Schiller, T. Krüger and O. Henrich, Soft Matter, 2018, 14, 9–26.
23. S. Karimnejad, A. A. Delouei, H. Basağaoğlu, M. Nazari, M. Shahmardan, G. Falcucci, M. Lauricella and S. Succi, Commun. Comput. Phys., 2022, 32, 899–950.
24. S. Succi, The Lattice Boltzmann Equation: For Complex States of Flowing Matter, Oxford University Press, 2018.
25. T. Krüger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva and E. M. Viggen, The Lattice Boltzmann Method, Springer International Publishing, 2017.
26. M. L. Manning, Phys. Rev. Lett., 2023, 130, 130002.
27. R. C. V. Coelho, N. A. M. Araújo and M. M. Telo da Gama, Soft Matter, 2022, 18, 7642–7653.
28. A. Beer and G. Ariel, Movement Ecology, 2019, 7, 1–17.
29. B. Dunweg and A. J. C. Ladd, Advanced Computer Simulation Approaches for Soft Matter Sciences III, Springer Berlin Heidelberg, 2009, pp. 89–166.
30. U. Frisch, B. Hasslacher and Y. Pomeau, Phys. Rev. Lett., 1986, 56, 1505–1508.
31. Y. H. Qian, D. D'Humières and P. Lallemand, Europhys. Lett., 1992, 17, 479–484.
32. X. He and L.-S. Luo, J. Stat. Phys., 1997, 88, 927–944.
33. A. K. Gunstensen, D. H. Rothman, S. Zaleski and G. Zanetti, Phys. Rev. A: At., Mol., Opt. Phys., 1991, 43, 4320–4327.
34. S. Chen and G. D. Doolen, Annu. Rev. Fluid Mech., 1998, 30, 329–364.
35. S. Succi, The Lattice Boltzmann Equation: For Fluid Dynamics and Beyond, Oxford Science Publications, 2001.
36. A. J. C. Ladd, J. Fluid Mech., 1994, 271, 311–339.
37. X. He, S. Chen and R. Zhang, J. Comput. Phys., 1999, 152, 642–663.
38. A. Cal, S. Succi, A. Cancelliere, R. Benzi and M. Gramignani, Phys. Rev. A: At., Mol., Opt. Phys., 1992, 45, 5771–5774.
39. X. He, S. Chen and G. D. Doolen, J. Comput. Phys., 1998, 146, 282–300.
40. S. A. Bawazeer, S. S. Baakeem and A. A. Mohamad, Arch. Comput. Methods Eng., 2021, 28, 4405–4423.
41. K. Sun, T. Wang, M. Jia and G. Xiao, Phys. A, 2012, 391, 3895–3907.
42. L. Chen, Q. Kang, Y. Mu, Y.-L. He and W.-Q. Tao, Int. J. Heat Mass Transfer, 2014, 76, 210–236.
43. D. Barthès-Biesel and J. M. Rallison, J. Fluid Mech., 1981, 113, 251.
44. H. Huang, M. C. Sukop and X.-Y. Lu, Multiphase lattice Boltzmann methods: theory and application, John Wiley and Sons, Inc., Chichester, West Sussex, 2015.
45. P. Asinari, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 78, 016701.
46. P. Lallemand and L.-S. Luo, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2000, 61, 6546–6562.
47. B. M. Boghosian, J. Yepez, P. V. Coveney and A. Wager, Proc. R. Soc. London, Ser. A, 2001, 457, 717–766.
48. A. Mazloomi M, S. S. Chikatamarla and I. V. Karlin, Phys. Rev. Lett., 2015, 114, 174502.
49. W. Lei, C. Xie, T. Wu, X. Wu and M. Wang, Sci. Rep., 2019, 9, 1453.
50. J. Zhang, P. C. Johnson and A. S. Popel, Phys. Biol., 2007, 4, 285–295.
51. W. R. Dodson and P. Dimitrakopoulos, J. Fluid Mech., 2009, 641, 263–296.
52. Y. Imai, T. Omori, Y. Shimogonya, T. Yamaguchi and T. Ishikawa, J. Biomech., 2016, 49, 2221–2228.
53. T. M. Fischer, Biophys. J., 2004, 86, 3304–3313.
54. T. Omori, T. Ishikawa, D. Barthès-Biesel, A.-V. Salsac, J. Walter, Y. Imai and T. Yamaguchi, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2011, 83, 041918.
55. D.-K. Sun and Z. Bo, Int. J. Heat Mass Transfer, 2015, 80, 139–149.
56. D. C. Bottino, J. Comput. Phys., 1998, 147, 86–113.
57. C. S. Peskin, Acta Numer., 2002, 11, 479–517.
58. J. M. Charrier, S. Shrivastava and R. Wu, J. Strain Anal. Eng. Des., 1989, 24, 55–74.
59. S. Shrivastava and J. Tang, J. Strain Anal. Eng. Des., 1993, 28, 31–51.
60. H. Goldstein, C. Poole and J. Safko, Classical Mechanics, Pearson, 3rd edn, 2002.
61. G. R. Lázaro, I. Pagonabarraga and A. Hernández-Machado, Eur. Phys. J. E: Soft Matter Biol. Phys., 2017, 40, 77.
62. T. Krüger, Computer Simulation Study of Collective Phenomena in Dense Suspensions of Red Blood Cells under Shear, Vieweg+Teubner Verlag, 2012.
63. D. Barthès-Biesel, A. Diaz and E. Dhenin, J. Fluid Mech., 2002, 460, 211–222.
64. E. Lac, D. Barthès-Biesel, N. A. Pelekasis and J. Tsamopoulos, J. Fluid Mech., 2004, 516, 303–334.
65. Z. Hashemi, M. Rahnama and S. Jafari, Sci. Iran., 2015, 22, 1877–1890.
66. S. Ramanujan and C. Pozrikidis, J. Fluid Mech., 1998, 361, 117–143.
67. C. Armstrong and Y. Peng, Phys. Rev. E, 2021, 103, 023309.
68. R. Skalak, A. Tozeren, R. P. Zarda and S. Chien, Biophys. J., 1973, 13, 245–264.
69. Y. Sui, Y. T. Chew, P. Roy, Y. P. Cheng and H. T. Low, Phys. Fluids, 2008, 20, 112106.
70. F. Pelusi, F. Guglietta, M. Sega, O. Aouane and J. Harting, Phys. Fluids, 2023, 35, 082126.
71. K. Vahidkhah, S. L. Diamond and P. Bagchi, Biophys. J., 2014, 106, 2529–2540.
72. J. Charrier, S. Shrivastava and R. Wu, J. Strain Anal. Eng. Des., 1987, 22, 115–125.
73. C. L. Armstrong, PhD thesis, Old Dominion University, 2021.
74. J. Li, G. Lykotrafitis, M. Dao and S. Suresh, Proc. Natl. Acad. Sci. U. S. A., 2007, 104, 4937–4942.
75. G. E. K. D. A. Fedosov and B. Caswell, Biophys. J., 2010, 98, 2215–2225.
76. K. Ichi Tsubota, S. Wada and T. Yamaguchi, Comput. Methods Programs Biomed., 2006, 83, 139–146.
77. C. Pozrikidis, Computational hydrodynamics of capsules and biological cells, CRC Press, 2010.
78. T. Wu and J. J. Feng, Biomicrofluidics, 2013, 7, 044115.
79. M. Nakamura, S. Bessho and S. Wada, Int. J. Numer. Method. Biomed. Eng., 2013, 29, 114–128.
80. S. Barns, M. A. Balanant, E. Sauret, R. Flower, S. Saha and Y. Gu, Biomed. Eng. Online, 2017, 16, 1–21.
81. M. Nakamura, S. Bessho and S. Wada, Int. J. Numer. Method. Biomed. Eng., 2014, 30, 42–54.
82. A. Guckenberger, M. P. Schraml, P. G. Chen, M. Leonetti and S. Gekle, Comput. Phys. Commun., 2016, 207, 1–23.
83. A. Guckenberger and S. Gekle, J. Phys.: Condens. Matter, 2017, 29, 203001.
84. X. Bian, S. Litvinov and P. Koumoutsakos, Comput. Methods Appl. Mech. Eng., 2020, 359, 112758.
85. W. Helfrich, Z. Naturforsch., 1973, 28, 693–703.
86. O.-Y. Zhong-can and W. Helfrich, Phys. Rev. A: At., Mol., Opt. Phys., 1989, 39, 5280–5288.
87. K. Sinha and M. D. Graham, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2015, 92, 042710.
88. M. Reuter, S. Biasotti, D. Giorgi, G. Patanè and M. Spagnuolo, Comput. Graph., 2009, 33, 381–390.
89. C. Zhang and T. Chen, Proceedings 2001 International Conference on Image Processing (Cat. No. 01CH37205), Thessaloniki, Greece, 2001, pp. 935–938.
90. Z. Gou, X. Ruan, F. Huang and X. Fu, Comput. Methods Biomech. Biomed. Eng., 2019, 22, 465–474.
91. C. W. Hirt, A. A. Amsden and J. L. Cook, J. Comput. Phys., 1974, 14, 227–253.
92. F.-B. Tian, H. Dai, H. Luo, J. F. Doyle and B. Rousseau, J. Comput. Phys., 2014, 258, 451–469.
93. C. Duprat and H. A. Shore, Fluid-structure interactions in low-Reynolds-number flows, Royal Society of Chemistry, 2016.
94. J. Tan, PhD thesis, Lehigh University, 2015.
95. C. S. Peskin and D. M. McQueen, J. Comput. Phys., 1989, 81, 372–405.
96. S. Kollmannsberger, S. Geller, A. Düster, J. Tölke, C. Sorger, M. Krafczyk and E. Rank, Int. J. Numer. Meth. Eng., 2009, 79, 817–845.
97. P. Ahlrichs and B. Dünweg, Int. J. Mod. Phys. C, 1998, 09, 1429–1438.
98. P. Ahlrichs and B. Dünweg, J. Chem. Phys., 1999, 111, 8225–8239.
99. C. W. Hsu and Y.-L. Chen, J. Chem. Phys., 2010, 133, 034906.
100. C. S. Peskin, J. Comput. Phys., 1972, 10, 252–271.
101. G. Herschlag and L. Miller, J. Theor. Biol., 2011, 285, 84–95.
102. J. Zhang, P. C. Johnson and A. S. Popel, J. Biomech., 2008, 41, 47–55.
103. Y. Liu and W. K. Liu, J. Comput. Phys., 2006, 220, 139–154.
104. L. M. Crowl and A. L. Fogelson, Int. J. Numer. Method. Biomed. Eng., 2010, 26, 471–487.
105. R. Mittal and G. Iaccarino, Annu. Rev. Fluid Mech., 2005, 37, 239–261.
106. T. Kajishima and K. Taira, Computational fluid dynamics: incompressible turbulent flows, Springer, 2016.
107. S. J. Shin, W.-X. Huang and H. J. Sung, Int. J. Numer. Methods Fluids, 2008, 58, 263–286.
108. X. Yang, X. Zhang, Z. Li and G.-W. He, J. Comput. Phys., 2009, 228, 7821–7836.
109. S. Ghosh and J. M. Stockie, Commun. Comput. Phys., 2015, 18, 380–416.
110. A. J. C. Ladd and R. Verberg, J. Stat. Phys., 2001, 104, 1191–1251.
111. A. J. C. Ladd, Phys. Fluids, 1997, 9, 491–499.
112. A. S. Nunes, R. C. V. Coelho, V. C. Braz, M. M. Telo da Gama and N. A. M. Araújo, Commun. Comput. Phys., 2023, 33, 1–21.
113. M. Gross, T. Krueger and F. Varnik, Soft Matter, 2014, 10, 4360.
114. Y. Kwon, Eng. Comput., 2006, 23, 860–875.
115. R. M. MacMeccan, J. R. Clausen, G. P. Neitzel and C. K. Aidun, J. Fluid Mech., 2009, 618, 13–39.
116. T. Krüger, F. Varnik and D. Raabe, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2009, 79, 046704.
117. B. Dünweg and A. J. C. Ladd, Advanced Computer Simulation Approaches for Soft Matter Sciences III, Springer Berlin Heidelberg, Berlin, Heidelberg, 2009.
118. P. Li and J. Zhang, Int. J. Numer. Method. Biomed. Eng., 2019, 35, e3200.
119. D. Taglienti, F. Guglietta and M. Sbragaglia, Phys. Rev. Fluids, 2023, 8, 013603.
120. F. F. Jackson, PhD thesis, University of Leeds, 2021.
121. O. Shardt, PhD thesis, University of Alberta, 2014.
122. F. F. Jackson, K. J. Kubiak, M. C. Wilson, M. Molinari and V. Stetsyuk, Langmuir, 2019, 35, 9564–9571.
123. D. Rothman and J. Keller, J. Stat. Phys., 1988, 52, 1119–1127.
124. A. Gunstensen, D. Rothman, S. Zaleski and G. Zanetti, Phys. Rev. A: At., Mol., Opt. Phys., 1991, 43, 4320.
125. D. Grunau, S. Chen and K. Eggert, Phys. Fluids A, 1993, 5, 2557–2562.
126. M. Latva-Kokko and D. Rothman, Phys. Rev., 2005, 71, 056702.
127. K. Connington and T. Lee, J. Mech. Sci. Technol., 2012, 26, 3857–3863.
128. X. Shan, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2006, 73, 047701.
129. S. Leclaire, M. Reggio and J.-Y. Trépanier, Appl. Math. Model., 2012, 36, 2237–2252.
130. H. Liu, A. Valocchi and Q. Kang, Phys. Rev., 2012, 85, 046309.
131. T. Reis and T. N. Phillips, J. Phys. A: Math. Theor., 2007, 40, 4033–4053.
132. S. Leclaire, A. Parmigiani, O. Malaspinas, B. Chopard and J. Latt, Phys. Rev. E, 2017, 95, 033306.
133. Y. Ba, H. Liu, Q. Li, Q. Kang and J. Sun, Phys. Rev. E, 2016, 94, 023310.
134. H. Liu, L. Wu, Y. Ba, G. Xi and Y. Zhang, J. Comput. Phys., 2016, 327, 873–893.
135. A. Montessori, M. Lauricella and S. Succi, Philos. Trans. R. Soc., A, 2019, 377, 20180149.
136. A. Montessori, M. Lauricella, N. Tirelli and S. Succi, J. Fluid Mech., 2019, 872, 327–347.
137. A. Montessori, A. Tiribocchi, F. Bonaccorso, M. Lauricella and S. Succi, Philos. Trans. R. Soc., A, 2020, 378, 20190406.
138. M. Swift, E. Orlandini, W. Osborn and J. Yeomans, Phys. Rev., 1996, 54, 5041–5052.
139. A. J. Briant and J. M. Yeomans, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2004, 69, 031603.
140. D. Marenduzzo, E. Orlandini, M. E. Cates and J. M. Yeomans, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 76, 031921.
141. R. C. V. Coelho, N. A. M. Araújo and M. M. Telo da Gama, Philos. Trans. R. Soc., A, 2021, 379, 20200394.
142. G. R. Lázaro, A. Hernández-Machado and I. Pagonabarraga, Soft Matter, 2014, 10, 7207–7217.
143. G. R. Lázaro, A. Hernández-Machado and I. Pagonabarraga, Soft Matter, 2014, 10, 7195–7206.
144. C. M. Pooley and K. Furtado, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 77, 046702.
145. A. Tiribocchi, M. Durve, M. Lauricella, A. Montessori, D. Marenduzzo and S. Succi, Nat. Commun., 2023, 14, 1096.
146. R. C. V. Coelho, H. R. J. C. Figueiredo and M. M. Telo da Gama, Phys. Rev. Res., 2023, 5, 033165.
147. M. Foglino, A. Morozov, O. Henrich and D. Marenduzzo, Phys. Rev. Lett., 2017, 119, 208002.
148. A. Tiribocchi, N. Stella, G. Gonnella and A. Lamura, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2009, 80, 026701.
149. A. Tiribocchi, A. Montessori, M. Lauricella, F. Bonaccorso, S. Succi, S. Aime, M. Milani and D. A. Weitz, Nat. Commun., 2021, 12, 82.
150. X. Shan and H. Chen, Phys. Rev., 1993, 47, 1815.
151. X. Shan and H. Chen, Phys. Rev., 1994, 49, 2941.
152. Q. Kang, D. Zhang and S. Chen, J. Fluid Mech., 2005, 545, 41–66.
153. J. Bao and L. Schaefer, Appl. Math. Model., 2013, 37, 1860–1871.
154. C. Stiles and Y. Xue, Comput. Fluids, 2016, 131, 81–90.
155. M. Sbragaglia, R. Benzi, L. Biferale, S. Succi, K. Sugiyama and F. Toschi, Phys. Rev., 2007, 75, 026702.
156. Z. Guo, C. Zheng and B. Shi, Phys. Rev., 2002, 65, 046308.
157. L. Fei, A. Scagliarini, A. Montessori, M. Lauricella, S. Succi and K. H. Luo, Phys. Rev. Fluids, 2018, 3, 104304.
158. Z.-H. Chai and T.-S. Zhao, Acta Mech. Sin., 2012, 28, 983–992.
159. Q. Li and K. Luo, Appl. Therm. Eng., 2014, 72, 56–61.
160. J. Huang, X. Yin and J. Killough, Phys. Rev. E, 2019, 100, 053304.
161. C. Peng, L. F. Ayala and O. M. Ayala, J. Comput. Phys., 2021, 429, 110018.
162. M. Sbragaglia, H. Chen, X. Shan and S. Succi, Europhys. Lett., 2009, 86, 24005.
163. M. Sbragaglia and X. Shan, Phys. Rev. E, 2011, 84, 036703.
164. R. Benzi, M. Sbragaglia, A. Scagliarini, P. Perlekar, M. Bernaschi, S. Succi and F. Toschi, Soft Matter, 2015, 11, 1271–1280.
165. S. S. Chikatamarla and I. V. Karlin, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2009, 79, 046701.
166. D. P. F. Silva, R. C. V. Coelho, M. M. T. da Gama and N. A. M. Araújo, Phys. Rev. E, 2023, 107, 035106.
167. R. C. V. Coelho, D. P. F. Silva, A. M. R. Maschio, M. M. Telo da Gama and N. A. M. Araújo, Phys. Fluids, 2023, 35, 013304.
168. A. Yiotis, J. Psihogios, M. Kainourgiakis, A. Papaioannou and A. Stubos, Colloids Surf., A, 2007, 300, 35–49.
169. G. Soligo, A. Roccon and A. Soldati, Meccanica, 2020, 55, 371–386.
170. C. D. Eggleton and A. S. Popel, Phys. Fluids, 1998, 10, 1834–1845.
171. M. Shapira and S. Haber, Int. J. Multiphase Flow, 1990, 16, 305–321.
172. A. Vananroye, P. Puyvelde and P. Moldenaers, J. Rheol., 2007, 51, 139–153.
173. G. Taylor, Proc. R. Soc. London, Ser. A, 1932, 138, 41–48.
174. G. Taylor, Proc. R. Soc. London, Ser. A, 1934, 146, 501–523.
175. S. Afkhami, P. Yue and Y. Renardy, Phys. Fluids, 2009, 21, 072106.
176. Z. Gou, F. Huang, X. Ruan and X. Fu, Commun. Comput. Phys., 2018, 24, 234–252.
177. A. Rahmat, J. Meng, D. R. Emerson, C.-Y. Wu, M. Barigou and A. Alexiadis, Microfluid. Nanofluid., 2021, 25, 1.
178. Y. Sui, Y. T. Chew, P. Roy, X. B. Chen and H. T. Low, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2007, 75, 066301.
179. M. Wouters, O. Aouane, T. Krüger and J. Harting, Phys. Rev. E, 2019, 100, 033309.
180. F. Guglietta, M. Behr, L. Biferale, G. Falcucci and M. Sbragaglia, Soft Matter, 2020, 16, 6191–6205.
181. D. Barthès-Biesel, J. Fluid Mech., 1980, 100, 831–853.
182. A. Ben Richou, A. Ambari, M. Lebey and J. Naciri, Chem. Eng. Sci., 2005, 60, 2535–2543.
183. O. H. Faxen, Proc. Roy. Irish Acad., 1946, 187, 1–13.
184. J. Happel and H. Brenner, Low Reynolds number hydrodynamics: with special applications to particulate media, Springer, 1st edn, 1983.
185. E. Montellà, B. Chareyre, S. Salager and A. Gens, MethodsX, 2020, 7, 101090.
186. F. Bonaccorso, M. Lauricella, A. Montessori, G. Amati, M. Bernaschi, F. Spiga, A. Tiribocchi and S. Succi, Comput. Phys. Commun., 2022, 277, 108380.
187. A. Tiribocchi, A. Montessori, G. Amati, M. Bernaschi, F. Bonaccorso, S. Orlandini, S. Succi and M. Lauricella, J. Chem. Phys., 2023, 158, 104101.
188. S. Succi, G. Amati, M. Bernaschi, G. Falcucci, M. Lauricella and A. Montessori, Comput. Fluids, 2019, 181, 107–115.
189. G. Falcucci, G. Amati, P. Fanelli, V. K. Krastev, G. Polverino, M. Porfiri and S. Succi, Nature, 2021, 595, 537–541.
190. J. Latt, O. Malaspinas, D. Kontaxakis, A. Parmigiani, D. Lagrava, F. Brogi, M. B. Belgacem, Y. Thorimbert, S. Leclaire, S. Li, F. Marson, J. Lemus, C. Kotsalos, R. Conradin, C. Coreixas, R. Petkantchin, F. Raynaud, J. Beny and B. Chopard, Comput. Math. Appl., 2021, 81, 334–350.
191. P. Coveney, R. Highfield and V. Ramakrishnan, Virtual You: How Building Your Digital Twin Will Revolutionize Medicine and Change Your Life, Princeton University Press, 2023.

---