2.1 The LB colour-gradient model for immiscible two-phase flows

The model we are using for two-phase flow simulations is a diffuse-interface LB colour-gradient model, where we introduce the interfacial tension force and Marangoni stress as well as the phase segregation in accordance with Liu et al. (Reference Liu, Zhang and Valocchi2012b ). In this model, distribution functions $f_{i}^{R}$ and $f_{i}^{B}$ are used to represent the red and blue fluids, and the total distribution function is defined as $f_{i}=f_{i}^{R}+f_{i}^{B}$ , where the subscript $i$ is the lattice velocity direction and ranges from 0 to $(n-1)$ for a given D $m$ Q $n$ lattice model. In this study, the D2Q9 and D3Q19 lattice models are used for 2D and 3D simulations respectively. Their corresponding lattice velocities $\boldsymbol{e}_{i}$ and weighting factors $w_{i}$ can be found in Mei et al. (Reference Mei, Yu, Shyy and Luo2002).

The present colour-gradient model consists of three steps, i.e. collision step, recolouring step and streaming step. First, the collision step reads as

(2.1) $$\begin{eqnarray}f_{i}^{\dagger }\left(\boldsymbol{x},t\right)=f_{i}\left(\boldsymbol{x},t\right)-\unicode[STIX]{x1D714}\left[f_{i}(\boldsymbol{x},t)-f_{i}^{eq}(\boldsymbol{x},t)\right]+\unicode[STIX]{x1D6F7}_{i},\end{eqnarray}$$

where $f_{i}\left(\boldsymbol{x},t\right)$ is the total distribution function in the $i$ th velocity direction at the position $\boldsymbol{x}$ and time $t$ , $f_{i}^{\dagger }$ and $f_{i}^{eq}$ are the post-collision and equilibrium distribution functions respectively, $\unicode[STIX]{x1D714}$ is the dimensionless relaxation parameter related to the fluid viscosity and $\unicode[STIX]{x1D6F7}_{i}$ is the perturbation term.

The equilibrium distribution function is obtained by a second-order Taylor expansion of the Maxwell–Boltzmann distribution with respect to the local fluid velocity $\boldsymbol{u}$ ,

(2.2) $$\begin{eqnarray}f_{i}^{eq}\left(\unicode[STIX]{x1D70C},\boldsymbol{u}\right)=\unicode[STIX]{x1D70C}w_{i}\left[1+\frac{\boldsymbol{e}_{i}\boldsymbol{\cdot }\boldsymbol{u}}{c_{s}^{2}}+\frac{(\boldsymbol{e}_{i}\boldsymbol{\cdot }\boldsymbol{u})^{2}}{2c_{s}^{4}}-\frac{\boldsymbol{u}^{2}}{2c_{s}^{2}}\right],\end{eqnarray}$$

where $\unicode[STIX]{x1D70C}=\unicode[STIX]{x1D70C}^{R}+\unicode[STIX]{x1D70C}^{B}$ is the total density, with the superscripts ‘ $R$ ’ and ‘ $B$ ’ referring to the red and blue fluids respectively, and $c_{s}=1/\sqrt{3}$ is the speed of sound.

The perturbation term contributes to the mixed interfacial regions and generates an interfacial force $\boldsymbol{F}_{s}$ . The perturbation term is given by (Guo, Zheng & Shi Reference Guo, Zheng and Shi2002a ; Liu et al. Reference Liu, Zhang and Valocchi2012b )

(2.3) $$\begin{eqnarray}\unicode[STIX]{x1D6F7}_{i}(\boldsymbol{x},t)=\left(1-\frac{\unicode[STIX]{x1D714}}{2}\right)w_{i}\left[\frac{\boldsymbol{e}_{i}-\boldsymbol{u}}{c_{s}^{2}}+\frac{(\boldsymbol{e}_{i}\boldsymbol{\cdot }\boldsymbol{u})\boldsymbol{e}_{i}}{c_{s}^{4}}\right]\boldsymbol{\cdot }\boldsymbol{F}_{s}(\boldsymbol{x},t)\unicode[STIX]{x1D6FF}_{t},\end{eqnarray}$$

where $\unicode[STIX]{x1D6FF}_{t}$ is the time step and the local fluid velocity is defined by $\unicode[STIX]{x1D70C}\boldsymbol{u}(\boldsymbol{x},t)=\sum _{i}f_{i}(\boldsymbol{x},t)\boldsymbol{e}_{i}+\boldsymbol{F}_{s}(\boldsymbol{x},t)\unicode[STIX]{x1D6FF}_{t}/2$ . Following Liu et al. (Reference Liu, Zhang and Valocchi2012b ), the concept of continuum surface force (CSF) is used to model the interfacial force with dynamic interfacial tension and Marangoni stress, which reads as

(2.4) $$\begin{eqnarray}\boldsymbol{F}_{s}(\boldsymbol{x},t)=-{\textstyle \frac{1}{2}}\unicode[STIX]{x1D70E}\unicode[STIX]{x1D705}\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}^{N}+{\textstyle \frac{1}{2}}|\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}^{N}|\unicode[STIX]{x1D735}_{s}\unicode[STIX]{x1D70E},\end{eqnarray}$$

where $\unicode[STIX]{x1D70C}^{N}=(\unicode[STIX]{x1D70C}^{R}-\unicode[STIX]{x1D70C}^{B})/\unicode[STIX]{x1D70C}$ is a colour indicator and is introduced to identify the location of the interface, $\unicode[STIX]{x1D70E}$ is the interfacial tension, $\unicode[STIX]{x1D735}_{s}=(\unicode[STIX]{x1D644}-\boldsymbol{n}\otimes \boldsymbol{n})\boldsymbol{\cdot }\unicode[STIX]{x1D735}$ is the surface gradient operator, $\unicode[STIX]{x1D644}$ is the second-order identity tensor, $\boldsymbol{n}$ is the interfacial unit normal vector, defined by $\boldsymbol{n}=-\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}^{N}/|\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}^{N}|$ and $\unicode[STIX]{x1D705}$ is the local interface curvature, related to $\boldsymbol{n}$ by

(2.5) $$\begin{eqnarray}\unicode[STIX]{x1D705}=-\unicode[STIX]{x1D735}_{s}\boldsymbol{\cdot }\boldsymbol{n}=-\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{n}.\end{eqnarray}$$

It should be noted that the first term on the right-hand side of (2.4) is the interfacial tension force (the normal force) and the second term is the Marangoni stress (the tangential force) that results from the non-uniform distribution of surfactants along the interface. It is worth emphasizing that neglect of the Marangoni stress term could lead to significant errors in simulation of interfacial flows with variable interfacial tension (Sui Reference Sui2014).

In the presence of surfactants, an equation of state (EOS) is required to relate the interfacial tension to the surfactant concentration. Often, a Langmuir EOS in the form

(2.6) $$\begin{eqnarray}\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D713})=\unicode[STIX]{x1D70E}_{0}[1+E_{0}\ln (1-\unicode[STIX]{x1D713}/\unicode[STIX]{x1D713}_{\infty })]\end{eqnarray}$$

is used, where $\unicode[STIX]{x1D713}$ is the surfactant concentration, $\unicode[STIX]{x1D713}_{\infty }$ is the surfactant concentration at the maximum packing, $\unicode[STIX]{x1D70E}_{0}$ is the interfacial tension of a clean interface (i.e. $\unicode[STIX]{x1D713}=0$ ), $E_{0}=RT\unicode[STIX]{x1D713}_{\infty }/\unicode[STIX]{x1D70E}_{0}$ is the elasticity number that measures the sensitivity of $\unicode[STIX]{x1D70E}$ to the variation of $\unicode[STIX]{x1D713}$ , $R$ is the constant of ideal gas and $T$ is the absolute temperature. For low surfactant concentration, the following linear approximation of the Langmuir equation can be used:

(2.7) $$\begin{eqnarray}\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D713})=\unicode[STIX]{x1D70E}_{0}\left(1-E_{0}\unicode[STIX]{x1D713}/\unicode[STIX]{x1D713}_{\infty }\right).\end{eqnarray}$$

As in Xu et al. (Reference Xu, Li, Lowengrub and Zhao2006), we quantify the initial average surfactant concentration $\unicode[STIX]{x1D713}_{0}$ through the surfactant coverage, which is defined as $x_{in}=\unicode[STIX]{x1D713}_{0}/\unicode[STIX]{x1D713}_{\infty }$ .

Using the Chapman–Enskog expansion, it can be proved that (2.1)–(2.3) exactly recover the Navier–Stokes equations for two-phase flow (Liu et al. Reference Liu, Zhang and Valocchi2012b ), where the pressure and the dynamic viscosity of the fluid mixture are defined by $p=\unicode[STIX]{x1D70C}c_{s}^{2}$  and $\unicode[STIX]{x1D707}=\left(1/\unicode[STIX]{x1D714}-1/2\right)\unicode[STIX]{x1D70C}c_{s}^{2}\unicode[STIX]{x1D6FF}_{t}$ . To allow for unequal viscosities of the two fluids, we determine $\unicode[STIX]{x1D707}$ by a harmonic mean (Liu et al. Reference Liu, Valocci, Werth, Kang and Oostrom2014),

(2.8) $$\begin{eqnarray}\frac{1}{\unicode[STIX]{x1D707}\left(\unicode[STIX]{x1D70C}^{N}\right)}=\frac{1+\unicode[STIX]{x1D70C}^{N}}{2\unicode[STIX]{x1D707}^{R}}+\frac{1-\unicode[STIX]{x1D70C}^{N}}{2\unicode[STIX]{x1D707}^{B}},\end{eqnarray}$$

where $\unicode[STIX]{x1D707}^{k}$ ( $k=R$ or $B$ ) is the dynamic viscosity of fluid $k$ .

Then, the recolouring (segregation) algorithm of Latva-Kokko & Rothman (Reference Latva-Kokko and Rothman2005) is used to promote phase segregation and ensure the immiscibility of both fluids. Following Latva-Kokko & Rothman (Reference Latva-Kokko and Rothman2005), the recoloured distribution functions of the red and blue fluids are

(2.9) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle f_{i}^{R\ddagger }(\boldsymbol{x},t)=\frac{\unicode[STIX]{x1D70C}^{R}}{\unicode[STIX]{x1D70C}}f_{i}^{\dagger }(\boldsymbol{x},t)+\unicode[STIX]{x1D6FD}w_{i}\frac{\unicode[STIX]{x1D70C}^{R}\unicode[STIX]{x1D70C}^{B}}{\unicode[STIX]{x1D70C}}\boldsymbol{e}_{i}\boldsymbol{\cdot }\boldsymbol{n},\\ \displaystyle f_{i}^{B\ddagger }(\boldsymbol{x},t)=\frac{\unicode[STIX]{x1D70C}^{B}}{\unicode[STIX]{x1D70C}}f_{i}^{\dagger }(\boldsymbol{x},t)-\unicode[STIX]{x1D6FD}w_{i}\frac{\unicode[STIX]{x1D70C}^{R}\unicode[STIX]{x1D70C}^{B}}{\unicode[STIX]{x1D70C}}\boldsymbol{e}_{i}\boldsymbol{\cdot }\boldsymbol{n},\end{array}\right\} & & \displaystyle\end{eqnarray}$$

where $f_{i}^{R\ddagger }$ and $f_{i}^{B\ddagger }$ are the recoloured distribution functions of the red and blue fluids respectively. Here, $\unicode[STIX]{x1D6FD}$ is the segregation parameter and is set to be 0.7 to maintain a narrow interface thickness (around 4–5 lattices) and keep spurious velocities low (Halliday et al. Reference Halliday, Hollis and Care2007; Gupta & Kumar Reference Gupta and Kumar2010). Moreover, Liu et al. (Reference Liu, Valocchi and Kang2012a ) indicated that this choice is necessary to reproduce the correct dynamical behaviour of droplets.

After the recolouring step, the red and blue distribution functions propagate to the neighbouring lattice nodes, known as the propagation or streaming step,

(2.10) $$\begin{eqnarray}f_{i}^{k}\left(\boldsymbol{x}+\boldsymbol{e}_{i}\unicode[STIX]{x1D6FF}_{t},t+\unicode[STIX]{x1D6FF}_{t}\right)=f_{i}^{k\ddagger }\left(\boldsymbol{x},t\right),\quad k=R~\text{or}~B,\end{eqnarray}$$

with the post-propagation distribution functions used to compute the densities of both fluids by $\unicode[STIX]{x1D70C}^{k}=\sum _{i}f_{i}^{k}$ .

2.2 Finite difference method for surfactant transport

With a sharp-interface representation, the transport of surfactant concentration is governed by a convection–diffusion equation with a source-like term to account for interfacial stretching and distortion (Stone & Leal Reference Stone and Leal1990),

(2.11) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D713}+\unicode[STIX]{x1D735}_{s}\boldsymbol{\cdot }(\unicode[STIX]{x1D713}\boldsymbol{u}_{s})=D_{s}\unicode[STIX]{x1D6FB}_{s}^{2}\unicode[STIX]{x1D713}-\unicode[STIX]{x1D705}\unicode[STIX]{x1D713}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n},\end{eqnarray}$$

where $\boldsymbol{u}_{s}=(\boldsymbol{I}-\boldsymbol{n}\otimes \boldsymbol{n})\boldsymbol{\cdot }\boldsymbol{u}$ is the tangential velocity and $D_{s}$ is the surfactant diffusivity. This equation has been commonly used to describe the surfactant concentration evolution in interface-tracking methods (Zhang et al. Reference Zhang, Eckmann and Ayyaswamy2006; Feigl et al. Reference Feigl, Megias-Alguacil, Fischer and Windhab2007; Jin & Stebe Reference Jin and Stebe2007; Lai et al. Reference Lai, Tseng and Huang2008; Muradoglu & Tryggvason Reference Muradoglu and Tryggvason2008) and is solved along a sharp interface represented by Lagrangian markers or particles. Assuming that the surfactant concentration $\unicode[STIX]{x1D713}$ could be extended off the interface $\unicode[STIX]{x1D6E4}$ , an equivalent Eulerian formulation was derived from (2.11) by Xu & Zhao (Reference Xu and Zhao2003),

(2.12) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D713}+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}-\unicode[STIX]{x1D713}(\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n})=D_{s}(\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D713}-\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\boldsymbol{\cdot }\boldsymbol{n}-\unicode[STIX]{x1D705}\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}).\end{eqnarray}$$

It should be noted that an additional Hamilton–Jacobian equation, similar to the one used for re-initialization in the level-set method, needs to be solved in order to extend the surfactant concentration on the interface to its neighbourhood (Xu & Zhao Reference Xu and Zhao2003; Xu et al. Reference Xu, Li, Lowengrub and Zhao2006, Reference Xu, Yang and Lowengrub2012). This extension enables one to solve the surfactant transport equation, i.e. (2.12), in an Eulerian framework by using FD or finite element methods.

By extending (2.12) to the general fluid domain $\unicode[STIX]{x1D6FA}$ with the aid of an auxiliary Dirac function, Teigen et al. (Reference Teigen, Song, Lowengrub and Voigt2011) recently proposed a surfactant evolution equation of diffuse-interface form,

(2.13) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}(\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D6E4}}\unicode[STIX]{x1D713})+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D6E4}}\unicode[STIX]{x1D713}\boldsymbol{u})=D_{s}\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D6E4}}\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}),\end{eqnarray}$$

where the Dirac function $\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D6E4}}$ is used to localize the surfactant concentration explicitly at the interface and should satisfy $\int _{\unicode[STIX]{x1D6E4}}\unicode[STIX]{x1D713}\,\text{d}\unicode[STIX]{x1D6E4}=\int _{\unicode[STIX]{x1D6FA}}\unicode[STIX]{x1D713}\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D6E4}}\,\text{d}\unicode[STIX]{x1D6FA}$ (Teigen et al. Reference Teigen, Song, Lowengrub and Voigt2011).

The surfactant evolution equation of diffuse-interface form has been combined with the phase-field method to simulate interfacial flows with surfactants. In spite of its great success and ongoing advances, the phase-field method is still criticized as it cannot conserve mass for each fluid, the choice of an optimal mobility remains an open question and small droplets or bubbles are prone to dissolving unphysically. By contrast, the LB colour-gradient model possesses many advantages in simulating multiphase and multicomponent flows, including strict mass conservation for each fluid, high numerical accuracy and good numerical stability for a broad range of viscosity ratios. In addition, the colour-gradient model is perfectly compatible with the surfactant evolution equation of diffuse-interface form because its interface is of diffuse nature and it adopts a Dirac function to model the interfacial force (Liu et al. Reference Liu, Zhang and Valocchi2012b ). To be consistent with the colour-gradient model, the Dirac function in (2.13) is taken as $\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D6E4}}=|\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}^{N}|/2$ . Thus, the surfactant evolution equation of diffuse-interface form can be further written as

(2.14) $$\begin{eqnarray}\unicode[STIX]{x2202}_{t}(|\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}^{N}|\unicode[STIX]{x1D713})+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D713}|\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}^{N}|\boldsymbol{u})=D_{s}\left[(\unicode[STIX]{x1D735}|\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}_{N}|)\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}+|\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}^{N}|\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D713}\right].\end{eqnarray}$$

It should be noted that (2.14), in combination with the colour-gradient model, is well suited to handling the interface breakup and coalescence without invoking any ad hoc criteria, and it does not need to solve an additional equation to extend the surfactant concentration off the interface, greatly simplifying the computation. More importantly, (2.14) can be easily extended to modelling of soluble surfactants (Teigen et al. Reference Teigen, Song, Lowengrub and Voigt2011).

In the present work, (2.14) is adopted to describe the surfactant transport, and it is solved by the FD method using a modified Crank–Nicholson scheme for the time discretization (Xu & Zhao Reference Xu and Zhao2003),

(2.15) $$\begin{eqnarray}\displaystyle & \displaystyle \left.\begin{array}{@{}c@{}}\displaystyle \frac{\unicode[STIX]{x1D713}^{t+\unicode[STIX]{x1D6FF}_{t}}|\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}^{N}|^{t+\unicode[STIX]{x1D6FF}_{t}}-\unicode[STIX]{x1D713}^{t}|\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}^{N}|^{t}}{\unicode[STIX]{x1D6FF}_{t}}=D_{s}|\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}^{N}|\frac{\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D713}^{t+\unicode[STIX]{x1D6FF}_{t}}+\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D713}^{t}}{2}+{\mathcal{S}}^{t},~~\text{for}~t=0,\\ \displaystyle \frac{\unicode[STIX]{x1D713}^{t+\unicode[STIX]{x1D6FF}_{t}}|\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}^{N}|^{t+\unicode[STIX]{x1D6FF}_{t}}-\unicode[STIX]{x1D713}^{t}|\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}^{N}|^{t}}{\unicode[STIX]{x0394}t}=D_{s}|\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}^{N}|\frac{\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D713}^{t+\unicode[STIX]{x1D6FF}_{t}}+\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D713}^{t}}{2}+\frac{3}{2}{\mathcal{S}}^{t}-\frac{1}{2}{\mathcal{S}}^{t-1},~~\text{for}~t>1,\end{array}\right\} & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

with ${\mathcal{S}}=-\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}(\unicode[STIX]{x1D713}|\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}^{N}|)+D_{s}(\unicode[STIX]{x1D735}|\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}^{N}|)\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}$ . In (2.15), all of the spatial derivatives are discretized using standard central difference schemes except for the convection term $\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}(\unicode[STIX]{x1D713}|\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}^{N}|)$ , which is discretized using a third-order weighted essentially non-oscillatory (WENO) scheme (Jiang & Shu Reference Jiang and Shu1996; Xu & Zhao Reference Xu and Zhao2003). The WENO scheme has the advantages that it can handle steep gradients well and fewer grid points are needed to achieve a high-order resolution.

The resulting linear system for $\unicode[STIX]{x1D713}^{t+\unicode[STIX]{x1D6FF}_{t}}$ from (2.15) is symmetric positive definite, similar to that from a transient heat conduction equation. It can be solved by various efficient and robust algorithms, e.g. the successive over relaxation (SOR) method or the conjugate gradient method (Saad Reference Saad2003). In our numerical computations, we simply use the SOR method with the relaxation factor fixed at 1.2.

3.1 Droplet deformation in an extensional flow

We first simulate the deformation of a surfactant-covered droplet in a 3D extensional flow, and compare the numerical results with those obtained from the BIM (Stone & Leal Reference Stone and Leal1990; Feigl et al. Reference Feigl, Megias-Alguacil, Fischer and Windhab2007). As sketched in figure 1, a spherical droplet (red fluid) of radius $R$ is initially suspended in an immiscible liquid (blue fluid) in an external extensional flow, where the far-field flow velocity is given by

(3.1) $$\begin{eqnarray}\boldsymbol{u}_{\infty }(\boldsymbol{x})=\dot{\unicode[STIX]{x1D6FE}}\left(\begin{array}{@{}ccc@{}}-0.5 & 0 & 0\\ 0 & -0.5 & 0\\ 0 & 0 & 1\end{array}\right)\boldsymbol{\cdot }\boldsymbol{x}.\end{eqnarray}$$

Here, $\dot{\unicode[STIX]{x1D6FE}}$ is the shear rate and $\boldsymbol{x}$ is the position vector. In the presence of surfactants, the droplet behaviour can be characterized by several important dimensionless numbers, including the Reynolds number ( $Re$ ), the capillary number ( $Ca$ ) and the surface Péclet number ( $Pe$ ), which are defined as follows:

(3.2a-c ) $$\begin{eqnarray}Re=\unicode[STIX]{x1D70C}^{B}\dot{\unicode[STIX]{x1D6FE}}R^{2}/\unicode[STIX]{x1D707}^{B},\quad Ca=\unicode[STIX]{x1D707}_{B}R\dot{\unicode[STIX]{x1D6FE}}/\unicode[STIX]{x1D70E}_{0},\quad Pe=R^{2}\dot{\unicode[STIX]{x1D6FE}}/D_{s}.\end{eqnarray}$$

In addition, for a diffuse-interface model, the Cahn number should be small enough to approach the sharp-interface limit, and it is defined by $Cn=\unicode[STIX]{x1D709}/R$ , where $\unicode[STIX]{x1D709}=1/(6k\unicode[STIX]{x1D6FD})$ is a measure of the diffuse-interface thickness, with $k=0.1504$ for the D2Q9 model and $k=0.134$ for the D3Q19 model (Riaud et al. Reference Riaud, Zhao, Wang, Cheng and Luo2014).

For creeping flows, i.e. $Re\ll 1$ and $Ca\ll 1$ , the droplet will eventually deform into an ellipsoidal shape, which is usually characterized by the deformation parameter $D_{f}$ and the dimensionless extensional length $L$ ,

(3.3a,b ) $$\begin{eqnarray}D_{f}=\frac{a-b}{a+b},\quad L=a/R,\end{eqnarray}$$

where $a$ and $b$ are the half-lengths of the major and minor axes of the droplet respectively.

The simulations are first run at $Re=0.05$ for a spherical droplet of $R=30$ lattices (the corresponding Cahn number is $Cn=1/(6k\unicode[STIX]{x1D6FD}R)=0.0592$ ), with the effective capillary number ( $Ca_{e}$ ) varying from 0.028 to 0.11. The effective capillary number is defined as

(3.4) $$\begin{eqnarray}Ca_{e}=\frac{Ca}{1-x_{in}E_{0}},\end{eqnarray}$$

and it is used here to keep consistency with Stone & Leal (Reference Stone and Leal1990). The linear EOS (2.6) is used with the elasticity number $E_{0}=0.1$ . Both fluids are assumed to have equal densities and viscosities, and the surface Péclet number is set as $Pe=10Ca$ . Due to symmetry, we only simulate one octant of the droplet, with a domain size of $100\times 100\times 160$ lattices. To obtain an extensional flow, we apply the non-equilibrium extrapolation boundary condition proposed by Guo, Zheng & Shi (Reference Guo, Zheng and Shi2002b ) at the far field. For all values of $Ca_{e}$ considered, the droplet is continuously elongated along the $z$ -direction until a steady shape is achieved.

Figure 2. The steady droplet shapes (represented in red) and streamlines at the $x=0$ plane for effective capillary numbers of (a) $Ca_{e}=0.028$ , (b) $Ca_{e}=0.057$ , (c) $Ca_{e}=0.083$ and (d) $Ca_{e}=0.11$ .

Figure 2 shows the steady droplet shapes and streamlines at the $x=0$ plane for various values of $Ca_{e}$ . It is seen that a single vortex is formed and fully filled inside the droplet in steady state. Moreover, the droplet deformation increases with increasing $Ca_{e}$ . Figure 3 plots the dimensionless surfactant concentration $\unicode[STIX]{x1D713}^{\ast }=\unicode[STIX]{x1D713}/\unicode[STIX]{x1D713}_{0}$ on the steady droplet interface (represented by $\unicode[STIX]{x1D70C}^{N}=0$ ) as a function of $z/R$ for different values of $Ca_{e}$ . It is observed that for each $Ca_{e}$ , the surfactant concentration reaches its highest value at the tip of the droplet and gradually decreases along the negative direction of the $z$ -axis.

Figure 3. The dimensionless surfactant concentration $\unicode[STIX]{x1D713}^{\ast }$ on the steady droplet interface ( $\unicode[STIX]{x1D70C}^{N}=0$ ) as a function of $z/R$ for various values of $Ca_{e}$ .

In addition, we quantify the deformation parameter $D_{f}$ in steady state, dimensionless extensional length $L$ , and the maximum ( $\unicode[STIX]{x1D713}_{max}^{\ast }$ ) and minimum ( $\unicode[STIX]{x1D713}_{min}^{\ast }$ ) values of the dimensionless surfactant concentration at $\unicode[STIX]{x1D70C}^{N}=0$ for various values of $Ca_{e}$ , which are plotted in figure 4. It is seen that our numerical results for $D_{f}$ , $L$ , $\unicode[STIX]{x1D713}_{max}^{\ast }$ and $\unicode[STIX]{x1D713}_{min}^{\ast }$ are all in good agreement with those obtained by Stone & Leal (Reference Stone and Leal1990) and Feigl et al. (Reference Feigl, Megias-Alguacil, Fischer and Windhab2007). This suggests that the present method can produce equally accurate results to the sharp-interface approach for a finite interface thickness.

Figure 4. Comparison of the simulated (a) deformation parameter $D_{f}$ and dimensionless extensional length $L$ , and (b) maximum ( $\unicode[STIX]{x1D713}_{max}^{\ast }$ ) and minimum ( $\unicode[STIX]{x1D713}_{min}^{\ast }$ ) values of the dimensionless surfactant concentration with those obtained by Stone & Leal (Reference Stone and Leal1990) and Feigl et al. (Reference Feigl, Megias-Alguacil, Fischer and Windhab2007).

3.2 Droplet deformation in a 2D shear flow

We then apply the hybrid method to simulate the droplet deformation in a 2D shear flow in the presence of surfactants, and compare the numerical results with those obtained from the level-set method (Xu et al. Reference Xu, Yang and Lowengrub2012). As illustrated in figure 5(a), a circular droplet is initially placed halfway between two parallel walls that are separated by a distance $H$ . Shear flow is introduced to this system by moving the walls with equal but opposite velocities $u_{w}$ , and the resulting shear rate is $\dot{\unicode[STIX]{x1D6FE}}=2u_{w}/H$ . The computational domain is set as $[-5R,5R]\times [-2R,2R]$ , where $R$ is the initial droplet radius and is taken to be 100 lattices. Periodic boundary conditions are applied in the horizontal direction, while the halfway bounce-back scheme is used on the upper and lower walls. The Langmuir EOS (2.7) is employed, with $\unicode[STIX]{x1D70E}_{0}=10^{-3}$ and $E_{0}=0.2$ . We run the simulations with $Ca=0.1$ , $Re=Pe=10$ and three different values of the surfactant coverage, i.e. $x_{in}=0$ , $0.3$ and $0.6$ (it should be noted that $x_{in}=0$ corresponds to the case of a clean droplet).

Figure 5. Illustration of (a) a 2D droplet and (b) a 3D droplet in a simple shear flow. In (b), the computational domain is taken as half of the full system, i.e. $[0,L]\times [0.5,0.5+W]\times [0,H]$ , and the droplet centre is located at $(L/2,0.5,H/2)$ .

Figure 6 shows the droplet shapes at $\unicode[STIX]{x1D70F}=\dot{\unicode[STIX]{x1D6FE}}t=9$ for three different values of $x_{in}$ , where the $x$ and $y$ coordinates are both normalized by the droplet radius $R$ . Increase in $x_{in}$ leads to a more prolate droplet and a decrease of the droplet inclination angle, defined as the angle between the major axis of the droplet and the horizontal ( $x$ ) axis. Moreover, the droplet shapes obtained by Xu et al. (Reference Xu, Yang and Lowengrub2012) are also plotted in this figure as a comparison, and good agreement is observed for all values of $x_{in}$ considered.

Figure 6. Droplet shapes at $\unicode[STIX]{x1D70F}=9$ for $x_{in}=0$ , $0.3$ and $0.6$ . Other parameters are taken as $Re=10$ , $Pe=10$ , $Ca=0.1$ and $E_{0}=0.2$ . The $x$ and $y$ coordinates are both normalized by the droplet radius $R$ .

Figure 7(a,b) shows the distributions of the dimensionless surfactant concentration $\unicode[STIX]{x1D713}^{\ast }$ and dimensionless interfacial tension $\unicode[STIX]{x1D70E}^{\ast }=\unicode[STIX]{x1D70E}/\unicode[STIX]{x1D70E}_{0}$ along the arclength $s$ , together with a comparison with those obtained by the level-set method (Xu et al. Reference Xu, Yang and Lowengrub2012). It should be noted that the arclength $s$ , measured counterclockwise from the right intersection point (i.e. the point ‘A’ in figure 5 a) between $y=0$ and the droplet interface, is normalized by the initial droplet radius $R$ . It is clearly seen that our numerical results agree well with those in Xu et al. (Reference Xu, Yang and Lowengrub2012). Specifically, for the surfactant-covered droplets, i.e. $x_{in}=0.3$ and $0.6$ , we find that $\unicode[STIX]{x1D713}^{\ast }$ is distributed non-uniformly along the interface and tends to accumulate at the droplet tips, resulting in non-uniform distribution of the interfacial tension where the lowest interfacial tension occurs at the droplet tips. In addition, increase in $x_{in}$ decreases the local interfacial tension (see figure 7 b) that resists the movement of the droplet, thereby promoting droplet deformation and reducing the inclination angle.

Figure 7. The distributions of (a) the dimensionless surfactant concentration $\unicode[STIX]{x1D713}^{\ast }$ , (b) the dimensionless interfacial tension $\unicode[STIX]{x1D70E}^{\ast }$ , (c) the signed dimensionless capillary force $F_{Ca}$ and (d) the signed dimensionless Marangoni force $F_{Ca}$ along the arclength $s$ of the interface at $\unicode[STIX]{x1D70F}=9$ for $x_{in}=0$ , $0.3$ and $0.6$ . Other parameters are taken as $Re=10$ , $Pe=10$ , $Ca=0.1$ and $E_{0}=0.2$ .

In figure 7(c,d), we plot the dimensionless capillary force $F_{Ca}$ and Marangoni force $F_{Ma}$ as a function of the arclength $s$ , where $F_{Ca}$ and $F_{Ma}$ are calculated by $\unicode[STIX]{x1D705}\unicode[STIX]{x1D70E}^{\ast }R$ and $\unicode[STIX]{x1D735}_{s}\unicode[STIX]{x1D70E}^{\ast }\boldsymbol{\cdot }\boldsymbol{\unicode[STIX]{x1D70F}}R$ with the unit tangential vector $\boldsymbol{\unicode[STIX]{x1D70F}}=(-n_{y},n_{x})$ . It can be observed that the distributions of $F_{Ca}$ and $F_{Ma}$ are in excellent agreement with the numerical results of Xu et al. (Reference Xu, Yang and Lowengrub2012). Along the arclength, $F_{Ca}$ reaches its maximum value at the droplet tips due to the high interface curvature, and increase of $x_{in}$ decreases the absolute value of the capillary force. The Marangoni force arises due to the non-uniform distribution of $\unicode[STIX]{x1D70E}^{\ast }$ . It acts in the tangential direction of the interface, and is directed from higher to lower $\unicode[STIX]{x1D713}^{\ast }$ regions, resisting the accumulation of surfactants (or, say, lowering the non-uniformity of the surfactant distribution). In figure 7(d), we find that $F_{Ma}$ is zero at the droplet tips due to symmetry, and takes its local maximum and minimum values near the tips because of the highest gradient of $\unicode[STIX]{x1D70E}^{\ast }$ in these regions. It is also noticed that the amplitude of $F_{Ma}$ increases with increasing $x_{in}$ from 0.3 to 0.6.

4.1 Droplet deformation

For a clean droplet in the Stokes flow regime, Taylor (Reference Taylor1934) derived a theoretical formula for the small deformation in the shear flow in terms of the viscosity ratio $\unicode[STIX]{x1D706}$ ( $\unicode[STIX]{x1D706}=\unicode[STIX]{x1D707}^{R}/\unicode[STIX]{x1D707}^{B}$ ) and the capillary number,

(4.1) $$\begin{eqnarray}D_{f,T}=\frac{19\unicode[STIX]{x1D706}+16}{16\unicode[STIX]{x1D706}+16}Ca.\end{eqnarray}$$

However, the prediction of Taylor does not consider the proximity of the droplet to the walls, i.e. the confinement ratio, defined as the ratio between the droplet diameter $2R$ and the wall separation $H$ ( $2R/H$ ). When the confinement ratio is higher than 0.4, the deformation cannot be predicted accurately by (4.1). Moreover, (4.1) is not able to describe the droplet deformation for very large viscosity ratios (Vananroye, Van Puyvelde & Moldenaers Reference Vananroye, Van Puyvelde and Moldenaers2007). More discussion about this equation can be found in several review papers (Cristini & Tan Reference Cristini and Tan2004; Guido & Greco Reference Guido and Greco2004; Puyvelde et al. Reference Puyvelde, Vananroye, Cardinaels and Moldenaers2008; Minale Reference Minale2010).

To address the influence of the confinement ratio, a number of theoretical or phenomenological models have been proposed to predict the deformation parameter. For instance, Shapira & Haber (Reference Shapira and Haber1990) provided a new expression that combines the Taylor deformation and an additional term accounting for the wall effects, and Maffettone & Minale (Reference Maffettone and Minale1998) proposed a phenomenological model to predict the droplet deformation under transient conditions. Later, Vananroye et al. (Reference Vananroye, Van Puyvelde and Moldenaers2007) developed the so-called MMSH model, in which the model of Shapira & Haber (Reference Shapira and Haber1990) is modified by replacing the Taylor model with the model of Maffettone & Minale (Reference Maffettone and Minale1998) as bulk reference. Vananroye et al. (Reference Vananroye, Van Puyvelde and Moldenaers2007) found that the droplet behaviour can be reasonably predicted by the MMSH model for a wide range of viscosity ratios and confinement ratios, and for $Ca\leqslant 0.3$ . The MMSH model is given by

(4.2) $$\begin{eqnarray}D_{f,MMSH}=D_{f,MM}\left[1+C_{s}\frac{1+2.5\unicode[STIX]{x1D706}}{1+\unicode[STIX]{x1D706}}\left(\frac{R}{H}\right)^{3}\right],\end{eqnarray}$$

with

(4.3) $$\begin{eqnarray}D_{f,MM}=\frac{\sqrt{m_{1}^{2}+Ca^{2}}-\sqrt{m_{1}^{2}+(1-m_{2}^{2})Ca^{2}}}{m_{2}Ca},\end{eqnarray}$$

where

(4.4a,b ) $$\begin{eqnarray}m_{1}=\frac{40(\unicode[STIX]{x1D706}+1)}{(2\unicode[STIX]{x1D706}+3)(19\unicode[STIX]{x1D706}+16)},\quad m_{2}=\frac{5}{2\unicode[STIX]{x1D706}+3}+\frac{3Ca^{2}}{2+6Ca^{2}}\end{eqnarray}$$

and $C_{s}$ is a shape parameter, which is taken as $5.6996$ for a droplet positioned halfway between two plates.

Figure 8. The time evolution of the deformation parameter $D_{f}$ for $x_{in}=0$ , 0.15 and 0.3 at $Ca=0.1$ , 0.2 and 0.3.

In our numerical simulation, the initial droplet radius is $R=25$ , and the domain height and width are $H=100$ and $W=50$ respectively, so that the confinement ratio is $2R/H=0.5$ . The simulations are run in a $200\times 50\times 100$ lattice domain for the capillary number ranging from 0.05 to 0.3 at $Re=0.1$ , which is small enough to meet the Stokes flow condition. Three different values of the surfactant coverage are considered, i.e. $x_{in}=0$ , $0.15$ and $0.3$ . Figure 8 shows the time evolution of the deformation parameter $D_{f}$ at $Pe=10$ for $Ca=0.1,0.2$ and $0.3$ . It is observed that $D_{f}$ first increases and then evolves to a steady value for each of the cases considered. For a fixed surfactant coverage $x_{in}$ , $D_{f}$ increases with $Ca$ , and for a constant $Ca$ , $D_{f}$ increases with $x_{in}$ . In addition, the droplet takes a longer time to achieve a steady shape for higher $Ca$ or $x_{in}$ .

In figure 9, we present the variation of the deformation parameter with $Ca$ for the cases of clean droplet, surfactant-laden droplet, and clean droplet with reduced interfacial tension $\unicode[STIX]{x1D70E}=\unicode[STIX]{x1D70E}_{0}/[1+E_{0}\ln (1-x_{in})]$ , and their corresponding deformation parameters are denoted as $D_{f0}$ , $D_{f2}$ and $D_{f1}$ for the convenience of subsequent discussion. Predictions from the MMSH model are also plotted as a comparison. It can be seen that our simulation results concerning the clean droplet agree well with the predictions from the MMSH model, suggesting that the present method is able to accurately simulate the deformation of a clean droplet in simple shear flows. The presence of the surfactants increases the droplet deformation for all values of $Ca$ . The increased deformation, i.e. $D_{f2}-D_{f0}$ , consists of two parts: (1) $\unicode[STIX]{x1D6FF}_{1}=D_{f1}-D_{f0}$ , which is attributed to a reduction of the interfacial tension caused by the average surfactant concentration $\unicode[STIX]{x1D713}_{0}$ , and (2) $\unicode[STIX]{x1D6FF}_{2}=D_{f2}-D_{f1}$ , which is attributed to the non-uniform effects from non-uniform capillary pressure on the droplet surface and Marangoni stresses along the interface. As can be clearly seen from the inset of figure 9, $\unicode[STIX]{x1D6FF}_{1}$ accounts for most of the effect of surfactants on the droplet deformation, and increases with $Ca$ . On the other hand, as $Ca$ increases, $\unicode[STIX]{x1D6FF}_{2}$ first increases and then decreases, with its maximum value occurring at $Ca=0.15$ . The reason may be that as $Ca$ increases, the dilution of the surfactants becomes increasingly significant, which weakens the non-uniform effects.

Figure 9. The deformation parameter $D_{f}$ as a function of $Ca$ for the clean droplet, the surfactant-laden droplet, and the clean droplet with the reduced interfacial tension $\unicode[STIX]{x1D70E}_{e}=\unicode[STIX]{x1D70E}_{0}/[1+E_{0}\ln (1-x_{in})]$ at (a) $x_{in}=0.15$ and (b) $x_{in}=0.3$ . The predictions from the MMSH model are also shown for comparison. The inset plots the increment of the deformation parameter $\unicode[STIX]{x1D6FF}$ (due to the presence of surfactants) and its two constituents $\unicode[STIX]{x1D6FF}_{1}$ and $\unicode[STIX]{x1D6FF}_{2}$ as a function of $Ca$ .

The surface Péclet number $Pe$ is an important parameter which affects the non-uniform distribution of surfactants, and here we take $Pe$ ranging from 0.1 to 10 to investigate the effect of $Pe$ on the droplet deformation. Different values of $Pe$ are achieved by varying $D_{s}$ while keeping other parameters unchanged. In figure 10, we plot the deformation parameter $D_{f}$ , and the $\unicode[STIX]{x1D713}_{min}$ and $\unicode[STIX]{x1D713}_{max}$ values of the surfactant concentration at $\unicode[STIX]{x1D70C}^{N}=0$ as a function of $Pe$ for $Ca=0.15$ and $x_{in}=0.3$ . It is observed that the difference between $\unicode[STIX]{x1D713}_{max}$ and $\unicode[STIX]{x1D713}_{min}$ increases with $Pe$ , which indicates an enhanced non-uniformity of the surfactant distribution. This non-uniformity leads to a larger $\unicode[STIX]{x1D6FF}_{2}$ , and thus an enhanced droplet deformation, as shown in figure 10(a).

Figure 10. (a) The deformation parameter $D_{f}$ and (b) the minimum and maximum values of the surfactant concentration on the interface ( $\unicode[STIX]{x1D70C}^{N}=0$ ) for various values of $Pe$ at $Ca=0.15$ and $x_{in}=0.3$ .

As pointed out in § 1, the Marangoni stress term was not explicitly considered by Farhat et al. (Reference Farhat, Celiker, Singh and Lee2011) in the modelling of interfacial flows with insoluble surfactants. To show the effect of neglecting the Marangoni stress on the droplet behaviour, we run the simulation with the Marangoni stress turned off for $Ca=0.1$ , $Re=0.1$ , $Pe=10$ and $x_{in}=0.15$ , and compare the simulated results with those obtained when the Marangoni stress is included. Figure 11 shows the final droplet shapes and surfactant concentration distributions on the interface ( $\unicode[STIX]{x1D70C}^{N}=0$ ) for the surfactant-laden cases with and without Marangoni stresses. It can be seen that neglect of the Marangoni stress leads to more surfactants distributed in the tip regions of the droplet, which cause a decrease in interfacial tension and thus a larger droplet deformation. Specifically, the deformation parameter is $D_{f}=0.148$ in the case with Marangoni stress, which is increased to $D_{f}=0.169$ when the Marangoni stress is turned off. However, it should be noted that, physically, the Marangoni stress will arise spontaneously as long as there exists an interfacial tension gradient (which can be caused by a concentration gradient or a temperature gradient), and neglect of the Marangoni stress in numerical modelling could lead to significant errors or even incorrect droplet behaviour when a spatially varying interfacial tension is considered (Sui Reference Sui2014).

Figure 11. Droplet shapes and surfactant concentration distributions at $Ca=0.1$ , $Re=0.1$ , $Pe=10$ and $x_{in}=0.15$ for the surfactant-laden cases (a) with and (b) without Marangoni stresses.

4.2 Droplet breakup

A large number of studies (Pathak & Migler Reference Pathak and Migler2003; Sibillo et al. Reference Sibillo, Pasquariello, Simeone, Cristini and Guido2006; Renardy Reference Renardy2007; Janssen et al. Reference Janssen, Vananroye, Van Puyvelde, Moldenaers and Anderson2010) have indicated that in the case of a clean droplet, there usually exists a critical capillary number $Ca_{cr}$ , above which the droplet continuously deforms without reaching a steady state and finally breaks up into daughter droplets. It has been experimentally and numerically recognized that $Ca_{cr}$ is a function of the confinement ratio and the viscosity ratio in the Stokes flow regime. Here, we investigate the critical capillary number of droplet breakup at a constant viscosity ratio of unity for varying confinement ratios and Reynolds numbers. The domain width is fixed at $4R$ with the droplet radius set as $R=30$ . The domain height varies with the confinement ratio, and the domain length is set to be at least $10R$ in order to avoid the occurrence that the droplet touches itself via the periodic boundary conditions before breakup.

4.2.1 Influence of the capillary number and confinement ratio

First, the critical capillary number for the breakup of a clean droplet is studied for $Re=0.1$ with confinement ratios ranging from 0.3 to 0.85. As previously done by Janssen et al. (Reference Janssen, Vananroye, Van Puyvelde, Moldenaers and Anderson2010), we define the critical capillary number ( $Ca_{cr}$ ) as the lowest capillary number found at which an initially spherical droplet breaks up. For each confinement ratio, the ‘lowest’ capillary number is found by increasing $Ca$ with an increment of $0.02$ or lower from a case where the droplet does not break up. Figure 12 illustrates the simulation results of our hybrid method along with the experimental and numerical results of Janssen et al. (Reference Janssen, Vananroye, Van Puyvelde, Moldenaers and Anderson2010). In the figure, we use discrete symbols without a ‘+’ to denote the binary breakup (see, e.g., figure 13 b) and discrete symbols with a ‘+’ to denote the ternary breakup, where the droplet splits into three daughter droplets (see, e.g., figure 13 c).

Figure 12. The critical capillary number of droplet breakup as a function of the confinement ratio in surfactant-free and surfactant-laden systems. Binary and ternary breakups are represented by the discrete symbols without a ‘+’ and the discrete symbols with a ‘+’ respectively. The experimental and BIM simulation results for clean droplets are taken from Janssen et al. (Reference Janssen, Vananroye, Van Puyvelde, Moldenaers and Anderson2010).

Figure 13. Snapshots of the droplet shape at $Re=0.1$ for (a) $Ca=0.35$ , $2R/H=0.5$ , (b) $Ca=0.415$ , $2R/H=0.7$ and (c) $Ca=0.46$ , $2R/H=0.8$ for the clean droplets.

As shown in figure 12, the simulated critical capillary number is varied from $0.415$ at $2R/H=0.3$ and $0.4$ , to $0.392$ at $2R/H=0.5$ , then to $0.395$ at $2R/H=0.6$ , next to $0.415$ at $2R/H=0.7$ , and finally to $0.46$ at $2R/H=0.8$ and $0.85$ , as highlighted by the dashed lines. This suggests that the critical capillary number first decreases and then increases with the confinement ratio, and the minimum critical capillary number is achieved at a confinement ratio of around 0.5. These findings agree well with previous experimental and numerical studies (Janssen et al. Reference Janssen, Vananroye, Van Puyvelde, Moldenaers and Anderson2010) in trend. Quantitatively speaking, the values of the critical capillary number obtained with the present method are a little higher than those obtained by Janssen et al. (Reference Janssen, Vananroye, Van Puyvelde, Moldenaers and Anderson2010) using the BIM; on the other hand, the present results are closer to the available experimental data in comparison with the numerical results of Janssen et al. (Reference Janssen, Vananroye, Van Puyvelde, Moldenaers and Anderson2010). In addition, it is clearly seen that a ternary breakup occurs instead of a binary breakup at higher confinement ratios, e.g. $2R/H=0.8$ (also see figure 13), consistent with previous experimental and numerical findings (Janssen et al. Reference Janssen, Vananroye, Van Puyvelde, Moldenaers and Anderson2010).

Then, we investigate the breakup of a surfactant-laden droplet for various capillary numbers and confinement ratios at $x_{in}=0.3$ , $Re=0.1$ and $Pe=10$ . The droplet and domain sizes are kept the same as those in the corresponding case of a clean droplet. Figure 12 shows the effect of the confinement ratio on the critical capillary number for droplet breakup in the surfactant-laden cases (represented by the dash-dotted lines). For low confinement ratios ( $2R/H<0.7$ ), the critical capillary number first decreases and then increases with increasing confinement ratio, and reaches its minimum value at a confinement ratio of around 0.5. Binary breakup always occurs at low confinement ratios, but ternary breakup is preferred at high confinement ratios. These findings are consistent with those in clean cases. Unlike in clean cases, however, the critical capillary number is almost independent of the confinement ratio for $2R/H\geqslant 0.7$ . It is worth noting that, to capture the exact variation of the critical capillary number with the confinement ratio, a large number of simulations are still required, which will be studied by exploiting massively parallel computing on graphics processing units (GPUs) in the near future. In addition, the presence of surfactants not only significantly decreases the critical capillary number for each confinement ratio but also may influence the mode of droplet breakup. For example, at a confinement ratio of 0.7, the breakup mode for a clean droplet is binary breakup, while ternary breakup occurs for a surfactant-laden droplet.

Figure 14. Snapshots of droplet shape at (a) $Ca=0.35$ , $2R/H=0.5$ and (b) $Ca=0.36$ , $2R/H=0.8$ for a surfactant-laden droplet. It should be noted that the droplet surface is coloured by the dimensionless surfactant concentration $\unicode[STIX]{x1D713}^{\ast }$ .

To further understand the role of surfactants in droplet breakup, we plot snapshots of a surfactant-laden droplet for $2R/H=0.5$ and $Ca=0.35$ in figure 14(a), where the droplet surface is coloured by the dimensionless surfactant concentration. In the stretching stage (see $\unicode[STIX]{x1D70F}=1$ and 10), $\unicode[STIX]{x1D713}_{max}^{\ast }$ increases but $\unicode[STIX]{x1D713}_{min}^{\ast }$ decreases, and the maximum surfactant concentration occurs near the droplet tips, leading to a more elongated droplet than the corresponding case of a clean droplet (see figure 13 a). After the initial stretching stage, a neck forms in the middle of the droplet and then gradually thins until the droplet finally breaks up into two equally sized daughter droplets. In the necking stage (see $\unicode[STIX]{x1D70F}=38$ , $54$ and $55$ ), $\unicode[STIX]{x1D713}_{max}^{\ast }$ decreases and $\unicode[STIX]{x1D713}_{min}^{\ast }$ almost remains constant, exhibiting a severe dilution of the surfactants, consistent with previous studies (Stone & Leal Reference Stone and Leal1990; Teigen et al. Reference Teigen, Song, Lowengrub and Voigt2011). Moreover, it is interesting to note that in the stretching stage, the minimum surfactant concentration occurs in the middle of the droplet where the neck later forms, but as soon as the neck forms, the minimum surfactant concentration $\unicode[STIX]{x1D713}_{min}^{\ast }$ moves away from the necking region where the interface curvature is relatively high. As indicated in figure 12, a high confinement ratio favours ternary breakup in both clean and surfactant-laden cases. Figure 14(b) shows snapshots of ternary breakup for a surfactant-laden droplet at $2R/H=0.8$ and $Ca=0.36$ . Clearly, the droplet stretches to a more elongated shape before two necks form a little away from the droplet centre. During the stretching and necking stages, $\unicode[STIX]{x1D713}_{max}^{\ast }$ and $\unicode[STIX]{x1D713}_{min}^{\ast }$ vary with the same trend as those at a confinement ratio of 0.5, but the dilution of surfactants in the necking stage is more severe as the droplet corresponds to a larger surface area. In the stretching stage, the minimum surfactant concentration still occurs in the middle of the droplet, but the necks do not form there; as soon as the necks form, the surfactant concentration in the necking region increases due to the increased interface curvature, which accelerates the necking process.

Figure 15. The critical capillary number of droplet breakup as a function of the confinement ratio for a clean droplet, a surfactant-laden droplet, and a surfactant-laden droplet with a reduced interfacial tension of $\unicode[STIX]{x1D70E}_{e}=\unicode[STIX]{x1D70E}_{0}/[1+E_{0}\ln (1-x_{in})]$ . Binary and ternary breakups are represented by discrete symbols without a ‘+’ and discrete symbols with a ‘+’ respectively. Here, $\unicode[STIX]{x1D6E5}$ is the decrement of the critical capillary number due to the presence of surfactants, and $\unicode[STIX]{x1D6E5}_{1}$ and $\unicode[STIX]{x1D6E5}_{2}$ are its constituents. Each inset plots the droplet shapes before and after the breakup.

As mentioned above, the presence of surfactants decreases the critical capillary number of droplet breakup for each confinement ratio, and the decrease in the critical capillary number is also attributed to two mechanisms: (1) the reduction of the interfacial tension caused by the average surfactant concentration $\unicode[STIX]{x1D713}_{0}$ and (2) the non-uniform effects associated with non-uniform capillary pressure and Marangoni stresses. The contributions from these two mechanisms are denoted as $\unicode[STIX]{x1D6E5}_{1}$ and $\unicode[STIX]{x1D6E5}_{2}$ , which are shown in figure 15. It should be noted in this figure that the dashed and dash-dotted lines are directly taken from figure 12, and the solid lines with inverted triangles are obtained by ruling out the reduction of the interfacial tension caused by $\unicode[STIX]{x1D713}_{0}$ in the surfactant-laden cases. For each of the confinement ratios considered, $\unicode[STIX]{x1D6E5}_{1}$ always has a positive contribution to the decrease in the critical capillary number. Moreover, we surprisingly find that the non-uniform effects ( $\unicode[STIX]{x1D6E5}_{2}$ ) increase $Ca_{cr}$ , thus retarding the droplet breakup, for low confinement ratios where binary breakup occurs, but decrease $Ca_{cr}$ , thus promoting the droplet breakup, for high confinement ratios where ternary breakup occurs. The difference at low and high confinement ratios is probably explained as follows. For low confinement ratios, it can be observed from, e.g., figure 14(a) that in the entire stretching stage, the minimum value of surfactant concentration is located in the middle of the droplet, where the interfacial tension is higher than the average value. Since the interfacial tension acts to stabilize the interface, the surfactants retard the droplet deformation in the middle region, thus against the formation of a neck. On the other hand, for high confinement ratios, it is noticed from, e.g., figure 14(b) that the necks are formed in the regions away from the droplet centre, where the surfactant concentration is relatively high. The resulting low interfacial tension favours the formation of necks, thus promoting the droplet breakup.

4.2.2 Influence of the Reynolds number

Another important parameter governing droplet breakup is the Reynolds number, and its influence on the critical capillary number is investigated for both clean and surfactant-laden droplets at a confinement ratio of 0.5. The Reynolds number is varied from 0.05 to 10. For surfactant-laden droplets, the surfactant coverage and Péclet number are still taken as $x_{in}=0.3$ and $Pe=10$ .

Figure 16. The critical capillary number of droplet breakup as a function of the Reynolds number at a confinement ratio of 0.5 for clean droplets (dashed lines with triangles), surfactant-laden droplets (dash-dotted lines with diamonds), and surfactant-laden droplets with a reduced interfacial tension of $\unicode[STIX]{x1D70E}_{e}$ (solid lines with inverted triangles). Binary and ternary breakups are represented by discrete symbols without a ‘+’ and discrete symbols with a ‘+’ respectively. Here, $\unicode[STIX]{x1D6E5}$ is the decrement of the critical capillary number due to the presence of surfactants, and $\unicode[STIX]{x1D6E5}_{1}$ and $\unicode[STIX]{x1D6E5}_{2}$ are its constituents. Each inset plots the droplet shapes before and after the breakup.

Figure 16 displays the simulation results of our method for clean and surfactant-laden droplets. Again, the critical capillary number is defined as the lowest capillary number found at which an initially spherical droplet breaks up. As the Reynolds number increases, the critical capillary number first remains almost unchanged until $Re=0.1$ and then decreases monotonically for the clean droplets, which is consistent with previous findings in a surfactant-free system (Renardy & Cristini Reference Renardy and Cristini2001). For the surfactant-laden droplets, the critical capillary number varies with $Re$ with the same trend as for the clean droplets. In either clean or surfactant-laden cases, binary breakup occurs for $Re\leqslant 2.5$ , whereas for $Re=10$ , ternary breakup occurs, and among the formed droplets, the middle one is very small. In addition, the presence of surfactants can decrease the value of $Ca_{cr}$ for each $Re$ , but the influence of surfactants on $Ca_{cr}$ is smaller in ternary breakup than in binary breakup. This is because in ternary breakup, the dominant convection sweeps more surfactants to both tips of the droplet, which leads to a lower surfactant concentration in the middle part of the droplet, thereby weakening the effect of surfactants on the breakup.

Figure 16 also plots $Ca_{cr}$ as a function of $Re$ for a surfactant-laden droplet with a reduced interfacial tension of $\unicode[STIX]{x1D70E}_{e}$ (represented by the solid lines with inverted triangles). For each $Re$ , $\unicode[STIX]{x1D6E5}_{1}$ , which accounts for the reduction of the interfacial tension by $\unicode[STIX]{x1D713}_{0}$ , always has a positive contribution to the decrease of $Ca_{cr}$ , consistent with the previous observation in figure 15. Depending on $Re$ , the non-uniform effects $\unicode[STIX]{x1D6E5}_{2}$ play different roles in influencing the droplet breakup. Specifically, the non-uniform effects prevent droplet breakup for low $Re$ , but promote breakup for high $Re$ . Moreover, we notice that the non-uniform effects on the droplet breakup are overall smaller at high $Re$ than at low $Re$ , which is attributed to the enhanced surfactant dilution at high $Re$ .