Modelling double emulsion formation in planar flow-focusing microchannels

TL;DR

This study uses a lattice Boltzmann model to systematically analyze double emulsion formation in flow-focusing microchannels, revealing how flow parameters, interfacial tension, and channel geometry influence droplet morphologies and formation regimes.

Contribution

It introduces a comprehensive 3D formation regime diagram and scaling laws, advancing understanding of emulsion behaviors in microfluidic devices.

Findings

Identified three main formation regimes: periodic two-step, one-step, and non-periodic.

Derived scaling laws for droplet size and size ratio in the one-step regime.

Demonstrated the impact of interfacial tension ratio and channel geometry on emulsion morphology.

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

Tables

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

Regime	Relevant experimental literature	Microfluidic device
1-2	Abate et al. (2011) figure 2	Two cross-junctions
	Kim et al. (2013) figure 3	Glass capillaries
	Nabavi et al. (2015a) figure 2	Glass capillaries
3	Nabavi et al. (2017b) figure 4	Glass capillaries
	Oh et al. (2006) figure 6	Assembled flow-focusing device
4	Kim et al. (2013) figure 2	Glass capillaries
5	Anna et al. (2003) figures 2-3	Single flow-focusing device
	Cubaud & Mason (2008) Sec V-A	Single cross-junction
6	Nabavi et al. (2017b) figure 4	Glass capillaries
	Shang et al. (2014) figure 2	Glass capillaries
	Nisisako et al. (2005) figure 2	Two T-junctions
8	Abate et al. (2011) figure 2	Two cross-junctions
	Utada et al. (2005) figure 1	Glass capillaries
9-10	Kim et al. (2013) figure 2	Glass capillaries
11	Wu & Gong (2013) figures 3	Assembled T-Y junctions

Equations

$$\mathcal{H}={\int }_{\Omega }{c}_s^2\rho \ln \rho dV+\mathcal{F}.$$

$$\mathcal{F}$$

$${\sigma }_{rg}$$

$${\sigma }_{rb}$$

$${\sigma }_{gb}$$

$$\varphi =C_1-C_2,\psi =C_3,$$

$${\partial }_t\varphi +\nabla \cdot (\varphi \mathbf{v})$$

$${\partial }_t\psi +\nabla \cdot (\psi \mathbf{v})$$

$${\partial }_t\rho +\nabla \cdot (\rho \mathbf{v})$$

$${\partial }_t(\rho \mathbf{v})+\nabla \cdot (\rho \mathbf{v}\mathbf{v})$$

$$E_{a}(C_{m})=\left\{\begin{array}[]{ll}\beta C_{m}^{2},&C_{m}<0\\[2.0pt] 0,&0\leq C_{m}\leq 1\\[2.0pt] \beta(C_{m}-1)^{2},&C_{m}>1,\end{array}\right.$$

$$f_i(\mathbf{r}+{\mathbf{e}}_i{\delta }_t,t+{\delta }_t)$$

$$g_i(\mathbf{r}+{\mathbf{e}}_i{\delta }_t,t+{\delta }_t)$$

$$h_i(\mathbf{r}+{\mathbf{e}}_i{\delta }_t,t+{\delta }_t)$$

$$\displaystyle M_{\phi}\ =\ \Gamma_{\phi}(\tau_{g}-\frac{\delta t}{2}),\$$

$$\displaystyle M_{\psi}\ =\ \Gamma_{\psi}(\tau_{h}-\frac{\delta t}{2}),\$$

$$f_i^{eq}$$

$$g_i^{eq}$$

$$h_i^{eq}$$

$$\rho =\sum _i{f}_i,\rho \mathbf{v}=\sum _i{f}_i{\mathbf{e}}_i+\frac{\mathbf{F}\delta t}{2},\varphi =\sum _i{g}_i,\psi =\sum _i{h}_i.$$

$$\displaystyle\left.{g_{i}}\vphantom{\Big{|}}\right|_{i=1,5,8}$$

$$\displaystyle\left.{h_{i}}\vphantom{\Big{|}}\right|_{i=1,5,8}$$

$${\chi }_i(L_x,y,t+\delta t)=\frac{{\chi }_i(L_x,y,t)+\zeta (L_x-1,y,t){\chi }_i(L_x-1,y,t+\delta t)}{1+\zeta (L_x-1,y,t)}.$$

$${\chi }^{\prime }(L_x+1,y,t+\delta t)=\frac{{\chi }^{\prime }(L_x+1,y,t)+\zeta (L_x,y,t){\chi }^{\prime }(L_x,y,t+\delta t)}{1+\zeta (L_x,y,t)}.$$

$$u_x(y)=\frac{-4u_{max}(y-y_1)(y-y_2)}{(y_2-y_1{)}^2},1\le y\le L_y-1,$$

$$C_1(x,y)$$

$$C_2(x,y)$$

$$C_3(x,y)$$

$$\frac{L_p}{w_1}=(\widetilde{ϵ}+\tilde{\gamma }\frac{Q_{dispersed}}{Q_{continuous}})C{a}_o^{\tilde{m}},$$

$$\frac{R_{emulsion}}{w_1}=(0.270+0.0526\frac{Q_m}{Q_i})C{a}_o^{-0.268},$$

$$\frac{R_{inner}}{R_{emulsion}}=0.904(\frac{Q_i}{Q_i+Q_m}{)}^{0.609}C{a}_m^{-0.060}C{a}_o^{0.030}.$$

$$\left\{\begin{array}[]{llll}C_{1}<0:&\mu_{\rho}=\beta\left(\frac{\rho+\phi-\psi}{2}\right),&\mu_{\phi}=\beta\left(\frac{\rho+\phi-\psi}{2}\right),&\mu_{\psi}=-\beta\left(\frac{\rho+\phi-\psi}{2}\right);\\[2.0pt] C_{2}<0:&\mu_{\rho}=\beta\left(\frac{\rho-\phi-\psi}{2}\right),&\mu_{\phi}=-\beta\left(\frac{\rho-\phi-\psi}{2}\right),&\mu_{\psi}=-\beta\left(\frac{\rho-\phi-\psi}{2}\right);\\[2.0pt] C_{3}<0:&\mu_{\rho}=0,&\mu_{\phi}=0,&\mu_{\psi}=2\beta\psi;\\[2.0pt] \end{array}\right.$$

$$\left\{\begin{array}[]{llll}C_{1}>1:&\mu_{\rho}=\beta\left(\frac{\rho+\phi-\psi}{2}-1\right),&\mu_{\phi}=\beta\left(\frac{\rho+\phi-\psi}{2}-1\right),&\mu_{\psi}=-\beta\left(\frac{\rho+\phi-\psi}{2}-1\right);\\[2.0pt] C_{2}>1:&\mu_{\rho}=\beta\left(\frac{\rho-\phi-\psi}{2}-1\right),&\mu_{\phi}=-\beta\left(\frac{\rho-\phi-\psi}{2}-1\right),&\mu_{\psi}=-\beta\left(\frac{\rho-\phi-\psi}{2}-1\right);\\[2.0pt] C_{3}>1:&\mu_{\rho}=0,&\mu_{\phi}=0,&\mu_{\psi}=2\beta(\psi-1);\\[2.0pt] \end{array}\right.$$

Full text

1 Introduction

Double emulsions are droplets with one other droplet inside. Their core-shell structure has attracted wide attentions in various fields (Vladisavljević et al., 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 technique 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 technique 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 (Azarmanesh et al., 2016; Lee & Weitz, 2008), 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 (Shang et al., 2017; Zhang et al., 2016).

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 et al., 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 multi-layer (Abate & Weitz, 2009) or multi-core emulsions (Nisisako et al., 2005; Nabavi et al., 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 (Pannacci et al., 2008; Abate et al., 2011; Seo et al., 2007), co-axial 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 analyzed 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 multi-core emulsions when they are formed using a series of T-junctions. The one-step formation regime is mostly encountered in co-axial 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 et al., 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 co-axial 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 behaviors, including using the volume of fluid (VOF) 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 (Guzowski et al., 2012; Pannacci et al., 2008; Chao et al., 2016; Nisisako & Hatsuzawa, 2010).

We have chosen to use the lattice Boltzmann method (LBM). LBM is highly favorable for the study of emulsion formation behaviors due to its simplicity in solving interface dynamics, including droplet break-ups 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 (Semprebon et al., 2016; Liang et al., 2016; Wöhrwag et al., 2018; Abadi et al., 2018), color-fluid model (Leclaire et al., 2013a, b; Yu et al., 2019b; Fu et al., 2016, 2017), 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 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 behaviors in planar hierarchical flow-focusing junctions. We focus on the two-dimensional (2D) case to reduce the computational time needed for parametric studies. The major physical difference in the flow dynamics between the 2D and the three-dimensional (3D) 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 2D system and a systematic 2D 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

[TABLE]

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

[TABLE]

It is constructed using a double-well potential form for the bulk free energy $E_{0}(C_{m})=(\kappa_{m}/2)C_{m}^{2}(1-C_{m})^{2}$ and a square gradient term for the interface region $E_{interface}(\nabla C_{m})=(\alpha^{2}\kappa_{m}/2)(\nabla C_{m})^{2}$. $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 $\rho_{m}$, which we have set to be 1.0 for simplicity. Thus the total density is related to the concentration fractions by defining $\rho=C_{1}\rho_{1}+C_{2}\rho_{2}+C_{3}\rho_{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 here could be denoted by $\mathsfbi{C}$ = $[C_{1}\quad C_{2}\quad C_{3}]$ = $[1\quad 0\quad 0]$, $[0\quad 1\quad 0]$ and $[0\quad 0\quad 1]$, respectively. $\alpha$ is the interface width parameter. The interfacial tension between red-green phases $\sigma_{rg}$, red-blue phases $\sigma_{rb}$, and green-blue phases $\sigma_{gb}$ are related to $\kappa$’s through

[TABLE]

To capture the interface dynamics, two order parameters $\phi$ and $\psi$ are introduced as

[TABLE]

and the original concentration fields can be reversely obtained from $\rho$, $\phi$ and $\psi$ via $C_{1}=(\rho+\phi-\psi)/2$, $C_{2}=(\rho-\phi-\psi)/2$, and $C_{3}=\psi$. 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:

[TABLE]

Here, $\boldsymbol{v}$ is the fluid velocity and $\eta$ is the dynamic viscosity. $M_{\phi}$ and $M_{\psi}$ are the mobility values for $\phi$ and $\psi$. The thermodynamic properties are related to the equations of motion via the chemical potential $\mu_{\phi}$, $\mu_{\psi}$ and the pressure tensor $\mathsfbi{P}$. The chemical potential is defined as the variational derivative of the free energy $\mu_{q}=\delta\mathcal{F}/\delta q$, where $q=\rho,\phi$ or $\psi$. The pressure tensor term in Eq. (8) is constructed by $\nabla\boldsymbol{\cdot}\mathsfbi{P}=\nabla(\rho c_{s}^{2})+\nabla\boldsymbol{\cdot}\mathsfbi{P_{0}}$ and $\nabla\boldsymbol{\cdot}\mathsfbi{P_{0}}=(\rho\nabla\mu_{\rho}+\phi\nabla\mu_{\phi}+\psi\nabla\mu_{\psi})$. The first term $\nabla(\rho 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}=-\nabla\boldsymbol{\cdot}\mathsfbi{P_{0}}$ term is treated as an external force term in the lattice Boltzmann implementation. The explicit expressions of $\mu_{\rho}$, $\mu_{\phi}$, $\mu_{\psi}$, and $\mathsfbi{P}$ 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 $\sigma_{rb}>\sigma_{gb}+\sigma_{rg}$. Given the relation between the $\sigma$’s and $\kappa_{m}$’s in Eq.(3), one can easily find that $\kappa_{2}$ should be negative while $\kappa_{1}$ and $\kappa_{3}$ are positive. In the free energy model, the negative $\kappa_{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 $\kappa_{m}$’s 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 LB 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 $\kappa_{m}$, here we introduce an additional free energy term given by

[TABLE]

where $\beta$ 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 Eq. (9), additional terms should be included in the chemical potentials accordingly, which are listed in Appendix A.

2.2 Lattice Boltzmann method

To solve Eq. (5)-(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 $\phi$ and $\psi$, 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

[TABLE]

Here, the force term is implemented through the exact difference method (Kupershtokh et al., 2009; Mazloomi et al., 2015), which is expressed as the last two terms enclosed in brackets in Eq. (10), with $\widetilde{\boldsymbol{u}}=\sum_{i}f_{i}\boldsymbol{e}_{i}/\rho$, i.e., the velocity without the force term, and $\delta\widetilde{\boldsymbol{u}}=\boldsymbol{F}\delta t/\rho$. The lattice time step $\delta_{t}$ is set to be 1.0. $\tau_{f}$ is the relaxation parameter given by $1/\tau_{f}=C_{1}/\tau_{1}+C_{2}/\tau_{2}+C_{3}/\tau_{3}$, where $\tau_{1,2,3}$ are related to the viscosity of each fluid by $\tau_{1,2,3}=3\eta_{r,g,b}/\rho+1/2$, respectively (Krüger et al., 2017). $\tau_{g}$ and $\tau_{h}$ are the relaxation parameters that are related to the mobility parameters $M_{\phi}$ and $M_{\psi}$ in the Cahn-Hilliard equations through

[TABLE]

where $\Gamma_{\phi}$ and $\Gamma_{\psi}$ are two tunable parameters. The mobility values $M_{\phi}$ and $M_{\psi}$ are relevant for the timescale 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 insure the reasonable damping of the convective term (Lim & Lam, 2013; Jacqmin, 1999). 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_{\phi}\kappa_{2})$. The absolute values of $Pe_{m}$ used is generally on the order of $o(10)-o(80)$, which is on similar magnitudes to those used in previous two-phase droplet behavior studies (Shardt et al., 2014; Zhou et al., 2010; Menech, 2006). Moreover, $M_{\psi}=M_{\phi}/3$ is considered to assign symmetric mobility for each concentration component (Semprebon et al., 2016).

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

[TABLE]

where the weight coefficient $\omega_{i}$ are given by $\omega_{0}=4/9$, $\omega_{1-4}=1/9$ and $\omega_{5-8}=1/36$. The macroscopic variables are related to the distribution functions through

[TABLE]

2.3 Boundary conditions

The boundary conditions involved in the present study contain: no-slip boundary, wetting boundary and the inlet-outlet boundary. No-slip boundary condition is used on the solid walls, which is realized by the half-way 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 LB ternary color-fluid model. Similarly for the free energy model used here, the macroscopic values of $\rho$, $\phi$ and $\psi$ 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 $\phi$ and $\psi$ are $\phi_{in}$ and $\psi_{in}$, respectively. The sum of the unknown distribution functions can be solved according to Eq. (22), and then $g_{1,5,8}$ and $h_{1,5,8}$ are allocated by their weight factors as

[TABLE]

For the outlet boundary, the convective boundary condition (CBC) (Lou et al., 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 $\chi_{i}$ = $f_{i}$, $g_{i}$ and $h_{i}$ at the outlet layer ($x=L_{x}$),

[TABLE]

The other is for the macroscopic quantities, such as $\chi^{\prime}$ = $\rho$, $\phi$, $\psi$ and $\mathsfbi{P_{0}}$ 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,

[TABLE]

Here, $\zeta$ is the characteristic velocity normal to the outlet boundary. For simplicity, we have explicitly computed $\chi_{i}$ and $\chi^{\prime}$ through $\zeta$ at time $t$. Three common choices for $\zeta$ in CBC’s 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

[TABLE]

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 right 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 Eq. (27) 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 centered 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 $\sigma$ 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 Eq. (27). 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}=X_{d}/D$ and $t^{*}=tu_{max}/D$, where $D$ is the droplet diameter. The $X^{*}_{d}$ 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}$ occurs when the droplet completely moves out of the channel.

It is seen in figure 2 (b1-b3) that the $X^{*}_{d}$ 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 CBC’s has little effect on the flow behaviors away from the outlet. Deviations in $X^{*}_{d}$ curves from the reference lines occur at around $t^{*}$ = 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}$ increases as $u_{max}$ decreases for the cases with CBC-AV and CBC-MV. When the $u_{max}$ is on 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 et al., 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.

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 $\sigma_{gb}/\sigma_{rg}$ and $\sigma_{rb}/\sigma_{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 $\sigma_{rb}/\sigma_{rg}>1+\sigma_{gb}/\sigma_{rg}$; (I-B) complete engulfing with the green droplet inside the red one for $\sigma_{gb}/\sigma_{rg}>1+\sigma_{rb}/\sigma_{rg}$; (II) non-engulfing, for $\sigma_{rb}/\sigma_{rg}+\sigma_{gb}/\sigma_{rg}<1$, where the red and green droplets tend to separate from each other; (III) partial engulfing (Janus droplet), for $|(\sigma_{rb}/\sigma_{rg})-(\sigma_{gb}/\sigma_{rg})|<1$ and $\sigma_{rb}/\sigma_{rg}+\sigma_{gb}/\sigma_{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)

[TABLE]

Periodic boundary conditions are used for all boundaries. To reproduce all the possible morphologies, simulations are performed at various groups of ($\sigma_{gb}/\sigma_{rg}$, $\sigma_{rb}/\sigma_{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); and (III) partial engulfing emulsions: (1.0, 1.0), (1.0, 1.5), (1.0, 0.5), (0.5, 1.0), (1.5, 1.0) and (100, 100). The interfacial tension $\sigma_{rg}$ is fixed at 0.005 except for the case with ($\sigma_{gb}/\sigma_{rg}$, $\sigma_{rb}/\sigma_{rg}$) = (100, 100), where $\sigma_{rg}$ = 0.00001 is used to reach the high interfacial tension ratio. The value of the coefficient $\beta$ 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}=\sigma_{jk}-\sigma_{ij}-\sigma_{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 $\beta$ in the additional free energy term introduced in Eq. (9), varying $\beta$ = 0, 0.0001, 0.001, 0.01, 0.1 and 1.0 for one typical double emulsion morphology at ($\sigma_{gb}/\sigma_{rg}$, $\sigma_{rb}/\sigma_{rg}$) = (1.4, 0.35). The obtained result at $\beta=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 $\beta=0.0001$ in figure 3 (c). For $\beta$ varying from 0.001 to 0.1, the complete engulfing morphology could be successfully reproduced and invisible difference is observed for different values of $\beta$. However, further increasing $\beta$ to 1.0 induces numerical instability, which indicates that the $\beta$ 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, $\beta=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 setup 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}$. 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.

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}=\rho_{i}u_{i}^{2}w_{1}/\sigma_{im}$; Capillary numbers (the ratios of viscous force to interfacial tension force) of the middle and outer fluids $Ca_{m}=\eta_{m}u_{m}/\sigma_{im}$, $Ca_{o}=\eta_{o}u_{o}/\sigma_{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 $\lambda_{im}=\eta_{i}/\eta_{m}$, $\lambda_{om}=\eta_{o}/\eta_{m}$; and interfacial tension ratios $\sigma_{io}/\sigma_{im}$ and $\sigma_{mo}/\sigma_{im}$. Here, $u$ is the average inlet velocity. However, in this work, we focus on double emulsion formation behaviors in the limit of small inertia (Nabavi et al., 2017b). As such, it is more appropriate to use $Ca_{i}=\eta_{i}u_{i}/\sigma_{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 $\lambda_{im}=\lambda_{om}=1.24$, and the interfacial tension ratio is fixed at $(\sigma_{mo}/\sigma_{im},\sigma_{io}/\sigma_{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 5 (b1)). 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 5 (b2)), while the decussate regime with one empty middle phase droplet is achieved by decreasing the inner flow rate to $Ca_{i}=0.008$ (figure 5 (b3)). 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 on 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 behaviors are still obtained. This suggests that, while $We_{i}$ can affect the resulting formation regimes, the rich flow behaviors with many different formation regimes are still present in the limit of small inertia (Wu et al., 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.

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$) is 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) and (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 behaviors.

To analyze 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 side of the red solid line, the inner phase breaks up in the jetting mode. The right 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, and the formation behavior 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 literatures 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 behaviors 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 analyzed in detail to help deepen the understanding of double emulsion formation behaviors.

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 behavior 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 behaviors 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; and (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.

(1) Size variations of double emulsions generated in the two-step regime

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 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 (a-i)-(a-III) shows the effect of $Ca_{i}$ on the two-step double emulsion formation behaviors. 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 co-axial 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 formaiton behaviors 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 behaviors. 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.

(2) Dripping-threading regime

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

(3) Bidisperse formation regime

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 behavior, the temporal evolution of the inner and middle thread tip locations are traced as denoted by $X_{i}$ and $X_{m}$ in figure 8 (a1) and (a2), respectively. The time $t$ and locations $X_{i,m}$ are normalized using $t^{*}=t(u_{m})_{max}/w_{1}$ and $X^{*}_{i,m}=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 neighboring 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).

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) and (a2) suggests that such bidisperse behaviors 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 et al., 2005; Utada et al., 2007; Garstecki et al., 2005).

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 easily 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 8 (b1)), but not in the middle phase fluid (figure 8 (b2)). 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 behavior could still happen, and it is named as the Bidisperse-Threading regime (Regime 7).

(4) Decussate regime with two empty droplets

Decussate regimes occupy a substantial proportion in the two-step formation region of figure 6. Among them, decussate regimes with one alternate empty droplet is 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 decusste regime with two empty droplets, where the two empty droplets are produced one by one. Our results show that by lowering $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.

Comparing figure 9 (a) and (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 distribution character 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.

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 sub-figure, the interface shapes are depicted by the solid lines. The leading viscous force component is displayed in the upper part, i.e., $\tau_{xy}=\eta(\p u_{y}/\p x+\p u_{x}/\p y)$ for a two-dimensional system, and it is normalized using $\tau^{*}_{xy}=\tau_{xy}w_{1}/\sigma_{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 figure 10 (a1) to (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 scalling laws of the present ternary system. To complete the scaling law involving the ternary flow characters, we first analyze the effects of governing flow parameter on the double emulsion sizes.

(1) Size variations of double emulsions generated in the one-step regime

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^{*}=A/(\pi 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.

The effect of $Ca_{m}$ on the size of each component and the corresponding snapshots are illustrated in figure 11 (b1) and (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 monotonously with $Ca_{m}$. These trends qualitatively agree with the size variation characters 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) and (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 monotonously. 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}$.

(2) Scaling laws for character sizes of the double emulsion

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; Tan et al., 2008; Christopher 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

[TABLE]

where the plug length is normalized by the inlet width $w_{1}$, and $\tilde{\epsilon}$, $\tilde{\gamma}$ 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 Eq. (31), 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 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 Eq. (31). 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 Eq. (31) for the ternary system. Firstly, the droplet length is only suitable for plug shape 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 et al., 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}/\pi}$ is used as the dependant 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 behaviors more complex. Based on the analysis of flow parameter effects shown in figure 11, $Q_{dispersed}/Q_{continuous}$ in Eq. (31) 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

[TABLE]

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 Eq. (32) 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 Eq. (32) can provide a general guidance for predicting double emulsion size.

Another size of interest is the ratio between the equivalent inner droplet radius $R_{inner}=\sqrt{A_{inner}/\pi}$ 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 co-axial 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 co-axial 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 fine tune the entire size.

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

[TABLE]

The way to obtain the values of the coefficients in Eq. (33) is the same to that used in Eq. (32). 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 Eq. (33) 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 Eq. (33) 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., $(\sigma_{mo}/\sigma_{im},\sigma_{io}/\sigma_{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, $\sigma_{im}$ is fixed at 0.005 except for the case at $(\sigma_{mo}/\sigma_{im},\sigma_{io}/\sigma_{im})$ = (100, 100), where $\sigma_{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 (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 $(\sigma_{mo}/\sigma_{im},\sigma_{io}/\sigma_{im})$ = (1.0, 2.2) (figure 13 (a1) and (b1)), the inverse engulfed double emulsion is captured in figure 13 (a2) and (b2) by reversing the interfacial tension ratios to $(\sigma_{mo}/\sigma_{im},\sigma_{io}/\sigma_{im})$ = (2.2, 1.0). With the inverse interfacial tension ratios, the inner phase fluid is more favored 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.

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 behaviors 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 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 co-axial glass capillaries, which could prevent the early formation of the three-phase contact line, and hence delay the interfacial tension effect.

At $(\sigma_{mo}/\sigma_{im},\sigma_{io}/\sigma_{im})$ = (0.48, 0.48), the red and green droplets tend to separate with 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 multi-core 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 figure 13 (a4-a6) and (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), (b4), and (b5), biconcave and biconvex emulsions are formed downstream. These structures are analogous to those experimentally fabricated by Nisisako et al. (2015) and Nie et al. (2006). Finally, for $(\sigma_{mo}/\sigma_{im},\sigma_{io}/\sigma_{im})$ = (100, 100) shown in figure 13 (a6) and (b6), $\sigma_{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 $\sigma_{im}$ easily leads to the parallel layered flow behaviors 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.

The effects of $w_{2}/w_{1}$, $w_{5}/w_{1}$ and $w_{7}/w_{1}$ on double emulsion formation behaviors are illustrated in figure 14. The (a) and (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 sub-column. The results for the original geometry used in previous sections are marked with an inverted triangle. In figure 14 (a1) and (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) and (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 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 14 (a2-ii)-(a2-iv)), but they both increase in the one-step formation regime (see figure 14 (b2-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 14 (a2-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) and (b3). The two-step formation regime shifts to the one-step formation regime at $w_{7}/w_{1}$ = 1.0 (figure 14 (a3-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 (Wu et al., 2017a; Utada et al., 2007). Even though the lengthening of the connection channel increases the flow resistance through it, the flow behaviors 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 behaviors are almost unchanged. Noteworthy, satellite droplets appear at $w_{7}/w_{1}\geq 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 behaviors 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}$. 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 behaviors 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 (Regime5), 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 (Nisisako et al., 2015; Nie et al., 2006), 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 et al. (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 2D 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 3D 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 2D 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 behaviors in flow-focusing channels are also related to the pressure buildup 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 Eqs.(32) and (33) 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 behaviors 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.

Acknowledgments

H. Liu, C. Zhang and N. Wang acknowledge financial supports from the National Key Research and Development Project of China (No. 2016YFB0200902), and the National Natural Science Foundation of China (Nos. 51876170, 51506168 and 51711530130). H. Kusumaatmaja acknowledges funding from EPSRC (No. EP/P007139/1). C. Semprebon acknowledges support from Northumbria University through the Vice-Chancellor’s Fellowship Programme, and funding from EPSRC (No. EP/S036857/1). N. Wang was supported by the China Scholarship Council for one year (2017-2018) at Durham University, UK.

Declaration of Interests

The authors report no conflict of interest.

Appendix A

The expressions for the additional terms in the chemical potentials due to the additional energy term in Eq. (9) are provided below:

[TABLE]

[TABLE]