Modelling ternary fluids in contact with elastic membranes

TL;DR

This paper introduces a thermodynamically consistent model for simulating ternary fluid interactions with elastic membranes, combining lattice Boltzmann and immersed boundary methods for versatile fluid-structure interaction analysis.

Contribution

A novel integrated modeling framework for ternary fluids and elastic membranes using lattice Boltzmann and immersed boundary methods.

Findings

Validated model with elastic capsule at fluid interface

Achieved Galilean invariance in simulations

Demonstrated versatility with deformable capsule capillary bridge

Abstract

We present a thermodynamically consistent model of a ternary fluid interacting with elastic membranes. Following a free-energy modelling approach and taking into account the thermodynamics laws, we derive the equations governing the ternary fluid flow and dynamics of the membranes. We also provide the numerical framework for simulating such fluid-structure interaction problems. It is based on the lattice Boltzmann method, employed for resolving the evolution equations of the ternary fluid in an Eulerian description, coupled to the immersed boundary method, allowing for the membrane equations of motion to be solved in a Lagrangian system. The configuration of an elastic capsule placed at a fluid-fluid interface is considered for validation purposes. Systematic simulations are performed for a detailed comparison with reference numerical results obtained by Surface Evolver, and the Galilean invariance of the proposed model is also proven. The proposed approach is versatile, and a wide range of geometries can be simulated. To demonstrate this, the problem of a capillary bridge formed between two deformable capsules is investigated here.

Equations

$$\frac{d\left(U+E\right)}{dt}=\frac{dW}{dt}+\frac{dQ}{dt},$$

$$F=U-T{S}_{sys}.$$

$$d{S}_{sur}=-\frac{dQ}{T}.$$

$$\frac{dS}{dt}=\frac{d{S}_{sys}}{dt}+\frac{d{S}_{sur}}{dt}=-\frac{1}{T}\frac{d\left(F+E\right)}{dt}+\frac{1}{T}\frac{dW}{dt}.$$

$$F_f={\int }_V{f}_{f}dV,E_f=\frac{1}{2}{\int }_V\rho |\mathbf{u}{|}^{2}dV,$$

$$f_f=f_b+f_{\nabla }.$$

$$f_{\nabla }=\frac{1}{2}\sum _{m,n=1}^N{c}_{mn}\nabla C_m\cdot \nabla C_n,$$

$$\rho =\sum _{m=1}^N{C}_m{M}_{w,m},$$

$$\frac{d{F}_f}{dt}={\int }_V\frac{\partial f_f}{\partial t}dV+{\int }_V\nabla \cdot \left(f_f\mathbf{u}\right)dV,$$

$$\frac{d{E}_f}{dt}$$

$$\frac{\partial C_m}{\partial t}+\nabla \cdot \left(C_m\mathbf{u}\right)+\nabla \cdot {\mathbf{J}}_m=0,$$

$$\frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho \mathbf{u}\right)+\sum _{m=1}^N{M}_{w,m}\nabla \cdot {\mathbf{J}}_m=0.$$

$$\frac{d{E}_f}{dt}$$

$$\frac{dW}{dt}={\int }_{\partial V}\left({\mathbf{F}}_{sur}\cdot \mathbf{u}\right)ds,$$

$$\frac{dW}{dt}$$

$$-\frac{d\left(F_f+E_f\right)}{dt}+\frac{dW}{dt}$$

$$p=\sum _{m=1}^N{C}_m{\mu }_m-f_f,$$

$${\mu }_m=\frac{\delta f_f}{\delta C_m}={\mu }_m^b-\underset{{\mu }_m^{\nabla }}{\underbrace{\sum _{n=1}^N\nabla \cdot \left(c_{mn}\nabla C_n\right)}},$$

$$p=\underset{p_b}{\underbrace{\sum _{m=1}^N{C}_m{\mu }_m^b-f_b}}-\sum _{m,n=1}^N{C}_m\nabla \cdot \left(c_{mn}\nabla C_n\right)-\frac{1}{2}\sum _{m,n=1}^N{c}_{mn}\nabla C_m\cdot \nabla C_n.$$

$$p_b=\underset{p_{b,f}}{\underbrace{\sum _{m=1}^N{C}_m{\mu }_m^{b,f}-f_{b,f}}}+\underset{p_{b,c}}{\underbrace{C_3\frac{\delta f_{b,c}\left(C_3,\mathcal{I}\right)}{\delta C_3}-f_{b,c}}},$$

$$\nabla p_b$$

$$\nabla p_{b,c}-C_3\nabla \left(\frac{\delta f_{b,c}}{\delta C_3}\right)=\nabla \left(C_3\frac{\delta f_{b,c}}{\delta C_3}-f_{b,c}\right)-C_3\nabla \left(\frac{\delta f_{b,c}}{\delta C_3}\right)=-\frac{\delta f_{b,c}}{\delta \mathcal{I}}\nabla \mathcal{I}.$$

$$\frac{\partial f_b}{\partial t}$$

$$\frac{\partial f_{\nabla }}{\partial t}$$

$$\nabla \cdot \left(f_{\nabla }\mathbf{u}\right)$$

$$\frac{\partial f_f}{\partial t}+\nabla \cdot \left(f_f\mathbf{u}\right)$$

$$\sum _{m,n=1}^N\left(\nabla C_m\right)\nabla \cdot \left(c_{mn}\nabla C_n\right)+\frac{1}{2}\sum _{m,n=1}^N\nabla \left(c_{mn}\nabla C_m\cdot \nabla C_n\right)=\sum _{m,n=1}^N\nabla \cdot \left(c_{mn}\nabla C_m\otimes \nabla C_n\right),$$

$$\frac{\partial f_f}{\partial t}+\nabla \cdot \left(f_f\mathbf{u}\right)$$

$$-\frac{d\left(F_f+E_f\right)}{dt}+\frac{dW}{dt}$$

$$F_m={\mathcal{E}}_s+{\mathcal{E}}_b,$$

Full text

I Introduction

Multi-phase/component flows in interaction with deformable, thin structures are encountered in a broad range of applications. One example of such flows, that is abundant in many biological systems, are single/multi-component flows in which soft particles, composed of an elastic shell enclosing one or more substances, are embedded. For instance, red blood cells are comprised of a thin lipid bilayer membrane containing a cytosol with viscosity about five times higher than that of the surrounding plasma (Freund, 2014). Likewise, viral capsids are formed of polymer DNA confined by a shell of protein molecules (Alberts et al., 2008). In general, capsules accommodating multiple aqueous solutions, such as polyethylene glycol and dextran, serve as model biological cells (Scott Long et al., 2005). Artificially fabricated elastic capsules (Yow & Routh, 2006) are also employed as container and delivery systems in many industrial applications. For example, polymeric multilayer capsules are engineered for drug-delivery purposes (De Cock et al., 2010), while liquid-core capsules are used in various processes of cosmetic (Martins et al., 2014) and food (Gharsallaoui et al., 2007) industries. Another case where a multicomponent flow interacting with soft particles occurs is in the self-assembly of colloidal aggregates into patchy particles, a phenomenon of particular importance for the successful design of bottom-up materials (Glotzer & Solomon, 2007; Pawar & Kretzschmar, 2010). These patchy particles can be formed of compartments, each containing different polymers, bonded together by a solvent. The wetting of liquid droplets surrounded by air/vapour on soft solids is another example of multiphase flows in contact with deformable structures, with relevance in biology, e.g. cell locomotion (Lo et al., 2000), medicine and engineering, e.g. in the identification of cancer cells (Discher et al., 2005), and in the development of smart coatings (Yao et al., 2013). Similar elastowetting problems also arise with fibrous materials, such as the coalescence of wet hair (Bico et al., 2004).

Even though these elastocapillary problems have been extensively studied both analytically and experimentally (Ivanovska et al., 2004; Buenemann & Lenz, 2008; Hegemann et al., 2018; Scott Long et al., 2008; Kusumaatmaja et al., 2009; Kusumaatmaja & Lipowsky, 2011; Meester et al., 2016; Marchand et al., 2012; Style & Dufresne, 2012; Style et al., 2013; Bico et al., 2018; Duprat et al., 2011; Singh et al., 2014; Sauret et al., 2015), the accurate computational modelling of such configurations is lagging behind. Few works have dealt with the simulation of elastic liquid-core capsules immersed in another fluid component (Ramanujan & Pozrikidis, 1998; Zhang et al., 2007; Kaoui et al., 2008). In these works, the ambient and enclosed fluids are assumed to have equal densities and different viscosities, and the fluid-fluid surface tension, if taken into account, takes the form of a force exerted on the membrane, rather than being an inherent quantity of the fluids. There has been an attempt by Lubbers et al. (2014) to numerically solve the equilibrium shape problem of liquid droplets on soft solids. The equilibrium shapes are obtained by minimizing the total elastocapillary energy of the vapour-liquid-substrate system; thus, no fluid dynamics properties of the droplets and surrounding air are available. Shao et al. (2013) have developed a numerical technique to study the contact line dynamics of a two-phase fluid in interaction with rigid, solid circular cylinders. Working towards that direction, Li et al. (2016) have proposed a numerical approach capable of simulating the interaction between a two-phase flow and elastic fibers. In their work, the wetting is controlled by empirical parameters entering the fluid solver through cohesion forces mimicking the role of surface tensions. Recently, Bueno et al. (2018) have developed a numerical framework for the simulation of binary fluids in contact with nonlinear hyperelastic solids to investigate the wetting of soft substrates and elastic micropillars.

The aim of the present work is to analytically derive a thermodynamically consistent model of multicomponent fluids in interaction with elastic membranes, and to provide a versatile numerical framework for the simulation of such fluid-structure interaction problems. We focus here on the case of a ternary fluid where one component is enclosed inside the membranes. The model can, however, be generalised to consider more fluid components contained in/surrounding the membranes. The availability of such computational method will allow us to systematically study a wide range of elastocapillary phenomena intractable to analytical solutions, in order to complement expanding experimental activities in this area.

We follow a top-down modelling approach for the ternary fluid using the lattice Boltzmann method (Swift et al., 1995, 1996; Briant et al., 2004; Briant & Yeomans, 2004; Lee & Lin, 2005; Mazloomi M. et al., 2015; Semprebon et al., 2016a; Wöhrwag et al., 2018); namely, the free energy of the fluid system is initially formulated, including the desired thermodynamics features, such as the immiscibility of the fluid mixture, the surface tensions between the different fluid components, and the energy due to the interaction of the elastic membranes with the confined component. Given the free energy, other thermodynamic properties, e.g. the chemical potentials of each component, and the pressure of the fluid mixture, as well as the macroscopic equations of motion of the multicomponent fluid can be derived. This technique is the contrary of the bottom-up modelling approach, where the macroscopic properties of the fluid arise from the microscopic interactions between the fluid particles. These interactions usually take the form of an interparticle potential (Shan & Chen, 1993; Chen et al., 2014). Krüger et al. (2017) have discussed in detail the advantages and limitations of each of these modelling approaches. The lattice Boltzmann model proposed by Semprebon et al. (2016a) has been in particular chosen here for resolving the flow in a uniform, Cartesian grid. This is a diffuse interface model, meaning that the fluid-fluid interfaces are spread over several lattices, rather than being tracked explicitly.

The capsules enclosing one fluid component are modelled as infinitely thin membranes, composed of a homogeneous and isotropic material, able to undergo stretching/compression and bending. The response of the membranes to the enclosing and surrounding fluid components is governed by constitutive equations for the strain and bending, whose formulations depend on the nature of the material. Here, we consider the simple case of a linear elastic material, for which the strain energy takes the form of Hooke’s law, and the bending energy is given by the squared mean curvature. However, the strain and bending energy formulations can be easily modified to model more realistic materials. For example, the neo-Hookean and Mooney-Rivlin laws are employed for hyperelastic rubber-like materials, while the model of Skalak et al. (1974) is widely used for red blood cells.

To resolve the interaction between the ternary fluid and the membranes, we adopt the immersed boundary method (IBM), initially proposed by Peskin (1977). While the equations of motion of the ternary fluid are solved in an Eulerian description, a Lagrangian coordinates system is used for the solution of the membranes’ constitutive laws. The membranes are discretised into points, referred to as Lagrangian markers, connected by linear tension/compression and torsion springs, whose coefficients can be related to the tensile and bending moduli of the desired material. The Lagrangian markers are allowed to move in order to follow the motion/deformation of the membranes’ boundary, without having to coincide with the underlying Eulerian fluid lattices. The desired strain and bending conditions of the membranes are then imposed on the Lagrangian markers. By transferring information between the two coordinates systems, the aforementioned Lagrangian conditions take the form of a forcing term, reproducing their effect on the Eulerian flow, in the lattice Boltzmann equation. The information transfer is achieved by interpolation and spreading operations.

The article is organised as follows. In §II, we derive a thermodynamically consistent model of a ternary fluid in contact with elastic membranes. The methods employed here for the numerical solution of the ternary fluid and its interacting membranes equations of motion are subsequently presented in §III.1 and §III.3. The Surface Evolver, an open-source software used for benchmarking purposes, is briefly discussed in §III.4. We then validate our model against Surface Evolver in §IV.1, considering an elastic capsule placed at a fluid-fluid interface as the benchmark configuration. We also prove the Galilean invariance of the equations governing the ternary fluid flow and the membranes’ dynamics. To demonstrate the versatility of our model, the problem of a capillary bridge formed between two elastic capsules is afterwards studied in §IV.2. It should be noted that although the simulations performed here are restricted to two dimensions, an extension to three-dimensional configurations is straightforward. Finally, the key contributions of the present work and its future perspectives are summarized in §V.

II Mathematical modelling of ternary fluids in interaction with elastic membranes

In this section, the equations of motion for the ternary fluid interacting with elastic membranes are derived. Following the rationale of Kou & Sun (2018), we firstly derive an entropy conservation equation for the ternary fluid and the membrane by making use of the first law of thermodynamics and splitting structure of the total entropy into contributions from the system of interest and its surroundings. To further reduce the entropy conservation equation, we then derive a transport equation of the free energy density of the fluid mixture. We subsequently formulate the free and kinetic energies of the elastic membrane. Taking into account the aforementioned and the second law of thermodynamics, we obtain the equations of motion of the ternary fluid and the elastic membrane. We then present briefly the free energy of the ternary fluid, which is based on the one proposed by Semprebon et al. (2016a) but modified accordingly in order to account for the interplay between the fluid and the membrane. By introducing auxiliary fields that allow for an easy numerical implementation of the equations of motion of the ternary fluid, the equivalent formulations of the continuity and Navier-Stokes equations are obtained along with the Cahn-Hilliard equations describing the evolution of these auxiliary fields.

II.1 Entropy equation

The first law of thermodynamics can be formulated as

[TABLE]

where $t$ is time, $U$ and $E$ are respectively the internal and kinetic energies, $W$ is the work done on the ternary fluid system by its surroundings, and $Q$ denotes the amount of heat supplied to the system for its temperature $T$ to be kept constant. The total entropy $S$ can be split into two parts: the entropy of the system $S_{sys}$, and the entropy of its surroundings $S_{sur}$. By definition, the total free energy of the system can be expressed as

[TABLE]

It is also known that the entropy of the system’s surroundings is related to the heat $Q$ by

[TABLE]

Taking into account (1)–(3), it can be shown that

[TABLE]

Here, the system of interest consists of the ternary fluid (denoted by subscript $f$) and the membrane (subscript $m$). Let us first concentrate on the terms related to the ternary fluid.

Given the simulation volume $V$, the free and kinetic energies of the ternary fluid can be defined as

[TABLE]

where $f_{f}$ is the free energy density of the ternary fluid, $\rho$ is the mass density of the fluid mixture, and $\mathbf{u}$ denotes the mass-averaged velocity. The fluid free energy density is the sum of two contributions: a free energy density of the bulk fluid $f_{b}$, and a local free energy density gradient contribution $f_{\nabla}$ allowing the existence of diffuse fluid-fluid interfaces

[TABLE]

The interfacial free energy density term is defined here as

[TABLE]

where $N$ denotes the number of components of the fluid mixture (thus, $N=3$ here), $c_{mn}$ is the cross influence parameter, and $C_{m}$ represents the concentration fraction of fluid $m$ $\left(m=1,2,3\right)$. The mass density is given by

[TABLE]

where $M_{w,m}$ denotes the weight of component $m$.

By applying the Reynolds transport and Gauss divergence theorems to (5), we can obtain

[TABLE]

[TABLE]

The mass transfer occurs through the convection and diffusion of each fluid component, when a velocity field is present. The mass conservation law for each component $m$ can then be written as

[TABLE]

where $\mathbf{J}_{m}=M_{m}\nabla\mu_{m}$ denotes the diffusion flux of component $m$, and $M_{m}$ represent the corresponding mobility parameters. The chemical potential $\mu_{m}$ is discussed in §II.2. By multiplying (11) with $M_{w,m}$, and summing over all fluid components $m$, the mass conservation equation of the fluid mixture can be obtained as

[TABLE]

By substituting (12) into (10), the latter can be rewritten as

[TABLE]

where $\frac{d\mathbf{u}}{dt}$ is the total derivative of $\mathbf{u}$, defined as $\frac{\partial\mathbf{u}}{\partial t}+\mathbf{u}\cdot\nabla\mathbf{u}$.

The rate of change of the work done by the force $\mathbf{F}_{sur}$ is given by

[TABLE]

where the integration takes place over the volume’s surface boundary $\partial V$. The force $\mathbf{F}_{sur}$ is related to the Cauchy stress tensor $\mathbf{\sigma}$ of the ternary fluid by $\mathbf{F}_{sur}=-\mathbf{\sigma}\cdot\mathbf{n}$, where $\mathbf{n}$ denotes the outward unit normal vector of $V$. Equation (14) can thus take the form

[TABLE]

By making use of (9), (13) and (15), the terms on the right-hand side of (4) related to the ternary fluid can be collected as follows

[TABLE]

II.2 Transport equation of the fluid free energy density

The pressure $p$ of the fluid mixture can be defined as

[TABLE]

where $\mu_{m}$ is the chemical potential of component $m$, given by

[TABLE]

where $\mu_{m}^{b}=\frac{\delta f_{b}}{\delta C_{m}}$ and $\mu_{m}^{\nabla}=\frac{\delta f_{\nabla}}{\delta C_{m}}$ stand for the bulk and interfacial chemical potentials of component $m$. Due to (18), equation (17) can be rewritten as

[TABLE]

The bulk pressure $p_{b}$ can be further divided into two parts: $p_{b,f}$ accounting for the contribution of the bulk free energy density owing to the chosen free-energy functional $f_{b,f}$ of the fluid mixture, and $p_{b,c}$ taking into account the contribution of the bulk free energy density due to a fluid-membrane coupling-energy functional $f_{b,c}$

[TABLE]

where $f_{b}=f_{b,f}+f_{b,c}$ and $\mu_{m}^{b,f}=\frac{\delta f_{b,f}\left(C_{m}\right)}{\delta C_{m}}$. The coupling-energy functional $f_{b,c}$ depends only on the component enclosed in the membrane, which is assumed to be fluid $3$ in this work, and an interfacial profile $\mathcal{I}$ across the elastic membrane for the ternary fluid. The latter will be discussed in more detail in §II.5. The gradient of the bulk pressure can be found as

[TABLE]

It can also be shown that

[TABLE]

Taking into account (7), (11), (21) and (22), the time derivative of the bulk $f_{b}$ and interfacial $f_{\nabla}$ free energy densities can be deduced as

[TABLE]

and

[TABLE]

The divergence of the interfacial free energy density in the presence of a velocity field is given by

[TABLE]

Making use of (6), (18), (19) and (23)–(II.2), the transport equation of the free energy density $f_{f}$ of the ternary fluid reads

[TABLE]

Taking into account that the following identity holds

[TABLE]

equation (26) becomes

[TABLE]

By substituting (27) into (II.1), the latter is reformulated as

[TABLE]

where $\mathbf{I}$ is the second-order identity tensor.

II.3 Free and kinetic energies of the elastic membrane

Let us focus now on the terms of (4) related to the membrane. The free energy $F_{m}$ of the membrane is the sum of two contributions: strain energy $\mathcal{E}_{s}$ and bending energy $\mathcal{E}_{b}$,

[TABLE]

where

[TABLE]

[TABLE]

The integration takes place over the surface $S$ of the elastic membrane, whose position in Eulerian coordinates at time $t$ is described by $\mathbf{X}=\mathbf{X}\left(\mathbf{s},t\right)$, where $\mathbf{s}$ denotes its Lagrangian coordinates. The parameters $\kappa_{s}$ and $\kappa_{b}$ are, respectively, the Young’s and bending moduli. The variable $\kappa=\partial\theta/\partial\mathbf{s}^{\prime}$, where $\theta=\theta\left(\mathbf{X}\right)$ is the tangential angle of an arc of length $\mathbf{s}^{\prime}=\mathbf{s}^{\prime}\left(\mathbf{X}\right)$, denotes the curvature of the surface $S$. The kinetic energy $E_{m}$ of the membrane is given by

[TABLE]

where $m$ is the mass of the membrane, here set to be equal to $1$.

II.4 Equations of motion

By combining equations (II.2) and (29)–(32), the entropy equation (4) can be rewritten as

[TABLE]

Here, we have used the natural boundary conditions

[TABLE]

where $\mathbf{n}_{\partial V}$ denotes the outward unit vector normal to the surface boundary $\partial V$ when integrating over the volume $\Omega$ of the fluid mixture.

The stress tensor $\mathbf{\sigma}$ of the ternary fluid can be considered as the sum of two contributions: a reversible stress tensor $\mathbf{\sigma}_{rev}$, and an irreversible one $\mathbf{\sigma}_{irrev}$

[TABLE]

In an ideal reversible process, there are no effects of viscosity and friction, the diffusion fluxes $\mathbf{J}_{m}$ are equal to [math], and the entropy is conserved. Taking these facts into consideration, equation (33) results in

[TABLE]

For an incompressible viscous Newtonian fluid, the irreversible stress tensor can be expressed as

[TABLE]

where $\eta$ is the dynamic viscosity.

According to the second law of thermodynamics, the total entropy cannot decrease over time. This, in combination with the non-negative nature of the third term on the right-hand side of (33), implies that

[TABLE]

In the present work, we consider that the weight of each fluid component is equal to $1$, namely $M_{w,m}=1,\ m=1,2,3$, and $\sum_{m=1}^{N}M_{w,m}\mathbf{J}_{m}=0$. As such, the continuity (12) and the Navier-Stokes (38) equations become

[TABLE]

and

[TABLE]

The equation of motion of the elastic membrane can be derived by considering the last four terms on the right-hand side of (33). Setting their sum to zero and writing the equation for a discretized membrane for concreteness,

[TABLE]

The variable $\mathbf{X}_{l}$ denotes the position of the $l^{th}$ Lagrangian marker (in Eulerian coordinates), with $l=1,\dots,N_{l}$. Taking the expressions in the square bracket to be zero for each Lagrangian marker, we obtain

[TABLE]

where $\mathbf{F}_{mem}$ denotes the total force exerted on the elastic membrane, $\mathbf{F}_{s,l}=-\frac{\partial\mathcal{E}_{s}}{\partial\mathbf{X}_{l}}$ and $\mathbf{F}_{b,l}=-\frac{\partial\mathcal{E}_{b}}{\partial\mathbf{X}_{l}}$ are the strain and bending forces, and $\mathbf{F}_{c,l}=-\int_{V}\frac{\delta f_{b,c}}{\delta\mathcal{I}}\frac{\partial\mathcal{I}}{\partial\mathbf{X}_{l}}\ dV$ is the coupling force. The detailed forms of the forces are provided in §III.2.

II.5 Free energy of the ternary fluid

As mentioned earlier, the free energy $F_{f}$ of the ternary fluid system considered here is the sum of two contributions: a Landau free-energy functional $\mathcal{E}_{f}$ allowing the coexistence of three fluid components, and a coupling-energy functional $\mathcal{E}_{c}$ taking into account the interaction between the membrane and its confined fluid component

[TABLE]

In this work, we have used the following forms for $\mathcal{E}_{f}$ and $\mathcal{E}_{c}$

[TABLE]

[TABLE]

The first term on the right-hand side of (44) corresponds to the ideal gas term, as it is commonly used in multicomponent free energy lattice Boltzmann models (Krüger et al., 2017).

Each concentration fraction $C_{m}$ has two bulk minima at $C_{m}=\left\{0,1\right\}$. For the ternary fluid system of interest, only the following three minimizers are considered:

[TABLE]

corresponding to the three bulk fluids in the ternary system. As our multiphase fluid model assumes diffuse fluid-fluid interfaces, a formulation of the interfacial profile for the concentration of each component $m$ is required. Here it takes the form

[TABLE]

where $d_{FI}$ measures the distance between the bulk fluid at position $\mathbf{x}$ and the fluid-fluid interface. Equation (47) ensures that $C_{m}\rightarrow 1$ for $d_{FI}\rightarrow\infty$, and $C_{m}\rightarrow 0$ for $d_{FI}\rightarrow-\infty$. The parameter $\alpha$ is proportional to the interface width, which here is the same for all three fluid-fluid interfaces. The fluid-fluid surface tensions are expressed as

[TABLE]

The coupling energy is formulated in such a way that for $\kappa_{c}>0$ there is a minimum when $C_{3}=\mathcal{I}$. The interfacial profile $\mathcal{I}$ across the elastic membrane for the ternary fluid is then defined as

[TABLE]

where $d\left(\mathbf{x},\mathbf{X}\right)$ denotes the distance between each bulk fluid component $m$ at position $\mathbf{x}$ and the membrane located at $\mathbf{X}$. This distance is assigned to be positive for the enclosed fluid component $3$, and negative for the surrounding fluid components $1$ and $2$. The width of the fluid-membrane interfaces is kept the same as the fluid-fluid interfaces one. It is obvious that with an increasing coupling coefficient $\kappa_{c}$, the interfacial profile for the concentration $C_{3}$ of the enclosed fluid component will be superimposed onto the profile of the fluid $3$-membrane interface. However, if $\kappa_{c}$ is too high, the coupling energy term will dominate over the rest of the free energy terms, which is undesirable.

By comparing (7) with the last term on the right-hand side of (44), it is obvious that

[TABLE]

where $\delta_{m,n}$ is the Kronecker delta. Taking into account (50), it can be shown that

[TABLE]

II.6 Coordinates Transformation

For convenience, we shall express the equations of motion of the ternary fluid in terms of the mass density $\rho$, and the auxiliary fields $\phi$ and $\psi$, given by

[TABLE]

Here we assume that all fluid components have the same density, and thus the mass density is set to be $\rho=1$. Taking the aforementioned into account, the Navier-Stokes (40) equations become

[TABLE]

where $\mu_{\rho}$, $\mu_{\phi}$ and $\mu_{\psi}$ are the chemical potentials coming from the chosen free energy density of the ternary fluid. The detailed formulations of the chemical potentials $\mu_{\rho}$, $\mu_{\phi}$ and $\mu_{\psi}$ in terms of the variables $\rho$, $\phi$ and $\psi$ are given in the appendix A. The last term on the right-hand side of (53) differentiates the present model from the ternary fluid model proposed by Semprebon et al. (2016a).

Equation (11) can be transformed into the Cahn-Hilliard equations governing the evolution of the order parameters $\phi$ and $\psi$ (Semprebon et al., 2016a)

[TABLE]

[TABLE]

where $M_{\phi}$ and $M_{\psi}$ are the mobility parameters, and $\mu_{\psi}^{\prime}=\mu_{\psi}+\frac{\delta f_{b,c}}{\delta\psi}$. Note that the corresponding diffusive term on the right-hand side of (55) in the ternary fluid model proposed by Semprebon et al. (2016a) does not include the contribution owing to the coupling-energy functional, $\frac{\delta f_{b,c}}{\delta\psi}$. Here we assume that the variables $C_{m}\left(m=1,2,3\right)$ have identical mobility parameters, resulting in $M_{\phi}=3M_{\psi}$.

III Numerical methods

In this section, the numerical techniques employed for the solution of the equations of motion of the ternary fluid and its interacting elastic membrane are discussed. The governing equations of the ternary fluid are solved numerically by means of a lattice Boltzmann method, as detailed in §III.1. The discretized forms of the strain, bending and coupling energies are presented in §III.2 along with the corresponding force formulations. The interaction between the ternary fluid and elastic membrane is solved by an immersed boundary method, described in §III.3, that is coupled to the lattice Boltzmann method following the algorithm presented in Krüger et al. (2011). Finally, we briefly report the equivalent energies implementation in Surface Evolver, a software used for benchmarking purposes, in §III.4.

To clarify notations, lower case letters are employed in the following for variables evaluated on the Eulerian lattices, while upper case ones refer to variables defined at the Lagrangian markers.

III.1 Lattice Boltzmann method

To solve the equations of motion of the ternary fluid (39) and (53)–(55), we employ the lattice Boltzmann method with three sets of distribution functions $f_{i}\left(\mathbf{x},t\right)$, $g_{i}\left(\mathbf{x},t\right)$, and $h_{i}\left(\mathbf{x},t\right)$, corresponding to the fluid density $\rho$ and the order parameters $\phi$ and $\psi$. The evolutions of the distribution functions are governed by the lattice Boltzmann equation, where the standard Bhatnagar-Gross-Krook (BGK) single relaxation time model (Bhatnagar et al., 1954) is used for the collision operator, and the exact difference scheme (Kupershtokh, 2004) is employed for the forcing term

[TABLE]

The variables $f_{i}\left(\mathbf{x},t\right)$, $g_{i}\left(\mathbf{x},t\right)$, and $h_{i}\left(\mathbf{x},t\right)$ refer to the distribution functions $f_{i}$, $g_{i}$, and $h_{i}$ at position $\mathbf{x}$ and time $t$ with velocity $\mathbf{c}_{i}$ along the $i^{th}$ lattice direction. The relaxation times $\tau$, $\tau_{\phi}$ and $\tau_{\psi}$ are linked to the dynamic viscosity $\eta$, and the mobility parameters $M_{\phi}$ and $M_{\psi}$ by

[TABLE]

where $c_{s}$ is the speed of sound, and $\Gamma_{\phi}$, $\Gamma_{\psi}$ are tunable parameters. The speed of sound is given by $c_{s}=\frac{1}{\sqrt{3}}c$, where $c=\frac{\delta x}{\delta t}$ is the lattice speed, and $\delta x$, $\delta t$ are the lattice spacing and time step, respectively. The variables $f_{i}^{eq}$, $g_{i}^{eq}$, and $h_{i}^{eq}$ denote the equilibrium distribution functions.

The latter are expressed as

[TABLE]

where $w_{i}$ are weight coefficients depending on the chosen lattice arrangement for the velocity discretisation, and $\mathbf{I}$ is the identity tensor. As mentioned earlier, the chemical potential $\mu_{\psi}^{\prime}$ is given by $\mu_{\psi}^{\prime}=\mu_{\psi}+\frac{\delta f_{b,c}}{\delta\psi}$, where it is obvious from (45) that $\frac{\delta f_{b,c}}{\delta\psi}=\kappa_{c}\left(\psi-\mathcal{I}\right)$. The formulations of the chemical potentials $\mu_{\phi}$ and $\mu_{\psi}$ are provided in the appendix A.

The macroscopic physical variables are defined as moments of the distribution functions

[TABLE]

The variable $\mathbf{u}$ represents the bare fluid velocity, and it is related to the actual fluid velocity $\mathbf{v}$ by

[TABLE]

where $\delta\mathbf{u}$ denotes the velocity correction given by

[TABLE]

The forcing term $\mathbf{f}$ can be considered as the sum of three contributions: a force $\mathbf{f}_{FE}$ taking into account the gradient terms on the right-hand side of (53), a force $\mathbf{f}_{IB}$ accounting for the interaction between the ternary fluid and elastic membrane, and a force $\mathbf{f}_{ext}$ allowing the existence of external forces

[TABLE]

The force $\mathbf{f}_{FE}$ can be written as

[TABLE]

where the pressure tensor $\mathbf{p}$ satisfies the condition

[TABLE]

The explicit form of $\mathbf{p}$ is provided in the appendix A. The form of the force $\mathbf{f}_{IB}$ is discussed in §III.3. In the present work, no external forces are considered, that is $\mathbf{f}_{ext}=\mathbf{0}$.

In the current normalisation, the lattice spacing is set equal to the time step, $\delta x=\delta t=1$, resulting in $c=1$ and $c_{s}=1/\sqrt{3}$. Both the relaxation times and the parameters $\Gamma_{\phi}$, $\Gamma_{\psi}$ are all considered to be equal to $1$, $\tau=\tau_{\phi}=\tau_{\psi}=1$ and $\Gamma_{\phi}=\Gamma_{\psi}=1$. Here, we employ the $D2Q9$ lattice arrangement, for which the lattice velocities are defined as

[TABLE]

and the weight coefficients are given by

[TABLE]

III.2 Membrane dynamics

The discretised strain and bending energies of the elastic membrane, given respectively by (30) and (31) in its continuous-space forms, are formulated as

[TABLE]

[TABLE]

The summation occurs over all the Lagrangian markers $l=1,\dots,N_{l}$. The variable $\mathbf{X}_{l}$ denotes the position of the $l^{th}$ Lagrangian marker (in Eulerian coordinates), and $\Delta s$ represents the initial distance between two consecutive Lagrangian markers. The membrane is initially discretized into Lagrangian markers such that $\Delta s=\delta x=1$. The tangential angle $\Delta\theta$ and arc length $\Delta s^{\prime}$ can be expressed as

[TABLE]

[TABLE]

Taking into account (78)–(81), the corresponding discretised strain and bending forces can then be found as

[TABLE]

and

[TABLE]

where

[TABLE]

[TABLE]

The $\delta_{m,l}$ is the Kronecker delta. The variables $\mathbf{F}_{s,l}$ and $\mathbf{F}_{b,l}$ denote the discretised strain and bending forces exerted on the $l^{th}$ Lagrangian marker.

The discretised coupling energy of the elastic membrane can be written as

[TABLE]

where $\sum_{\mathbf{x}}$ implies summation over all the Eulerian lattices $\mathbf{x}$, and $d$ is the domain dimensionality. The interfacial profile $\mathcal{I}$ is given by (49). The corresponding discretised coupling force can then be obtained as

[TABLE]

The variable $\mathbf{F}_{c,l}$ represents the discretised coupling force exerted on the $l^{th}$ Lagrangian marker.

The discretised total force acting on the elastic membrane can thus be written as

[TABLE]

In the computational implementation, we exert the forces $\mathbf{F}_{s,l}$, $\mathbf{F}_{b,l}$ and $\mathbf{F}_{c,l}$ on each Lagrangian marker $l$.

III.3 Immersed boundary method

To reproduce the effect of $\mathbf{F}_{mem,l}$ on the Eulerian fluid flow, denoted here by $\mathbf{f}_{IB}$, a spreading operation is used

[TABLE]

The term $\delta_{h}$ denotes the discretised Dirac delta function, and the following formulation proposed by Peskin (2002) is chosen here to perform the convolution in the spreading and interpolation operations

[TABLE]

If the variable r is a vector, $\mathbf{r}=\left(r_{x},r_{y},r_{z}\right)$, then the multidimensional $\delta_{h}$ is given by $\delta_{h}\left(\mathbf{r}\right)=\frac{1}{\delta x^{3}}\delta_{h}\left(\frac{r_{x}}{\delta x}\right)\delta_{h}\left(\frac{r_{y}}{\delta x}\right)\delta_{h}\left(\frac{r_{z}}{\delta x}\right)$. Once the force $\mathbf{f}_{IB}$ is computed, the density $\rho$, order parameters $\phi$ and $\psi$, and velocity $\mathbf{v}$ fields can be obtained at the next time step $t+\delta t$ by solving (56). To calculate the forces $\mathbf{F}_{mem,l}$ at $t+\delta t$, the position of the Lagrangian markers $\mathbf{X}_{l},\ l=1,\dots,N_{l}$ at $t+\delta t$ needs to be known. For this reason, the known $\mathbf{v}\left(\mathbf{x},t+\delta t\right)$ is interpolated at the Lagrangian markers as

[TABLE]

The updated position of the Lagrangian markers can then be found by Euler’s rule

[TABLE]

III.4 Surface Evolver

In the absence of closed-form solutions for non-trivial elastocapillary problems, we have benchmarked the proposed dynamic simulation method against a finite element approach. In particular, the shapes of elastic capsules placed at a fluid-fluid interface in mechanical equilibrium have been modelled. For the finite element approach we employ Surface Evolver, an open-source software developed by Brakke (1992), which has been extensively used to model the equilibrium shapes of liquid interfaces and capillary forces (Grzybowski et al., 2001; Collicott & Weislogel, 2004; Semprebon et al., 2016b), detect energy barriers of morphological transitions (Bi et al., 2015; Semprebon et al., 2018), and model deformations of elastic membranes (Manyuhina et al., 2010).

In Surface Evolver the interfaces are discretized by triangulated meshes constituted of vertices, edges and facets. Surface Evolver provides a large library of pre-implemented energy functionals, including surface tension, elastic strain and bending energies, from which the forces acting on each vertex are automatically computed within the internal routines. Local and global constraints can also be implemented in Surface Evolver, such as to maintain volume/area conservation and to keep certain vertices fixed. Configurations in mechanical equilibrium correspond to minima of the total energy, obtained through a conjugate gradient descent method.

For our benchmarks, we employ a 2D model in Surface Evolver, where the elastic capsules are initialized with exactly the same geometry as in the proposed coupled lattice Boltzmann-immersed boundary method (LB-IBM). The same free energy of the elastic membrane is implemented using the provided scripting language to formulate the strain and bending energies as in (78) and (79). In particular, the functions edge_length and sqcurve_string_marked are employed for the calculation of the strain and bending energies. The main difference with the LB-IBM is in the definition of the surface tensions; here, the corresponding energy is simply accounted for by adding a term proportional to the total length of each fluid or membrane interface, multiplied by a constant parameter matching the surface tension arising from the diffuse interface in the lattice Boltzmann method. No coupling energy has been considered, as both the elastic energies and surface tension are provided by the discrete representation of the capsule and the fluid-fluid interface. The constraint of conservation of the total capsule area has also been considered. Finally, the position of the vertices at the end of the computational domain has been kept fixed throughout the minimization procedure, as an equivalent to implementing periodic boundary conditions in the LB-IBM.

IV Results

The proposed fluid-structure solver is validated in §IV.1. The steady configuration of an elastic capsule positioned at a fluid-fluid interface is chosen for this purpose. We perform a thorough comparison with the reference results of Surface Evolver for different cases of surface tension ratios, and various combinations of Young’s and bending moduli. We also establish the Galilean invariance of the governing equations of the ternary fluid and elastic structure. Finally, a more complex configuration is considered in §IV.2 to demonstrate the capabilities of the proposed model.

IV.1 Elastic capsule at fluid-fluid interface

To benchmark the proposed model, the configuration of an elastic capsule placed at a fluid-fluid interface, as shown in figure 1$\left(a\right)$, is studied. An initially circular capsule of radius $R=20$ is located at the center of a computational domain of dimensions $12R\times 6R$. The capsule relaxes to a mechanical equilibrium shape, which depends on the balance of the elastic strain and bending forces and the surface tensions $\gamma_{12}$, $\gamma_{13}$ and $\gamma_{23}$. We can define a dimensionless deformation parameter, the Taylor deformation $D$, given by $D=\left(L-B\right)/\left(L+B\right)$, where $L$ and $B$ are, respectively, the major and minor axes of the final elliptical capsule shape. Periodic boundary conditions are considered at all the domain boundaries. Simulations are performed for several combinations of Young’s and bending moduli, $\kappa_{s}=\left\{10^{-3},10^{-2},10^{-1}\right\}$ and $\kappa_{b}=\left\{10^{-4},10^{-2}\right\}$, and various surface tensions. The parameter $\alpha$ in (48) is set to be $\alpha=2$. The coupling coefficient is taken as $\kappa_{c}=10^{-2}$, which was found to be the threshold before the domination of the coupling energy over the rest of the free energy terms occurs. This value of the coupling coefficient also ensured that the interfacial profile for the concentration $C_{3}$ of the enclosed fluid component is superimposed onto the profile of the fluid 3-membrane interface.

We first consider the case where $\gamma_{13}=\gamma_{23}$, resulting in a capsule shape that is symmetrical along the domain centreline. We will refer thereafter to this case as the symmetric case. For the surface tensions $\gamma_{13}$ and $\gamma_{23}$ to be equal, the values of the parameters $\kappa_{1}$ and $\kappa_{2}$ in (48) have to be the same, taken here as $\kappa_{1}=\kappa_{2}=8\cdot 10^{-3}$. The parameter $\kappa_{3}$ ranges from $5\cdot 10^{-4}$ to $10^{-2}$ with a step size of $5\cdot 10^{-4}$. The corresponding surface tension ratio $\gamma_{12}/\gamma_{13}=\gamma_{12}/\gamma_{23}$ varies from $0.88\bar{8}$ to $1.882$. Figure 2$\left(a\right)$ shows the Taylor deformations obtained by our simulations, whose results are represented by lines, compared to the ones measured by Surface Evolver, denoted by dot symbols. As expected, the variations in the Taylor deformations found for a particular combination of $\kappa_{s}$ and $\kappa_{b}$ become more apparent with a decrease in the Young’s modulus. By lowering the latter, the results for the elastic capsule tend to the ones obtained for the pure liquid lens configuration, depicted by filled square symbols. It is also obvious that the effect of the bending coefficient on the capsule’s shape is negligible for a given Young’s modulus.

Our results generally agree well with those of Surface Evolver, with a typical relative error in $D$ of $5.6\%$. The discrepancies in the Taylor deformation between the two numerical methods increase for highly deformable capsules, that is $\kappa_{s}\leq 10^{-3}$, subjected to high surface tension ratios. We find that this is due to spurious forces at the three-phase contact points. An increase in the interface width $\alpha$ and the computational domain size could help in the damping of these spurious effects, with an increase, however, in the corresponding computational cost. In any case, the relative errors in $D$ between the results of our simulations and Surface Evolver are less than $10\%$, a reasonable limit given the significant differences in the two solution techniques.

To illustrate the effect of these discrepancies in the Taylor deformation on the mechanical equilibrium shape of the elastic capsule, the latter is plotted for different combinations of $\kappa_{s}$ and $\kappa_{b}$ at the highest surface tension ratio studied here, that is $\gamma_{12}/\gamma_{13}=1.882$, as shown in figure 3$\left(a\right)$. For stiff and moderately deformable capsules, corresponding to $\kappa_{s}=10^{-1}$ and $\kappa_{s}=10^{-2}$, an excellent agreement is observed between our simulation results and those of Surface Evolver. For highly deformable capsules, a slightly elongated shape is obtained by the reference software compared to the one found by our fluid-structure solver. The relative error in $L$ is measured to be $5.75\%$. The capsule width agrees, however, perfectly.

We subsequently investigate the case where $\gamma_{13}\neq\gamma_{23}$, resulting in the capsule to be more immersed in one of the surrounding fluid phases. This case will be mentioned in the following as the asymmetric case. We assign the parameters $\kappa_{1}$ and $\kappa_{2}$ to be $\kappa_{1}=8\cdot 10^{-3}$ and $\kappa_{2}=2\cdot 10^{-3}$. The parameter $\kappa_{3}$ varies from $5\cdot 10^{-4}$ to $8\cdot 10^{-3}$ in steps of $5\cdot 10^{-4}$, with the corresponding surface tension ratio $\gamma_{13}/\gamma_{23}$ ranging from $1.60$ to $3.40$. The dependence of the measured Taylor deformations on the aforementioned surface tension ratio is presented for the LB-IBM and Surface Evolver results in figure 2$\left(b\right)$. Similar observations to the symmetric case can be made on the trends of $D$ for different Young’s and bending moduli. The plots of the mechanical equilibrium capsule shapes are depicted for different combinations of $\kappa_{s}$ and $\kappa_{b}$ at the highest surface tension ratio achieved here, namely $\gamma_{13}/\gamma_{23}=3.40$, in figure 3$\left(b\right)$. We can demonstrate again an excellent agreement between the results of our algorithm and Surface Evolver.

We finally prove the Galilean invariance of the proposed model. To do so, simulations are performed in an inertial frame of reference, and their results are compared to the ones obtained by the previous stationary simulations. At $t=0$, all the fluid components are given a constant horizontal velocity $U_{x,0}$. Due to the fluid-structure interaction, the initially circular capsule travels in the same direction with an equal velocity. Three different velocity values are tested here, $U_{x,0}=\left\{10^{-4},10^{-3},10^{-2}\right\}$. The Galilean invariance is checked in both the symmetric and asymmetric cases for $\kappa_{s}=10^{-3}$ and $\kappa_{b}=10^{-2}$ at the highest surface tension ratios examined respectively here. The shapes of the elastic capsule, after the latter reaches mechanical equilibrium in the inertial reference frame, found for the different velocities are compared to each other and to the ones of the corresponding stationary simulations. These comparisons can be seen in figure 4. In both cases, the results for the capsule shapes in the moving and stationary frames are perfectly superimposed, indicating that the proposed model is Galilean invariant.

IV.2 Capillary bridge between two elastic capsules

To show the capabilities of our model, the configuration of a capillary bridge formed between two deformable capsules, as shown in figure 1$\left(b\right)$, is now investigated. Two initially circular capsules of radius $R=20$ are placed at $\left(x_{c_{1}},y_{c_{1}}\right)=\left(9R/2,3R\right)$ and $\left(x_{c_{2}},y_{c_{2}}\right)=\left(15R/2,3R\right)$ of a computational domain of dimensions $12R\times 6R$. The capillary bridge, composed of the fluid component $2$, is initialized as a rectangular area of dimensions $2S\times H_{b}=R\times 31R/10$ located at the center of the computational domain. Both the capillary bridge and elastic capsules are surrounded by the fluid component $1$. Due to the presence of the surface tensions $\gamma_{12}$, $\gamma_{13}$ and $\gamma_{23}$ as well as the elastic strain and bending forces, the capsules relax to a deformed mechanical equilibrium shape depending on the balance of the aforementioned forces. This mechanical equilibrium shape depends also on the initial distance between the two capsules, $2S$, and the volume of the capillary bridge. To confine the parameter space of the current study, all the simulations are performed only for a dimensionless initial distance $S^{\prime}=2S/R=1$, and a relative area $A_{rel}=A_{b}/A_{c}\approx 1.02$, where $A_{b}$ and $A_{c}$ are the areas of the capillary bridge and elastic capsules. Three different surface tension ratios are tested here, $\gamma_{12}/\gamma_{13}\approx\left\{0.57,0.67,0.80\right\}$. The respective parameters in (48) are set to be $\alpha=2$, $\kappa_{1}=\left\{1\cdot 10^{-3},3\cdot 10^{-3},9\cdot 10^{-3}\right\}$, $\kappa_{2}=3\cdot 10^{-3}$ and $\kappa_{3}=6\cdot 10^{-3}$. The simulations are performed for various Young’s moduli, $\kappa_{s}=\left\{10^{-4},10^{-3},10^{-2},10^{-1},10^{0}\right\}$, and constant bending and coupling coefficients, $\kappa_{b}=10^{-2}$ and $\kappa_{c}=10^{-2}$. Periodic boundary conditions are applied at all the domain boundaries. At the converged state, the dimensionless aspect ratio of the elastic capsules can be defined as $L_{c}/H_{c}$.

The effect of Young’s modulus on the aspect ratio of the capsules is presented in figure 5 for the different surface tension ratios. For low $\gamma_{12}/\gamma_{13}$, the capsules seem to take similar mechanical equilibrium shapes independently of their degree of elasticity. For moderate and high $\gamma_{12}/\gamma_{13}$, the aspect ratio changes significantly between stiff and moderately deformable $\left(\kappa_{s}=10^{-2}\right)$ capsules, while it reaches a plateau for highly deformable $\left(\kappa_{s}\leq 10^{-3}\right)$ capsules. It can also be noted that the surface tension ratio affects considerably the capsules’ aspect ratio for a given Young’s modulus. This can be clearly seen in figure 6$\left(a\right)$, where the mechanical equilibrium shapes of the highly deformable $\left(\kappa_{s}=10^{-3}\right)$ capsules are plotted for the different $\gamma_{12}/\gamma_{13}$. At $\gamma_{12}/\gamma_{13}\approx 0.57$, the capsules retain an almost circular shape, having only the part of their surface coming in contact with the capillary bridge slightly compressed. As the surface tension ratio increases, the capsules move towards each other, taking a semi-circular shape and causing the formation of a narrower and higher capillary bridge. Despite the fact that Young’s modulus has a notable effect on the aspect ratio of the elastic capsules at high $\gamma_{12}/\gamma_{13}$, the corresponding variations in their shapes are small for different $\kappa_{s}$, as shown in figure 6$\left(b\right)$. The capsule becomes slightly shorter and wider with an increase in $\kappa_{s}$.

Finally, the transient shapes of the highly deformable capsules at a surface tension ratio $\gamma_{12}/\gamma_{13}\approx 0.80$ are shown in figure 7. For clarity, we present only the results of the left-hand side capsule; the transient shapes of the right-hand side capsule are symmetric to the ones presented across the vertical centreline. As mentioned previously, the elastic capsules have initially, at $t_{0}=0$, a circular shape. It is worth noting that the capsules take quickly, already at $t=t_{1}$, a semi-circular shape similar to the mechanical equilibrium one. As the time passes, the capsules become narrower and more elongated along the $x-$ and $y-$axis, respectively. Minimal changes in the capsules shapes can be observed between $t_{5}$, which corresponds to half of the total simulation time, and $t_{f}$, where the capsules are considered to have reached the mechanical equilibrium shape.

V Conclusions

In this work, we have presented a novel simulation technique for the coupling of a multicomponent fluid flow with deformable, infinitely thin structures. For the sake of simplicity, we have considered the case where elastic membranes enclosing a fluid component $3$ are immersed in a two-component/phase flow consisted of fluids $1$ and $2$. For this case, we have introduced a modified formulation of the free energy of the ternary fluid, taking into consideration its interaction with the elastic membranes, to the original one proposed by Semprebon et al. (2016a). Having as a starting point an equation for the conservation of the entropy, we have made use of the splitting structure of the latter and thermodynamic laws in order to derive the equations governing the evolution of the ternary fluid. Following a similar reasoning, the equation describing the response of the structure to the surrounding fluids has been deduced.

The macroscopic equations of motion of the fluid-structure system are solved here by a monolithic numerical approach. This approach consists of: a mesoscopic lattice Boltzmann method for resolving the ternary fluid flow in an Eulerian description, a finite difference method to evolve the membranes equations of motion in a Lagrangian framework, and an immersed boundary method to couple the Eulerian and Lagrangian solvers. The fluid and structure solvers are coupled through a forcing source term in the lattice Boltzmann equation, acting as a feedback of the structure’s response on the flow.

We have subsequently validated our computational algorithm against Surface Evolver, an open-source software capable of modelling steady liquid surfaces problems employing an energy minimization approach. The configuration of an elastic capsule placed at a fluid-fluid interface was considered as the benchmark test. We have compared in detail the equilibrium shapes of the capsule, and its corresponding deformation parameters for different scenarios of the surface tensions and combinations of the capsule’s Young’s and bending moduli. An overall good agreement has been observed between our results and the reference ones. We have also demonstrated the Galilean invariance of our model equations. Finally, our algorithm has been applied to a more complex configuration, that is the capillary bridge formed between two elastic capsules. Although an extensive investigation of the parameter space was beyond the scope of the present work, it should be noted that this configuration is a particularly rich phenomenon, where the criteria for the bridge rupture and the case of unequal capsules are worth being studied in the future.

We have assumed here that all fluid components have the same density. By modifying the Landau free-energy functional $\mathcal{E}_{f}$ of (44) in an appropriate way (Wöhrwag et al., 2018) and following a rationale similar to the one presented in §II, our model could be extended to the case where the components of the ternary fluid mixture have different densities. The model could also be generalised so as to include more fluid components enclosed in/surrounding the elastic membranes, allowing us to tackle a wider range of applications, for example capsules containing multiple aqueous phases (Scott Long et al., 2008; Kusumaatmaja et al., 2009). Also, extensions to three-dimensional configurations and other lattice Boltzmann models are straightforward. In addition, the formulations of the strain and bending energies could be readily modified to consider materials obeying different constitutive laws, such as non-linear hyperelastic materials or biological membranes. Finally, the structure solver could be adapted in order to simulate open surfaces, encountered for instance in the wetting of a soft substrate (Marchand et al., 2012; Style & Dufresne, 2012).

Acknowledgements

We would like to thank Professor Ken Brakke for his assistance in implementing our free energy model of the elastic membranes in Surface Evolver. MP and HK acknowledge EPSRC (EP/P007139/1) for funding. ACMS is grateful to the EPSRC Centre for Doctoral Training in Soft Matter and Functional Interfaces (EP/L015536/1) for financial support. TK acknowledges the award of a Chancellor’s Fellowship from the University of Edinburgh. C. Semprebon acknowledges support from Northumbria University through the Vice-Chancellor’s Fellowship Programme.

Appendix A

The chemical potentials $\mu_{\rho}$, $\mu_{\phi}$ and $\mu_{\psi}$ are expressed as

[TABLE]

The pressure tensor $\mathbf{p}$ is given by

[TABLE]

where the mixing coefficients take the form

[TABLE]

The bulk pressure term $p_{b}$ is given by

[TABLE]

We have ignored terms that provide constant pressure contribution throughout the whole system.