Modelling double emulsion formation in planar flow-focusing microchannels

Abstract

Double emulsion formation in a hierarchical flow-focusing channel is systematically investigated using a free-energy ternary lattice Boltzmann model. A three-dimensional formation regime diagram is constructed based on the capillary numbers of the inner ($Ca_{i}$), middle ($Ca_{m}$) and outer ($Ca_{o}$) phase fluids. The results show that the formation diagram can be classified into periodic two-step region, periodic one-step region, and non-periodic region. By varying $Ca_{i}$ and $Ca_{m}$ in the two-step formation region, different morphologies are obtained, including the regular double emulsions, decussate regimes with one or two alternate empty droplets, and structures with multiple inner droplets contained in the continuous middle phase thread. Bidisperse behaviours are also frequently encountered in the two-step formation region. In the periodic one-step formation region, scaling laws are proposed for the double emulsion size and for the size ratio between the inner droplet and the overall double emulsion. Furthermore, we show that the interfacial tension ratio can greatly change the morphologies of the obtained emulsion droplets, and the channel geometry plays an important role in changing the formation regimes and the double emulsion sizes. In particular, narrowing the side inlets or the distance between the two side inlets promotes the conversion from the two-step formation regime to the one-step formation regime.

1 Introduction

Double emulsions are droplets with one other droplet inside. Their core–shell structure has attracted wide attention in various fields (Vladisavljević, Nuumani & Nabavi 2017). In pharmaceuticals, one common technique is to use double emulsion for drug encapsulation of highly hydrophilic molecules. It solves the low encapsulation efficiency problem faced in single emulsion techniques due to the quick partitioning of the hydrophilic molecules into the external aqueous phase (Iqbal et al. 2015). Double emulsions are also suitable containers for performing chemical and biochemical reactions under well-defined conditions. Compared to bulk reactions, the greatly reduced volume needed in double emulsion techniques is beneficial for high throughput screening assays (Chang et al. 2018). Furthermore, double emulsions can be used as templates for the synthesis of more complex colloidosomes (Lee & Weitz 2008; Azarmanesh, Farhadi & Azizian 2016), as well as for controlled release of the inner contents (Datta et al. 2014). To ensure the successful applications of double emulsions, one of the key issues is to provide precise control over the double emulsion structure, size and monodispersity at a sufficient production rate (Zhang et al. 2016; Shang, Cheng & Zhao 2017).

Traditional double emulsion fabrication methods, such as the bulk emulsification and the membrane emulsification methods (Vladisavljević et al. 2017), are attractive to many industries (e.g. food and cosmetic) where scalability for large production is important (Varka et al. 2012). However, these techniques have poor size and monodispersity control (Silva, Rodríguez-Abreu & Vilanova 2016), which makes them inadequate for applications requiring high precision, such as in medical, pharmaceutical and material industries. The emergence of microfluidic technology (Utada et al. 2005; Whitesides 2006) opens up a new horizon. It provides more detailed control over the operating conditions (Vladisavljević et al. 2017) and offers great flexibility in designing multilayer (Abate & Weitz 2009) or multicore emulsions (Nisisako, Okushima & Torii 2005; Nabavi, Vladisavljević & Manović 2017b). So far, the microfluidic devices for producing double emulsions can be roughly classified into a series of two T-junctions (Okushima et al. 2004), two flow-focusing junctions (Seo et al. 2007; Pannacci et al. 2008; Abate, Thiele & Weitz 2011), coaxial capillaries (Utada et al. 2005; Nabavi et al. 2017b) and the possible combinations and variations of the aforementioned geometries (Nisisako & Hatsuzawa 2016; Zhu & Wang 2016).

The understanding of double emulsion formation dynamics are crucial for microfluidic control and equipment optimization. Double emulsions are commonly generated either in a two-step or one-step formation regime, depending on whether the inner part of the double emulsion is sheared simultaneously with the middle layer in the outer fluid (Utada et al. 2005). Due to the distinct flow details in the two-step and one-step formation regimes, the influence of flow conditions, fluid properties and geometrical parameters on each regime should be analysed respectively. For the two-step formation regime, Okushima et al. (2004) have systematically showed the effect of flow rates on the double emulsion sizes and the number of internal droplet for multicore emulsions when they are formed using a series of T-junctions. The one-step formation regime is mostly encountered in coaxial microcapillary devices. Experimental studies have been carried out on the effect of flow rates (Lee & Weitz 2008; Kim et al. 2013) and geometrical settings (Nabavi et al. 2017a). Scaling laws have also been developed for the emulsion size predictions (Utada et al. 2005; Chang et al. 2009).

Complementary to experiments, numerical studies on double emulsion formation dynamics in microfluidic channels have also garnered strong interest. For instance, great efforts were made to elucidate the effects of flow rates, interfacial tension ratios, geometry (Chen, Wu & Lin 2015; Nabavi et al. 2015b) and viscosities (Fu et al. 2016) on the double emulsion properties and the flow regime predictions (Nabavi et al. 2017a) for coaxial flow-focusing capillary devices. Simulations are particularly advantageous for providing accurate flow details and for allowing each relevant parameter in the system to be varied systematically. In the literature, a number of ternary fluid models have been successfully developed and applied in the study of multiple emulsion formation behaviours, including using the volume of fluid method (Chen et al. 2015; Nabavi et al. 2015b, 2017a; Azarmanesh et al. 2016), the front-tracking method (Vu et al. 2013), the free energy finite element method (Park & Anderson 2012) and the lattice Boltzmann method (Fu et al. 2016, 2017).

In this work, our focus is on the planar flow-focusing cross-junctions. They are promising for integration with other devices and they can be parallelized to raise the production rate of the emulsion droplets, while still ensuring good size control (Lee et al. 2016). Furthermore, in contrast to other microfluidic geometries, systematic parametric study is rarely reported on planar flow-focusing devices. Several works, such as Abate et al. (2011) and Azarmanesh et al. (2016), briefly discussed the possible conversion between the two-step and one-step formation regimes and the variation of shell thickness. However, it remains unclear in which flow rate regions monodisperse double emulsions are produced, and correspondingly, how the droplet sizes can be varied in those regions. It is likely that the droplet sizes have different dependencies on the flow rates for the two-step and one-step formation processes. There are also open questions on the role of channel geometry in the formation regime conversion, and on the effects of interfacial tension ratio in determining the morphologies of the emulsion droplets, including the possibility of complete, partial and non-engulfing shapes (Pannacci et al. 2008; Nisisako & Hatsuzawa 2010; Guzowski et al. 2012; Chao, Mak & Shum 2016).

We have chosen to use the lattice Boltzmann method (LBM). The lattice Boltzmann method is highly favourable for the study of emulsion formation behaviours due to its simplicity in solving interface dynamics, including droplet breakups and coalescences, as well as its ability to deal with complex geometries, and its high efficiency in parallel computation (Krüger et al. 2017). So far, three types of ternary lattice Boltzmann models have been developed, including the free-energy model (Liang, Shi & Chai 2016; Semprebon, Krüger & Kusumaatmaja 2016; Abadi, Fakhari & Rahimian 2018; Wöhrwag et al. 2018), colour-fluid model (Leclaire et al. 2013a; Leclaire, Reggio & Trépanier 2013b; Fu et al. 2016, 2017; Yu et al. 2019b), and the Shan–Chen type models (Bao & Schaefer 2013; Wei et al. 2018).

Here we improve on the free-energy lattice Boltzmann model developed by Semprebon et al. (2016). A major progress is that our model allows a wider range of interfacial tension ratios, such that all possible biphasic emulsion morphologies can be captured (Guzowski et al. 2012), including complete and non-engulfing shapes. The model developed by Semprebon et al. (2016) only allows partial engulfing shapes. Coupling the free-energy model with the advantages of the lattice Boltzmann method, we conduct a systematic study on the dynamics of double emulsion formation behaviours in planar hierarchical flow-focusing junctions. We focus on the two-dimensional (2-D) case to reduce the computational time needed for parametric studies. The major physical difference in the flow dynamics between the 2-D and the three-dimensional (3-D) systems lies in the lack of an additional Laplace pressure induced by the out-of-plane curvature (Chen & Deng 2017). Such contribution can accelerate the droplet pinch-off process (Hoang et al. 2013), especially at large Weber number. However, we believe most of the fundamental flow physics are still involved in the 2-D system and a systematic 2-D study can still capture some of the key qualitative trends in the formation regimes and emulsion sizes as function of the flow rates of each fluid phase.

The paper is organized as follows. In § 2, we describe the improved ternary free-energy model, the lattice Boltzmann method and the boundary conditions involved. In § 3, we validate the model and boundary conditions by Poiseuille flow, moving droplet and static emulsion morphology tests. In § 4, our systematic parametric study allows us to construct a flow regime diagram, where we describe a wide range of formation regimes, including the periodic two-step and one-step double emulsion formation regimes, decussate regime, bidisperse regime and even the continuous structure with multiple inner droplets. Scaling laws are also proposed for the double emulsions produced in the one-step formation regime, and the effects of the interfacial tension ratios and the geometrical parameters are investigated. Finally, we summarize our main findings and forecast prospects for future work in § 5.

2 Numerical method

2.1 Free-energy model

The present model is developed based on the equal-density ternary free-energy lattice Boltzmann model proposed by Semprebon et al. (2016). The system is described by the free-energy functional

(2.1)$$\begin{eqnarray}{\mathcal{H}}=\int _{\unicode[STIX]{x1D6FA}}c_{s}^{2}\unicode[STIX]{x1D70C}\ln \unicode[STIX]{x1D70C}\,\text{d}V+{\mathcal{F}}.\end{eqnarray}$$

The first term is the standard ideal gas term in the lattice Boltzmann method with $c_{s}=1/\sqrt{3}$ and $\unicode[STIX]{x1D70C}$ the total density. Here, $\unicode[STIX]{x1D6FA}$ is the system volume. The first term does not affect the phase behaviour. To realise three fluid components, the second term ${\mathcal{F}}$ is introduced and it is given by

(2.2)$$\begin{eqnarray}\displaystyle {\mathcal{F}} & = & \displaystyle \mathop{\sum }_{m=1}^{3}\int _{\unicode[STIX]{x1D6FA}}(E_{0}(C_{m})+E_{interface}(\unicode[STIX]{x1D735}C_{m}))\,\text{d}V\nonumber\\ \displaystyle & = & \displaystyle \mathop{\sum }_{m=1}^{3}\int _{\unicode[STIX]{x1D6FA}}\left[\frac{\unicode[STIX]{x1D705}_{m}}{2}C_{m}^{2}(1-C_{m})^{2}+\frac{\unicode[STIX]{x1D6FC}^{2}\unicode[STIX]{x1D705}_{m}}{2}(\unicode[STIX]{x1D735}C_{m})^{2}\right]\,\text{d}V.\end{eqnarray}$$

It is constructed using a double-well potential form for the bulk free energy $E_{0}(C_{m})=(\unicode[STIX]{x1D705}_{m}/2)C_{m}^{2}(1-C_{m})^{2}$ and a square gradient term for the interface region $E_{interface}(\unicode[STIX]{x1D735}C_{m})=(\unicode[STIX]{x1D6FC}^{2}\unicode[STIX]{x1D705}_{m}/2)(\unicode[STIX]{x1D735}C_{m})^{2}$. Here, $C_{m}$ ($m=1,2,3$) are the concentration fractions with two minimum values at $C_{m}=0$ and $C_{m}=1$ for each component $m$. In the current model, all components have the same density $\unicode[STIX]{x1D70C}_{m}$, which we have set to be 1.0 for simplicity. Thus the total density is related to the concentration fractions by defining $\unicode[STIX]{x1D70C}=C_{1}\unicode[STIX]{x1D70C}_{1}+C_{2}\unicode[STIX]{x1D70C}_{2}+C_{3}\unicode[STIX]{x1D70C}_{3}=C_{1}+C_{2}+C_{3}$, which is equal to 1.0 in this model. Three physically meaningful bulk phases termed red, green and blue could be denoted by $\unicode[STIX]{x1D63E}=[C_{1}\;\;C_{2}\;\;C_{3}]=[1\;\;0\;\;0]$, $[0\;\;1\;\;0]$ and $[0\;\;0\;\;1]$, respectively. Here, $\unicode[STIX]{x1D6FC}$ is the interface width parameter. The interfacial tension between red–green phases $\unicode[STIX]{x1D70E}_{rg}$, red–blue phases $\unicode[STIX]{x1D70E}_{rb}$ and green–blue phases $\unicode[STIX]{x1D70E}_{gb}$ are related to different $\unicode[STIX]{x1D705}$ through

(2.3)$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x1D70E}_{rg}=\frac{\unicode[STIX]{x1D6FC}}{6}(\unicode[STIX]{x1D705}_{1}+\unicode[STIX]{x1D705}_{2}),\\ \displaystyle \unicode[STIX]{x1D70E}_{rb}=\frac{\unicode[STIX]{x1D6FC}}{6}(\unicode[STIX]{x1D705}_{1}+\unicode[STIX]{x1D705}_{3}),\\ \displaystyle \unicode[STIX]{x1D70E}_{gb}=\frac{\unicode[STIX]{x1D6FC}}{6}(\unicode[STIX]{x1D705}_{2}+\unicode[STIX]{x1D705}_{3}).\end{array}\right\}\end{eqnarray}$$

To capture the interface dynamics, two order parameters $\unicode[STIX]{x1D719}$ and $\unicode[STIX]{x1D713}$ are introduced as

(2.4a,b)$$\begin{eqnarray}\unicode[STIX]{x1D719}=C_{1}-C_{2},\quad \unicode[STIX]{x1D713}=C_{3},\end{eqnarray}$$

and the original concentration fields can be reversely obtained from $\unicode[STIX]{x1D70C}$, $\unicode[STIX]{x1D719}$ and $\unicode[STIX]{x1D713}$ via $C_{1}=(\unicode[STIX]{x1D70C}+\unicode[STIX]{x1D719}-\unicode[STIX]{x1D713})/2$, $C_{2}=(\unicode[STIX]{x1D70C}-\unicode[STIX]{x1D719}-\unicode[STIX]{x1D713})/2$ and $C_{3}=\unicode[STIX]{x1D713}$. The order parameters and the hydrodynamics of the ternary fluid system are governed by two Cahn–Hilliard equations, the continuity and the Navier–Stokes equations

(2.5)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D719}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D719}\boldsymbol{v})=M_{\unicode[STIX]{x1D719}}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D707}_{\unicode[STIX]{x1D719}}, & \displaystyle\end{eqnarray}$$

(2.6)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D713}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D713}\boldsymbol{v})=M_{\unicode[STIX]{x1D713}}\unicode[STIX]{x1D6FB}^{2}\unicode[STIX]{x1D707}_{\unicode[STIX]{x1D713}}, & \displaystyle\end{eqnarray}$$

(2.7)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D70C}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}\boldsymbol{v})=0, & \displaystyle\end{eqnarray}$$

(2.8)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{t}(\unicode[STIX]{x1D70C}\boldsymbol{v})+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}\boldsymbol{v}\boldsymbol{v})=\unicode[STIX]{x1D735}\boldsymbol{\cdot }[\unicode[STIX]{x1D702}(\unicode[STIX]{x1D735}\boldsymbol{v}+\unicode[STIX]{x1D735}\boldsymbol{v}^{\text{T}})]-\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D64B}. & \displaystyle\end{eqnarray}$$

Here, $\boldsymbol{v}$ is the fluid velocity and $\unicode[STIX]{x1D702}$ is the dynamic viscosity. The mobility values for $\unicode[STIX]{x1D719}$ and $\unicode[STIX]{x1D713}$ are $M_{\unicode[STIX]{x1D719}}$ and $M_{\unicode[STIX]{x1D713}}$. The thermodynamic properties are related to the equations of motion via the chemical potential $\unicode[STIX]{x1D707}_{\unicode[STIX]{x1D719}}$, $\unicode[STIX]{x1D707}_{\unicode[STIX]{x1D713}}$ and the pressure tensor $\unicode[STIX]{x1D64B}$. The chemical potential is defined as the variational derivative of the free energy $\unicode[STIX]{x1D707}_{q}=\unicode[STIX]{x1D6FF}{\mathcal{F}}/\unicode[STIX]{x1D6FF}q$, where $q=\unicode[STIX]{x1D70C},\unicode[STIX]{x1D719}$ or $\unicode[STIX]{x1D713}$. The pressure tensor term in (2.8) is constructed by $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D64B}=\unicode[STIX]{x1D735}(\unicode[STIX]{x1D70C}c_{s}^{2})+\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D64B}_{\unicode[STIX]{x1D7EC}}$ and $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D64B}_{\unicode[STIX]{x1D7EC}}=(\unicode[STIX]{x1D70C}\unicode[STIX]{x1D735}\unicode[STIX]{x1D707}_{\unicode[STIX]{x1D70C}}+\unicode[STIX]{x1D719}\unicode[STIX]{x1D735}\unicode[STIX]{x1D707}_{\unicode[STIX]{x1D719}}+\unicode[STIX]{x1D713}\unicode[STIX]{x1D735}\unicode[STIX]{x1D707}_{\unicode[STIX]{x1D713}})$. The first term $\unicode[STIX]{x1D735}(\unicode[STIX]{x1D70C}c_{s}^{2})$ is the standard ideal gas term in LBM and it is simply incorporated in the equilibrium distribution function (Briant & Yeomans 2004; Zhang 2011). The $\boldsymbol{F}=-\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D64B}_{\unicode[STIX]{x1D7EC}}$ term is treated as an external force term in the lattice Boltzmann implementation. The explicit expressions of $\unicode[STIX]{x1D707}_{\unicode[STIX]{x1D70C}}$, $\unicode[STIX]{x1D707}_{\unicode[STIX]{x1D719}}$, $\unicode[STIX]{x1D707}_{\unicode[STIX]{x1D713}}$ and $\unicode[STIX]{x1D64B}$ are given in Semprebon et al. (2016).

Consider now a case where a red droplet is completely engulfed by a green one and they are submerged in a blue phase fluid at thermodynamic equilibrium. According to the theoretical analysis of Guzowski et al. (2012), the interfacial tensions should satisfy $\unicode[STIX]{x1D70E}_{rb}>\unicode[STIX]{x1D70E}_{gb}+\unicode[STIX]{x1D70E}_{rg}$. Given the relation between the $\unicode[STIX]{x1D70E}$ and $\unicode[STIX]{x1D705}_{m}$ values in (2.3), one can easily find that $\unicode[STIX]{x1D705}_{2}$ should be negative while $\unicode[STIX]{x1D705}_{1}$ and $\unicode[STIX]{x1D705}_{3}$ are positive. In the free-energy model, the negative $\unicode[STIX]{x1D705}_{2}$ will invert the shape of the bulk free-energy profile $E_{0}(C_{m})$: the two minimum values at 0 and 1.0 become two maximum values as shown by the blue solid line in figure 1. As such, setting one of the $\unicode[STIX]{x1D705}_{m}$ to be negative often leads to incorrect results or even numerical instability as the concentration value deviates significantly from [0, 1.0]. A similar situation has been encountered in other lattice Boltzmann models, and a simple remedy has been proposed by introducing an additional free-energy term (see Lee & Liu (2010) and Abadi et al. (2018)). Inspired by these works, to solve the problem induced by negative $\unicode[STIX]{x1D705}_{m}$, here we introduce an additional free-energy term given by

(2.9)$$\begin{eqnarray}E_{a}(C_{m})=\left\{\begin{array}{@{}ll@{}}\displaystyle \unicode[STIX]{x1D6FD}C_{m}^{2}, & C_{m}<0,\\ \displaystyle 0, & 0\leqslant C_{m}\leqslant 1,\\ \displaystyle \unicode[STIX]{x1D6FD}(C_{m}-1)^{2}, & C_{m}>1,\end{array}\right.\end{eqnarray}$$

where $\unicode[STIX]{x1D6FD}$ is an adjustable positive parameter controlling the slope of the energy profile $E_{0}+E_{a}$ as depicted by the red dashed line in figure 1. Since we add a new free-energy term in (2.9), additional terms should be included in the chemical potentials accordingly, which are listed in appendix A.

Figure 1. Illustration of the bulk free-energy profile without (blue solid line) and with (red dashed line) the additional free-energy term $E_{a}$, given by (2.9), for a negative $\unicode[STIX]{x1D705}_{m}$.

Figure 2. (a) Illustration of the moving droplet set-up; (b1–b3) time histories of $X_{d}$ for $u_{max}=1.5\times 10^{-3}$, $3.0\times 10^{-4}$ and $7.5\times 10^{-5}$ with different convective boundary conditions; (c1–c3) typical snapshots of the droplet shape and velocity vectors when the droplet passes through the outlet layer with each CBC boundary at $u_{max}$ corresponding to (b1–b3), respectively. The $X_{d}$ and time values are normalized using $X_{d}^{\ast }=X_{d}/D$ and $t^{\ast }=tu_{max}/D$.

2.2 Lattice Boltzmann method

To solve (2.5)–(2.8) using the lattice Boltzmann method, three distribution functions are introduced: $f_{i}(\boldsymbol{r},t)$ for the density, and $g_{i}(\boldsymbol{r},t)$ and $h_{i}(\boldsymbol{r},t)$ for the order parameters $\unicode[STIX]{x1D719}$ and $\unicode[STIX]{x1D713}$, respectively. The distribution functions are discretized in space $\boldsymbol{r}$ and time $t$, according to a set of lattice velocity vectors $\boldsymbol{e}_{i}$. In the D2Q9 discrete scheme (two-dimension with nine discrete velocities) used here, the lattice velocities are given as $\boldsymbol{e}_{0}=(0,0)$, $\boldsymbol{e}_{1,3}=(\pm 1,0)$, $\boldsymbol{e}_{2,4}=(0,\pm 1)$, $\boldsymbol{e}_{5,7}=(\pm 1,\pm 1)$ and $\boldsymbol{e}_{6,8}=(\mp 1,\pm 1)$, as shown in figure 2(a). The time evolution of the distribution functions includes the collision and streaming steps, which can be written as

(2.10)$$\begin{eqnarray}\displaystyle f_{i}(\boldsymbol{r}+\boldsymbol{e}_{\boldsymbol{i}}\unicode[STIX]{x1D739}_{\boldsymbol{t}},t+\unicode[STIX]{x1D6FF}_{t}) & = & \displaystyle f_{i}(\boldsymbol{r},t)-\frac{1}{\unicode[STIX]{x1D70F}_{f}}[f_{i}(\boldsymbol{r},t)-f_{i}^{eq}(\unicode[STIX]{x1D70C},\widetilde{\boldsymbol{u}})]\nonumber\\ \displaystyle & & \displaystyle +\,[f_{i}^{eq}(\unicode[STIX]{x1D70C},\widetilde{\boldsymbol{u}}+\unicode[STIX]{x1D6FF}\widetilde{\boldsymbol{u}})-f_{i}^{eq}(\unicode[STIX]{x1D70C},\widetilde{\boldsymbol{u}})],\end{eqnarray}$$

(2.11)$$\begin{eqnarray}\displaystyle & \displaystyle g_{i}(\boldsymbol{r}+\boldsymbol{e}_{\boldsymbol{i}}\unicode[STIX]{x1D739}_{\boldsymbol{t}},t+\unicode[STIX]{x1D6FF}_{t})=g_{i}(\boldsymbol{r},t)-\frac{1}{\unicode[STIX]{x1D70F}_{g}}[g_{i}(\boldsymbol{r},t)-g_{i}^{eq}(\boldsymbol{r},t)], & \displaystyle\end{eqnarray}$$

(2.12)$$\begin{eqnarray}\displaystyle & \displaystyle h_{i}(\boldsymbol{r}+\boldsymbol{e}_{\boldsymbol{i}}\unicode[STIX]{x1D739}_{\boldsymbol{t}},t+\unicode[STIX]{x1D6FF}_{t})=h_{i}(\boldsymbol{r},t)-\frac{1}{\unicode[STIX]{x1D70F}_{h}}[h_{i}(\boldsymbol{r},t)-h_{i}^{eq}(\boldsymbol{r},t)]. & \displaystyle\end{eqnarray}$$

Here, the force term is implemented through the exact difference method (Kupershtokh, Medvedev & Karpov 2009; Mazloomi, Chikatamarla & Karlin 2015), which is expressed as the last two terms enclosed in brackets in (2.10), with $\widetilde{\boldsymbol{u}}=\sum _{i}f_{i}\boldsymbol{e}_{i}/\unicode[STIX]{x1D70C}$, i.e. the velocity without the force term, and $\unicode[STIX]{x1D6FF}\widetilde{\boldsymbol{u}}=\boldsymbol{F}\unicode[STIX]{x1D6FF}t/\unicode[STIX]{x1D70C}$. The lattice time step $\unicode[STIX]{x1D6FF}_{t}$ is set to be 1.0. Here, $\unicode[STIX]{x1D70F}_{f}$ is the relaxation parameter given by $1/\unicode[STIX]{x1D70F}_{f}=C_{1}/\unicode[STIX]{x1D70F}_{1}+C_{2}/\unicode[STIX]{x1D70F}_{2}+C_{3}/\unicode[STIX]{x1D70F}_{3}$, where $\unicode[STIX]{x1D70F}_{1,2,3}$ are related to the viscosity of each fluid by $\unicode[STIX]{x1D70F}_{1,2,3}=3\unicode[STIX]{x1D702}_{r,g,b}/\unicode[STIX]{x1D70C}+1/2$, respectively (Krüger et al. 2017); $\unicode[STIX]{x1D70F}_{g}$ and $\unicode[STIX]{x1D70F}_{h}$ are the relaxation parameters that are related to the mobility parameters $M_{\unicode[STIX]{x1D719}}$ and $M_{\unicode[STIX]{x1D713}}$ in the Cahn–Hilliard equations through

(2.13)$$\begin{eqnarray}\displaystyle & \displaystyle M_{\unicode[STIX]{x1D719}}=\unicode[STIX]{x1D6E4}_{\unicode[STIX]{x1D719}}\left(\unicode[STIX]{x1D70F}_{g}-\frac{\unicode[STIX]{x1D6FF}t}{2}\right), & \displaystyle\end{eqnarray}$$

(2.14)$$\begin{eqnarray}\displaystyle & \displaystyle M_{\unicode[STIX]{x1D713}}=\unicode[STIX]{x1D6E4}_{\unicode[STIX]{x1D713}}\left(\unicode[STIX]{x1D70F}_{h}-\frac{\unicode[STIX]{x1D6FF}t}{2}\right), & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6E4}_{\unicode[STIX]{x1D719}}$ and $\unicode[STIX]{x1D6E4}_{\unicode[STIX]{x1D713}}$ are two tunable parameters. The mobility values $M_{\unicode[STIX]{x1D719}}$ and $M_{\unicode[STIX]{x1D713}}$ are relevant for the time scale of Cahn–Hilliard diffusion and the relaxation time of the interface. Generally, the mobility values should be sufficiently large to retain the interfacial thickness, but small enough to ensure the reasonable damping of the convective term (Jacqmin 1999; Lim & Lam 2013). At present it remains an open problem to assign mobility values in numerical studies. Indeed most papers use comparison with experiments to set the values, and one common solution is to use mobility related dimensionless parameters, e.g. the Peclet number ($Pe$). In our microfluidic study, a characteristic $Pe$ is defined based on the middle phase fluid as $Pe_{m}=(u_{m}w_{1})/(M_{\unicode[STIX]{x1D719}}\unicode[STIX]{x1D705}_{2})$. The absolute values of $Pe_{m}$ used is generally of the order of $O(10){-}O(80)$, which is of similar magnitude to those used in previous two-phase droplet behaviour studies (Menech 2006; Zhou et al. 2010; Shardt, Mitra & Derksen 2014). Moreover, $M_{\unicode[STIX]{x1D713}}=M_{\unicode[STIX]{x1D719}}/3$ is considered to assign symmetric mobility for each concentration component (Semprebon et al. 2016).

Here, $f_{i}^{eq}$, $g_{i}^{eq}$ and $h_{i}^{eq}$ are the local equilibrium distribution functions, which are given by

(2.15)$$\begin{eqnarray}\displaystyle & \displaystyle f_{i}^{eq}=\unicode[STIX]{x1D714}_{i}\unicode[STIX]{x1D70C}[1+3\boldsymbol{e}_{i}\boldsymbol{\cdot }\widetilde{\boldsymbol{u}}+{\textstyle \frac{9}{2}}(\boldsymbol{e}_{i}\boldsymbol{\cdot }\widetilde{\boldsymbol{u}})^{2}-{\textstyle \frac{3}{2}}\widetilde{\boldsymbol{u}}^{2}], & \displaystyle\end{eqnarray}$$

(2.16)$$\begin{eqnarray}\displaystyle & \displaystyle g_{i}^{eq}=\left\{\begin{array}{@{}ll@{}}\displaystyle \unicode[STIX]{x1D714}_{i}\left[3\unicode[STIX]{x1D6E4}_{\unicode[STIX]{x1D719}}\unicode[STIX]{x1D707}_{\unicode[STIX]{x1D719}}+3\unicode[STIX]{x1D719}\boldsymbol{e}_{i}\boldsymbol{\cdot }\boldsymbol{v}+\frac{9\unicode[STIX]{x1D719}}{2}(\boldsymbol{e}_{i}\boldsymbol{\cdot }\boldsymbol{v})^{2}-\frac{3\unicode[STIX]{x1D719}}{2}\boldsymbol{v}^{2}\right], & i=1{-}8,\\ \displaystyle \unicode[STIX]{x1D719}-\mathop{\sum }_{i=1}^{8}g_{i}^{eq}, & i=0,\end{array}\right. & \displaystyle\end{eqnarray}$$

(2.17)$$\begin{eqnarray}\displaystyle & \displaystyle h_{i}^{eq}=\left\{\begin{array}{@{}ll@{}}\displaystyle \unicode[STIX]{x1D714}_{i}\left[3\unicode[STIX]{x1D6E4}_{\unicode[STIX]{x1D713}}\unicode[STIX]{x1D707}_{\unicode[STIX]{x1D713}}+3\unicode[STIX]{x1D713}\boldsymbol{e}_{i}\boldsymbol{\cdot }\boldsymbol{v}+\frac{9\unicode[STIX]{x1D713}}{2}(\boldsymbol{e}_{i}\boldsymbol{\cdot }\boldsymbol{v})^{2}-\frac{3\unicode[STIX]{x1D713}}{2}\boldsymbol{v}^{2}\right], & i=1{-}8,\\ \displaystyle \unicode[STIX]{x1D713}-\mathop{\sum }_{i=1}^{8}h_{i}^{eq}, & i=0,\end{array}\right. & \displaystyle\end{eqnarray}$$

where the weight coefficients $\unicode[STIX]{x1D714}_{i}$ are given by $\unicode[STIX]{x1D714}_{0}=4/9$, $\unicode[STIX]{x1D714}_{1-4}=1/9$ and $\unicode[STIX]{x1D714}_{5-8}=1/36$. The macroscopic variables are related to the distribution functions through

(2.18a-d)$$\begin{eqnarray}\unicode[STIX]{x1D70C}=\mathop{\sum }_{i}f_{i},\quad \unicode[STIX]{x1D70C}\boldsymbol{v}=\mathop{\sum }_{i}f_{i}\boldsymbol{e}_{i}+\frac{\boldsymbol{F}\unicode[STIX]{x1D6FF}t}{2},\quad \unicode[STIX]{x1D719}=\mathop{\sum }_{i}g_{i},\quad \unicode[STIX]{x1D713}=\mathop{\sum }_{i}h_{i}.\end{eqnarray}$$

2.3 Boundary conditions

The boundary conditions involved in the present study contain the no-slip boundary, the wetting boundary and the inlet–outlet boundary. The no-slip boundary condition is used on the solid walls, which is realized by the halfway bounceback rule (Ladd 1994). The solid walls should have a preferential affinity with the continuous phase fluid to generate stable droplets/emulsions (Abate et al. 2011). Fu et al. (2016) successfully implemented the wetting boundary condition by setting a fictive density on the walls in a lattice Boltzmann ternary colour-fluid model. Similarly for the free-energy model used here, the macroscopic values of $\unicode[STIX]{x1D70C}$, $\unicode[STIX]{x1D719}$ and $\unicode[STIX]{x1D713}$ on the walls are designated to be the same as those of the continuous phase fluid that is assumed to completely wet the walls. For the velocity inlet, the Zou–He velocity boundary condition (Zou & He 1997) is applied to solve the unknown density distribution functions of $f_{i}$. To obtain the unknown $g_{i}$ and $h_{i}$ values at the inlet, the method used by Hao & Cheng (2009) and Liu & Zhang (2011) is adopted. Take figure 2(a) for instance, given an inlet boundary with the inflow direction pointing to the right, $g_{1,5,8}$ and $h_{1,5,8}$ are unknown after the streaming step. We assume that one pure single fluid exists at the inlet, where the prescribed values of $\unicode[STIX]{x1D719}$ and $\unicode[STIX]{x1D713}$ are $\unicode[STIX]{x1D719}_{in}$ and $\unicode[STIX]{x1D713}_{in}$, respectively. The sum of the unknown distribution functions can be solved according to (2.18), and then $g_{1,5,8}$ and $h_{1,5,8}$ are allocated by their weight factors as

(2.19a,b)$$\begin{eqnarray}\left.g_{i}\vphantom{\Big|}\right|_{i=1,5,8}=\frac{g^{\ast }\unicode[STIX]{x1D714}_{i}}{\unicode[STIX]{x1D714}_{1}+\unicode[STIX]{x1D714}_{5}+\unicode[STIX]{x1D714}_{8}},\quad g^{\ast }=g_{1}+g_{5}+g_{8}=\unicode[STIX]{x1D719}_{in}-\mathop{\sum }_{i}\left.g_{i}\vphantom{\Big|}\right|_{i=0,2,3,4,6,7},\end{eqnarray}$$

(2.20a,b)$$\begin{eqnarray}\left.h_{i}\vphantom{\Big|}\right|_{i=1,5,8}=\frac{h^{\ast }\unicode[STIX]{x1D714}_{i}}{\unicode[STIX]{x1D714}_{1}+\unicode[STIX]{x1D714}_{5}+\unicode[STIX]{x1D714}_{8}},\quad h^{\ast }=h_{1}+h_{5}+h_{8}=\unicode[STIX]{x1D713}_{in}-\mathop{\sum }_{i}\left.h_{i}\vphantom{\Big|}\right|_{i=0,2,3,4,6,7}.\end{eqnarray}$$

For the outlet boundary, the convective boundary condition (CBC) (Lou, Guo & Shi 2013; Chen & Deng 2017) is used for its good performance in multicomponent flow simulations. In the present model, the CBC is harnessed in two aspects. One is for the unknown distribution functions $\unicode[STIX]{x1D712}_{i}=f_{i}$, $g_{i}$ and $h_{i}$ at the outlet layer ($x=L_{x}$):

(2.21)$$\begin{eqnarray}\unicode[STIX]{x1D712}_{i}(L_{x},y,t+\unicode[STIX]{x1D6FF}t)=\frac{\unicode[STIX]{x1D712}_{i}(L_{x},y,t)+\unicode[STIX]{x1D701}(L_{x}-1,y,t)\unicode[STIX]{x1D712}_{i}(L_{x}-1,y,t+\unicode[STIX]{x1D6FF}t)}{1+\unicode[STIX]{x1D701}(L_{x}-1,y,t)}.\end{eqnarray}$$

The other is for the macroscopic quantities, such as $\unicode[STIX]{x1D712}^{\prime }=\unicode[STIX]{x1D70C}$, $\unicode[STIX]{x1D719}$, $\unicode[STIX]{x1D713}$ and $\unicode[STIX]{x1D64B}_{\unicode[STIX]{x1D7EC}}$ on the ghost layer right outside the outlet, i.e. $x=L_{x}+1$, which is needed to compute the derivative terms at the outlet fluid layer:

(2.22)$$\begin{eqnarray}\unicode[STIX]{x1D712}^{\prime }(L_{x}+1,y,t+\unicode[STIX]{x1D6FF}t)=\frac{\unicode[STIX]{x1D712}^{\prime }(L_{x}+1,y,t)+\unicode[STIX]{x1D701}(L_{x},y,t)\unicode[STIX]{x1D712}^{\prime }(L_{x},y,t+\unicode[STIX]{x1D6FF}t)}{1+\unicode[STIX]{x1D701}(L_{x},y,t)}.\end{eqnarray}$$

Here, $\unicode[STIX]{x1D701}$ is the characteristic velocity normal to the outlet boundary. For simplicity, we have explicitly computed $\unicode[STIX]{x1D712}_{i}$ and $\unicode[STIX]{x1D712}^{\prime }$ through $\unicode[STIX]{x1D701}$ at time $t$. Three common choices for $\unicode[STIX]{x1D701}$ in convective boundary conditions are the average velocity (CBC-AV), local velocity (CBC-LV) and the maximum velocity (CBC-MV) (Lou et al. 2013).

3 Model validation

3.1 Convective outlet boundary conditions

In this section, the performance of the CBC in the present model is tested by simulating a single-phase Poiseuille flow and a Poiseuille flow with a moving droplet. In the single-phase Poiseuille flow settings, a fluid with viscosity of 0.167 flows in the $x$ direction with a maximum velocity of $u_{max}=0.0015$ in a computational domain of $L_{x}\times L_{y}=99\times 39$. No-slip boundaries are used both on the top and bottom walls. The Zou–He velocity inlet is applied with a parabolic velocity distribution given as

(3.1)$$\begin{eqnarray}u_{x}(y)=\frac{-4u_{max}(y-y_{1})(y-y_{2})}{(y_{2}-y_{1})^{2}},\quad 1\leqslant y\leqslant L_{y}-1,\end{eqnarray}$$

where $y_{1}=0.5$ and $y_{2}=L_{y}-0.5$ are the locations of the bottom and top walls. All three options of the CBC mentioned above are implemented at the outlet, and their accuracy is quantified using the relative velocity error computed by $E_{u}=\sqrt{\sum ((u_{x})_{ana}-(u_{x})_{simu})^{2}/\sum ((u_{x})_{ana}^{2})}$, where $(u_{x})_{ana}$ is the analytical velocity given by (3.1) and $(u_{x})_{simu}$ denotes the simulated velocity. The obtained values of $E_{u}$ under CBC-AV, CBC-LV and CBC-MV conditions are $1.449\times 10^{-4}$, $1.151\times 10^{-4}$ and $2.254\times 10^{-4}$ in the middle of the channel, i.e. $x=49$, and $1.454\times 10^{-4}$, $5.4\times 10^{-3}$ and $3.144\times 10^{-4}$ at the outlet layer. It is seen that all three outlet boundaries give satisfactory results for flow far away from the outlet. However, the accuracy at the outlet layer varies: the CBC-AV provides the highest accuracy, CBC-MV is slightly lower and CBC-LV shows the poorest performance.

In the moving droplet test, a droplet with radius $R=20$ is centred at $(60,49.5)$ in a channel of $L_{x}\times L_{y}=199\times 99$, as illustrated in figure 2(a). The two fluid phases have the same viscosity of 0.167 and their interfacial tension $\unicode[STIX]{x1D70E}$ is 0.005. All the boundary conditions are the same as those in the single-phase Poiseuille flow simulations. The whole fluid domain is initialized with a uniform parabolic velocity profile as given by (3.1). Three different values of $u_{max}$ are tested, i.e. $u_{max}=1.5\times 10^{-3}$, $3.0\times 10^{-4}$ and $7.5\times 10^{-5}$. To make a quantitative comparison, the time history of the distance $X_{d}$ measured from the inlet to the leftmost point of the droplet is recorded and shown in figure 2(b1–b3). The $X_{d}$ and time are normalized using $X_{d}^{\ast }=X_{d}/D$ and $t^{\ast }=tu_{max}/D$, where $D$ is the droplet diameter. The $X_{d}^{\ast }$ curve of the droplet moving in a longer channel ($L_{x}\times L_{y}=399\times 99$) computed with CBC-AV is used as the reference result for each flow condition. Note in figure 2(b1–b3) that the sharp decrease of $X_{d}^{\ast }$ occurs when the droplet completely moves out of the channel.

It is seen in figure 2(b1–b3) that the $X_{d}^{\ast }$ increases linearly with time and agrees with the reference line before the droplet interface touches the outlet boundary for each of the tested flow conditions. The option of the convective boundary conditions has little effect on the flow behaviours away from the outlet. Deviations in $X_{d}^{\ast }$ curves from the reference lines occur at around $t^{\ast }=4$ when the droplet passes through the outlet. Compared to the reference lines, the case with CBC-AV slightly lags behind, and the case with CBC-MV moves a bit faster. Also, the case with CBC-LV gives the most accurate results for moderate characteristic velocities, as illustrated in figure 2(b1,b2). The deviation in $X_{d}^{\ast }$ increases as $u_{max}$ decreases for the cases with CBC-AV and CBC-MV. When the $u_{max}$ is of the same order of magnitude as the spurious velocities of the present model, i.e. $u_{max}=7.5\times 10^{-5}$ in (b3), numerical instability arises for the case with CBC-LV, whereas the cases with CBC-AV and CBC-MV show better robustness. Due to the low velocity often encountered in double emulsion generation, the robustness of the outlet boundary at low velocities is of great significance. On the other hand, for low velocity cases shown in (c2,c3), the velocity in regions close to the walls is less affected for the case with CBC-AV than that with CBC-MV. The momentum deficit or surplus around the outlet region could be attributed to the momentum imbalance at the outlet, which is not fully ensured by the CBC when the external force term is solved in the potential form (Li, Jia & Liu 2017). The resulting velocity profile is also affected by the form of the characteristic velocity used in the CBC. Considering all the above tests, CBC-AV generally shows better performance and it is therefore used in the following studies.

In addition, since we find the flow behaviours are unaffected away from the outlet, we always use channel length which is much larger compared to the typical emulsion droplet, in order to minimise any undesirable effect from the outlet boundary condition.

Figure 3. (a) Morphology diagram for two equal-sized droplets with $\unicode[STIX]{x1D6FD}=0.001$; (b) emulsion shape with $\unicode[STIX]{x1D6FD}=0$ for $(\unicode[STIX]{x1D70E}_{gb}/\unicode[STIX]{x1D70E}_{rg},\unicode[STIX]{x1D70E}_{rb}/\unicode[STIX]{x1D70E}_{rg})=(1.4,0.35)$; (c) emulsion shape with $\unicode[STIX]{x1D6FD}=0.0001$ for $(\unicode[STIX]{x1D70E}_{gb}/\unicode[STIX]{x1D70E}_{rg},\unicode[STIX]{x1D70E}_{rb}/\unicode[STIX]{x1D70E}_{rg})=(1.4,0.35)$.

3.2 Morphology diagram

Since the interfacial tension relations are crucial in determining ternary emulsion morphologies (Pannacci et al. 2008; Guzowski et al. 2012), another validation test is conducted to show the capability of the current model in simulating a wide range of interfacial tension ratios. Following the theoretical analysis of Guzowski et al. (2012), two equal-sized red and green droplets are initially put next to each other and dispersed in the outer blue fluid. Three typical thermodynamic equilibrium morphologies can be obtained depending on the interfacial tension ratios of $\unicode[STIX]{x1D70E}_{gb}/\unicode[STIX]{x1D70E}_{rg}$ and $\unicode[STIX]{x1D70E}_{rb}/\unicode[STIX]{x1D70E}_{rg}$, as divided by the solid lines shown in figure 3(a): (I-A), complete engulfing with the red droplet inside the green one for $\unicode[STIX]{x1D70E}_{rb}/\unicode[STIX]{x1D70E}_{rg}>1+\unicode[STIX]{x1D70E}_{gb}/\unicode[STIX]{x1D70E}_{rg}$; (I-B), complete engulfing with the green droplet inside the red one for $\unicode[STIX]{x1D70E}_{gb}/\unicode[STIX]{x1D70E}_{rg}>1+\unicode[STIX]{x1D70E}_{rb}/\unicode[STIX]{x1D70E}_{rg}$; (II), non-engulfing, for $\unicode[STIX]{x1D70E}_{rb}/\unicode[STIX]{x1D70E}_{rg}+\unicode[STIX]{x1D70E}_{gb}/\unicode[STIX]{x1D70E}_{rg}<1$, where the red and green droplets tend to separate from each other; (III), partial engulfing (Janus droplet), for $|(\unicode[STIX]{x1D70E}_{rb}/\unicode[STIX]{x1D70E}_{rg})-(\unicode[STIX]{x1D70E}_{gb}/\unicode[STIX]{x1D70E}_{rg})|<1$ and $\unicode[STIX]{x1D70E}_{rb}/\unicode[STIX]{x1D70E}_{rg}+\unicode[STIX]{x1D70E}_{gb}/\unicode[STIX]{x1D70E}_{rg}>1$, where the interfacial tensions satisfy a Neumann triangle.

In our numerical test, the red and green droplets are both initialized with radius $R=60$ surrounded by the blue fluid in a domain of $L_{x}\times L_{y}=399\times 399$. All the fluid viscosities are 0.167. The initial concentration fractions for three fluids are given by (Yu et al. 2019b)

(3.2)$$\begin{eqnarray}\displaystyle & \displaystyle C_{1}(x,y)=0.5+0.5\tanh \left[\frac{R-\sqrt{(x-199.5)^{2}+(y-199.5-R)^{2}}}{2\unicode[STIX]{x1D6FC}}\right], & \displaystyle\end{eqnarray}$$

(3.3)$$\begin{eqnarray}\displaystyle & \displaystyle C_{2}(x,y)=0.5+0.5\tanh \left[\frac{R-\sqrt{(x-199.5)^{2}+(y-199.5+R)^{2}}}{2\unicode[STIX]{x1D6FC}}\right], & \displaystyle\end{eqnarray}$$

(3.4)$$\begin{eqnarray}\displaystyle & \displaystyle C_{3}(x,y)=1-C_{1}(x,y)-C_{2}(x,y). & \displaystyle\end{eqnarray}$$

Periodic boundary conditions are used for all boundaries. To reproduce all the possible morphologies, simulations are performed at various groups of ($\unicode[STIX]{x1D70E}_{gb}/\unicode[STIX]{x1D70E}_{rg}$, $\unicode[STIX]{x1D70E}_{rb}/\unicode[STIX]{x1D70E}_{rg}$): (I-A), complete engulfing with red droplet inside – (0.3, 1.35), (0.75, 1.8), (1.0, 2.05); (I-B), complete engulfing with green droplet inside – (1.4, 0.35), (1.75, 0.7), (2.0, 0.95); (II), non-engulfing – (0.48, 0.48), (0.25, 0.72), (0.75, 0.23); (III), partial engulfing emulsions – (1.0, 1.0), (1.0, 1.5), (1.0, 0.5), (0.5, 1.0), (1.5, 1.0), (100, 100). The interfacial tension $\unicode[STIX]{x1D70E}_{rg}$ is fixed at 0.005 except for the case with $(\unicode[STIX]{x1D70E}_{gb}/\unicode[STIX]{x1D70E}_{rg},\unicode[STIX]{x1D70E}_{rb}/\unicode[STIX]{x1D70E}_{rg})=(100,100)$, where $\unicode[STIX]{x1D70E}_{rg}=0.00001$ is used to reach the high interfacial tension ratio. The value of the coefficient $\unicode[STIX]{x1D6FD}$ is set to be 0.001 for the additional free-energy term. The simulated equilibrium morphologies are shown by the insets in figure 3(a). Good agreements with theoretical morphologies are achieved for all types of emulsions. Moreover, Pannacci et al. (2008) experimentally investigated the equilibrium states of compound emulsions. Their results are presented as a function of the spreading coefficients, i.e. $S_{i}=\unicode[STIX]{x1D70E}_{jk}-\unicode[STIX]{x1D70E}_{ij}-\unicode[STIX]{x1D70E}_{ik}$ with $i,j,k=r,g,b$, respectively. By converting the values of the interfacial tension ratios tested in figure 3 to spreading coefficients, our numerically obtained emulsion morphologies are also consistent with their experimental observations.

It is worth noting that we have investigated the optimal value of the coefficient $\unicode[STIX]{x1D6FD}$ in the additional free-energy term introduced in (2.9), varying $\unicode[STIX]{x1D6FD}=0$, 0.0001, 0.001, 0.01, 0.1 and 1.0 for one typical double emulsion morphology at ($\unicode[STIX]{x1D70E}_{gb}/\unicode[STIX]{x1D70E}_{rg}$, $\unicode[STIX]{x1D70E}_{rb}/\unicode[STIX]{x1D70E}_{rg}$) = (1.4, 0.35). The obtained result at $\unicode[STIX]{x1D6FD}=0$ (corresponding to the model without the additional term) is shown in figure 3(b). As highlighted by the dashed squares in figure 3(b), two unphysical light blue regions caused by negative $C_{1}$ appear around the three-phase contact line and lead to incorrect result. The incorrect region is also observed for $\unicode[STIX]{x1D6FD}=0.0001$ in figure 3(c). For $\unicode[STIX]{x1D6FD}$ varying from 0.001 to 0.1, the complete engulfing morphology could be successfully reproduced and invisible difference is observed for different values of $\unicode[STIX]{x1D6FD}$. However, further increasing $\unicode[STIX]{x1D6FD}$ to 1.0 induces numerical instability, which indicates that the $\unicode[STIX]{x1D6FD}$ value cannot be large enough to dominate the double-well potential terms. Meanwhile, for the partial engulfing cases, correct morphologies could be captured even without the additional term, and they are generally unaffected by a small additional term. Based on the above findings, $\unicode[STIX]{x1D6FD}=0.001$ will be used in subsequent simulations.

4 Results and discussion

4.1 Previously observed formation regimes and grid independence test

The two-dimensional set-up of the hierarchical flow-focusing device is illustrated in figure 4. The inner red fluid is injected through the leftmost inlet with a width of $w_{1}$, and the middle green and outer blue fluids are injected by two vertical side inlets with widths of $w_{2}$ and $w_{3}$, respectively. All the inlet widths are set equal in this section, i.e. $w_{1}=w_{2}=w_{3}$. The channel connecting the two side inlets has a width of $w_{4}=w_{1}$, and the main channel width is $w_{5}=1.6w_{1}$. The length of the first inlet is $w_{6}=2w_{1}$, and the distance between the two side inlets is $w_{7}=3w_{1}$. Considering the symmetry of the flow problem in the $y$ direction, only a half of the geometry is simulated and the domain size is $L_{x}\times L_{y}=30w_{1}\times 2w_{1}$. The Zou–He velocity inlet boundary condition (Zou & He 1997) is used for all the inlets, and the CBC-AV is applied for the outlet. In addition to the no-slip boundary condition, the wetting boundary condition is also imposed on the solid surfaces, where the first and second junctions are fully wetted by the middle and outer phase fluids, respectively.

Figure 4. Illustration of the geometry and boundary settings of the planar hierarchical flow-focusing device in this work.

In the following, the subscripts $i$, $m$ and $o$ are used to denote the inner, middle and outer fluids. The dimensionless parameters that characterize the double emulsion formation process are typically defined as follows (Abate et al. 2011; Azarmanesh et al. 2016): the Weber number (the ratio of inertia force to interfacial tension force) of the inner fluid $We_{i}=\unicode[STIX]{x1D70C}_{i}u_{i}^{2}w_{1}/\unicode[STIX]{x1D70E}_{im}$; capillary numbers (the ratios of viscous force to interfacial tension force) of the middle and outer fluids $Ca_{m}=\unicode[STIX]{x1D702}_{m}u_{m}/\unicode[STIX]{x1D70E}_{im}$, $Ca_{o}=\unicode[STIX]{x1D702}_{o}u_{o}/\unicode[STIX]{x1D70E}_{mo}$; flow rate ratios $Q_{i}/Q_{m}=u_{i}/(2u_{m})$, $Q_{o}/Q_{m}=(2u_{o})/(2u_{m})=u_{o}/u_{m}$; viscosity ratios $\unicode[STIX]{x1D706}_{im}=\unicode[STIX]{x1D702}_{i}/\unicode[STIX]{x1D702}_{m}$, $\unicode[STIX]{x1D706}_{om}=\unicode[STIX]{x1D702}_{o}/\unicode[STIX]{x1D702}_{m}$; interfacial tension ratios $\unicode[STIX]{x1D70E}_{io}/\unicode[STIX]{x1D70E}_{im}$ and $\unicode[STIX]{x1D70E}_{mo}/\unicode[STIX]{x1D70E}_{im}$. Here, $u$ is the average inlet velocity. However, in this work, we focus on double emulsion formation behaviours in the limit of small inertia (Nabavi et al. 2017b). As such, it is more appropriate to use $Ca_{i}=\unicode[STIX]{x1D702}_{i}u_{i}/\unicode[STIX]{x1D70E}_{im}$ instead of $We_{i}$ for the inner fluid. We will change the values of $Ca_{i}$, $Ca_{m}$ and $Ca_{o}$ by adjusting the flow rate of each phase fluid and investigate their roles in formation regime conversions and double emulsion sizes. The viscosity ratios are kept at $\unicode[STIX]{x1D706}_{im}=\unicode[STIX]{x1D706}_{om}=1.24$, and the interfacial tension ratio is fixed at $(\unicode[STIX]{x1D70E}_{mo}/\unicode[STIX]{x1D70E}_{im},\unicode[STIX]{x1D70E}_{io}/\unicode[STIX]{x1D70E}_{im})=(1.0,2.2)$.

Four basic flow regimes identified in the double emulsion preparation (Abate et al. 2011; Azarmanesh et al. 2016) are shown in figure 5(a1–a4): (a1) the two-step formation regime, (a2) one-step formation regime, (a3) decussate regime with one empty droplet, and (a4) decussate regime with two empty droplets. Our simulations are able to successfully reproduce all four regimes. Specially, the two-step formation regime is obtained at $Ca_{i}=0.012$, $Ca_{m}=0.011$ and $Ca_{o}=0.035$ (figure 5b1). With the same $Ca_{m}$ and $Ca_{o}$ values, the one-step formation regime is observed by increasing the inner flow rate to $Ca_{i}=0.018$ (figure 5b2), while the decussate regime with one empty middle phase droplet is achieved by decreasing the inner flow rate to $Ca_{i}=0.008$ (figure 5b3). Moreover, if the decussate regime happens at higher $Ca_{o}$ values, e.g. $Ca_{o}=0.065$ with $Ca_{i}=0.005$ and $Ca_{m}=0.011$, two empty alternate middle phase droplets are found, as shown by figure 5(b4). The corresponding $We_{i}$ values to the $Ca_{i}$ values used in figure 5 are generally of the order of $O(10^{-3})$ to $O(10^{-2})$, which are considerably lower than those used in previous studies ($O(1)$) (Abate et al. 2011; Azarmanesh et al. 2016). However, we note that similar two-step and one-step formation behaviours are still obtained. This suggests that, while $We_{i}$ can affect the resulting formation regimes, the rich flow behaviours with many different formation regimes are still present in the limit of small inertia (Wu, Cao & Sundén 2017b). Thus, we shall focus here on the limit of small $We_{i}$ to understand the interplay between viscous and interfacial tension forces. For this reason, it is reasonable to use $Ca_{i}$ for the inner phase fluid in our study, which also highlights the importance of the flow rate ratios in determining the resulting formation regimes.

Figure 5. The present work can qualitatively reproduce common flow regimes previously reported in experimental (Abate et al. 2011) (a1,a2) and simulation (Azarmanesh et al. 2016) (a3–a4) results. Parameters for the present work are (b1) $Ca_{i}=0.012$, $Ca_{m}=0.011$ and $Ca_{o}=0.035$; (b2) $Ca_{i}=0.018$, $Ca_{m}=0.011$ and $Ca_{o}=0.035$; (b3) $Ca_{i}=0.008$, $Ca_{m}=0.011$ and $Ca_{o}=0.035$; (b4) $Ca_{i}=0.008$, $Ca_{m}=0.011$ and $Ca_{o}=0.065$. Movies for the cases shown in panels (b1–b4) are provided online as supplementary materials available at https://doi.org/10.1017/jfm.2020.299.

A grid independence test is conducted for the two-step formation regime mentioned in figure 5(b1). Four different grid resolutions are tested, i.e. $w_{1}=40$, 50, 80 and 100. To make a quantitative comparison, the results from the highest grid resolution ($w_{1}=100$) are used as a reference. The relative errors ($E_{w_{1}}=|X_{w_{1}}-X_{w_{1}=100}|/X_{w_{1}=100}$) of the entire emulsion size, pinch-off location and generation period are calculated, and their maximum values are recorded for each grid resolution. The maximum relative errors for $w_{1}=40$, 50 and 80 are $7.18\,\%$, $4.66\,\%$ and $0.83\,\%$, respectively. This suggests that increasing grid resolution from $w_{1}=50$ to 100 leads to the relative error less than $5\,\%$, and thus an inlet width of $w_{1}=50$ is used for the following studies, as a good balance between computational accuracy and cost.

4.2 Effect of flow rates

In the formation of double emulsions, it is known that two-step, one-step and decussate formation regimes can be obtained by varying $Ca_{i}$, $Ca_{m}$ and $Ca_{o}$ values. However, the dependence of each formation regime on $Ca_{i}$, $Ca_{m}$ and $Ca_{o}$ values has not been systematically studied, and how the obtained emulsion sizes vary is not very clear. In figure 6, a three-dimensional phase diagram is constructed to illustrate how the formation regimes are influenced by $Ca_{i}$, $Ca_{m}$ and $Ca_{o}$. The ranges for these influencing parameters are $Ca_{i}=(0.008,0.01,0.012,0.014,0.016,0.018,0.02,0.022,0.025,0.028,0.03)$, $Ca_{m}=(0.005,0.011,0.015,0.02,0.025,0.03)$ and $Ca_{o}=(0.025,0.035,0.05,0.065)$. It is seen that more formation regimes are obtained besides those reported in figure 5.

To differentiate these regimes, each regime is represented by a unique symbol. The nomenclature for each regime is generally a combination of the breakup modes of the inner and middle phases. To distinguish the dripping and jetting breakup modes, we use the pinch-off location. It is considered as jetting if the distance between the pinch-off location of the inner (or middle) fluid and the downstream edge of the middle (or outer) fluid side inlet is longer than $3w_{1}$ (Utada et al. 2005). Otherwise, it is categorized as dripping. According to our nomenclature, the periodic two-step and one-step formation regimes shown in figure 5(a1,a2) are therefore called as dripping–dripping (regime 1) and jetting–dripping (regime 8) regimes in figure 6, which distinguishes them from other irregular two-step or one-step formation behaviours.

Figure 6. Flow regimes as a function of $Ca_{i}$, $Ca_{m}$ and $Ca_{o}$. Each symbol represents a unique formation pattern. Movies are provided online as supplementary materials for regimes 1–10.

Table 1. Relevant experimental literature for the formation regimes shown by figure 6.

To analyse these formation regimes, we classify the formation regimes on each $Ca_{o}$ plane. Firstly, the formation regimes are divided into two regions by the red solid line according to the breakup mode of the inner phase fluid. The inner phase fluid breaks up in the dripping mode on the left region of the red solid line. All the points in this region are periodic and they could be further subdivided into regimes 1 to 7. On the right-hand side of the red solid line, the inner phase breaks up in the jetting mode. The right-hand region can be further divided into two subsections by the dashed blue line based on the breakup mode of the middle phase fluid. Below the dashed line, the middle phase fluid breaks up in the dripping mode and we obtain the periodic one-step double emulsion formation regime (regime 8). Above the dashed line, the middle phase fluid also breaks up in the jetting mode, but the formation behaviour loses the periodicity. For instance, the inner and middle phases are pinched off together irregularly (regime 9), or multiple inner droplets of different sizes are randomly encapsulated in the middle phase droplet (regime 10). In an extreme case at $Ca_{i}=0.03$, $Ca_{m}=0.005$ and $Ca_{o}=0.065$, parallel layered flow is observed (regime 11).

To put the formation regimes obtained in figure 6 into experimental context, relevant experimental literature to each formation regime are listed in table 1. Regimes 1, 2, 8, 11 have been reported in planar microfluidic devices, while regimes 1–4, 8–10 have been observed in capillary microfluidic devices. However, the bidisperse behaviours observed in regimes 5–7 have only been reported in two-phase experiments so far, some of which are listed to regime 5 in table 1 for reference. A few experiments also present a multiple emulsion formation regime similar to regime 6 but with two equal-sized inner droplets. Some of these studies are listed to regime 6 in table 1. In all, figure 6 establishes the connection among different formation regimes.

In the following sections, we focus on the two-step (regimes 1–7) and one-step (regime 8) periodic regions. The effects of flow parameters on the conversion of formation regimes and emulsion sizes will be analysed in detail to help deepen the understanding of double emulsion formation behaviours.

4.2.1 Two-step periodic region

In the two-step formation region of figure 6, two types of periodic double emulsion formation regimes are observed, i.e. the dripping–dripping regime (regime 1) and the dripping–jetting regime (regime 2). Regime 1 is limited to a small range of governing flow parameters owing to the strict criterion in pinch-off locations. On the other hand, regime 2 occupies a relatively wider region, and the applicable range of $Ca_{m}$ for regime 2 shrinks to higher values as $Ca_{o}$ increases, due to the appearance of decussate regime at lower $Ca_{m}$. Moreover, the shape of the red solid lines varies little with $Ca_{o}$ over the entire range considered. It indicates that the breakup behaviour of the inner phase fluid is mainly determined by $Ca_{i}$ and $Ca_{m}$.

To clarify the effects of $Ca_{i}$, $Ca_{m}$ and $Ca_{o}$ on the two-step formation regimes, we illustrate the typical formation behaviours as a function of the $Ca_{i}$, $Ca_{m}$ and $Ca_{o}$, respectively, in figure 7. The parameters are (a) $Ca_{i}=0.008$, 0.012, 0.014 and 0.016 at $Ca_{m}=0.015$, $Ca_{o}=0.025$; (b) $Ca_{m}=0.005$, 0.015, 0.02 and 0.03 at $Ca_{i}=0.008$, $Ca_{o}=0.025$; (c) $Ca_{o}=0.025$, 0.035, 0.05 and 0.065 at $Ca_{i}=0.008$, $Ca_{m}=0.015$. Based on figure 7, we will discuss the size variations of double emulsion generated in the two-step formation regime. Moreover, new insights on other typical formation regimes obtained in the ternary system will also be discussed.

For the typical two-step formation regime shown in figure 7, the dripping/squeezing breakup mode of the inner phase fluid in the ternary system is similar to that which happens in a binary system. It is generally attributed to the action of the leading viscous force and the squeezing effect that overcome the interfacial tension force (Fu et al. 2012; Cubaud & Mason 2008; Yu et al. 2019a). For the breakup of the middle phase fluid, it is subject to both the viscous force of the outer phase fluid and the resulting flow from the generated inner droplets. The expansion of the middle phase front in the main channel also leads to accumulated upstream pressure in the outer phase fluid. All these factors assist in the breakup of the middle phase fluid.

Figure 7. Dynamics of two-step formation regimes as a function of (a-i–a-iv) $Ca_{i}=0.008$, 0.012, 0.014 and 0.016 at $Ca_{m}=0.015$ and $Ca_{o}=0.025$; (b-i–b-iv) $Ca_{m}=0.005$, 0.015, 0.02 and 0.03 at $Ca_{i}=0.008$ and $Ca_{o}=0.025$; (c-i–c-iv) $Ca_{o}=0.025$, 0.035, 0.05 and 0.065 at $Ca_{i}=0.008$ and $Ca_{m}=0.015$.

Figure 7(a-i–a-iii) shows the effect of $Ca_{i}$ on the two-step double emulsion formation behaviours. With increasing $Ca_{i}$, the inner droplet size decreases and its formation frequency increases. This trend has also been numerically observed by Fu et al. (2016) for double emulsions generated at high flow rates of the middle phase fluid in a two-dimensional simplified coaxial device. The increased formation frequency of the inner phase droplet shortens the time to accumulate the upstream pressure in the outer phase fluid and actuates the pinch-off of the middle phase front. Thus, the size of the middle layer of the entire double emulsion also decreases.

The effect of $Ca_{m}$ on two-step formation behaviours are given in figure 7(b). From figure 7(b-i) to (b-iii), it is noticed that the inner droplet size decreases while the middle part of the double emulsion increases. As the intermediate layer, the middle phase fluid has dual effects. With increasing $Ca_{m}$, the increased viscous force of the middle phase fluid exerted on the inner phase fluid leads to the size reduction of the inner droplet. Meanwhile, the increased middle flow rate decreases its velocity difference to that of the outer phase, which effectively lowers the outer shear stress and extends the time for the middle part of the double emulsion to grow larger. Our results show that the increase in the middle part size is usually more significant than the decrease in the inner droplet size. Thus, the entire double emulsion size increases.

Figure 7(c) displays the effects of $Ca_{o}$ on double emulsion formation behaviours. As seen in figure 7(c-i–c-ii), increasing $Ca_{o}$ does not change the double emulsion size in any obvious way before the formation regime changes, but the distance between the generated double emulsions gets larger due to the increased outer flow rate. Further increasing $Ca_{o}$ to 0.05, decussate regime with one alternate empty droplet appears and the shell thickness of the double emulsion is greatly reduced. When $Ca_{o}$ reaches 0.065, double emulsions with two alternate empty droplets are captured.

For the periodic two-step region shown in figure 6, increasing $Ca_{i}$ or $Ca_{m}$ would both lead to the dripping–threading regime, where the inner droplet is produced periodically in the continuous middle phase thread. Two examples are given in figure 7(a-iv,b-iv). We would like to compare the differences between the dripping–threading morphologies obtained by adjusting $Ca_{i}$ and $Ca_{m}$, respectively. The inner droplets are produced in small sizes in both cases, while the formation frequency is higher at large $Ca_{i}$ than that obtained at large $Ca_{m}$. As a result, the capillary perturbations on the middle phase fluid in case (a-iv) is counteracted by the high formation frequency of the inner phase droplets and the obtained dripping–threading regime is very stable. Nabavi et al. (2017b) experimentally reported a similar stable dripping–threading regime also by increasing the inner flow rate in a capillary device. Unlike in case (a-iv), it is seen in (b-iv) that some necking regions develop at the middle phase fluid as it flows downstream, which would possibly lead to the middle phase breakup somewhere more downstream. Such unstable regimes as shown in figure 7(b-iii–b-iv) remind us of the varicose shape reported in binary experiments (Cubaud & Mason 2008). The narrow main channel limits the expansion of the middle phase front to form an emulsion, and the following embryonic emulsion shape begins to grow before the front one is pinched off.

In view of applications, the compound structure generated in the dripping–threading regime is capable of producing bundles of microcapsules that are promising for storing, handling and arrayed assay of small volumes (Oh et al. 2006). To remove the dripping–threading regime, we can increase $Ca_{o}$ to produce regular double emulsions.

In figure 7(a-ii), size variations are observed in the generated double emulsion sequence in the main channel: a smaller double emulsion is followed by a larger one, and this pattern repeats itself. To reveal the periodicity of this behaviour, the temporal evolution of the inner and middle thread tip locations are traced as denoted by $X_{i}$ and $X_{m}$ in figures 8(a1) and 8(a2), respectively. The time $t$ and locations $X_{i,m}$ are normalized using $t^{\ast }=t(u_{m})_{max}/w_{1}$ and $X_{i,m}^{\ast }=X_{i,m}/w_{1}$, where $(u_{m})_{max}=0.0015$ is the maximum flow rate of the middle phase fluid used in the current study. After the double emulsions are produced regularly, the points corresponding to the pinch-off moment and location in each formation period are marked by the superimposed round circles for the inner phase fluid in figure 8(a1) and diamonds for the middle phase fluid in figure 8(a2). Clearly, periodic fluctuations in pinch-off locations and formation periods are observed in both the inner and middle phases between every two neighbouring droplets, which is consistent with the variation in emulsion sizes observed in figure 7(a-ii). This flow pattern is named as the in-mid-bidisperse regime (regime 5).

Figure 8. Temporal evolutions of the thread tip locations of the inner ($X_{i}^{\ast }$) and middle ($X_{m}^{\ast }$) phases obtained for: (a1,a2), $Ca_{i}=0.012$, $Ca_{m}=0.015$ and $Ca_{o}=0.025$; and (b1,b2), $Ca_{i}=0.012$, $Ca_{m}=0.02$ and $Ca_{o}=0.025$. The time and location are normalized using $t^{\ast }=t(u_{m})_{max}/w_{1}$ and $X_{i,m}^{\ast }=X_{i,m}/w_{1}$, where $(u_{m})_{max}=0.0015$ is the maximum flow rate of the middle phase used in the current study. The superimposed empty round circles in (a1,b1) and diamond symbols in (a2,b2) mark the periodic points that correspond to each pinch-off moment and location of the inner and middle phases, respectively. The inset snapshots in (a1) and (b1) show the corresponding flow behaviours in each case.

Regime 5 is frequently observed for $Ca_{i}$ between 0.012 and 0.014 with $Ca_{m}$ approximately from 0.015 to 0.03 on each $Ca_{o}$ plane. The earlier occurrence time of the inner phase bidispersity observed in figure 8(a1,a2) suggests that such bidisperse behaviours mostly originate from the inner phase fluid and then propagate to the middle phase fluid. Since the breakup mode of the inner phase fluid is rarely affected by $Ca_{o}$, the reason for the inner phase bidispersity should be similar to that in a binary system. It is normally attributed to the oscillations in the amount of residual liquid on the entrance side after the previous droplet is pinched off (Coullet, Mahadevan & Riera 2005; Garstecki, Fuerstman & Whitesides 2005; Utada et al. 2007).

Noteworthy, the influence of inner phase bidispersity brings richer dynamics in the present ternary system depending on $Ca_{m}$ and $Ca_{o}$ values. For cases like the one shown in figure 7(a-ii), the middle phase fluid is easy to be pinched off and it follows the bidisperse breakup frequencies of the inner phase fluid (regime 5). However, for a flow condition with a high $Ca_{m}$ and a low $Ca_{o}$, the thicker middle phase fluid could extend its pinch-off time and engulf every two inner phase droplets inside, as shown by one typical case at $Ca_{i}=0.012$, $Ca_{m}=0.02$ and $Ca_{o}=0.025$ in figure 8(b1,b2). As such, the variation in formation frequency only happens in the inner phase fluid (figure 8b1), but not in the middle phase fluid (figure 8b2). It is named as the in-bidisperse regime (regime 6). In addition, even if the middle phase fluid forms a continuous thread, e.g. at $Ca_{i}=0.014$, $Ca_{m}=0.03$ and $Ca_{o}=0.025$, the inner bidisperse behaviour could still happen, and it is named as the bidisperse–threading regime (regime 7).

Decussate regimes occupy a substantial proportion in the two-step formation region of figure 6. Among them, decussate regimes with one alternate empty droplet are commonly observed while decussate regimes with two empty droplets mainly happen at high $Ca_{o}$ values. Figure 7(c-iv) gives one example of the decussate regime with two empty droplets, and the formation process of the two empty droplets is shown in figure 9(a). It is seen that a long section of the middle phase fluid is pinched off entirely, and then it breaks up into two daughter droplets during the retraction process of the stretched structure when flowing downstream. Azarmanesh et al. (2016) numerically reported another type of formation process for decussate regime with two empty droplets, where the two empty droplets are produced one by one. Our results show that by lowering the $Ca_{m}$ of the case shown in figure 9(a) to 0.011 in figure 9(b), the formation process reported by Azarmanesh et al. (2016) is reproduced in our work.

Figure 9. Decussate regimes with two empty droplets: (a) $Ca_{i}=0.008$, $Ca_{m}=0.015$ and $Ca_{o}=0.065$; (b) $Ca_{i}=0.008$, $Ca_{m}=0.011$ and $Ca_{o}=0.065$.

Comparing figure 9(a,b), the only difference lies in $Ca_{m}$. A lower $Ca_{m}$ signifies a higher velocity difference between the middle and the outer phases, which leads to a stronger viscous force exerted on the middle phase fluid and contributes to the early pinch-off of the middle phase front around the bulb neck. With regard to figure 9(a) at a higher $Ca_{m}$, the middle phase front is not pinched off until the entrance of the inner droplet that prevents the continuous injection of the middle phase fluid to its thread tip.

Decussate regimes are also of practical significance. For instance, if the downstream channel is connected to an expansion channel, the empty droplet can catch up with the double emulsion droplet ahead and merge to form a large double emulsion with thicker middle layer (Azarmanesh et al. 2016). Moreover, the empty droplet and the double emulsion droplet can be viewed as two distinct inner components to produce more complex functional multiple emulsions (Nisisako et al. 2005).

4.2.2 One-step periodic region

Even though both the two-step (regimes 1–2) and one-step (regime 8) formation regimes can be used for producing double emulsions, the one-step regime is advantageous over the two-step regime in several aspects, such as better robustness in wetting conditions, producing double emulsions with thinner middle layer and emulsifying non-Newtonian fluids (Abate et al. 2011). Moreover, from the point of view of the applicable parameter ranges shown in figure 6, a wide range of one-step formation points distributed continuously, different from the characteristic distribution of the periodic two-step double emulsion formation regimes with interference from other flow regimes. Therefore, the periodic one-step double emulsion region has more statistical significance over the periodic two-step region. It enables us to investigate the one-step formation mechanism more quantitatively and construct possible scaling laws for the double emulsion sizes.

Figure 10. Time evolution of the interface dynamics at $Ca_{i}=0.018$, $Ca_{m}=0.011$ and $Ca_{o}=0.05$. In each subfigure, the interface shapes are denoted by the solid lines. The distribution of the normalized viscous force component $\unicode[STIX]{x1D70F}_{xy}^{\ast }=\unicode[STIX]{x1D70F}_{xy}w_{1}/\unicode[STIX]{x1D70E}_{im}$ is shown in the upper part, while the streamlines are shown in the lower part. The time is normalized using $t^{\ast }=t(u_{m})_{max}/w_{1}$.

In figure 6, the applicable range of $Ca_{i}$ for regime 8 increases with $Ca_{o}$ and decreases with $Ca_{m}$. The latter trend is consistent with the experimental observation of Kim et al. (2013) in the study of periodic one-step double emulsion using a capillary device. To better understand the one-step double emulsion formation process, the typical temporal evolution of the interface dynamics at $Ca_{i}=0.018$, $Ca_{m}=0.011$ and $Ca_{o}=0.05$ is shown in figure 10. In each subfigure, the interface shapes are depicted by the solid lines. The leading viscous force component is displayed in the upper part, i.e. $\unicode[STIX]{x1D70F}_{xy}=\unicode[STIX]{x1D702}(\unicode[STIX]{x2202}u_{y}/\unicode[STIX]{x2202}x+\unicode[STIX]{x2202}u_{x}/\unicode[STIX]{x2202}y)$ for a two-dimensional system, and it is normalized using $\unicode[STIX]{x1D70F}_{xy}^{\ast }=\unicode[STIX]{x1D70F}_{xy}w_{1}/\unicode[STIX]{x1D70E}_{im}$. The streamlines are shown in the lower part. Figure 10(a1) corresponds to the moment just after a previous double emulsion is pinched off, where a strong shear stress region is activated to resist the retraction of the highly deformed middle–outer interface. During the evolution from figures 10(a1) to 10(a3), the middle phase thread tip approximately recovers to a semicircular shape under the effect of interfacial tension (Utada et al. 2007; Fu et al. 2012). In the meantime, the highest shear stress is lowered, and a more evenly distributed high shear stress region is formed along the inner–middle interface. The inflation of the compound inner and middle thread tip partially blocks the inflow of the outer phase fluid. Then, the outer fluid squeezes back the expanded compound thread tip and stretches it downstream. An obvious neck region is formed in figure 10(a4) and it keeps shrinking until the pinch-off happens in the inner phase fluid as shown in figure 10(a5). It is seen that a higher positive and a lower negative shear stress regions are induced immediately near the newly pinched inner thread tip and the generated inner droplet, respectively. The weakly connected middle thin thread is pinched off just after the configuration in figure 10(a6).

Based on the analysis of figure 10, the double emulsion formation process in the one-step regime can be approximately viewed as the sum of a partial blocking period and a squeezing period, which is analogous to that of a single droplet formation process in squeezing or dripping regime in binary flow-focusing systems (Cubaud & Mason 2008; Liu & Zhang 2011; Fu et al. 2012). Thus, the scaling laws proposed in binary flow-focusing systems could be used as the basis for the construction of scaling laws of the present ternary system. To complete the scaling law involving the ternary fluid flows, we first analyse the effects of governing flow parameter on the double emulsion sizes.

As shown in figure 11(a1), the effect of $Ca_{i}$ is investigated for $Ca_{i}=0.016$, 0.018, 0.02 and 0.022 at $Ca_{m}=0.011$ and $Ca_{o}=0.05$. The areas of the entire double emulsion $A_{emulsion}$, the inner part $A_{inner}$ and the middle part $A_{middle}$ are measured after the double emulsion is produced periodically. The area quantities are normalized using $A^{\ast }=A/(\unicode[STIX]{x03C0}w_{1}^{2})$. Compared to the effect of increasing $Ca_{i}$ in the two-step double emulsion formation regime given in figure 7(a-i–a-iii), the inner droplet size decreases in the two-step formation regime, while it varies little except for an initial minor increase in the one-step formation regime. Since the inner droplet size varies little, the time needed for the inner phase fluid to breakup is shortened with the increase of $Ca_{i}$ (Utada et al. 2007). This leads to the size reduction in the middle part and the entire double emulsion size, which is qualitatively similar to the effect of $Ca_{i}$ observed in two-step formation regimes. By investigating all the periodic one-step data shown in figure 6, the variations in double emulsion sizes caused by $Ca_{i}$ are qualitatively the same for other $Ca_{m}$ and $Ca_{o}$ conditions.

Figure 11. Inner part, middle part and the entire double emulsion size variations in the periodic one-step formation regime (regime 8) by changing: (a1) $Ca_{i}=0.016$, 0.018, 0.02 and 0.022 at $Ca_{m}=0.011$ and $Ca_{o}=0.05$; (b1) $Ca_{m}=0.005$, 0.011 and 0.015 at$Ca_{i}=0.018$ and $Ca_{o}=0.05$; (c1) $Ca_{o}=0.025$, 0.035, 0.05 and 0.065 at $Ca_{i}=0.018$ and $Ca_{m}=0.011$. The snapshots shown in (a2-i–a2-iv), (b2-i–b2-iii) and (c2-i–c2-iv) correspond to the flow conditions mentioned in (a1–c1).

The effect of $Ca_{m}$ on the size of each component and the corresponding snapshots are illustrated in figure 11(b1,b2) at $Ca_{m}=0.005$, 0.011 and 0.015 with $Ca_{i}=0.018$ and $Ca_{o}=0.05$. As $Ca_{m}$ increases, the inner droplet size decreases and the middle part increases. The middle part size always increases faster than the decrease of the inner droplet size. Thus, the entire double emulsion size increases monotonically with $Ca_{m}$. These trends qualitatively agree with the characteristic size variation obtained in the two-step regimes as shown in figure 7(b-i–b-iii). It indicates the same effects of $Ca_{m}$ on both formation regimes. We further verify that varying $Ca_{i}$ and $Ca_{o}$ conditions in figure 6 does not change the effects of $Ca_{m}$.

We have learned the effects of $Ca_{o}$ on two-step formation regimes in figure 7(c): the inner droplet size is almost independent of $Ca_{o}$, but the breakup frequency of the middle phase increases with increasing $Ca_{o}$, which could further lead to the decussate regime. However, a different effect of $Ca_{o}$ is expected in the one-step formation regime since the inner and middle phase fluids are emulsified simultaneously. In figure 11(c1,c2), $Ca_{o}$ is increased from 0.025, 0.035, 0.05 to 0.065 at $Ca_{i}=0.018$ and $Ca_{m}=0.011$. As $Ca_{o}$ increases, identical variation trends occur to the inner part, middle part and the entire double emulsion sizes: the sizes consistently increase slightly at the very beginning and then decrease monotonically. For other $Ca_{i}$ and $Ca_{m}$ values investigated in figure 6, the initial increase in sizes is not common with increasing $Ca_{o}$, but the decreasing trend is always obtained due to the enhanced viscous force at larger $Ca_{o}$. Therefore, for the purpose of constructing the scaling law on the double emulsion sizes, the occasional increasing trend is neglected, and we will assume the size has a decreasing trend with increasing $Ca_{o}$.

To construct a phenomenological scaling law for the size of double emulsion produced in the one-step regime, we take inspiration from the scaling laws developed for the size of single droplet generated in squeezing regime within a single cross-junction. Several researchers have contributed to the development of droplet size scaling laws in such binary systems (Garstecki et al. 2006; Christopher et al. 2008; Tan et al. 2008; Liu & Zhang 2011). Specifically, Liu & Zhang (2011) developed a scaling law for the length $L_{p}$ of the obtained plug shape droplet, which is given by

(4.1)$$\begin{eqnarray}\frac{L_{p}}{w_{1}}=\left(\tilde{\unicode[STIX]{x1D716}}+\tilde{\unicode[STIX]{x1D6FE}}\frac{Q_{dispersed}}{Q_{continuous}}\right)Ca_{o}^{\tilde{m}},\end{eqnarray}$$

where the plug length is normalized by the inlet width $w_{1}$, and $\tilde{\unicode[STIX]{x1D716}}$, $\tilde{\unicode[STIX]{x1D6FE}}$ and $\tilde{m}$ are fitting parameters. $Q_{dispersed}/Q_{continuous}$ is the flow rate ratio between the dispersed droplet phase and the continuous carrier phase.

In (4.1), the contributions of the blocking and squeezing processes for the size of the obtained droplet are described by the first and second terms in the bracket, respectively. It also includes the power-law dependence of the droplet size on the outer phase capillary number as pointed out by Christopher et al. (2008). Moreover, the work of Liu & Zhang (2011) showed that for a droplet produced at different width or height conditions, the fitting parameters will be affected, but the variation of droplet size still obeys the generalized expression of (4.1). The good agreement with available results justifies the validity of this scaling law (Liu & Zhang 2011), and it is therefore used as the basis to construct a size scaling law for the double emulsion.

Besides the similarities, we would like to highlight the differences between the binary and ternary systems so as to extend (4.1) for the ternary system. Firstly, the droplet length is only suitable for plug-shaped droplets whose diameter is wider than the channel width (Garstecki et al. 2006; Liu & Zhang 2011). The volume values are more general quantities for the ellipsoid-like double emulsions (Steegmans, Schron & Boom 2009; Chang et al. 2009; Fu et al. 2016). The volume quantities are also more convenient to measure the size of each part of the double emulsion than the length quantities. Hence, the areas of each part of the double emulsion are monitored as discussed in the above subsection. The equivalent radius of the double emulsion area defined by $R_{emulsion}=\sqrt{A_{emulsion}/\unicode[STIX]{x03C0}}$ is used as the dependent variable of the scaling law. Secondly, the dispersed phase in the one-step formation regime is made up of both the inner and middle phases for the continuous outer phase fluid. The independent control over the inner and middle flow properties make the flow behaviours more complex. Based on the analysis of flow parameter effects shown in figure 11, $Q_{dispersed}/Q_{continuous}$ in (4.1) is replaced by $Q_{m}/Q_{i}$ to incorporate the positive effect of $Ca_{m}$ and the negative effect of $Ca_{i}$ on the entire double emulsion size.

Thus, the scaling law for $R_{emulsion}$ is constructed as

(4.2)$$\begin{eqnarray}\frac{R_{emulsion}}{w_{1}}=\left(0.270+0.0526\frac{Q_{m}}{Q_{i}}\right)Ca_{o}^{-0.268},\end{eqnarray}$$

where the parameter values 0.270, 0.0526 and $-0.268$ are obtained by fitting all the investigated periodic one-step data shown in figure 6 with the principle of minimum residual norm. To test the obtained scaling law, the values of the double emulsion radius $(R_{emulsion}/w_{1})_{pred}$ computed from (4.2) are plotted against the simulated radius values $(R_{emulsion}/w_{1})_{simu}$ in figure 12(a). The line of parity is plotted as a reference, and the closer the scattered data points are to the line of parity, the better the agreement is between the scaling law and the simulated results. It is seen that most of the points scatter around the line of parity, and the simple formula of (4.2) can provide a general guidance for predicting double emulsion size.

Figure 12. The parity plots of (a) the normalized entire double emulsion radius $(R_{emulsion}/w_{1})_{pred}$ computed from (4.2) and the simulated values of $(R_{emulsion}/w_{1})_{simu}$; (b) the radius ratio of the inner part to the entire double emulsion $(R_{inner}/R_{emulsion})_{pred}$ computed from (4.3) and the simulated values of $(R_{inner}/R_{emulsion})_{simu}$. The points in both plots are based on all periodic one-step flow conditions (regime 8) obtained in figure 6. The legend table shows that the flow conditions of all feasible $Ca_{i}$ at each $Ca_{m}$ and $Ca_{o}$ combination are represented by symbols of the same type, and the values of $Ca_{m}$ and $Ca_{o}$ are differentiated through the symbol colours and shapes, respectively. The inset in subfigure (a) shows the snapshot of one typical periodic one-step formation regime at $Ca_{i}=0.02$, $Ca_{m}=0.011$ and $Ca_{o}=0.05$.

Another size of interest is the ratio between the equivalent inner droplet radius $R_{inner}=\sqrt{A_{inner}/\unicode[STIX]{x03C0}}$ and the entire double emulsion radius: $R_{inner}/R_{emulsion}$. Chang et al. (2009) experimentally proposed a scaling law for the double emulsion generated in coaxial capillaries. The inner droplet and the entire double emulsion are viewed to have the same formation time before being pinched off together in the dripping mode. According to the mass conservation law, $R_{inner}/R_{emulsion}$ is predicted by $R_{inner}/R_{emulsion}=(Q_{i}/(Q_{i}+Q_{m}))^{n}$, and the power-law exponent $n$ is $1/3$ and $1/2$ for three- and two-dimensional studies, respectively. Recently, Fu et al. (2016) numerically confirmed this relation in their two-dimensional study using a coaxial capillary device. However, the inner phase fluid actually breaks up slightly earlier than that of the middle phase fluid, especially in the current planar hierarchical flow-focusing device (see figure 10). The difference in formation time between the two phases is also observed to be moderately affected by $Ca_{m}$ and $Ca_{o}$. To consider their effects, two power-law relations are assumed between $R_{inner}/R_{emulsion}$ and $Ca_{m}$ and $Ca_{o}$, respectively. A scale factor is also added to finely tune the entire size.

Based on the above analysis, the scaling law for $R_{inner}/R_{emulsion}$ is constructed as

(4.3)$$\begin{eqnarray}\displaystyle \frac{R_{inner}}{R_{emulsion}}=0.904\left(\frac{Q_{i}}{Q_{i}+Q_{m}}\right)^{0.609}Ca_{m}^{-0.060}Ca_{o}^{0.030}. & & \displaystyle\end{eqnarray}$$

The way to obtain the values of the coefficients in (4.3) is the same as that used in (4.2). The fitted power-law exponent of $Q_{i}/(Q_{i}+Q_{m})$ is 0.609, which is close to 0.5 mentioned in the work of Chang et al. (2009) and Fu et al. (2016). The difference can be attributed to the inconsistency in the breakup time of the inner and middle phases. Nevertheless, the difference in the formation time is small, which is also reflected by the scale factor 0.904 that is close to 1.0, and the near zero power-law exponents of $Ca_{m}^{-0.060}$ and $Ca_{o}^{0.030}$. Similar to figure 12(a), the parity plot for the computed values of $(R_{inner}/R_{emulsion})_{pred}$ using (4.3) and the measured values $(R_{inner}/R_{emulsion})_{simu}$ are shown in figure 12(b). The good agreement between the scattered points and the parity line justifies the validity of the scaling law of (4.3) for the $R_{inner}/R_{emulsion}$ values.

4.3 Effect of interfacial tension ratio

In figure 3, we show that a variation in the interfacial tension ratio could result in distinct equilibrium morphologies of two droplets of different fluids. To elucidate the role of interfacial tension ratios on the emulsion structure in different double emulsion formation processes, six groups of interfacial tension ratios that cover different regions of figure 3 are investigated, i.e. $(\unicode[STIX]{x1D70E}_{mo}/\unicode[STIX]{x1D70E}_{im},\unicode[STIX]{x1D70E}_{io}/\unicode[STIX]{x1D70E}_{im})=(1.0,2.2),(2.2,1.0),(0.48,0.48),(1.0,0.5),(1.0,1.5)$ and (100, 100) under two flow conditions for periodic two-step (regime 1) and one-step (regime 8) formation regimes. The flow parameters for the two-step and one-step formation regimes are given at $Ca_{i}=0.012$ and $Ca_{i}=0.02$, respectively, with $Ca_{m}=0.011$ and $Ca_{o}=0.035$. The corresponding flow rate ratios are $Q_{i}:Q_{m}:Q_{o}=0.171:0.390:1$ and $0.286:0.390:1$. To obtain different interfacial tension ratios, $\unicode[STIX]{x1D70E}_{im}$ is fixed at 0.005 except for the case at $(\unicode[STIX]{x1D70E}_{mo}/\unicode[STIX]{x1D70E}_{im},\unicode[STIX]{x1D70E}_{io}/\unicode[STIX]{x1D70E}_{im})=(100,100)$, where $\unicode[STIX]{x1D70E}_{im}=0.00001$ is used, similar to those used in figure 3. Figure 13 illustrates the snapshots of the (a) two-step and (b) one-step flow rates for each interfacial tension ratio group. Note that the first column before the panel (a) series shows the corresponding static equilibrium morphology of each interfacial tension ratio group as shown in figure 3. Relevant experimental works are marked next to the related snapshots.

It is seen in figure 13 that the formation details and the emulsion morphologies are greatly affected by the interfacial tension ratios in both formation regimes. Firstly, compared to the double emulsions obtained at $(\unicode[STIX]{x1D70E}_{mo}/\unicode[STIX]{x1D70E}_{im},\unicode[STIX]{x1D70E}_{io}/\unicode[STIX]{x1D70E}_{im})=(1.0,2.2)$ (figure 13a1,b1), the inverse engulfed double emulsion is captured in figure 13(a2,b2) by reversing the interfacial tension ratios to $(\unicode[STIX]{x1D70E}_{mo}/\unicode[STIX]{x1D70E}_{im},\unicode[STIX]{x1D70E}_{io}/\unicode[STIX]{x1D70E}_{im})=(2.2,1.0)$. With the inverse interfacial tension ratios, the inner phase fluid is more favoured to the outer phase fluid and tends to engulf the middle phase droplet to lower the system’s interfacial energy. In the two-step formation regime shown in figure 13(a2), as the individually generated inner droplet approaches the second cross-junction, it is getting closer to the middle–outer interface. Once the inner droplet touches the middle–outer phase interface, the attraction between the inner and outer phases would prompt the pinch-off of the middle phase layer between them and actuate the formation of the middle phase droplet. Afterwards, the inner droplet itself becomes a bridge connecting the newly formed middle phase droplet and the remaining middle phase front. Soon it breaks into two parts under the viscous force of the outer fluid. The inner phase portion adhered to the middle phase droplet evolves to wrap the middle phase droplet and the inverse double emulsion morphology is finally formed. In the one-step formation regime of figure 13(b2), the inverse double emulsion is also obtained. However, the formation details are different due to the continuous supply of the inner phase fluid in the jetting mode. A string of small middle phase droplets are formed and connected by the inner phase fluid. The compound thread tip is then emulsified by the outer fluid for every two front middle phase droplets. The detached two middle phase droplets covered by the inner phase fluid soon merge with each other and produce a pure double emulsion.

Figure 13. Snapshots of emulsion formation behaviours under the effect of interfacial tension ratios in (a) two-step and (b) one-step formation regimes. The interfacial tension ratios for cases (a1,b1–a6,b6) are $(\unicode[STIX]{x1D70E}_{mo}/\unicode[STIX]{x1D70E}_{im},\unicode[STIX]{x1D70E}_{io}/\unicode[STIX]{x1D70E}_{im})=(1.0,2.2)$, (2.2, 1.0), (0.48, 0.48), (1.0, 0.5), (1.0, 1.5) and (100, 100). The first column before the panel (a) series shows the corresponding static equilibrium morphology of two equal-sized droplets at each interfacial tension ratio group. Relevant experimental works are put next to the related snapshots. The white square marked in panel (a1) highlights the region that the inner droplet is about to touch the middle–outer interface.

Chao et al. (2016) experimentally captured the conversion from an initial double emulsion to its inverse structure using a glass-based capillary microfluidic device. Using the terminology of our work, an intermediate red-in-green-in-blue double emulsion is initially produced in their work, and the thermodynamic equilibrium green-in-red-in-blue configuration is only obtained after the external flow is stopped. However, in our work, the final configuration is formed directly without the intermediate red-in-green-in-blue configuration. This implies that the moment for interfacial tension dominating over the hydrodynamic effects in the formation behaviours is earlier in our simulations than that in the experimental work of Chao et al. (2016). This could be explained by the experimental findings of Pannacci et al. (2008). They pointed out that it is necessary for the inner droplet to touch the inner boundary of its host to evolve to thermodynamic equilibrium under the capillary forces. In other words, the sooner the three-phase contact line is formed, the faster the interfacial tension starts to dominate. For instance, if we look into the formation details in figure 13(a2), there should be an instantaneous moment, like that highlighted in the square region in figure 13(a1), where the inner droplet is approaching the middle–outer interface due to the squeezing of the outer fluid. It allows the capillary force to act earlier. Regarding the experimental work of Chao et al. (2016), a relatively thick middle layer surrounds the inner phase orifice in the coaxial glass capillaries, which could prevent the early formation of the three-phase contact line, and hence delay the interfacial tension effect.

At $(\unicode[STIX]{x1D70E}_{mo}/\unicode[STIX]{x1D70E}_{im},\unicode[STIX]{x1D70E}_{io}/\unicode[STIX]{x1D70E}_{im})=(0.48,0.48)$, the red and green droplets tend to separate from each other at thermodynamic equilibrium. In the two-step formation regime shown in figure 13(a3), the inner and middle phase droplets are successively formed and flow downstream without touching each other in the outer fluid, consistent with their static equilibrium morphologies. These alternately generated single droplets of two phase fluids have possible applications in being the source materials for producing multicore emulsions (Nisisako et al. 2005). In figure 13(b3), a more complex multiple emulsion is obtained in the one-step formation regime: an inner phase droplet is seized by two middle phase droplets on both sides in the flow direction. The contact length between the components of the multiple emulsion is decreasing when flowing downstream, but the components do not completely separate from each other in the finite computational domain. It can be attributed to two possible reasons. The first one is that the sequence structure results from a transient double emulsion rather than separately produced like in the two-step formation regime. Thus, it takes longer for the sequence structure to evolve to its thermodynamic equilibrium. Secondly, once the middle phase thread tip is pinched off, the lateral outer phase fluid rapidly fills the pinch-off region. Consequently, the most upstream component in the sequence is more accelerated and the hydrodynamic effects keep the three components staying next to each other. The complete separation of the components could be expected after the inflow pumps are stopped.

For the three cases shown in figures 13(a4–a6) and 13(b4–b6), since the interfacial tensions in each case satisfy the Neumann triangle relation, the partial engulfing (Janus) emulsion should be achieved at thermodynamic equilibrium. Zhang et al. (2017) experimentally captured the transformation from the core–shell structure to the Janus droplets based on prefabricated double emulsions. Here, our results in figure 13(a4) show that the Janus droplet could be produced directly in the two-step formation regime within the same device for producing double emulsions. For figures 13(a5), 13(b4) and 13(b5), biconcave and biconvex emulsions are formed downstream. These structures are analogous to those experimentally fabricated by Nisisako, Ando & Hatsuzawa (2015) and Nie et al. (2006). Finally, for $(\unicode[STIX]{x1D70E}_{mo}/\unicode[STIX]{x1D70E}_{im},\unicode[STIX]{x1D70E}_{io}/\unicode[STIX]{x1D70E}_{im})=(100,100)$ shown in figure 13(a6,b6), $\unicode[STIX]{x1D70E}_{im}$ is so small that the inner and middle phase fluids can be approximately viewed as the same fluid, and the high $Ca_{i}$ induced by small $\unicode[STIX]{x1D70E}_{im}$ easily leads to the parallel layered flow behaviours for both the two-step and one-step flow conditions.

4.4 Effect of geometry

Geometrical parameters in microfluidics are usually the key factors in single or double emulsion preparations (Liu & Zhang 2011; Nabavi et al. 2015b; Wu et al. 2017a). In this section, we focus on the effect of the geometrical parameters in changing the double emulsion formation regimes and the obtained double emulsion sizes. For the geometry shown in figure 4, six normalized geometrical parameters can be defined as $w_{2}/w_{1}$, $w_{3}/w_{1}$, $w_{4}/w_{1}$, $w_{5}/w_{1}$, $w_{6}/w_{1}$ and $w_{7}/w_{1}$. Among them, the inlet length $w_{6}/w_{1}$ can be neglected, since the fully developed velocity distribution is always provided at the inlet, and the inner phase flow profile varies little before it reaches the middle phase inlet junction. Then, for simplicity, we make two assumptions to reduce the governing geometrical parameters, i.e. the side inlets for the middle and outer phase fluids have equal widths ($w_{3}=w_{2}$), and the width of the channel connected the side inlets is set equal to that of the inner phase inlet ($w_{4}=w_{1}$). Therefore, the main geometrical factors are reduced to the side inlet width ($w_{2}/w_{1}$), the main channel width ($w_{5}/w_{1}$) and the distance between the side inlets of the middle and outer phase fluids ($w_{7}/w_{1}$). Those geometrical factors are all investigated at two flow rates that lead to two-step and one-step formation regimes, respectively, for the original geometry. Different from the flow conditions used in the interfacial tension effect section, two closer $Ca_{i}$ values of 0.014 and 0.016 are used in this section at $Ca_{m}=0.011$ and $Ca_{o}=0.035$, to show the geometrical effect more obviously in changing the formation regimes.

Figure 14. Snapshots of double emulsion formation behaviours under the effect of geometrical parameters using the flow conditions that lead to (a) two-step and (b) one-step regimes in the original geometry. The results for the original geometry are marked with an inverted triangle. (a1,b1) $w_{2}/w_{1}$ ranges from 0.8, 1.0, 1.2 to 1.4 at $w_{5}/w_{1}=1.6$ and $w_{7}/w_{1}=3.0$; (a2,b2) $w_{5}/w_{1}$ ranges from 1.0, 1.4, 1.6, 1.8 to 2.0 at $w_{2}/w_{1}=1.0$ and $w_{7}/w_{1}=3.0$; (a3,b3) $w_{7}/w_{1}$ ranges from 1.0, 2.0, 3.0, 4.0 to 5.0 at $w_{2}/w_{1}=1.0$ and $w_{5}/w_{1}=1.6$.

The effects of $w_{2}/w_{1}$, $w_{5}/w_{1}$ and $w_{7}/w_{1}$ on double emulsion formation behaviours are illustrated in figure 14. The panel (a) and panel (b) series correspond to the two-step and one-step flow rate conditions, and each parameter of concern increases from top to bottom in each subcolumn. The results for the original geometry used in previous sections are marked with an inverted triangle. In figure 14(a1,b1), $w_{2}/w_{1}$ is increased from 0.8, 1.0, 1.2 to 1.4 at $w_{5}/w_{1}=1.6$ and $w_{7}/w_{1}=3.0$. It is seen that the breakup mode of the inner phase apparently changes from the jetting mode to the dripping mode with increasing $w_{2}/w_{1}$ at both flow rates. A larger $w_{2}/w_{1}$ is required to induce the inner breakup mode transition at higher $Ca_{i}$ values. Increasing the side inlet width increases the viscous force of the side-injected fluids to overcome the unaltered interfacial tension force, which leads to the breakup mode transition of the inner phase fluid. Figure 14(a1-ii–a1-iv) and figure 14(b1-i–b1-iii) illustrate the effect of increasing $w_{2}/w_{1}$ on emulsion sizes in the two-step and the one-step formation regimes. The size of the middle part increases in both formation regimes. However, the inner droplet size varies little in the two-step formation regime but decreases in the one-step formation regime.

The effect of the main channel width is studied for $w_{5}/w_{1}=1.0$, 1.4, 1.6, 1.8 and 2.0 at $w_{2}/w_{1}=1.0$ and $w_{7}/w_{1}=3.0$ as shown in figure 14(a2,b2). It is seen that decreasing $w_{5}/w_{1}$ does not affect the formation regime of the inner phase fluid, but it could increase the breakup frequency of the middle phase fluid and induce the decussate regime, as observed in figure 14(a2-i,b2-i). For the flow-focusing geometry, all three inflow fluids converge to the main channel. Thus, narrowing the width of the main channel ($w_{5}$) increases the fluid velocity in the axial central region of channel, which creates a larger velocity gradient in the direction perpendicular to the main flow. During the expansion of the middle phase thread tip, it is subject to a higher shear stress, and as such the middle phase is more likely to break up. With increasing $w_{5}/w_{1}$, the inner part and the entire double emulsion size vary little in the two-step formation regime (see figure 14a2-ii–a2-iv), but they both increase in the one-step formation regime (see figure 14b2-ii–b2-v). It indicates that the main channel width has a more obvious effect on the size of double emulsions generated in the one-step regime. Additionally, emulsions with two inner droplets are regularly obtained in the two-step formation regime at a wider collection channel, i.e. $w_{5}/w_{1}=2.0$ (see figure 14a2-v), similar to those experimentally captured in a double cross-junction device (Deng et al. 2011) and capillary devices (Nabavi et al. 2017b; Levenstein et al. 2016).

At last, the distance between the two side inlets is investigated at $w_{7}/w_{1}=1.0$, 2.0, 3.0, 4.0 and 5.0, $w_{2}/w_{1}=1.0$ and $w_{5}/w_{1}=1.6$, as shown in figure 14(a3,b3). The two-step formation regime shifts to the one-step formation regime at $w_{7}/w_{1}=1.0$ (figure 14a3-i), where the inner phase front reaches the second junction before it breaks up in the dripping mode. However, more generally, the breakup modes and double emulsion sizes vary little with $w_{7}/w_{1}$ in both formation regimes, similar to the findings in binary systems using flow-focusing type geometries (Utada et al. 2007; Wu et al. 2017a). Even though the lengthening of the connection channel increases the flow resistance through it, the flow behaviours inside vary little due to the slightly affected viscous force. As such, the velocities of the inner and middle phases are almost unaffected when they flow into the outer phase junction. Therefore, the overall flow behaviours are almost unchanged. It is noteworthy that satellite droplets appear when $w_{7}/w_{1}\geqslant 4.0$ in the two-step flow regime, due to the highly stretched middle phase thread tip during the emulsification process. It suggests that narrowing the distance between the side inlets could be a possible solution to avoid satellite droplets in producing double emulsions.

5 Conclusions

In this work, a two-dimensional ternary free-energy lattice Boltzmann model is developed and used to systematically study the double emulsion formation behaviours in a planar hierarchical flow-focusing channel under variations of the flow rate, interfacial tension ratio and geometrical settings.

The periodic two-step, one-step and decussate double emulsion formation regimes previously reported in the literature are qualitatively reproduced. A three-dimensional phase diagram is then constructed to show the distribution of each formation regime governed by $Ca_{i}$, $Ca_{m}$ and $Ca_{o}$ values. Depending on the breakup mode of the inner and middle phases, three distinct domains are classified as the periodic two-step, periodic one-step and non-periodic regions. The range for the periodic two-step region is almost unaffected by $Ca_{o}$, and it can be subdivided into seven formation regimes according to the pinch-off locations and the uniqueness of formation frequencies. Among them, periodic double emulsions are produced in regime 1 and 2. In these two regimes, the entire double emulsion size decreases with $Ca_{i}$, increases with $Ca_{m}$, and varies little with $Ca_{o}$. The dripping–threading regime (regime 3) occurs when the middle phase fluid forms a continuous protective layer and carries multiple inner droplets. Decussate regimes (regime 4) with one or two alternate empty droplets are both obtained. Noteworthy, the two empty droplets in the decussate regime could be produced either in a one-by-one sequence, or by breaking an initially formed large empty droplet into two daughter droplets. The bidisperse behaviours in double emulsion size and formation frequency are captured in a certain range of $Ca_{i}$ values in the two-step formation regime. The bidispersity could exist simultaneously for both the inner and middle phase fluids (regime 5), or only occur to the inner phase fluid (regime 6 and 7). In the periodic one-step region for double emulsions (regime 8), the entire double emulsion size is found to decrease with $Ca_{i}$ and $Ca_{o}$, but increases with $Ca_{m}$. Compared to the two-step formation regime, $Ca_{o}$ has a more obvious effect on the size of double emulsions formed in the one-step regime. Based on the one-step data (regime 8), two empirical scaling laws are constructed for the size of the entire double emulsion and the proportion of the inner droplet. The good predictions of both scaling laws justify that the one-step formation process of double emulsions can be analogously viewed as a sum of a blocking period and a squeezing period.

Another contribution of this work is that the presented free-energy model is capable of dealing with a wide range of interfacial tension ratios, and provides accurate results for predicting complete engulfing double emulsions, partial engulfing Janus droplets and non-engulfing separate droplets. In particular, it was necessary to include an additional free-energy term to capture the complete engulfing double emulsions. In the current microfluidic device, a variation in the interfacial tension ratios leads to distinct emulsion morphologies, including the inverse engulfing double emulsions (Chao et al. 2016), non-engulfing single droplets (Nisisako et al. 2005), Janus droplets (Zhang et al. 2017), biconcave and biconvex emulsions (Nie et al. 2006; Nisisako et al. 2015) and even parallel flows.

Regarding channel geometrical parameters, the breakup mode of the inner phase fluid is changed from dripping to jetting by decreasing the side inlet width $w_{2}/w_{1}$, or by narrowing the distance between the two phase side inlets $w_{7}/w_{1}$. This leads to the conversion from the two-step formation regime to the one-step formation regime. The main channel width $w_{5}/w_{1}$ should not be too small in order to avoid the decussate regime. Moreover, narrowing $w_{7}/w_{1}$ is a possible solution to get rid of the satellite droplets for double emulsions generated in the two-step regime. The entire double emulsion size increases with $w_{2}/w_{1}$, but is rarely affected by $w_{5}/w_{1}$ or $w_{7}/w_{1}$ in the two-step formation regime. For the one-step formation regime, the double emulsion size increases with $w_{2}/w_{1}$ and $w_{5}/w_{1}$, but is independent of $w_{7}/w_{1}$.

We would like to point out that the above work is carried out in a two-dimensional scheme. The present ternary free-energy model could be directly extended to three dimensions. The main differences lie in the spatial and velocity discretization schemes, which we have resolved, e.g. in Sadullah, Semprebon & Kusumaatmaja (2018) and Wöhrwag et al. (2018). Based on the fundamental knowledge achieved in the present work, a three-dimensional study is the next step. Indeed, we mostly obtain the jetting regime of the inner phase fluid at $We_{i}\sim 1.0$ with the current 2-D ternary LBM, in contrast to the dripping regime obtained by the reference studies. We believe the main reason for this difference could be attributed to the 3-D effects. Since the dispersed thread grows faster at larger $We_{i}$, a stronger capillary effect is needed to promote droplet formation in the dripping regime (Utada et al. 2005; Fu et al. 2012). It means that the Laplace pressure contribution induced by the out-of-plane curvature plays an important role to promote droplet breakup at larger $We_{i}$ (Hoang et al. 2013), which is absent in 2-D simulations. In addition, wall confinement in the out-of-plane direction can also become important. For example, Azarmanesh et al. (2016) pointed out that the double emulsion formation behaviours in flow-focusing channels are also related to the pressure build-up at the upstream of the inner phase fluid, which is influenced by both the in-plane and out-of-plane wall confinements.

It will be interesting to generalize the scaling laws presented here to three dimensions, and to compare them against experimental observations. We expect the forms of the scaling laws in (4.2) and (4.3) to remain the same but different exponents will be obtained in three dimensions. Furthermore, equal density fluids are used at present. Our newly developed high-density ternary free-energy model (Wöhrwag et al. 2018) could be applied to investigate double emulsion formation behaviours with other fluid types, where density difference is an important factor. It is also worth extending the current ternary free-energy model to deal with multiple emulsions with more components $(N>3)$, or introducing variable interfacial tensions governed by the surfactants (Liu et al. 2018) to study more complex fluid systems.