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. (Reference Semprebon, Krüger and Kusumaatmaja2016). 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 Reference Briant and Yeomans2004; Zhang Reference Zhang2011). 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. (Reference Semprebon, Krüger and Kusumaatmaja2016).

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. (Reference Guzowski, Korczyk, Jakiela and Garstecki2012), 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 (Reference Lee and Liu2010) and Abadi et al. (Reference Abadi, Fakhari and Rahimian2018)). 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 Reference Kupershtokh, Medvedev and Karpov2009; Mazloomi, Chikatamarla & Karlin Reference Mazloomi, Chikatamarla and Karlin2015), 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. Reference Krüger, Kusumaatmaja, Kuzmin, Shardt, Silva and Viggen2017); $\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 Reference Jacqmin1999; Lim & Lam Reference Lim and Lam2013). 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 Reference Menech2006; Zhou et al. Reference Zhou, Yue, Feng, Ollivier-Gooch and Hu2010; Shardt, Mitra & Derksen Reference Shardt, Mitra and Derksen2014). Moreover, $M_{\unicode[STIX]{x1D713}}=M_{\unicode[STIX]{x1D719}}/3$ is considered to assign symmetric mobility for each concentration component (Semprebon et al. Reference Semprebon, Krüger and Kusumaatmaja2016).

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 Reference Ladd1994). The solid walls should have a preferential affinity with the continuous phase fluid to generate stable droplets/emulsions (Abate et al. Reference Abate, Thiele and Weitz2011). Fu et al. (Reference Fu, Zhao, Bai, Jin and Cheng2016) 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 Reference Zou and He1997) 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 (Reference Hao and Cheng2009) and Liu & Zhang (Reference Liu and Zhang2011) 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 Reference Lou, Guo and Shi2013; Chen & Deng Reference Chen and Deng2017) 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. Reference Lou, Guo and Shi2013).

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 Reference Li, Jia and Liu2017). 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. Reference Pannacci, Bruus, Bartolo, Etchart, Lockhart, Hennequin, Willaime and Tabeling2008; Guzowski et al. Reference Guzowski, Korczyk, Jakiela and Garstecki2012), 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. (Reference Guzowski, Korczyk, Jakiela and Garstecki2012), 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. Reference Yu, Liu, Liang and Zhang2019b)

(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. (Reference Pannacci, Bruus, Bartolo, Etchart, Lockhart, Hennequin, Willaime and Tabeling2008) 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.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 Reference Zou and He1997) 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. Reference Abate, Thiele and Weitz2011; Azarmanesh et al. Reference Azarmanesh, Farhadi and Azizian2016): 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. Reference Nabavi, Vladisavljević and Manović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. Reference Abate, Thiele and Weitz2011; Azarmanesh et al. Reference Azarmanesh, Farhadi and Azizian2016) 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. Reference Abate, Thiele and Weitz2011; Azarmanesh et al. Reference Azarmanesh, Farhadi and Azizian2016). 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 Reference Wu, Cao and Sundén2017b). 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. Reference Abate, Thiele and Weitz2011) (a1,a2) and simulation (Azarmanesh et al. Reference Azarmanesh, Farhadi and Azizian2016) (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. Reference Utada, Lorenceau, Link, Kaplan, Stone and Weitz2005). 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. Reference Fu, Wu, Ma and Li2012; Cubaud & Mason Reference Cubaud and Mason2008; Yu et al. Reference Yu, Liu, Zhao and Chen2019a). 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. (Reference Fu, Zhao, Bai, Jin and Cheng2016) 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. (Reference Nabavi, Vladisavljević and Manović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 Reference Cubaud and Mason2008). 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. Reference Oh, Kim, Baek, Seong and Lee2006). 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 Reference Coullet, Mahadevan and Riera2005; Garstecki, Fuerstman & Whitesides Reference Garstecki, Fuerstman and Whitesides2005; Utada et al. Reference Utada, Fernandez-Nieves, Stone and Weitz2007).

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. (Reference Azarmanesh, Farhadi and Azizian2016) 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. (Reference Azarmanesh, Farhadi and Azizian2016) 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. Reference Azarmanesh, Farhadi and Azizian2016). 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. Reference Nisisako, Okushima and Torii2005).

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. Reference Abate, Thiele and Weitz2011). 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. (Reference Kim, Kim, Kim, Han and Weitz2013) 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. Reference Utada, Fernandez-Nieves, Stone and Weitz2007; Fu et al. Reference Fu, Wu, Ma and Li2012). 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 Reference Cubaud and Mason2008; Liu & Zhang Reference Liu and Zhang2011; Fu et al. Reference Fu, Wu, Ma and Li2012). 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. Reference Utada, Fernandez-Nieves, Stone and Weitz2007). 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. Reference Garstecki, Fuerstman, Stone and Whitesides2006; Christopher et al. Reference Christopher, Noharuddin, Taylor and Anna2008; Tan et al. Reference Tan, Xu, Li and Luo2008; Liu & Zhang Reference Liu and Zhang2011). Specifically, Liu & Zhang (Reference Liu and Zhang2011) 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. (Reference Christopher, Noharuddin, Taylor and Anna2008). Moreover, the work of Liu & Zhang (Reference Liu and Zhang2011) 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 Reference Liu and Zhang2011), 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. Reference Garstecki, Fuerstman, Stone and Whitesides2006; Liu & Zhang Reference Liu and Zhang2011). The volume values are more general quantities for the ellipsoid-like double emulsions (Steegmans, Schron & Boom Reference Steegmans, Schron and Boom2009; Chang et al. Reference Chang, Serra, Bouquey, Prat and Hadziioannou2009; Fu et al. Reference Fu, Zhao, Bai, Jin and Cheng2016). 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. (Reference Chang, Serra, Bouquey, Prat and Hadziioannou2009) 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. (Reference Fu, Zhao, Bai, Jin and Cheng2016) 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. (Reference Chang, Serra, Bouquey, Prat and Hadziioannou2009) and Fu et al. (Reference Fu, Zhao, Bai, Jin and Cheng2016). 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. (Reference Chao, Mak and Shum2016) 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. (Reference Chao, Mak and Shum2016). This could be explained by the experimental findings of Pannacci et al. (Reference Pannacci, Bruus, Bartolo, Etchart, Lockhart, Hennequin, Willaime and Tabeling2008). 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. (Reference Chao, Mak and Shum2016), 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. Reference Nisisako, Okushima and Torii2005). 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. (Reference Zhang, Ge, Geng, Luo, Chen and Xu2017) 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 (Reference Nisisako, Ando and Hatsuzawa2015) and Nie et al. (Reference Nie, Li, Seo, Xu and Kumacheva2006). 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 Reference Liu and Zhang2011; Nabavi et al. Reference Nabavi, Vladisavljević, Gu and Ekanem2015b; Wu et al. Reference Wu, Liu, Zhao and Chen2017a). 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. Reference Deng, Meng, Xie, Ju, Mou, Wang and Chu2011) and capillary devices (Nabavi et al. Reference Nabavi, Vladisavljević and Manović2017b; Levenstein et al. Reference Levenstein, Bawazer, Nally, Marchant, Gong, Meldrum and Kapur2016).

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. Reference Utada, Fernandez-Nieves, Stone and Weitz2007; Wu et al. Reference Wu, Liu, Zhao and Chen2017a). 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.