Hostname: page-component-7c8c6479df-94d59 Total loading time: 0 Render date: 2024-03-18T09:46:03.340Z Has data issue: false hasContentIssue false

A hybrid lattice Boltzmann and finite difference method for droplet dynamics with insoluble surfactants

Published online by Cambridge University Press:  21 December 2017

Haihu Liu*
Affiliation:
School of Energy and Power Engineering, Xi’an Jiaotong University, 28 West Xianning Road, Xi’an  710049, China
Yan Ba
Affiliation:
School of Astronautics, Northwestern Polytechnical University, 127 West Youyi Road, Xi’an 710072, China
Lei Wu
Affiliation:
James Weir Fluids Laboratory, Department of Mechanical and Aerospace Engineering, University of Strathclyde, Glasgow G1 1XJ, UK
Zhen Li
Affiliation:
School of Energy and Power Engineering, Xi’an Jiaotong University, 28 West Xianning Road, Xi’an  710049, China
Guang Xi
Affiliation:
School of Energy and Power Engineering, Xi’an Jiaotong University, 28 West Xianning Road, Xi’an  710049, China
Yonghao Zhang
Affiliation:
James Weir Fluids Laboratory, Department of Mechanical and Aerospace Engineering, University of Strathclyde, Glasgow G1 1XJ, UK
*
Email address for correspondence: haihu.liu@xjtu.edu.cn

Abstract

Droplet dynamics in microfluidic applications is significantly influenced by surfactants. It remains a research challenge to model and simulate droplet behaviour including deformation, breakup and coalescence, especially in the confined microfluidic environment. Here, we propose a hybrid method to simulate interfacial flows with insoluble surfactants. The immiscible two-phase flow is solved by an improved lattice Boltzmann colour-gradient model which incorporates a Marangoni stress resulting from non-uniform interfacial tension, while the convection–diffusion equation which describes the evolution of surfactant concentration in the entire fluid domain is solved by a finite difference method. The lattice Boltzmann and finite difference simulations are coupled through an equation of state, which describes how surfactant concentration influences interfacial tension. Our method is first validated for the surfactant-laden droplet deformation in a three-dimensional (3D) extensional flow and a 2D shear flow, and then applied to investigate the effect of surfactants on droplet dynamics in a 3D shear flow. Numerical results show that, at low capillary numbers, surfactants increase droplet deformation, due to reduced interfacial tension by the average surfactant concentration, and non-uniform effects from non-uniform capillary pressure and Marangoni stresses. The role of surfactants on the critical capillary number ($Ca_{cr}$) of droplet breakup is investigated for various confinements (defined as the ratio of droplet diameter to wall separation) and Reynolds numbers. For clean droplets, $Ca_{cr}$ first decreases and then increases with confinement, and the minimum value of $Ca_{cr}$ is reached at a confinement of 0.5; for surfactant-laden droplets, $Ca_{cr}$ exhibits the same variation in trend for confinements lower than 0.7, but, for higher confinements, $Ca_{cr}$ is almost a constant. The presence of surfactants decreases $Ca_{cr}$ for each confinement, and the decrease is also attributed to the reduction in average interfacial tension and non-uniform effects, which are found to prevent droplet breakup at low confinements but promote breakup at high confinements. In either clean or surfactant-laden cases, $Ca_{cr}$ first remains almost unchanged and then decreases with increasing Reynolds number, and a higher confinement or Reynolds number favours ternary breakup. Finally, we study the collision of two equal-sized droplets in a shear flow in both surfactant-free and surfactant-contaminated systems with the same effective capillary numbers. It is identified that the non-uniform effects in the near-contact interfacial region immobilize the interfaces when two droplets are approaching each other and thus inhibit their coalescence.

Type
JFM Papers
Creative Commons
Creative Common License - CCCreative Common License - BY
This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted re-use, distribution, and reproduction in any medium, provided the original work is properly cited.
Copyright
© 2017 Cambridge University Press

1 Introduction

Emulsions, which consist of droplets of one fluid dispersed in another fluid, are common in our daily lives and can be found in the production processes of food, personal care products and pharmaceuticals (Sjöblom Reference Sjöblom2005; Seemann et al. Reference Seemann, Brinkmann, Pfohl and Herminghaus2011). Since emulsion properties such as the stability, rheology and particle morphology are heavily dependent on the droplet size distribution, it is of great importance to control the droplet deformation, breakup and coalescence during emulsification. Once the emulsification device and the compositions of emulsions have been determined, the control of droplet dynamics can be realized through the addition of surfactants (Kobayashi et al. Reference Kobayashi, Nakajima, Chun, Kikuchi and Fujita2002; Kobayashi, Nakajima & Mukataka Reference Kobayashi, Nakajima and Mukataka2003; Kobayashi, Mukataka & Nakajima Reference Kobayashi, Mukataka and Nakajima2005; Van der Graaf et al. Reference Van der Graaf, Steegmans, Van der Sman, Schroën and Boom2005; Kobayashi et al. Reference Kobayashi, Takano, Maeda, Wada, Uemura and Nakajima2008; Chen et al. Reference Chen, Liu, Xu and Luo2015; Josephides & Sajjadi Reference Josephides and Sajjadi2015).

In the investigation of droplet dynamics, numerical simulation is playing an increasingly important role, as it can provide easier access to quantities such as droplet deformation and orientation as well as the velocity, pressure and surfactant concentration distributions inside and outside the droplet, which are difficult to measure experimentally. However, the efficient and accurate computational modelling of droplet dynamics with surfactants in a multiphase system remains a challenging task. The surfactant concentration at the interface alters the interfacial tension locally. The non-uniform surfactant concentration creates non-uniform interfacial tension forces and Marangoni stresses along the interface, which in turn affect the flow field in a complicated manner. Meanwhile, the flow field will influence the surfactant distribution. The interaction between surfactants and the flow field is highly nonlinear. Numerically, the surfactant evolution equation must be solved together with the hydrodynamic equations on moving and deforming interfaces that may undergo topological changes such as breakup and coalescence. In addition, for soluble surfactants, they are not only able to convect and diffuse along the interface, but also to be adsorbed to or desorbed from the interface. The exchange of surfactants between the interface and the bulk phase further increases the computational complexity.

Numerical methods to simulate interfacial flows with surfactants based on conventional Navier–Stokes solvers can roughly be divided into two categories: interface-tracking and interface-capturing methods. Interface-tracking methods, including mainly the boundary integral method (BIM) (Stone Reference Stone1990; Stone & Leal Reference Stone and Leal1990; Milliken & Leal Reference Milliken and Leal1994; Stone Reference Stone1994; Eggleton, Pawar & Stebe Reference Eggleton, Pawar and Stebe1999; Johnson & Borhan Reference Johnson and Borhan2000; Eggleton, Tsai & Stebe Reference Eggleton, Tsai and Stebe2001; Kruijt-Stegeman, van de Vosse & Meijer Reference Kruijt-Stegeman, van de Vosse and Meijer2004; Bazhlekov, Anderson & Meijer Reference Bazhlekov, Anderson and Meijer2006; Feigl et al. Reference Feigl, Megias-Alguacil, Fischer and Windhab2007), the front-tracking method (Zhang, Eckmann & Ayyaswamy Reference Zhang, Eckmann and Ayyaswamy2006; Muradoglu & Tryggvason Reference Muradoglu and Tryggvason2008) and the immersed boundary method (Jin & Stebe Reference Jin and Stebe2007; Lai, Tseng & Huang Reference Lai, Tseng and Huang2008, Reference Lai, Tseng and Huang2010), use the Lagrangian approach to explicitly represent the interface. The BIM employs a surface mesh to track the interface. It was first proposed by Stone & Leal (Reference Stone and Leal1990) to simulate the droplet deformation with insoluble surfactants, in which the surfactant concentration was obtained by solving a time-dependent convective–diffusion equation (Stone Reference Stone1990). This method was extended to 3D flows by Bazhlekov et al. (Reference Bazhlekov, Anderson and Meijer2006) and Feigl et al. (Reference Feigl, Megias-Alguacil, Fischer and Windhab2007), and to soluble surfactants by Milliken & Leal (Reference Milliken and Leal1994). Eggleton et al. (Reference Eggleton, Pawar and Stebe1999) found that the surfactants accumulate at the droplet tips and may drive the interfacial tension to near zero, and thus they improved the BIM by accounting for the maximum limit of surfactant concentration. To date, the BIM has become the most successful method for studying the deformation of surfactant-laden droplets. Another popular interface-tracking method is the front-tracking method (Tryggvason et al. Reference Tryggvason, Bunner, Esmaeeli, Juric, Al-Rawahi, Tauber, Han, Nas and Jan2001), in which a fixed grid is used to determine the fluid flow while the interface is tracked using a lower-dimensional grid; this method for soluble surfactants was further developed by Zhang et al. (Reference Zhang, Eckmann and Ayyaswamy2006) and Muradoglu & Tryggvason (Reference Muradoglu and Tryggvason2008). The immersed boundary method for interfacial flows with insoluble surfactants was proposed by Lai et al. (Reference Lai, Tseng and Huang2008), in which a mixture of Eulerian variables in the fluid domain and Lagrangian variables represented by a set of discrete Lagrangian markers is used, and this method was later extended to contact line simulations (Lai et al. Reference Lai, Tseng and Huang2010). Although interface-tracking methods have been widely used to simulate the dynamics of surfactant-laden droplets, they have difficulty in dealing with droplets undergoing topological changes such as breakup and coalescence.

Interface-capturing methods use an indicator function to implicitly represent the interface in an Eulerian grid. This simplifies gridding, discretization and handling of topological changes. Commonly used interface-capturing methods are the volume-of-fluid method (James & Lowengrub Reference James and Lowengrub2004; Alke & Bothe Reference Alke and Bothe2009), the level-set method (Ceniceros Reference Ceniceros2003; Deng, Ito & Li Reference Deng, Ito and Li2003; Xu & Zhao Reference Xu and Zhao2003; Xu et al. Reference Xu, Li, Lowengrub and Zhao2006; Xu, Yang & Lowengrub Reference Xu, Yang and Lowengrub2012) and the phase-field method (Teigen et al. Reference Teigen, Li, Lowengrub, Wang and Voigt2009, Reference Teigen, Song, Lowengrub and Voigt2011). James & Lowengrub (Reference James and Lowengrub2004) proposed an axisymmetric volume-of-fluid method for insoluble surfactants on a moving fluid interface, in which they tracked the surfactant mass and surface area independently instead of directly solving the surfactant concentration. This model was extended to 3D flows with soluble surfactants by Alke & Bothe (Reference Alke and Bothe2009). Based on an Eulerian formulation of the surfactant concentration equation, Xu & Zhao (Reference Xu and Zhao2003) developed a level-set method for interfacial Stokes flows with insoluble surfactants, in which the surfactant concentration was extended to a neighbourhood of the interface through an additional Hamilton–Jacobian equation. This method was improved with the use of an immersed interface method for the Stokes flow field (Xu et al. Reference Xu, Li, Lowengrub and Zhao2006) and a continuum surface force model for the stress jump condition (Xu et al. Reference Xu, Yang and Lowengrub2012), and was recently extended to model contact line dynamics (Xu & Ren Reference Xu and Ren2014). Teigen et al. (Reference Teigen, Song, Lowengrub and Voigt2011) developed a phase-field method to simulate soluble surfactants, in which the interface structure was resolved via a free energy functional based on the Cahn–Hilliard theory. With the aid of a delta function, the surfactant evolution equations in the bulk phases and on the interface were extended to the entire fluid domain by incorporating additional terms to approximate the physical boundary condition at the interface. The resulting surfactant evolution equations are known as the surfactant evolution equations of diffuse-interface form in order to distinguish from their equivalence to the sharp-interface form (Xu & Zhao Reference Xu and Zhao2003). The surfactant evolution equation of diffuse-interface form does not need additional procedures to extend the surfactant concentration to a neighbourhood of the interface, greatly simplifying the modelling of surfactant dynamics. Despite their great success in simulating droplet breakup and coalescence, the volume-of-fluid and level-set methods require either sophisticated interface reconstruction algorithms or unphysical re-initialization processes to represent the interface and surfactant concentration, and the phase-field method yields an interface thickness far greater than its actual value, which may lead to unphysical dissolution of small droplets (Zhang & Wang Reference Zhang and Wang2010) and mobility-dependent numerical results (van der Sman & van der Graaf Reference van der Sman and van der Graaf2008; Magaletti et al. Reference Magaletti, Picano, Chinappi, Marino and Casciola2013). It still remains an open question for the phase-field method to choose an appropriate interface thickness so that the numerical difficulties become amenable while the correct physics can be captured.

Recently, the lattice Boltzmann (LB) method has shown great potential for modelling interfacial interactions. It is a pseudomolecular method tracking the evolution of the distribution function of an assembly of molecules and built upon microscopic models and mesoscopic kinetic equations. Its mesoscopic nature can provide many advantages of molecular dynamics, making the LB method particularly useful for simulation of multiphase multicomponent flows. A number of multiphase multicomponent models have been proposed, which can be classified into four main categories: the colour-gradient model (Gunstensen et al. Reference Gunstensen, Rothman, Zaleski and Zanetti1991; Halliday et al. Reference Halliday, Law, Care and Hollis2006; Liu, Valocchi & Kang Reference Liu, Valocchi and Kang2012a ), the phase-field-based model (Swift, Osborn & Yeomans Reference Swift, Osborn and Yeomans1995; Pooley, Kusumaatmaja & Yeomans Reference Pooley, Kusumaatmaja and Yeomans2008; Wang et al. Reference Wang, Shu, Huang and Teo2015; Kusumaatmaja, Hemingway & Fielding Reference Kusumaatmaja, Hemingway and Fielding2016), the interparticle-potential model (Shan & Chen Reference Shan and Chen1993; Benzi et al. Reference Benzi, Biferale, Sbragaglia, Succi and Toschi2006; Sbragaglia et al. Reference Sbragaglia, Benzi, Biferale, Succi, Sugiyama and Toschi2007; Falcucci et al. Reference Falcucci, Ubertini, Bella, Maio and Palpacelli2010a ; Falcucci, Ubertini & Succi Reference Falcucci, Ubertini and Succi2010b ; Falcucci et al. Reference Falcucci, Jannelli, Ubertini and Succi2013; Montessori et al. Reference Montessori, Falcucci, Rocca, Ansumali and Succi2015) and the mean-field theory model (He & Doolen Reference He and Doolen2002). Among these models, the colour-gradient model has been extensively used to simulate immiscible multiphase flow problems because of its strengths such as low spurious velocities, high numerical accuracy, strict mass conservation for each fluid and good numerical stability for a broad range of fluid properties. Recently, it was extended to model thermocapillary flows (Liu, Zhang & Valocchi Reference Liu, Zhang and Valocchi2012b ; Liu & Zhang Reference Liu and Zhang2015) and contact line dynamics with contact angle hysteresis (Ba et al. Reference Ba, Liu, Sun and Zheng2013; Liu et al. Reference Liu, Ju, Wang, Xi and Zhang2015). The colour-gradient model has recently been attempted to simulate interfacial flows with insoluble surfactants by Farhat et al. (Reference Farhat, Celiker, Singh and Lee2011). However, their model is not viable because of the following limitations. (1) The model solves a surfactant evolution equation of sharp-interface form (Stone & Leal Reference Stone and Leal1990) over a diffuse interface, and incompatibility would arise as a volume-distribution surfactant concentration should be used instead of its surface counterpart from a diffuse-interface point of view. (2) It solves the surfactant evolution equation in the interfacial region only rather than in the entire fluid domain, and thus initialization techniques for the surfactant concentration at the emerging interfacial nodes need to be developed, which would inevitably increase the modelling complexity and errors. (3) It does not consider the Marangoni stress, which are unphysical and could have important effects on the droplet behaviour. In addition, this model may not be suitable for extension to more complicated systems, such as soluble surfactant solutions.

To model and simulate droplet dynamics with insoluble surfactants, a hybrid method is proposed with a recently improved LB colour-gradient model to describe the immiscible two-phase flow, and a finite difference (FD) method to solve a convection–diffusion equation of diffuse-interface form, originally designed in the context of the phase-field method, for surfactant transport. The interfacial tension force and the Marangoni stress are both modelled using the concept of a continuum surface force. An equation of state that describes how the interfacial tension is influenced by the surfactant concentration will couple the LB and FD simulations dynamically. The capability and accuracy of this hybrid method will be assessed by simulating the deformation of a surfactant-laden droplet in a 3D extensional flow and in a 2D shear flow. It will then be used to investigate the effect of surfactants on the deformation and breakup of a single droplet in a 3D shear flow, where the critical conditions of droplet breakup are explored in a surfactant-laden system with varying confinement ratios and Reynolds numbers for the first time. Finally, we will study the collision of two equal-sized droplets in a shear flow.

2 Numerical method

We develop a hybrid LB–FD method to simulate the interfacial dynamics with insoluble surfactants. This method uses a recently improved LB colour-gradient model (Liu et al. Reference Liu, Zhang and Valocchi2012b ) to simulate immiscible two-phase flows. Compared with the commonly used LB colour-gradient models (Lishchuk, Care & Halliday Reference Lishchuk, Care and Halliday2003; Halliday, Hollis & Care Reference Halliday, Hollis and Care2007; Liu et al. Reference Liu, Valocchi and Kang2012a ; Ba et al. Reference Ba, Liu, Sun and Zheng2013), an additional interfacial Marangoni force that arises from the interfacial tension gradient due to the non-uniform distribution of surfactants is incorporated into the perturbation term. To be compatible with the diffuse nature of the two-phase interface obtained by the colour-gradient model, a convection–diffusion equation of diffuse-interface form to describe the transport of the surfactant concentration is introduced and solved by the FD method. The use of an evolution equation of diffuse-interface form not only allows us to solve the surfactant concentration in the entire fluid domain without additional extension or initialization procedures, but also offers great simplicity in dealing with topological changes such as the interface rupturing and merging. However, it should be noted that the present method is limited to binary fluids with similar densities and to relatively low Reynolds number, which could be improved by using the multiple-relaxation-time model developed in Ba et al. (Reference Ba, Liu, Li, Kang and Sun2016).

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 Numerical validation

The hybrid LB–FD method is validated against the available literature data for a surfactant-laden droplet subject to a 3D extensional flow and a 2D shear flow.

Figure 1. Illustration of a droplet in a 3D extensional flow.

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 Deformation and breakup of a droplet in a 3D shear flow

Having validated the hybrid LB–FD method, we investigate the influence of surfactants on the droplet deformation and breakup in a 3D shear flow. As depicted in figure 5(b), a 3D droplet (red fluid) with radius $R$ is initially placed in the centre of two plates, both parallel to the $x$ $y$ plane at distance $H$ apart. The two plates move at equal speeds of $u_{w}$ but in opposite directions, thereby creating a linear shear flow with a shear rate of $\dot{\unicode[STIX]{x1D6FE}}=2u_{w}/H$ . Periodic boundary conditions are applied in the $x$ -direction, while halfway bounce-back boundary conditions are imposed on the top and bottom walls with constant velocities of $(u_{w},0,0)$ and $(-u_{w},0,0)$ , respectively. Due to symmetry, only half of the flow domain in the $y$ -direction is considered. Assuming that the symmetric boundary is positioned at $y=0.5$ , we place an extra grid layer at $y=0$ to facilitate the implementation of a symmetric boundary condition. For the distribution functions, symmetry is achieved by setting $f_{i^{\ast }}^{k\ddagger }(x,0,z)=f_{i}^{k\ddagger }(x,1,z)$ before the propagation step, where $\boldsymbol{e}_{i^{\ast }}=\boldsymbol{e}_{i}-2(\boldsymbol{e}_{i}\boldsymbol{\cdot }\bar{\boldsymbol{n}})\bar{\boldsymbol{n}}$ , with $\bar{\boldsymbol{n}}=(0,1,0)$ . Besides, the symmetry requires that $\unicode[STIX]{x1D719}(x,0,z)=\unicode[STIX]{x1D719}(x,1,z)$ for a scalar $\unicode[STIX]{x1D719}$ and $v_{x}(x,0,z)=v_{x}(x,1,z)$ , $v_{y}(x,0,z)=-v_{y}(x,1,z)$ , $v_{z}(x,0,z)=v_{z}(x,1,z)$ for a vector $\boldsymbol{v}$ . Throughout this section, both fluids are of equal density and equal viscosity, with their values depending on the specific flow conditions. The Langmuir equation (2.6) is employed to relate the interfacial tension to the surfactant concentration, in which the interfacial tension of a clean interface is taken as $\unicode[STIX]{x1D70E}_{0}=10^{-3}$ and the elasticity number is $E_{0}=0.5$ .

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

5 Collision of two droplets in a simple shear flow

One of the important functions of surfactant in microfluidic applications is to inhibit the coalescence of droplets (Link et al. Reference Link, Anna, Weitz and Stone2004; Baret Reference Baret2012; Krebs, Schroën & Boom Reference Krebs, Schroën and Boom2012). To demonstrate the capability of the hybrid method in dealing with interface merging and to explore how a non-uniform distribution of surfactants affects droplet coalescence, we consider the collision of two equal-sized droplets subject to a simple shear flow in both surfactant-free and surfactant-contaminated systems (see figure 17).

Figure 17. Illustration of the collision of two droplets in a 2D shear flow. It should be noted that the point ‘A’ is the left intersection point of the right droplet and the horizontal centre line, and the points ‘B’ and ‘C’ are the near tip and the far tip of the right droplet.

The computational domain is taken as $720\times 480$ , and two circular droplets with equal radius $R=80$ are initially placed at (200, 280) and (540, 200). In the surfactant-contaminated system, a surfactant coverage of $x_{in}=0.3$ is considered. Other dimensionless parameters are selected as $Re=0.4$ , $E_{0}=0.4$ and $Pe=10$ . It is known that a reduction of the interfacial tension, and thus an increase of the effective capillary number, due to the presence of surfactants can inhibit droplet coalescence (Shardt, Derksen & Mitra Reference Shardt, Derksen and Mitra2013). To rule out the effect of average surfactant concentration in reducing the interfacial tension, the effective capillary number

(5.1) $$\begin{eqnarray}Ca_{e}=\frac{Ca}{1+E_{0}\ln (1-x_{in})}=0.1\end{eqnarray}$$

is used in both the surfactant-free and surfactant-contaminated systems.

Figure 18. Collision of two droplets in simple shear flow for clean (the dashed lines) and surfactant-laden cases at $Re=0.4$ and $Ca_{e}=0.1$ . In the surfactant-laden case, $Pe=10$ , $x_{in}=0.3$ and $E_{0}=0.4$ , and the interface is coloured by the dimensionless surfactant concentration $\unicode[STIX]{x1D713}^{\ast }$ .

Figure 18 shows the evolution of the droplet shape for the clean droplets and surfactant-contaminated droplets at times ( $\unicode[STIX]{x1D70F}=\dot{\unicode[STIX]{x1D6FE}}t$ ) of 0, 2.3, 4.6, 5.6, 6.4 and 8.1. In this figure, we plot the contour of $\unicode[STIX]{x1D70C}^{N}=0$ for the clean droplets and the surfactant concentration distributions for the surfactant-contaminated droplets. It is seen that the clean droplets first move towards each other (figure 18 a,b), then the film between the droplets becomes very thin (figure 18 c), and finally the two droplets merge together as a result of the interfacial tension (figure 18 df). On the contrary, for surfactant-contaminated droplets, the thin film between the droplets maintains a certain thickness, and the droplets slide over each other without coalescence. This clearly shows that the non-uniform effects from non-uniform capillary pressure and Marangoni stresses act as an additional repulsive force to prevent droplet coalescence, which is consistent with previous experimental and numerical findings (Hu, Pine & Leal Reference Hu, Pine and Leal2000; Dai & Leal Reference Dai and Leal2008; Liu & Zhang Reference Liu and Zhang2010).

Figure 19. 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_{Ma}$ along the arclength $s$ for $Re=0.4$ , $Pe=10$ , $Ca_{e}=0.1$ , $E_{0}=0.4$ , $x_{in}=0.3$ at dimensionless times of $\unicode[STIX]{x1D70F}=0.58$ , 2.3, 4.6 and 5.8. It should be noted that the near tip represents the one first encountered as $s$ increases.

As shown in figure 18, for the surfactant-contaminated droplets, surfactants are swept towards droplet tips at the beginning, but the maximum values of the surfactant concentration at the two tips are unequal due to the asymmetry of the shear flow with respect to the droplet centre and the effect from the neighbouring droplet. As the two droplets move close to each other, the symmetry of the surfactant concentration distribution is further broken, and we can observe that the surfactant concentration is greatly enhanced at the near tip (point ‘B’ in figure 17), but remains almost unchanged at the far tip (point ‘C’ in figure 17). The variation of the surfactant concentration can be observed more clearly in figure 19, which plots the distributions of the dimensionless surfactant concentration $\unicode[STIX]{x1D713}^{\ast }$ , the dimensionless interfacial tension $\unicode[STIX]{x1D70E}^{\ast }$ , and the signed magnitudes of the dimensionless capillary force $F_{Ca}$ and Marangoni force $F_{Ma}$ along the arclength $s$ at times of $\unicode[STIX]{x1D70F}=0.58$ , 2.3, 4.6 and 5.8. Here, $s$ is measured from the left intersection point (point ‘A’ in figure 17) between the droplet interface and the horizontal centreline, and increases anticlockwise along the interface. As shown in figure 19(a), four stages can be identified in the evolution of the surfactant concentration: (1) in the early stage, e.g. $\unicode[STIX]{x1D70F}=0.58$ , the surfactants predominantly accumulate at the droplet tips due to the outer shear flow, producing nearly equal maximum surfactant concentration; (2) as $\unicode[STIX]{x1D70F}$ increases from 0.58 to 4.6, $\unicode[STIX]{x1D713}^{\ast }$ at the far tip almost holds a constant value while $\unicode[STIX]{x1D713}^{\ast }$ at the near tip increases; (3) when the gap between the two droplets decreases to a minimum value at $\unicode[STIX]{x1D70F}\simeq 4.6$ , the maximum $\unicode[STIX]{x1D713}^{\ast }$ is achieved at the far tip; and (4) as the droplets slide over each other (see, e.g., $\unicode[STIX]{x1D70F}=5.8$ ), $\unicode[STIX]{x1D713}^{\ast }$ at the far tip decreases, but $\unicode[STIX]{x1D713}^{\ast }$ at the near tip increases.

Figure 20. The velocity vectors and the streamlines in the vicinity of the gap at $\unicode[STIX]{x1D70F}=4.6$ for $Re=0.4$ , $Pe=10$ , $x_{in}=0.3$ , $E_{0}=0.4$ and $Ca_{e}=0.1$ : (a) clean droplets; (b) surfactant-laden droplets.

The distribution of $\unicode[STIX]{x1D713}^{\ast }$ determines the interfacial tension, thus influencing the outcome of the droplet collision. From figure 19(b), we see that the dimensionless interfacial tension $\unicode[STIX]{x1D70E}^{\ast }$ is distributed unevenly on the droplet interface, and its non-uniformity increases as the droplets come close to each other, which leads to increased non-uniform effects associated with non-uniform capillary forces and Marangoni stresses (see figure 19 d). The non-uniform capillary forces and Marangoni stresses are able to reduce the interface mobility during the droplet–droplet interaction (an indicator is the broader gap between the droplets in the presence of surfactants, which can be seen from the comparison between figures 20 a and 20 b), thus hindering the coalescence of droplets. Interestingly, we find that at $s=0$ , as the droplets approach each other, $\unicode[STIX]{x1D713}^{\ast }$ does not change much, but $F_{Ca}$ increases drastically, which is attributed to the increase in the local interface curvature. In figure 19(d), it is clearly shown that $F_{Ma}$ is fixed at 0 at both tips, and the maximum or minimum values of $F_{Ma}$ always occur at locations near the far tip. Finally, it is worth noting that the addition of surfactants can lead to considerable variation of the velocity vectors and the streamlines in the vicinity of the gap, which can be seen from figure 20. For clean droplets, a single vortex is formed with its centre in the gap, whereas for surfactant-contaminated droplets, a counter-rotating vortex pair is observed, with each vortex centred near the gap inside its corresponding droplet.

6 Conclusions

In summary, a hybrid LB–FD method has been proposed to simulate interfacial flows with insoluble surfactants. The hybrid method solves the immiscible two-phase flows by using the LB colour-gradient model, in which the Marangoni stress is incorporated into the perturbation term to account for the tangential stress created by the non-uniform interfacial tension, and describes the transport of surfactant concentration through a convection–diffusion equation of diffuse-interface form, which is solved in the entire fluid domain by the FD method. The LB and FD methods are coupled by an EOS that relates the interfacial tension to the surfactant concentration and can be selected arbitrarily. The accuracy of the hybrid method has been validated for the deformation of a surfactant-laden droplet in an 3D extensional flow and a 2D shear flow.

The hybrid method has been used to study the effect of surfactants on droplet deformation and breakup in a three-dimensional shear flow. It has been found that the presence of surfactants increases droplet deformation due to (1) a reduction of the interfacial tension caused by the average surfactant concentration and (2) the non-uniform effects from non-uniform capillary pressure and Marangoni stresses. The former accounts for most of the effect of the surfactants on the droplet deformation. Increase of the surface Péclet number can increase the non-uniformity of the surfactant distribution, thus leading to enhanced droplet deformation. Regarding the droplet breakup, we have focused on the critical capillary number ( $Ca_{cr}$ ) of droplet breakup for various confinement ratios ( $2R/H$ ) and Reynolds numbers. For clean droplets, $Ca_{cr}$ first increases and then decreases with the confinement ratio, and the minimum $Ca_{cr}$ is achieved at a confinement ratio of around 0.5. For surfactant-laden droplets, $Ca_{cr}$ varies with the confinement ratio with the same trend as for clean droplets for $2R/H<0.7$ , but for higher confinement ratios, $Ca_{cr}$ remains almost constant. The presence of surfactants decreases $Ca_{cr}$ for each confinement ratio, and the decrease in $Ca_{cr}$ is also attributed to the reduction in average interfacial tension and the non-uniform effects. It has been identified that the non-uniform effects prevent droplet breakup at low confinement ratios where binary breakup occurs, but promote breakup at high confinement ratios where ternary breakup occurs. In either clean or surfactant-laden cases, $Ca_{cr}$ first remains almost unchanged until $Re=0.1$ , and then decreases monotonically with increasing Reynolds number. The mode of droplet breakup depends on the confinement ratio and Reynolds number: a higher confinement ratio or Reynolds number favours ternary breakup.

We have also simulated the collision of two equal-sized droplets subject to a simple shear flow in both surfactant-free and surfactant-contaminated systems with equal effective capillary numbers. It is observed that the non-uniform effects associated with non-uniform capillary forces and Marangoni stresses can immobilize the interfaces during the droplet–droplet approach, thus acting as a repulsive force to prevent droplet coalescence.

The present study not only provides a simple, accurate and reliable numerical method for the simulation of interfacial flows with insoluble surfactants, but also improves our understanding of droplet dynamics in confined shear flows with the presence of surfactants.

Acknowledgements

This work is financially supported by the National Natural Science Foundation of China (no. 51506168 and 51711530130), the National Key Research and Development Program of China (no. 2016YFB0200902), the China Postdoctoral Science Foundation (no. 2016M590943), and the Thousand Youth Talents Program for Distinguished Young Scholars. This work is also supported by the UK’s Engineering and Physical Sciences Research Council (EPSRC) under grant EP/L00030X/1.

References

Alke, A. & Bothe, D. 2009 3D numerical modeling of soluble surfactant at fluidic interfaces based on the volume-of-fluid method. Fluid Dyn. Mater. Process. 5 (4), 345372.Google Scholar
Ba, Y., Liu, H., Li, Q., Kang, Q. & Sun, J. 2016 Multiple-relaxation-time color-gradient lattice Boltzmann model for simulating two-phase flows with high density ratio. Phys. Rev. E 94, 023310.Google Scholar
Ba, Y., Liu, H., Sun, J. & Zheng, R. 2013 Color-gradient lattice Boltzmann model for simulating droplet motion with contact-angle hysteresis. Phys. Rev. E 88 (4), 043306.Google Scholar
Baret, J.-C. 2012 Surfactants in droplet-based microfluidics. Lab on a Chip 12 (3), 422433.10.1039/C1LC20582JGoogle Scholar
Bazhlekov, I. B., Anderson, P. D. & Meijer, H. E. 2006 Numerical investigation of the effect of insoluble surfactants on drop deformation and breakup in simple shear flow. J. Colloid Interface Sci. 298 (1), 369394.10.1016/j.jcis.2005.12.017Google Scholar
Benzi, R., Biferale, L., Sbragaglia, M., Succi, S. & Toschi, F. 2006 Mesoscopic modeling of a two-phase flow in the presence of boundaries: the contact angle. Phys. Rev. E 74, 021509.Google Scholar
Ceniceros, H. D. 2003 The effects of surfactants on the formation and evolution of capillary waves. Phys. Fluids 15 (1), 245256.10.1063/1.1528940Google Scholar
Chen, Y., Liu, G.-T., Xu, J.-H. & Luo, G.-S. 2015 The dynamic mass transfer of surfactants upon droplet formation in coaxial microfluidic devices. Chem. Engng Sci. 132, 18.10.1016/j.ces.2015.04.006Google Scholar
Cristini, V. & Tan, Y.-C. 2004 Theory and numerical simulation of droplet dynamics in complex flows: a review. Lab on a Chip 4, 257264.10.1039/B403226HGoogle Scholar
Dai, B. & Leal, L. G. 2008 The mechanism of surfactant effects on drop coalescence. Phys. Fluids 20, 040802.10.1063/1.2911700Google Scholar
Deng, S., Ito, K. & Li, Z. 2003 Three-dimensional elliptic solvers for interface problems and applications. J. Comput. Phys. 184 (1), 215243.10.1016/S0021-9991(02)00028-1Google Scholar
Eggleton, C. D., Pawar, Y. P. & Stebe, K. J. 1999 Insoluble surfactants on a drop in an extensional flow: a generalization of the stagnated surface limit to deforming interfaces. J. Fluid Mech. 385, 7999.10.1017/S0022112098004054Google Scholar
Eggleton, C. D., Tsai, T.-M. & Stebe, K. J. 2001 Tip streaming from a drop in the presence of surfactants. Phys. Rev. Lett. 87 (4), 048302.10.1103/PhysRevLett.87.048302Google Scholar
Falcucci, G., Jannelli, E., Ubertini, S. & Succi, S. 2013 Direct numerical evidence of stress-induced cavitation. J. Fluid Mech. 728, 362375.10.1017/jfm.2013.271Google Scholar
Falcucci, G., Ubertini, S., Bella, G., Maio, A. D. & Palpacelli, S. 2010a Lattice Boltzmann modeling of diesel spray formation and break-up. SAE Intl J. Fuels Lubr. 3 (1), 582593.10.4271/2010-01-1130Google Scholar
Falcucci, G., Ubertini, S. & Succi, S. 2010b Lattice Boltzmann simulations of phase-separating flows at large density ratios: the case of doubly-attractive pseudo-potentials. Soft Matt. 6, 43574365.10.1039/c002974bGoogle Scholar
Farhat, H., Celiker, F., Singh, T. & Lee, J. S. 2011 A hybrid lattice Boltzmann model for surfactant-covered droplets. Soft Matt. 7 (5), 19681985.10.1039/c0sm00569jGoogle Scholar
Feigl, K., Megias-Alguacil, D., Fischer, P. & Windhab, E. J. 2007 Simulation and experiments of droplet deformation and orientation in simple shear flow with surfactants. Chem. Engng Sci. 62 (12), 32423258.10.1016/j.ces.2007.02.008Google Scholar
Guido, S. & Greco, F. 2004 Dynamics of a liquid drop in a flowing immiscible liquid. Rheol. Rev. 2004, 99142.Google Scholar
Gunstensen, A. K., Rothman, D. H., Zaleski, S. & Zanetti, G. 1991 Lattice Boltzmann model of immiscible fluids. Phys. Rev. A 43 (8), 43204327.10.1103/PhysRevA.43.4320Google Scholar
Guo, Z.-L., Zheng, C.-G. & Shi, B.-C. 2002b Non-equilibrium extrapolation method for velocity and pressure boundary conditions in the lattice Boltzmann method. Chin. Phys. 11 (4), 366.Google Scholar
Guo, Z., Zheng, C. & Shi, B. 2002a Discrete lattice effects on the forcing term in the lattice Boltzmann method. Phys. Rev. E 65, 110126.Google Scholar
Gupta, A. & Kumar, R. 2010 Effect of geometry on droplet formation in the squeezing regime in a microfluidic T-junction. Microfluid Nanofluid 8, 799812.10.1007/s10404-009-0513-7Google Scholar
Halliday, I., Hollis, A. P. & Care, C. M. 2007 Lattice Boltzmann algorithm for continuum multicomponent flow. Phys. Rev. E 76, 026708.Google Scholar
Halliday, I., Law, R., Care, C. M. & Hollis, A. 2006 Improved simulation of drop dynamics in a shear flow at low Reynolds and capillary number. Phys. Rev. E 73 (5), 056708.Google Scholar
He, X. & Doolen, G. D. 2002 Thermodynamic foundations of kinetic theory and lattice Boltzmann models for multiphase flows. J. Stat. Phys. 107 (1), 309328.10.1023/A:1014527108336Google Scholar
Hu, Y. T., Pine, D. J. & Leal, L. G. 2000 Drop deformation, breakup, and coalescence with compatibilizer. Phys. Fluids 12 (3), 484489.10.1063/1.870254Google Scholar
James, A. J. & Lowengrub, J. 2004 A surfactant-conserving volume-of-fluid method for interfacial flows with insoluble surfactant. J. Comput. Phys. 201 (2), 685722.10.1016/j.jcp.2004.06.013Google Scholar
Janssen, P. J. A., Vananroye, A., Van Puyvelde, P., Moldenaers, P. & Anderson, P. D. 2010 Generalized behavior of the breakup of viscous drops in confinements. J. Rheol. 54 (5), 10471060.10.1122/1.3473924Google Scholar
Jiang, G.-S. & Shu, C.-W. 1996 Efficient implementation of weighted ENO schemes. J. Comput. Phys. 126 (1), 202228.10.1006/jcph.1996.0130Google Scholar
Jin, F. & Stebe, K. J. 2007 The effects of a diffusion controlled surfactant on a viscous drop injected into a viscous medium. Phys. Fluids 19 (11), 112103.10.1063/1.2775055Google Scholar
Johnson, R. A. & Borhan, A. 2000 Stability of the shape of a surfactant-laden drop translating at low Reynolds number. Phys. Fluids 12 (4), 773784.10.1063/1.870334Google Scholar
Josephides, D. N. & Sajjadi, S. 2015 Increased drop formation frequency via reduction of surfactant interactions in flow-focusing microfluidic devices. Langmuir 31 (3), 12181224.10.1021/la504299rGoogle Scholar
Kobayashi, I., Mukataka, S. & Nakajima, M. 2005 Effects of type and physical properties of oil phase on oil-in-water emulsion droplet formation in straight-through microchannel emulsification, experimental and CFD studies. Langmuir 21 (13), 57225730.10.1021/la050039nGoogle Scholar
Kobayashi, I., Nakajima, M., Chun, K., Kikuchi, Y. & Fujita, H. 2002 Silicon array of elongated through-holes for monodisperse emulsion droplets. AIChE J. 48 (8), 16391644.10.1002/aic.690480807Google Scholar
Kobayashi, I., Nakajima, M. & Mukataka, S. 2003 Preparation characteristics of oil-in-water emulsions using differently charged surfactants in straight-through microchannel emulsification. Colloids Surf. A 229 (1), 3341.10.1016/j.colsurfa.2003.08.005Google Scholar
Kobayashi, I., Takano, T., Maeda, R., Wada, Y., Uemura, K. & Nakajima, M. 2008 Straight-through microchannel devices for generating monodisperse emulsion droplets several microns in size. Microfluid Nanofluid 4 (3), 167177.10.1007/s10404-007-0167-2Google Scholar
Krebs, T., Schroën, K. & Boom, R. 2012 Coalescence dynamics of surfactant-stabilized emulsions studied with microfluidics. Soft Matt. 8, 1065010657.10.1039/c2sm26122gGoogle Scholar
Kruijt-Stegeman, Y. W., van de Vosse, F. N. & Meijer, H. E. H. 2004 Droplet behavior in the presence of insoluble surfactants. Phys. Fluids 16 (8), 27852796.10.1063/1.1756168Google Scholar
Kusumaatmaja, H., Hemingway, E. J. & Fielding, S. M. 2016 Moving contact line dynamics: from diffuse to sharp interfaces. J. Fluid Mech. 788, 209227.10.1017/jfm.2015.697Google Scholar
Lai, M.-C., Tseng, Y.-H. & Huang, H. 2008 An immersed boundary method for interfacial flows with insoluble surfactant. J. Comput. Phys. 227 (15), 72797293.10.1016/j.jcp.2008.04.014Google Scholar
Lai, M. C., Tseng, Y. H. & Huang, H. 2010 Numerical simulation of moving contact lines with surfactant by immersed boundary method. Commun. Comput. Phys. 8 (4), 735757.10.4208/cicp.281009.120210aGoogle Scholar
Latva-Kokko, M. & Rothman, D. H. 2005 Diffusion properties of gradient-based lattice Boltzmann models of immiscible fluids. Phys. Rev. E 71 (5), 056702.Google Scholar
Link, D. R., Anna, S. L., Weitz, D. A. & Stone, H. A. 2004 Geometrically mediated breakup of drops in microfluidic devices. Phys. Rev. Lett. 92, 054503.10.1103/PhysRevLett.92.054503Google Scholar
Lishchuk, S. V., Care, C. M. & Halliday, I. 2003 Lattice Boltzmann algorithm for surface tension with greatly reduced microcurrents. Phys. Rev. E 67 (3), 036701.Google Scholar
Liu, H., Ju, Y., Wang, N., Xi, G. & Zhang, Y. 2015 Lattice Boltzmann modeling of contact angle and its hysteresis in two-phase flow with large viscosity difference. Phys. Rev. E 92, 033306.Google Scholar
Liu, H., Valocchi, A. J. & Kang, Q. 2012a Three-dimensional lattice Boltzmann model for immiscible two-phase flow simulations. Phys. Rev. E 85 (4), 046309.Google Scholar
Liu, H., Valocci, A. J., Werth, C., Kang, Q. & Oostrom, M. 2014 Pore-scale simulation of liquid CO2 displacement of water using a two-phase lattice Boltzmann model. Adv. Water Resour. 73, 144158.10.1016/j.advwatres.2014.07.010Google Scholar
Liu, H. & Zhang, Y. 2010 Phase-field modeling droplet dynamics with soluble surfactants. J. Comput. Phys. 229 (24), 91669187.10.1016/j.jcp.2010.08.031Google Scholar
Liu, H. & Zhang, Y. 2015 Modeling thermocapillary migration of a microfluidic droplet on a solid surface. J. Comput. Phys. 280, 3753.10.1016/j.jcp.2014.09.015Google Scholar
Liu, H., Zhang, Y. & Valocchi, A. J. 2012b Modeling and simulation of thermocapillary flows using lattice Boltzmann method. J. Comput. Phys. 231 (12), 44334453.10.1016/j.jcp.2012.02.015Google Scholar
Maffettone, P. L. & Minale, M. 1998 Equation of change for ellipsoidal drops in viscous flow. J. Non-Newtonian Fluid Mech. 78 (2), 227241.10.1016/S0377-0257(98)00065-2Google Scholar
Magaletti, F., Picano, F., Chinappi, M., Marino, L. & Casciola, C. M. 2013 The sharp-interface limit of the Cahn–Hilliard/Navier–Stokes model for binary fluids. J. Fluid Mech. 714, 95126.10.1017/jfm.2012.461Google Scholar
Mei, R., Yu, D., Shyy, W. & Luo, L.-S. 2002 Force evaluation in the lattice Boltzmann method involving curved geometry. Phys. Rev. E 65, 041203.Google Scholar
Milliken, W. J. & Leal, L. G. 1994 The influence of surfactant on the deformation and breakup of a viscous drop: the effect of surfactant solubility. J. Colloid Interface Sci. 166 (2), 275285.10.1006/jcis.1994.1296Google Scholar
Minale, M. 2010 Models for the deformation of a single ellipsoidal drop: a review. Rheol. Acta 49 (8), 789806.10.1007/s00397-010-0442-0Google Scholar
Montessori, A., Falcucci, G., Rocca, M. L., Ansumali, S. & Succi, S. 2015 Three-dimensional lattice pseudo-potentials for multiphase flow simulations at high density ratios. J. Stat. Phys. 161 (6), 14041419.10.1007/s10955-015-1318-6Google Scholar
Muradoglu, M. & Tryggvason, G. 2008 A front-tracking method for computation of interfacial flows with soluble surfactants. J. Comput. Phys. 227 (4), 22382262.10.1016/j.jcp.2007.10.003Google Scholar
Pathak, J. A. & Migler, K. B. 2003 Droplet-string deformation and stability during microconfined shear flow. Langmuir 19 (21), 86678674.10.1021/la0346907Google Scholar
Pooley, C. M., Kusumaatmaja, H. & Yeomans, J. M. 2008 Contact line dynamics in binary lattice Boltzmann simulations. Phys. Rev. E 78, 056709.Google Scholar
Puyvelde, P. V., Vananroye, A., Cardinaels, R. & Moldenaers, P. 2008 Review on morphology development of immiscible blends in confined shear flow. Polymer 49 (25), 53635372.10.1016/j.polymer.2008.08.055Google Scholar
Renardy, Y. 2007 The effects of confinement and inertia on the production of droplets. Rheol. Acta 46 (4), 521529.10.1007/s00397-006-0150-yGoogle Scholar
Renardy, Y. Y. & Cristini, V. 2001 Effect of inertia on drop breakup under shear. Phys. Fluids 13 (1), 713.10.1063/1.1331321Google Scholar
Riaud, A., Zhao, S., Wang, K., Cheng, Y. & Luo, G. 2014 Lattice-Boltzmann method for the simulation of multiphase mass transfer and reaction of dilute species. Phys. Rev. E 89, 053308.Google Scholar
Saad, Y. 2003 Iterative Methods for Sparse Linear Systems. Society for Industrial and Applied Mathematics.10.1137/1.9780898718003Google Scholar
Sbragaglia, M., Benzi, R., Biferale, L., Succi, S., Sugiyama, K. & Toschi, F. 2007 Generalized lattice Boltzmann method with multirange pseudopotential. Phys. Rev. E 75, 026702.Google Scholar
Seemann, R., Brinkmann, M., Pfohl, T. & Herminghaus, S. 2011 Droplet based microfluidics. Rep. Prog. Phys. 75 (1), 016601.Google Scholar
Shan, X. & Chen, H. 1993 Lattice Boltzmann model for simulating flows with multiple phases and components. Phys. Rev. E 47 (3), 18151820.Google Scholar
Shapira, M. & Haber, S. 1990 Low Reynolds number motion of a droplet in shear flow including wall effects. Intl J. Multiphase Flow 16 (2), 305321.10.1016/0301-9322(90)90061-MGoogle Scholar
Shardt, O., Derksen, J. J. & Mitra, S. K. 2013 Simulations of droplet coalescence in simple shear flow. Langmuir 29 (21), 62016212.10.1021/la304919pGoogle Scholar
Sibillo, V., Pasquariello, G., Simeone, M., Cristini, V. & Guido, S. 2006 Drop deformation in microconfined shear flow. Phys. Rev. Lett. 97 (5), 054502.10.1103/PhysRevLett.97.054502Google Scholar
Sjöblom, J. 2005 Emulsions and Emulsion Stability. Surfactant Science Series, vol. 132. CRC Press.10.1201/9781420028089Google Scholar
van der Sman, R. G. M. & van der Graaf, S. 2008 Emulsion droplet deformation and breakup with lattice Boltzmann model. Comput. Phys. Commun. 178, 492504.10.1016/j.cpc.2007.11.009Google Scholar
Stone, H. A. 1990 A simple derivation of the time-dependent convective–diffusion equation for surfactant transport along a deforming interface. Phys. Fluids A 2 (1), 111112.10.1063/1.857686Google Scholar
Stone, H. A. 1994 Dynamics of drop deformation and breakup in viscous fluids. Annu. Rev. Fluid Mech. 26 (1), 65102.10.1146/annurev.fl.26.010194.000433Google Scholar
Stone, H. A. & Leal, L. G. 1990 The effects of surfactants on drop deformation and breakup. J. Fluid Mech. 220, 161186.10.1017/S0022112090003226Google Scholar
Sui, Y. 2014 Moving towards the cold region or the hot region? Thermocapillary migration of a droplet attached on a horizontal substrate. Phys. Fluids 26, 092102.10.1063/1.4894077Google Scholar
Swift, M. R., Osborn, W R. & Yeomans, J. M. 1995 Lattice Boltzmann simulation of nonideal fluids. Phys. Rev. Lett. 75 (5), 830.10.1103/PhysRevLett.75.830Google Scholar
Taylor, G. I. 1934 The formation of emulsions in definable fields of flow. Proc. R. Soc. Lond. A 146 (858), 501523.Google Scholar
Teigen, K. E., Li, X., Lowengrub, J., Wang, F. & Voigt, A. 2009 A diffuse-interface approach for modeling transport, diffusion and adsorption/desorption of material quantities on a deformable interface. Commun. Math. Sci. 4 (7), 1009.Google Scholar
Teigen, K. E., Song, P., Lowengrub, J. & Voigt, A. 2011 A diffuse-interface method for two-phase flows with soluble surfactants. J. Comput. Phys. 230 (2), 375393.10.1016/j.jcp.2010.09.020Google Scholar
Tryggvason, G., Bunner, B., Esmaeeli, A., Juric, D., Al-Rawahi, N., Tauber, W., Han, J., Nas, S. & Jan, Y.-J. 2001 A front-tracking method for the computations of multiphase flow. J. Comput. Phys. 169 (2), 708759.10.1006/jcph.2001.6726Google Scholar
Vananroye, A., Van Puyvelde, P. & Moldenaers, P. 2007 Effect of confinement on the steady-state behavior of single droplets during shear flow. J. Rheol. 51 (1), 139153.10.1122/1.2399089Google Scholar
Van der Graaf, S., Steegmans, M. L. J., Van der Sman, R. G. M., Schroën, C. G. P. H. & Boom, R. M. 2005 Droplet formation in a T-shaped microchannel junction: a model system for membrane emulsification. Colloids Surf. A 266 (1), 106116.10.1016/j.colsurfa.2005.06.019Google Scholar
Wang, Y., Shu, C., Huang, H. B. & Teo, C. J. 2015 Multiphase lattice Boltzmann flux solver for incompressible multiphase flows with large density ratio. J. Comput. Phys. 280, 404423.10.1016/j.jcp.2014.09.035Google Scholar
Xu, J.-J., Li, Z., Lowengrub, J. & Zhao, H. 2006 A level-set method for interfacial flows with surfactant. J. Comput. Phys. 212 (2), 590616.10.1016/j.jcp.2005.07.016Google Scholar
Xu, J.-J. & Ren, W. 2014 A level-set method for two-phase flows with moving contact line and insoluble surfactant. J. Comput. Phys. 263, 7190.10.1016/j.jcp.2014.01.012Google Scholar
Xu, J.-J., Yang, Y. & Lowengrub, J. 2012 A level-set continuum method for two-phase flows with insoluble surfactant. J. Comput. Phys. 231 (17), 58975909.10.1016/j.jcp.2012.05.014Google Scholar
Xu, J.-J. & Zhao, H.-K. 2003 An Eulerian formulation for solving partial differential equations along a moving interface. J. Sci. Comput. 19 (1–3), 573594.10.1023/A:1025336916176Google Scholar
Zhang, J., Eckmann, D. M. & Ayyaswamy, P. S. 2006 A front tracking method for a deformable intravascular bubble in a tube with soluble surfactant transport. J. Comput. Phys. 214 (1), 366396.10.1016/j.jcp.2005.09.016Google Scholar
Zhang, T. & Wang, Q. 2010 Cahn–Hilliard vs singular Cahn–Hilliard equations in phase field modeling. Commun. Comput. Phys. 7 (2), 362382.Google Scholar
Figure 0

Figure 1. Illustration of a droplet in a 3D extensional flow.

Figure 1

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

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}$.

Figure 3

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 (1990) and Feigl et al. (2007).

Figure 4

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 5

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 6

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

Figure 7

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.

Figure 8

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

Figure 9

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

Figure 10

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.

Figure 11

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

Figure 12

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.

Figure 13

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 }$.

Figure 14

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.

Figure 15

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

Figure 17. Illustration of the collision of two droplets in a 2D shear flow. It should be noted that the point ‘A’ is the left intersection point of the right droplet and the horizontal centre line, and the points ‘B’ and ‘C’ are the near tip and the far tip of the right droplet.

Figure 17

Figure 18. Collision of two droplets in simple shear flow for clean (the dashed lines) and surfactant-laden cases at $Re=0.4$ and $Ca_{e}=0.1$. In the surfactant-laden case, $Pe=10$, $x_{in}=0.3$ and $E_{0}=0.4$, and the interface is coloured by the dimensionless surfactant concentration $\unicode[STIX]{x1D713}^{\ast }$.

Figure 18

Figure 19. 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_{Ma}$ along the arclength $s$ for $Re=0.4$, $Pe=10$, $Ca_{e}=0.1$, $E_{0}=0.4$, $x_{in}=0.3$ at dimensionless times of $\unicode[STIX]{x1D70F}=0.58$, 2.3, 4.6 and 5.8. It should be noted that the near tip represents the one first encountered as $s$ increases.

Figure 19

Figure 20. The velocity vectors and the streamlines in the vicinity of the gap at $\unicode[STIX]{x1D70F}=4.6$ for $Re=0.4$, $Pe=10$, $x_{in}=0.3$, $E_{0}=0.4$ and $Ca_{e}=0.1$: (a) clean droplets; (b) surfactant-laden droplets.