*Article* **Dynamic Behaviours of Monodisperse Double Emulsion Formation in a Tri-Axial Capillary Device**

**Yuchen Dai , Haotian Cha, Nhat-Khuong Nguyen, Lingxi Ouyang, Fariba Galogahi, Ajeet Singh Yadav, Hongjie An , Jun Zhang , Chin Hong Ooi and Nam-Trung Nguyen \***

> Queensland Micro-Nanotechnology Centre, Griffith University, Nathan, QLD 4111, Australia **\*** Correspondence: nam-trung.nguyen@griffith.edu.au

**Abstract:** We investigated experimentally, analytically, and numerically the formation process of double emulsion formations under a dripping regime in a tri-axial co-flow capillary device. The results show that mismatches of core and shell droplets under a given flow condition can be captured both experimentally and numerically. We propose a semi-analytical model using the match ratio between the pinch-off length of the shell droplet and the product of the core growth rate and its pinchoff time. The mismatch issue can be avoided if the match ratio is lower than unity. We considered a model with the wall effect to predict the size of the matched double emulsion. The model shows slight deviations with experimental data if the Reynolds number of the continuous phase is lower than 0.06 but asymptotically approaches good agreement if the Reynolds number increases from 0.06 to 0.14. The numerical simulation generally agrees with the experiments under various flow conditions.

**Keywords:** double emulsion; core–shell droplet; microfluidics; tri-axial capillary; computational fluid dynamics; dripping regime

#### **1. Introduction**

Double emulsions core–shell droplets are referred to as dispersed droplets containing a smaller one as the core. Double emulsions are highly desired for applications in drug delivery, food science, controlled release of substances, etc., because they can encapsulate a cargo in the core. Due to its significance, various designs for the formation of double emulsions have been extensively studied. Double emulsions are performed either simultaneously or successively and categorised accordingly into one-step and two-step processes. The former is simple and robust [1], while the latter is more efficient and controllable [2]. The one-step approach is the focus of our present work. One-step formation of double emulsions has been achieved with microfluidic devices. The working principle relies on the so-called Rayleigh-Plateau instability, where a liquid jet becomes unstable with its wavelength larger than the circumference. The jet consequently breaks up into segments with minimum surface area, forming a spherical shape [3,4]. Utada et al. [5] connected a glass cylindrical capillary with a square glass one and forced all liquid phases through an exit downstream. Varying the dimension of the outlet orifice can adjust the size of the droplets. The team also employed the breakup mechanism of a single emulsion to explain the size distribution of the double ones and defined an effective capillary number to classify the formation process as dripping and jetting regimes. Nie et al. [6] developed a co-flow microfluidic device in polyurethane elastomer to form double emulsions. The sizes of the core and the outer droplets were tuned with the flow rates of the three liquid phases. In general, double emulsion droplets with diameters from ten to hundreds of micrometres can be formed [1,5,7,8]. However, larger core–shell droplets with diameters ranging from hundreds of micrometres to several millimetres are also required for applications such as cell manipulations [9,10], macrocapsules [11], and microreactors [12]. Producing double emulsion droplets in this size range using typical microfluidic devices requires a significantly low flow rate of the outer phase, usually resulting in a drop accumulation issue

**Citation:** Dai, Y.; Cha, H.; Nguyen, N.-K.; Ouyang, L.; Galogahi, F.; Yadav, A.S.; An, H.; Zhang, J.; Ooi, C.H.; Nguyen, N.-T. Dynamic Behaviours of Monodisperse Double Emulsion Formation in a Tri-Axial Capillary Device. *Micromachines* **2022**, *13*, 1877. https://doi.org/10.3390/ mi13111877

Academic Editor: Pingan Zhu

Received: 13 October 2022 Accepted: 28 October 2022 Published: 31 October 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

because the gravity or buoyancy effect dominates the inertial drag. To prevent this problem, capillary devices were oriented vertically and scaled up to the so-called "millifluidics" [13]. Note that the term "microfluidics" is still valid for droplet smaller than a millimetre even if the device dimension is larger than that [14,15]. Shao et al. [15] experimentally investigated double emulsions in a dual-coaxial capillary with a 2.2 mm diameter tube as the outer channel. Adjusting the relative position of the ends of the inner and middle capillaries can switch the operation between two-step and one-step process. Their results demonstrated that the dynamic effect of the inner phase on the size of shell droplet is insignificant in the two-step set up, but drastically affecting it in the one-step process. Additionally, the team found that the mismatch in formation frequencies of core and shell droplets in the two-step process can be avoided with the one-step device. Schmit et al. [16] utilised a pendant drop method to produce millimetre-size double emulsions in a capillary tube with a 4 mm outer diameter. Their experimental results quantitatively demonstrated the effects of the inner and middle phase flow rates on the number as well as the size of the core in each droplet.

Though it is ultimately essential to obtain experimental data on double emulsions, numerical simulations are also of great importance for understanding and optimising the operation parameters for droplet formations. The velocity and pressure fields can be readily visualised, which are not always accessible with experimental methods. Zhou et al. [17] simulated double emulsions in a capillary device using a diffuse-interface framework. The calculation domain was simplified as a 2-dimensional axisymmetric geometry. The interfacial thickness and position were determined by a phase-field variable. By applying a similar method, Park and Anderson [18] successfully predicted the dripping and jetting regimes as well as their transition for the capillary device reported by Utada et al. [5]. Vu et al. [19] numerically investigated double emulsions in a capillary device using the front-tracking method proposed by Tryggvason et al. [20]. The results showed that the jetting mode can be promoted by increasing the Reynolds number (*Re*) and the Webber number (*We*). Fu et al. [21] established a model based on the ternary Lattice Boltzmann method to simulate the one-step process. Herrada et al. [22] carried out numerical simulations to validate their linear stability analysis on double emulsions in a capillary device. The team adopted the volume-of-fluid (VOF) method to simulate multiphases and interfaces. The effects of viscosity ratio, flow rate ratio, and interfacial tension ratio on the one-step process were also quantitatively demonstrated using the VOF method [23,24]. Azarmanesh et al. [25] simulated both one-step and two-step processes using the VOF method and found that the Capillary number (*Ca*) of the inner phase can significantly affect the outer droplet size in the two-step process, while the *Ca* of the outer phase plays a significant role in one-step counterpart. More recently, Yang et al. [26] numerically investigated the deformation behaviours of one-step double emulsions using VOF method. Examining the streamlines in the channel, the team discovered that the core moving forward relative to the shell was caused by the large vortex passing through the core–shell interface.

The dripping regime is well known to produce a monodisperse droplet size distribution while the jetting regime leads to a polydisperse counterpart. For most applications, a monodisperse emulsion is preferred. As a result, a simple force balance might be conducted to elucidate the formation mechanism and to predict the size distribution for various flow conditions. In general, the forces exerted on a double emulsion are the same as those on a single emulsion. These forces are kinetic, drag, interfacial tension, buoyancy, and Laplace pressure forces [27]. Note that the Laplace pressure was commonly neglected in symmetric and axisymmetric flows but has to be considered in a T-junction configuration due to the significant difference between the head and tail curvatures [28,29]. In addition, the drag force may have various forms while the others remain the same. The deviations were mainly caused by the different wall effect corrections when extending the Stokes drag of unbounded flows to confined counterparts [30,31]. Adjusting the continuous phase flow rate can tune the size of the droplets as its changes the drag force. Temperature can affect the surface tension and the viscosity and thus can be used to tune the droplet size [32,33]. Furthermore, to connect the model of the single droplet with the core–shell

droplet, an equivalence of droplets volume ratio and input flow rates ratio was usually introduced [13,34]. Note that this hypothesis is valid if each droplet contains a single core.

Although efforts have been dedicated to double emulsion formation with microfluidics, most reported works were mainly based on either experiments or CFD simulation. Only a few can provide a relatively simple analytical model. Often, results from experiments, theory, and numerical simulation are inconsistent, preventing effective optimisation of devices and operation parameters for the generation of double emulsion. To our best knowledge, no past studies concurrently examined experiments, CFD simulations, and analytical models under the same condition. We propose in the present study a simple experimental setup, a vertically oriented tri-axial co-flow capillary device, to generate core–shell droplets in a one-step process. We investigate the mechanism of the double emulsion formation both experimentally and numerically. Furthermore, we develop and validate a simple but effective analytical model for predicting the droplet size distribution with the wall effect and according to various flow conditions.

#### **2. Experimental Materials and Methods**

#### *2.1. Design and Fabrication of the Capillary Device*

The tri-axial needles were purchased from Ramé-Hart instrument Co. (Succasunna, NJ, USA). The device structure and dimensions are detailed in Figure 1A and Table 1, with an inner needle tip extension ∆*z* = 0.1 mm. A straight circular channel was connected with the outlet of the needles. The outlet channel was made of poly-dimethylsiloxane (PDMS) moulded on a straight wire with the same diameter as that of the needles. The inlets of the needles are individually connected with three syringes via tubing. Figure 1B,C show the schematic diagram and the images of the experimental setup. *Micromachines* **2022**, *13*, x FOR PEER REVIEW 4 of 15

**Figure 1.** The capillary device used in this study: (**A**) the tri-axial needles; (**B**) the schematic diagram; (**C**) the photo of the experiment setup. **Figure 1.** The capillary device used in this study: (**A**) the tri-axial needles; (**B**) the schematic diagram; (**C**) the photo of the experiment setup.

**Needle/Phase Inner Diameter (µm) Outer Diameter (µm) Thickness (µm)** 1 178 356 89

Fluorinated oil (HFE, Novec 7500 3 M, Merck, Darmstadt, Germany) was used as the core dispersed phase, whose density and dynamic viscosity are 1.61 g/mL and 1.31 mPa·s, respectively [8]. The shell dispersed phase was a polymer consisting of 0.06 g ethyl-4(dimethylamino) benzoate (Merck, Darmstadt, Germany), 0.05 g camphorquinone (Merck), and 10 g trimethylolpropane trimethacrylate (TMPTMA, Merck). The corresponding density and dynamic viscosity are 1.07 g/mL and 42 mPa·s, respectively [35]. The continuous phase was 80% glycerol dissolved in 20% DI water, with 0.1% Tween 20, and its density and dynamic viscosity are 1.21 g/mL and 75.42 mPa·s, respectively [36]. The laboratory

The interfacial tension σ between each pair of two liquid phases were measured using the reverse pendant drop method with an optical tensiometer (Theta Flex from Biolin Scientific, Gothenburg, Sweden). Figure 2A reports the mean contact angles of each pair of two phases, while Figure 2B,C show the images of the pendant drops. The interfacial tensions of the inner and outer interfaces were measured for 10 s for each case with 33 frames per second, and the measured time-averaged mean interfacial coefficients are

**Table 1.** The dimensions of the tri-axial needles.

temperature was maintained at 23 °C.

3.45 and 8.48 mN/m, respectively.

*2.2. Materials*


**Table 1.** The dimensions of the tri-axial needles.

#### *2.2. Materials*

Fluorinated oil (HFE, Novec 7500 3 M, Merck, Darmstadt, Germany) was used as the core dispersed phase, whose density and dynamic viscosity are 1.61 g/mL and 1.31 mPa·s, respectively [8]. The shell dispersed phase was a polymer consisting of 0.06 g ethyl-4(dimethylamino) benzoate (Merck, Darmstadt, Germany), 0.05 g camphorquinone (Merck), and 10 g trimethylolpropane trimethacrylate (TMPTMA, Merck). The corresponding density and dynamic viscosity are 1.07 g/mL and 42 mPa·s, respectively [35]. The continuous phase was 80% glycerol dissolved in 20% DI water, with 0.1% Tween 20, and its density and dynamic viscosity are 1.21 g/mL and 75.42 mPa·s, respectively [36]. The laboratory temperature was maintained at 23 ◦C.

The interfacial tension σ*ij* between each pair of two liquid phases were measured using the reverse pendant drop method with an optical tensiometer (Theta Flex from Biolin Scientific, Gothenburg, Sweden). Figure 2A reports the mean contact angles of each pair of two phases, while Figure 2B,C show the images of the pendant drops. The interfacial tensions of the inner and outer interfaces were measured for 10 s for each case with 33 frames per second, and the measured time-averaged mean interfacial coefficients are 3.45 and 8.48 mN/m, respectively. *Micromachines* **2022**, *13*, x FOR PEER REVIEW 5 of 15

**Figure 2.** The interfacial tension coefficients measurement: (**A**) the mean contact angles (error bars represent the standard deviations and *n* = 7); (**B**) middle-phase reversed pendant drop in inner phase; (**C**) Middle phase reversed pendant drop in outer phase. **Figure 2.** The interfacial tension coefficients measurement: (**A**) the mean contact angles (error bars represent the standard deviations and *n* = 7); (**B**) middle-phase reversed pendant drop in inner phase; (**C**) Middle phase reversed pendant drop in outer phase.

#### *2.3. Experimental Setup and Data Analysis. 2.3. Experimental Setup and Data Analysis*

analyse the captured videos.

and momentum equations are expressed as:

(⃗)

above are related to the volume fraction of each phase

every single computational cell were calculated as follows:

**3. Numerical Methodology** *3.1. Governing Equations*

as:

the constraint:

The three phases were delivered into the capillary device at specific flow rates using syringe pumps (MNT-SPM-100, Welland, Austrlia). The formation and detachment of the droplets from the needle tips under various flow conditions were recorded using a high-speed camera (Photron FASTCAM SA3, San Diego, CA, USA, attached with a 25 mm F2.8 Ultra Macro 2.5-5.0X lens) at 250 frames per second. The open-source software ImageJ 1.53m (National Institutes of Health USA, Bethesda, MD, USA) was adopted to The three phases were delivered into the capillary device at specific flow rates using syringe pumps (MNT-SPM-100, Welland, Austrlia). The formation and detachment of the droplets from the needle tips under various flow conditions were recorded using a high-speed camera (Photron FASTCAM SA3, San Diego, CA, USA, attached with a 25 mm F2.8 Ultra Macro 2.5-5.0X lens) at 250 frames per second. The open-source software ImageJ 1.53 m (National Institutes of Health USA, Bethesda, MD, USA) was adopted to analyse the captured videos.

<sup>+</sup> <sup>∇</sup> <sup>∙</sup> (⃗) <sup>=</sup> <sup>0</sup> (1)

<sup>+</sup> <sup>∇</sup> <sup>∙</sup> (⃗) <sup>=</sup> <sup>0</sup> (4)

)] + ⃗ + ⃗ (2)

. The density and viscosity in

(3)

where *ρ*, *t*, *v*, *p*, *µ*, , and denote density, time, velocity, pressure, dynamic viscosity, gravitational acceleration, and a source term, respectively. The Navier–Stokes equations

3

=1

The volume of fluid (VOF) method was adopted to model the fluid–fluid interfaces

The sum of volume fractions of three phases within each computational cell follows

(, ) = ∑

( , )

<sup>+</sup> <sup>∇</sup> <sup>∙</sup> (⃗⃗) <sup>=</sup> −∇p <sup>+</sup> <sup>∇</sup> <sup>∙</sup> [(∇⃗ <sup>+</sup> ∇⃗

All phases were modelled as incompressible fluids and laminar flows. The continuity

#### **3. Numerical Methodology**

#### *3.1. Governing Equations*

All phases were modelled as incompressible fluids and laminar flows. The continuity and momentum equations are expressed as:

$$\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \vec{v}) = 0 \tag{1}$$

$$\frac{\partial(\rho \stackrel{\rightarrow}{\boldsymbol{\upsilon}})}{\partial t} + \nabla \cdot (\rho \stackrel{\rightarrow}{\boldsymbol{\upsilon}} \stackrel{\rightarrow}{\boldsymbol{\upsilon}}) = -\nabla \mathbf{p} + \nabla \cdot [\mu (\nabla \stackrel{\rightarrow}{\boldsymbol{\upsilon}} + \nabla \stackrel{\rightarrow}{\boldsymbol{\upsilon}}^T)] + \rho \stackrel{\rightarrow}{\boldsymbol{\varsigma}} + \stackrel{\rightarrow}{F} \tag{2}$$

where *ρ*, *t*, *v*, *p*, *µ*, *g*, and *F* denote density, time, velocity, pressure, dynamic viscosity, gravitational acceleration, and a source term, respectively. The Navier–Stokes equations above are related to the volume fraction of each phase *α<sup>i</sup>* . The density and viscosity in every single computational cell were calculated as follows:

$$\sigma(\rho\_\prime \,\mu) = \sum\_{i=1}^3 a\_i(\rho\_{i\prime} \,\mu\_i) \tag{3}$$

The volume of fluid (VOF) method was adopted to model the fluid–fluid interfaces as:

$$\frac{\partial \boldsymbol{\alpha}\_{i}}{\partial t} + \nabla \cdot (\boldsymbol{\alpha}\_{i} \stackrel{\rightarrow}{\boldsymbol{\upsilon}}) = \mathbf{0} \tag{4}$$

The sum of volume fractions of three phases within each computational cell follows the constraint:

$$\sum\_{i=1}^{3} \alpha\_i = 1\tag{5}$$

The continuous surface force (CSF) model developed by Brackbill et al. [37] was used to model the surface tension via the source term → *F* in the momentum equation as:

$$\vec{F} = \sum\_{pairs(i,j),\ i$$

where *κ* stands for the interface curvature determined by:

$$\kappa\_{\dot{i}} = \nabla \cdot \mathfrak{h}\_{\dot{i}} \tag{7}$$

and ˆ*n* is defined as:

$$
\hat{\mathfrak{m}}\_{i} = \frac{\nabla \mathfrak{a}\_{i}}{|\nabla \mathfrak{a}\_{i}|} \tag{8}
$$

#### *3.2. Computational Domain and Boundary Conditions*

As the droplets are formed at the needle outlet, the corresponding circular cross-section channel downstream was selected to be the computational domain. The centreline was regarded as the axis due to the axisymmetric geometry. Boundary conditions at the entrance, exit, and walls were set to be velocity inlet, pressure outlet, and no-slip wall, respectively, Figure 3. The contact angles between each interface and the needle wall were determined based on experimental results in Figure 2B,C because the tri-axial needle and the hook needle used in the interfacial tension measurement are both made of stainless steel. The hydrodynamic entrance length for each needle inlet was estimated by *lent*,*<sup>i</sup>* = 0.0575*ReiD<sup>i</sup>* assuming a laminar pipe flow [38]. The corresponding order of magnitude in this study was calculated as *lent*,1 ∼ *O*(1 µm), *lent*,2 ∼ *O*(0.1 µm), *lent*,3 ∼ *O*(10 µm), respectively. To ensure each flow inlet is fully developed before reaching the expansion region, the inlet lengths for inner, middle, and outer needles in the simulation were fixed, respectively as 300, 200, and 200 µm, much longer than the hydrodynamic entrance lengths. Note that the

inner needle inlet length was determined by considering its extension length. Structural meshes were generated, and a grid independence analysis was conducted with three cell sizes: 15 µm, 10 µm, and 5 µm with corresponding grid numbers of 23,652 (coarse), 51,785 (medium), and 204,428 (fine), respectively. The results showed that the medium and fine meshes, under the same input conditions, can model the interfaces with insignificant differences. Since the fine mesh drastically increases the computational resources, the medium one was selected for this study. *Micromachines* **2022**, *13*, x FOR PEER REVIEW 7 of 15

**Figure 3.** Schematic illustration of the computational domain, structural mesh, and forces exerting on a growing double emulsion droplet. **Figure 3.** Schematic illustration of the computational domain, structural mesh, and forces exerting on a growing double emulsion droplet.

We performed simulations using the commercially available CFD software package ANSYS 2022. The SIMPLE algorithm was utilised for the pressure-velocity coupling. A PRESTO scheme was selected to calculate the exact value of the pressure term. We implemented the Geo-Reconstruct algorithm to acquire the interface interpolation. The governing equations were discretised by a second-order implicit scheme for the spatial terms and transient formulation. The time step was determined by the Courant–Fredrichs–Lewy number lower than 0.25. We performed simulations using the commercially available CFD software package ANSYS 2022. The SIMPLE algorithm was utilised for the pressure-velocity coupling. A PRESTO scheme was selected to calculate the exact value of the pressure term. We implemented the Geo-Reconstruct algorithm to acquire the interface interpolation. The governing equations were discretised by a second-order implicit scheme for the spatial terms and transient formulation. The time step was determined by the Courant–Fredrichs– Lewy number lower than 0.25.

#### **4. Analytical Model 4. Analytical Model**

regime (*Ca*<sup>3</sup> < 1, where <sup>3</sup> =

Both experiments and CFD simulations can capture the droplet formation in detail, but building empirical relations based on these methods usually requires relatively intensive case studies. Since we focus on the dripping regime, a simpler analytical model is more convenient for elucidating the droplet breakup mechanism and approximately predicting the droplet size. In the dripping regime, several forces are exerted on the growing droplets along the axial direction: viscous drag and kinetic forces pull them downstream while surface tension and buoyancy forces act as opposite counterparts, Both experiments and CFD simulations can capture the droplet formation in detail, but building empirical relations based on these methods usually requires relatively intensive case studies. Since we focus on the dripping regime, a simpler analytical model is more convenient for elucidating the droplet breakup mechanism and approximately predicting the droplet size. In the dripping regime, several forces are exerted on the growing droplets along the axial direction: viscous drag *F<sup>D</sup>* and kinetic *F<sup>K</sup>* forces pull them downstream while surface tension *F<sup>σ</sup>* and buoyancy *F<sup>B</sup>* forces act as opposite counterparts, Figure 3.

Figure 3. The assumptions for the analytical model are (i) the continuous flow is in the Stokes regime (*Re*<sup>3</sup> << 1, where <sup>3</sup> = 3ℎ̅<sup>3</sup> 3 ); (ii) the droplet formation is in the dripping The assumptions for the analytical model are (i) the continuous flow is in the Stokes regime (*Re*<sup>3</sup> << 1, where *Re*<sup>3</sup> = *ρ*3*rchannelv*<sup>3</sup> *µ*3 ); (ii) the droplet formation is in the dripping regime (*Ca*<sup>3</sup> < 1, where *Ca*<sup>3</sup> = *µ*3*v*<sup>3</sup> *σ*23 ); (iii) the shape of droplets is spherical; (iv) laminar and

and incompressible flow are also presumed as per ansatzes in CFD simulations. The forces

= 11<sup>1</sup> + 22<sup>2</sup>

); (iii) the shape of droplets is spherical; (iv) laminar

= 1σ<sup>12</sup> + 2σ<sup>23</sup> (10)

(9)

3̅<sup>3</sup> <sup>23</sup>

incompressible flow are also presumed as per ansatzes in CFD simulations. The forces involved in the formation process are [13,27,39]:

$$F\_K = \rho\_1 Q\_1 v\_1 + \rho\_2 Q\_2 v\_2 \tag{9}$$

$$F\_{\sigma} = \pi d\_{o1} \sigma\_{12} + \pi d\_{o2} \sigma\_{23} \tag{10}$$

$$F\_B = \frac{4}{3}\pi \left(\frac{d\_{drop}}{2}\right)^3 \rho \text{g} \tag{11}$$

where *Q* and *d* represent volume flow rate input and diameter, respectively; subscripts *o* and *drop* denote the outer and the whole droplet. As for the Stokes drag force, we employed the model with the wall effect considered since the diameter of the droplets formed in our device has a similar magnitude as that of the channel [31,40]:

$$F\_D = 3\pi\mu\_3 d\_{drop} \left(\overline{v}\_3 K\_1 - \overline{v}\_{drop} K\_2\right) \tag{12}$$

where *v*<sup>3</sup> and *vdrop* are the mean velocities of the continuous phase and the whole droplet:

$$
\overline{v}\_3 = \frac{4Q\_3}{\pi d\_{channel}^2} \tag{13}
$$

$$F\_D = 3\pi\mu\_3 d\_{drop} \left(\overline{v}\_3 K\_1 - \overline{v}\_{drop} K\_2\right) \tag{14}$$

with *K*<sup>1</sup> and *K*<sup>2</sup> being the wall effect correction coefficients as:

$$\begin{array}{c} \text{K}\_{1} = 1/\left[1 - 2.10443\lambda + 2.08877\lambda^{3} - 0.94813\lambda^{5} - 1.372\lambda^{6} + 3.87\lambda^{8} \\ \quad - 4.19\lambda^{10} + \mathcal{O}(\lambda^{11})\right] \end{array} \tag{15}$$

$$K\_2 = K\_1 \left[ 1 - \frac{2}{3} \lambda - 0.1628 \lambda^3 - 0.4059 \lambda^7 + 0.5236 \lambda^9 + 1.51 \lambda^{10} + \mathcal{O}(\lambda^{11}) \right] \tag{16}$$

where *λ* is the dimensionless droplet diameter as follows:

$$
\lambda = \frac{d\_{drop}}{d\_{channel}}\tag{17}
$$

Note that Equations (15) and (16) have been proven to show good agreements with the exact solution with *λ* ≤ 0.8 [31,40]. For larger *λ*, we may need higher order expansions or the exact solution of Haberman and Sayre [41], which is too long to be expressed here.

By scaling *ddrop* ∼ *do*2, we can obtain the orders of magnitude:

$$F\_D \sim \mathcal{O}(10^{-5} \,\text{N}), \; F\_\text{K} \sim \mathcal{O}(10^{-9} \,\text{N}), \; F\_\sigma \sim \mathcal{O}(10^{-5} \,\text{N}), \; F\_\text{B} \sim \mathcal{O}(10^{-7} \,\text{N})\tag{18}$$

Thus, the force balance equation in the current study can be reasonably simplified as:

$$F\_D = F\_\sigma \tag{19}$$

Now the diameter of the whole droplet *ddrop* can be obtained by numerically solving the algebraic system above. To estimate the size of the core, we employed [13,34]:

$$\left(\frac{d\_{core}}{2}\right)^3 = \frac{Q\_1}{Q\_1 + Q\_2} \left(\frac{d\_{drop}}{2}\right)^3 \tag{20}$$

#### **5. Results and Discussion**

We first examined the formation process with the continuous phase flow rate *Q*<sup>3</sup> of 250, 500, and 1000 µL/min, while the inner and middle phase flow rates *Q*<sup>1</sup> and *Q*<sup>2</sup> were fixed at 3 µL/min and 30 µL/min, respectively. Figure 4 presents the double emulsion formation

in the case of *Q*<sup>1</sup> = 3 µL/min, *Q*<sup>2</sup> = 30 µL/min, *Q*<sup>3</sup> = 250 µL/min over a formation period. The Experimental (top row) and the corresponding simulation results (bottom row) are obtained at nine time steps, with an excellent agreement in the droplet formation pattern. However, the total formation time showed a 16% relative deviation. We observed that the core droplet breaks up before the shell does. A satellite droplet also forms after the break-up of the shell droplet. Over the entire formation process, double emulsions are formed in a regular pattern, with each shell droplet containing a core. In addition, as the whole droplet flows downstream, the core droplet moves downwards quicker than the shell one due to the higher density. This phenomenon is slightly overpredicted in the simulation; see the second image in Figure 4. This might also be caused by the fact that the needles are not precisely tri-axial due to the manufacturing error, evidently by the left skewness in the experimental results. *Micromachines* **2022**, *13*, x FOR PEER REVIEW 9 of 15 whole droplet flows downstream, the core droplet moves downwards quicker than the shell one due to the higher density. This phenomenon is slightly overpredicted in the simulation; see the second image in Figure 4. This might also be caused by the fact that the needles are not precisely tri-axial due to the manufacturing error, evidently by the left skewness in the experimental results.

**Figure 4.** Droplet formation with regular double emulsions at *Q*<sup>1</sup> = 3 µL/min, *Q*<sup>2</sup> = 30 µL/min, *Q*<sup>3</sup> = 250 µL/min. **Figure 4.** Droplet formation with regular double emulsions at *Q*<sup>1</sup> = 3 µL/min, *Q*<sup>2</sup> = 30 µL/min, *Q*<sup>3</sup> = 250 µL/min.

Furthermore, the other two cases with a continuous phase flow rate *Q*3 of 500 and 1000 µL/min, while the inner *Q*1 and middle *Q*2 phase flow rates were fixed at 3 µL/min and 30 µL/min, respectively, Figure 5. Figure 5A shows that with *Q*<sup>3</sup> = 500 µL/min, double and single emulsions occur alternately, meaning that only one of every two droplets contains a core. Although the time period still shows a 14% deviation between the experiment and simulation, the flow patterns and droplet sizes were captured with remarkable agreement. Worse conditions can be found with *Q*<sup>3</sup> is 1000 µL/min, Figure 5B, where only one in three droplets contains a core in the experiment, and one in four droplets contains a core in the simulation. Now the flow patterns show a difference within a period. Since it is necessary to avoid the mismatch in applications, the deviation on the mismatched frequencies is not discussed in detail. These unmatched phenomena can be explained by the mismatch of pinch-off locations. The pinch-off location of the middle phase is stretched downwards due to the increased capillary number of the continuous phase, while the pinch-off location of the inner phase remains near the needle tip. Additionally, we also observed the negligible difference between the size of the double emulsion and the single droplet in both experiments and simulations, meaning that the inner phase flow rate only has a slight effect on the whole droplet size. This observation agrees with previous studies [24,26]. Furthermore, the other two cases with a continuous phase flow rate *Q*<sup>3</sup> of 500 and 1000 µL/min, while the inner *Q*<sup>1</sup> and middle *Q*<sup>2</sup> phase flow rates were fixed at 3 µL/min and 30 µL/min, respectively, Figure 5. Figure 5A shows that with *Q*<sup>3</sup> = 500 µL/min, double and single emulsions occur alternately, meaning that only one of every two droplets contains a core. Although the time period still shows a 14% deviation between the experiment and simulation, the flow patterns and droplet sizes were captured with remarkable agreement. Worse conditions can be found with *Q*<sup>3</sup> is 1000 µL/min, Figure 5B, where only one in three droplets contains a core in the experiment, and one in four droplets contains a core in the simulation. Now the flow patterns show a difference within a period. Since it is necessary to avoid the mismatch in applications, the deviation on the mismatched frequencies is not discussed in detail. These unmatched phenomena can be explained by the mismatch of pinch-off locations. The pinch-off location of the middle phase is stretched downwards due to the increased capillary number of the continuous phase, while the pinch-off location of the inner phase remains near the needle tip. Additionally, we also observed the negligible difference between the size of the double emulsion and the single droplet in both experiments and simulations, meaning that the inner phase flow rate only has a slight effect on the whole droplet size. This observation agrees with previous studies [24,26].

**Figure 5.** Droplet formation with mismatch at *Q*<sup>1</sup> = 3 µL/min, *Q*<sup>2</sup> = 30 µL/min and: (**A**) *Q*<sup>3</sup> = 500 µL/min; (**B**) *Q*<sup>3</sup> = 1000 µL/min. **Figure 5.** Droplet formation with mismatch at *Q*<sup>1</sup> = 3 µL/min, *Q*<sup>2</sup> = 30 µL/min and: (**A**) *Q*<sup>3</sup> = 500 µL/min; (**B**) *Q*<sup>3</sup> = 1000 µL/min.

The mismatch phenomenon can be further illustrated along with the velocity streamlines in the case with *Q*<sup>1</sup> = 3 µL/min, *Q*<sup>2</sup> = 30 µL/min, *Q*<sup>3</sup> = 500 µL/min in Figure 6. The blue lines represent the interface between the inner and middle phases, while the red ones stand for the interface between the middle and outer phases. Due to the abrupt expansion of the channel, the velocity profile drastically changes in the vicinity of the needle tips and then develops downstream. Two pairs of vortices, which in fact are two vortex rings due to the axisymmetric geometry, can be observed near the needle tips as a result of the sudden change in the velocity profile, Figure 6A,B. The lower vortex ring formed inside the shell droplet pushes the growing core droplet head backward. The upper vortex ring also acts as an obstacle due to the reverse flow direction along the centreline. As both core and shell droplets grow, Figure 6A–C, the size of the lower vortex ring gradually decreases to null by the shear force of the continuous phase. In contrast, the upper vortex ring increases since the core droplet does not grow beyond the pinch-off location of the shell droplet. After the breakup of the shell droplet, the upper vortex ring crosses the interface between the inner and middle phases, Figure 6D, and is then squeezed by the middle-outer interface and the inner-middle interface. This process forms a new vortex ring from the inner The mismatch phenomenon can be further illustrated along with the velocity streamlines in the case with *Q*<sup>1</sup> = 3 µL/min, *Q*<sup>2</sup> = 30 µL/min, *Q*<sup>3</sup> = 500 µL/min in Figure 6. The blue lines represent the interface between the inner and middle phases, while the red ones stand for the interface between the middle and outer phases. Due to the abrupt expansion of the channel, the velocity profile drastically changes in the vicinity of the needle tips and then develops downstream. Two pairs of vortices, which in fact are two vortex rings due to the axisymmetric geometry, can be observed near the needle tips as a result of the sudden change in the velocity profile, Figure 6A,B. The lower vortex ring formed inside the shell droplet pushes the growing core droplet head backward. The upper vortex ring also acts as an obstacle due to the reverse flow direction along the centreline. As both core and shell droplets grow, Figure 6A–C, the size of the lower vortex ring gradually decreases to null by the shear force of the continuous phase. In contrast, the upper vortex ring increases since the core droplet does not grow beyond the pinch-off location of the shell droplet. After the breakup of the shell droplet, the upper vortex ring crosses the interface between the inner and middle phases, Figure 6D, and is then squeezed by the middle-outer interface and the inner-middle interface. This process forms a new vortex ring from the inner needle

needle tip, Figure 6E. After the mismatch, the growing core droplet accumulated inside

where ∆1

velocity estimated as:

tip, Figure 6E. After the mismatch, the growing core droplet accumulated inside the next shell one due to the fixed flow rate input, Figure 6F. Again, as both droplets grow, the size of the lower vortex ring decreases, but the size of the upper one decreases as well in this period, Figure 6F–H. This is because the accumulated growing core droplet is larger, and the vortex ring inside does not attach to the inner needle tip even before the core breaks up. Consequently, the vortex ring is stretched longer and narrower. Thus, more inertia assists the core droplet heading downstream through the space between the vortex ring and the inner-middle interface, Figure 6H. Finally, the core droplet grows beyond the pinch-off location of the shell droplet and breaks up along with the outer emulsion, Figure 6I,J. The whole process then repeats periodically. A similar explanation can be used for the mismatch in the case with *Q*<sup>1</sup> = 3 µL/min, *Q*<sup>2</sup> = 30 µL/min, *Q*<sup>3</sup> = 1000 µL/min. The difference is that the shell droplet pinches off two or three times before the core droplet accumulates to grow beyond the pinch-off position. in this period, Figure 6F–H. This is because the accumulated growing core droplet is larger, and the vortex ring inside does not attach to the inner needle tip even before the core breaks up. Consequently, the vortex ring is stretched longer and narrower. Thus, more inertia assists the core droplet heading downstream through the space between the vortex ring and the inner-middle interface, Figure 6H. Finally, the core droplet grows beyond the pinch-off location of the shell droplet and breaks up along with the outer emulsion, Figure 6I,J. The whole process then repeats periodically. A similar explanation can be used for the mismatch in the case with *Q*<sup>1</sup> = 3 µL/min, *Q*<sup>2</sup> = 30 µL/min, *Q*<sup>3</sup> = 1000 µL/min. The difference is that the shell droplet pinches off two or three times before the core droplet accumulates to grow beyond the pinch-off position.

*Micromachines* **2022**, *13*, x FOR PEER REVIEW 11 of 15

**Figure 6.** Velocity streamlines near the needle tip at *Q*<sup>1</sup> = 3 µL/min, *Q*<sup>2</sup> = 30 µL/min, *Q*<sup>3</sup> = 500 µL/min: (**A**–**E**) when the mismatch occurs; (**F**–**J**) when the core matches with the shell droplet. **Figure 6.** Velocity streamlines near the needle tip at *Q*<sup>1</sup> = 3 µL/min, *Q*<sup>2</sup> = 30 µL/min, *Q*<sup>3</sup> = 500 µL/min: (**A**–**E**) when the mismatch occurs; (**F**–**J**) when the core matches with the shell droplet.

To solve the mismatch issue, the core droplet needs to grow beyond the pinch-off location of the shell droplet within each pinch-off period, or 2−∆1 1 ≤ 2 . Hence, we define a match ratio of: To solve the mismatch issue, the core droplet needs to grow beyond the pinch-off location of the shell droplet within each pinch-off period, or *<sup>l</sup>p*2−∆*ztip*<sup>1</sup> *vg*<sup>1</sup> ≤ *tp*2. Hence, we define a match ratio of:

$$\mathfrak{F} = \frac{l\_{p2} - \Delta z\_{tip1}}{v\_{g1} t\_{p2}} \tag{21}$$

1 2 is the inner needle tip elevation difference; 1 is the core droplet growth where ∆*ztip*<sup>1</sup> is the inner needle tip elevation difference; *vg*<sup>1</sup> is the core droplet growth velocity estimated as:

$$v\_{\mathcal{S}^1} \approx \frac{Q\_1}{\pi \left(\frac{d\_{\mathcal{S}^1}}{2}\right)^2} \tag{22}$$

1 ≈ ( 1 2 ) <sup>2</sup> (22) and the pinch-off time of the shell *tp*<sup>2</sup> calculated based on the theoretical model as:

1

2

<sup>1</sup> + <sup>2</sup>

$$t\_{p2} = \frac{\frac{4}{3}\pi \left(\frac{d\_{drop}}{2}\right)^3}{Q\_1 + Q\_2} \tag{23}$$

(23)

, cannot

4 3 (

2 =

=

However, to our best knowledge, the pinch-off length of the shell droplet, *lp*2, cannot be approximated through a simple model. Thus, we measured it and the pinch-off time experimentally based on the averaged value of the upper and lower pinch-off lengths, Figure 7A. Error bars were calculated from standard deviations of the mean (*n* = 3). We examined *Q*<sup>3</sup> ranging from 200 to 1200 µL/min with a 100-µL/min interval with the inner phase flow rate *Q*<sup>1</sup> of 3, 6, 12 µL/min, and the middle phase flow rate *Q*<sup>2</sup> of 30, 15, 7.5 µL/min, respectively. Note that the corresponding continuous phase Reynolds number, *Re*3, and Capillary number, *Ca*3, range from 0.02 to 0.14 and from 0.01 to 0.09, respectively. Figure 7B,C show the measured averaged pinch-off length and time against the continuous phase flow rate at different flow rates of the inner and middle phase, respectively. The standard deviations of the mean were calculated from three samples for each case [42].

Figure 8 **Figure 7.** Pinch-off parameters: (**A**) pinch-off length measurement; (**B**) pinch-off length versus the continuous phase flow rate; (**C**) pinch-off time versus the continuous phase flow rate; (**D**) match ratio versus the continuous phase Reynolds number.

Figure 8 With increasing *Q*3, the pinch-off length and time of the shell droplet generally increase and decrease, respectively. We also noticed that, when a mismatch occurs, the standard deviation of the pinch-off length drastically increases, while the counterpart of the pinch-off time shows a relatively stable pattern. Additionally, when the sum of *Q*<sup>1</sup> and *Q*<sup>2</sup> increases, the pinch-off time for each case decreases and gradually approaches closely as *Q*<sup>3</sup> increases. In contrast, as the sum of *Q*<sup>1</sup> and *Q*<sup>2</sup> increases, the pinch-off length shows a reduction pattern before mismatches occur but a relatively chaotic counterpart due to the mismatches, leading to abrupt increases for different cases. The mismatch issue does not show up in the case of *Q*<sup>1</sup> = 12 µL/min and *Q*<sup>2</sup> = 30 µL/min within the *Q*<sup>3</sup> range in this study, and thus the corresponding line has a negligible change. Equations (21)–(23) indicate that the mismatch issue can be theoretically avoided once the match ratio, *ξ*, is lower than unity. Increasing *Q*<sup>1</sup> or decreasing *Q*<sup>2</sup> can reduce *ξ*. Increasing ∆*ztip*<sup>1</sup> can also decrease *ξ*, but this option was not examined due to the limitation of the tri-axial needles used in this study. Figure 7D shows the operation map with the match ratio versus the continuous phase Reynolds number for various Reynolds numbers of inner and middle phases. The operation map shows the occurrence of mismatches in experiments. We found a very close threshold at *ξ* = 1, beyond which the mismatch occurs. Additionally, the match ratio

1

generally decreases with an increase in the inner-middle phase flow rate ratio. However, even at the same inner-middle phase flow rate ratio *Q*1/*Q*<sup>2</sup> = 0.4, the one with a higher total flow rate shows a better performance than the other. Hence, it would be better to increase the inner phase flow rate rather than decrease the middle one.

After resolving the mismatch issue by increasing the inner phase flow rate, we compared the size distribution among experiments, simulations, and theory. The diameters of the droplets in the experiment and simulation were determined through ImageJ by measuring the vertical and horizontal lengths of at least three whole droplets for each case and then averaged as *d<sup>i</sup>* = *ai*+*b<sup>i</sup>* 2 , as shown in Figure 8A. Error bars were again calculated from standard deviations of the mean. Figure 8B shows the effect of the continuous phase Reynolds number on the matched droplet size normalized by the channel diameter. As observed, the analytical result on the whole droplet shows a slight deviation from the experimental one when *Re*<sup>3</sup> ranges from 0.02 to 0.06 and asymptotically approaches good agreement with *Re*<sup>3</sup> increasing from 0.06 to 0.14. The analytical result of the core droplet diameter shows a better prediction against the experimental data within the range of *Re*<sup>3</sup> in our current study. This can be explained by the approximated wall effect correction model since increasing *Re*<sup>3</sup> decreases the dimensionless droplet diameter *λ*, leading to less influence of the wall effect; see Equations (15) and (16). Furthermore, the reduction in the whole droplet diameter also reduces the deformation effect, which results in a better agreement with the assumption of a spherical droplet. In contrast, the CFD simulation results generally overlap with experimental data, indicating that the current numerical model is valid to capture the double emulsions in the dripping regime through the tri-axial capillary device. Figure 8

Figure 8 **Figure 8.** Matched droplets: (**A**) core and shell droplet diameter measurement; (**B**) size distribution comparisons among experiments, simulations, and theory.

#### **6. Conclusions**

1 We experimentally, numerically, and theoretically investigated the formation process of double emulsion in a vertically arranged tri-axial co-flow capillary device. The volume of fluid (VOF) method and continuous surface force (CSF) model were used for the numerical simulation to catch the details of fluid flow. An analytical model based on a simple force balance was built to estimate the droplet size distribution in accordance with the continuous phase Reynolds number. We found that mismatches of core and shell droplets under certain flow conditions can be captured both experimentally and numerically. To quantify the mismatch, we proposed a semi-theoretical model by matching the pinchoff length of the shell droplets with the product of the growth rate of the core and the pinch-off time of the shell. The mismatch issue is expected to be avoided if the match ratio is lower than unity, which was validated with experimental data. Regarding the reduction in the match ratio, we found that increasing the inner phase flow rate shows better performance than reducing the middle one. Considering the wall effect, the analytical model for predicting the size of matched double emulsions showed slight deviations from experiments if the continuous phase Reynolds number is lower than 0.06. However, the behaviour asymptotically approaches good agreement in the Reynolds number range of

0.06 to 0.14. The numerical simulation generally agreed with the experimental data under the investigated flow conditions.

**Author Contributions:** Conceptualisation, N.-T.N. and Y.D.; methodology, Y.D. and J.Z.; software, Y.D. and H.C.; material, Y.D. and F.G.; formal analysis, Y.D. and N.-K.N.; investigation, Y.D. and L.O.; data curation, Y.D., A.S.Y. and J.Z.; writing—original draft preparation, Y.D.; writing—review and editing, J.Z., H.A. and N.-T.N.; visualisation, Y.D.; supervision, N.-T.N.; project administration, N.-T.N.; funding acquisition, N.-T.N. and C.H.O. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Discovery Project (Grant No. DP220100261), ARC DECRA fellowship (Grant No. DE210100692), and ARC Future Fellowships (FT180100361). This work was partly performed at the Queensland node-Griffith-of the Australian National Fabrication Facility, a company established under the National Collaborative Research Infrastructure Strategy to provide nano and microfabrication facilities for Australian researchers.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**

