Author's Accepted Manuscript

Double emulsion production in glass capillary microfluidic device: Parametric investigation of droplet generation behaviour

Seyed Ali Nabavi, Goran T. Vladisavljevic, Sai Gu , Ekanem E. Ekanem

www.elsevier.com/locate/ces

PII: S0009-2509(15)00173-6

DOI: http://dx.doi.org/10.1016/j.ces.2015.03.004

Reference: CES12212

To appear in: Chemical Engineering Science

Received date: 9 December 2014

Revised date: 15 February 2015

Accepted date: 6 March 2015

Cite this article as: Seyed Ali Nabavi, Goran T. Vladisavljevic, Sai Gu, Ekanem E. Ekanem, Double emulsion production in glass capillary microfluidic device: Parametric investigation of droplet generation behaviour, Chemical Engineering Science, http://dx.doi.org/10.1016/jxes.2015.03.004

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Double emulsion production in glass capillary microfluidic device: Parametric investigation of droplet generation behaviour

Seyed Ali Nabavia, Goran T. Vladisavljevicb*, Sai Gua'*, Ekanem E. Ekanemb

a Offshore, Process and Energy Engineering Department, Cranfield University, Cranfield MK43 0AL, United Kingdom.

b Department of Chemical Engineering, Loughborough University, Loughborough LE11 3TU, United Kingdom.

* Corresponding author at: School of Engineering, Cranfield University, Cranfield, MK43 0AL, United Kingdom. Tel.: +441234 755277; fax: +441234 754685. E-mail address: s.gu@cranfield.ac.uk (S. Gu).

Corresponding author. Fax: +44 1509223923.E-mail address: g.vladisavljevic@lboro.ac.uk (G.T. Vladisavljevic)

Abstract

A three-phase axisymmetric numerical model based on Volume of Fluid-Continuum Surface Force (VOF-CSF) model was developed to perform parametric analysis of compound droplet production in three-phase glass capillary devices that combine co-flow and countercurrent flow focusing. The model predicted successfully generation of core-shell and multi-cored double emulsion droplets in dripping and jetting (narrowing and widening) regime and was used to investigate the effects of phase flow rates, fluid properties, and geometry on the size, morphology, and production rate of droplets. As the outer fluid flow rate increased, the size of compound droplets was reduced until a dripping-to-jetting transition occurred. By increasing the

middle fluid flow rate, the size of compound droplets increased, which led to a widening jetting regime. The jetting was supressed by increasing the orifice size in the collection capillary or increasing the interfacial tension at the outer interface up to 0.06 N/m. The experimental and simulation results can be used to encapsulate CO2 solvents within gas-permeable microcapsules.

Keywords: Dripping regime; jetting regime; core-shell droplets; flow focusing; solvent encapsulation; glass capillary microfluidics.

1. Introduction

Double water-in-oil-in-water emulsions are three-phase dispersions composed of inner aqueous droplets dispersed in larger oil droplets, which are themselves dispersed in another aqueous phase (Utada et al., 2005). The presence of intermediate (middle) fluid as a protective shell or semipermeable barrier which separates the inner aqueous phase from the outer one, makes double emulsions suitable in a wide range of applications including food industry (Edris and Bergnstahl, 2001; Muschiolik, 2007), cosmetics (Gallarate et al., 1999), controlled delivery (Kim and Park, 2004; Yamaguchi et al., 2002), and encapsulation (Aines et al., 2013; Martinez et al., 2012). Microcapsules for encapsulation of CO2 solvents within a CO2 permeable polymer are a promising alternative to conventional liquid CO2 capture materials such as monoethanolamine (Aines et al., 2013). The most challenging step in the production of tailor-made CO2 solvent capsules is to fabricate monodispersed core/shell template droplets with controllable shell thickness. The thinner the shell and the higher the surface area of the innermost droplet, the

higher the permeability to CO2 and the CO2 loading capacity, but the mechanical stability of capsules over multiple absorption-desorption cycles is also important. In order to achieve active control over the capsule production, understanding the underlying physics behind double emulsion formation in microfluidic devices is obligatory.

Conventional emulsification devices such as high-shear homogenizers and colloid mill (King and Keswani, 1994; Maa and Hsu, 1996) are based on high shear mixing of immiscible liquids. These devices, however, suffer from poor droplet size reproducibility and are not suitable for generation of core/shell droplets. Recently, microfluidic emulsification devices have attracted much attention due to their unprecedented level of control over droplet size and morphology. The most common microfluidic strategies for drop generation are flow focusing (Chen et al., 2009; Pannacci et al., 2008; Seo et al., 2007), T-junction (Nisisako et al., 2005; Okushima et al., 2004) and co-flowing (Cramer et al., 2004; Herrada et al., 2010; Perro et al., 2011) . Utada et al. (2005; 2007) developed a glass capillary device by combining co-flow and flow focusing in coaxial glass capillaries. In microfluidic devices, droplets can be produced in two main regimes, dripping and jetting, depending on the balance between interfacial, viscous, inertial and gravity forces (Lagus and Edd, 2013). The experimental studies of drop generation in glass capillary devices were focused on of the effects of fluid flow rates and geometry on the size and morphology of compound droplets (Chang and Su, 2008; Chang et al., 2009; Chen et al., 2012; Chu et al., 2007; Erb et al., 2011; Lee and Weitz, 2008; Martinez et al., 2012; Shirk et al., 2013) and dripping-to-jetting transition (Utada et al., 2005).

Although considerable experimental research has been devoted to double emulsion formation in microfluidic devices, there are only a few numerical parametric studies on microfluidic compound droplet production (Herrada et al., 2010; Park and Anderson, 2012; Radev and

Tchavdarov, 1988; Suryo et al., 2006; Vu et al., 2013; Zhou et al., 2006). Zhou et al. (2006) have studied the formation of compound droplets in a planar (2D) flow focusing microfluidic device using diffuse-interface method. However, in their setup all inlet streams are delivered from the same direction, while in 3D glass capillary device, the outer fluid flows counter-currently to the inner and middle fluid. Park and Anderson (2012) developed a two dimensional axisymmetric numerical model based on a ternary diffuse-interface method to predict dripping-to-jetting transition in glass capillary device and the effect of outer fluid flow rate on the size of double emulsion droplets. However, the model was not validated and a simplified setup was assumed with all inlet streams flowing co-currently. Vu et al. (2013) have used the front-tracking finite difference method to study the effect of surface tension, inlet velocities and radius ratio on the size of compound droplets formed in an axisymmetric co-flow device under jetting and dripping regime. However, the effect of viscosities and densities of middle and outer fluid was not investigated and again, a co-flowing outer fluid was considered.

There are many approaches used for simulation of free surface flow, which can be classified into interface tracking methods, such as immersed boundary method (Shin and Juric, 2002; Tryggvason et al., 2001), finite element method (Notz et al., 2001; Wilkes et al., 1999), boundary integral method (Cristini et al., 2001, 1998), and interface capturing methods including Level set method (Osher and Fedkiw, 2001), constrained-interpolation-profile method (Yabe et al., 2001), and volume of fluid method (VOF) (Hirt and Nichols, 1981). The VOF method coupled with CSF (continuum surface force) (Brackbill et al., 1992) approach was found to be very successful in predicting droplet formation and break up (Chen et al., 2013b; Vladisavljevic et al., 2014; Zhang, 1999), droplet deformation (Chen et al., 2013a) and dripping to jetting transitions (Chen et al., 2013b; Herrada et al., 2008).

To the best of our knowledge, this is the first numerical simulation study of double emulsion formation in an axisymmetric three-phase glass capillary device developed by Utada et al. (2005). We have used a two-dimensional axisymmetric VOF approach (Hirt and Nichols, 1981) coupled with continuum surface force (Brackbill et al., 1992) to simulate generation of core/shell droplets encapsulating CO2 solvents. Overall 13 parameters (8 physical parameters of the fluids, 3 fluid flow rates and 2 geometrical parameters of the device) have been systematically varied to investigate their effects on the droplet generation behaviour.

2. Governing equations and numerical method 2.1 Governing equations

The governing equations are those for the conservation of mass (continuity), Eq. (1), and momentum, Eq. (2).

*+v 2 at

\ / + V 2

=-VP+V2

where U is the velocity vector, P is the pressure, and t, ^ and p are the time, dynamic viscosity and density, respectively. A source term is denoted by Fb and includes two body forces, the gravitational force and interfacial tension force, F^. Since the length-scale is in the order of

micro-scale, the gravitational acceleration is negligible and Fb = F^

In the VOF model, a momentum equation, Eq. (2), was solved for all phases, and tracking of the interface was achieved by solving a continuity equation for volume fraction, f', of one or more phases. Eq. (3).

— + U.Vf' = 0 (3)

Eq. (3) explains the advection of fluids through the cells to gather the required information for reconstruction of the interface. The portion of each cell filled with fluids is determined by volume fraction in each cell where

the cell is filled with fluid 1

the interface of fluids exists in the cell the cell is filled with fluid 2

f' = 0 0 < f' < 1 f'= 1 <

The dynamic viscosity and density in the momentum equation, Eq. (2) were calculated as follows:

I = f| +(1 - f) |i2 (5)

p = f pi +(1 - f' ) p 2 (6)

The subscripts 1 and 2 present the first and second existing phases in the cell when the cell is filled with two phases. Continuum surface force (CSF) method was utilized to calculate the interfacial surface force term at the free surface in Eq. (2), where:

Fa= oKVf (7)

where o is the interfacial tension and k is the local curvature of the interface calculated as follows:

k = V-n (8)

where n is the unit normal defined as:

n = ^ (9)

2.2 Numerical method

The governing equations were solved by utilizing an unsteady pressure-based segregated algorithm established in finite-volume based commercial software Ansys Fluent v. 14.0. The discretised moment equation was approximated using both second order upwind and QUICK scheme in order to get the most accurate solution. Since there was no discrepancy between the solutions obtained by the two schemes, second order upwind scheme was chosen to carry on the simulations. The interpolation of pressure term was achieved by PRESTO scheme which directly calculates the pressure term on the faces. Although PRESTO scheme is computationally costly, it results in direct calculation of pressure term on faces and avoids interpolation errors. The pressure-velocity coupling was achieved by SIMPLE scheme. The interface interpolation was accomplished by Geo-Reconstruct algorithm. A variable time step method using Courant number, Co, was utilized in order to reduce the computational cost. The prescribed Co in the present simulations was 0.35. Fig. 1 is a schematic diagram of the three-phase glass capillary model used in this study along with the generated mesh. The red line shows the computational

domain and since the model is axisymmetric, the results are based on three-dimensional solutions. Grid dependency study was performed by constructing five meshes with a resolution of 8 ^m, 4 ¡¡m, 3 ¡¡m, 2 ¡¡m and 1 ¡¡m. No difference in behaviour was found for the resolutions of 1-4 ¡¡m. Considering computational cost and the large number of simulations required for this study, it was very time-consuming to employ very fine mesh for the whole computational domain. Therefore, a very fine mesh, 2 ¡¡m, was used inside the collection tube and around injection nozzle, where droplet formation occurs, while a coarser mesh was used for the rest of the domain. The boundary conditions are summarised in Table 1.

2.3 Dimensionless numbers

Non-dimensional numbers used in this study are Weber numbers of inner phase, We1, and that of compound jet, We2, and capillary numbers of inner, Ca1, middle, Ca2 and outer phase, Ca3:

2 2 We1 = We2 = P 2uj*~ja

°12 °23

Ca, = ^L Ca2 = ^L Ca3 = Ml

012 0 23 0 23

where p , o, ^, and u are density, interfacial tension, viscosity, and average velocity, respectively, D is a diameter, and subscripts 1, 2, 3, jet, and N stand for inner, middle and outer phase, respectively, compound jet, and nozzle. The values of ujet, u2, u3, and Djet were directly

measured in the vicinity of exit orifice from the numerical data, while u1 was calculated as

u1 = 4Q1 / (nDN). Since the difference in density between the inner and middle phase is very

small, the outer fluid density is used for estimation of We2 .

3. Experimental Method

3.1 Material

The inner phase was composed of 95 wt% deionized water and 5 wt% glycerol (Fisher scientific, UK). The middle phase was a mixture of 2 wt% xiameter® rsn- 0744 resin (UNIVAR, UK) and 98 wt% polydimethyl siloxane (PDMS) fluid (Dow Corning 200/ 10c S fluid, VWR, UK). The outer phase was a mixture composed of 58 wt% deionized water, 40 wt% glycerol and 2 wt% polyvinyl alcohol (PVA, Mw = 13,000-23,000 g mol-1, 87-89% hydrolysed, Sigma-Aldrich, UK). Prior to each use, the outer phase was kept in a quiescent state for 10 min to allow entrained air bubbles to rise to the surface. Table 2 lists the density and viscosity of all phases, measured using a pycnometer and capillary viscometer, respectively. The interfacial tension at the inner and outer interface of 22.6 and 7.3 mN/m, respectively, was measured using a Kruss DSA-100 pendant drop tensiometer.

3.2 Device fabrication

Round capillaries were supplied from Intracel (Royston, UK) with an inner diameter of 0.58 mm, an outer diameter of 1 mm, and a length of 150 mm. A Flaming/Brown micropipette puller (P-97, Sutter Instrument Co., Linton Instrumentation, Norfolk, UK) was used to break each round capillary into two identical parts, each with a tapered end. The pulled capillary was mounted onto a microforge microscope (Narishige MF-830) to observe the shape of the tip while it was sandpapered to the final size. The tip of both capillaries was then treated with octadecyltrimethoxysilane and 2-[methoxy(polyethyleneoxy)propyl]trimethoxysilane to make the glass surface hydrophobic and hydrophilic, respectively.

The injection capillary was inserted into a square capillary (AIT Glass, Rockaway, US) with an inner width of 1 mm. The two capillaries were axially centred on a microscope slide and fix in a position using two-part epoxy glue (5-Minute Epoxy® Devcon). The collection tube was then inserted into the square capillary, placed at a desirable distance from the injection capillary, axially centred and fix to the microscope slide. Three hypodermic needles with polypropylene hubs (BD Precisionglide® 20G, Sigma-Aldrich, UK) were attached to the capillaries to introduce the three phases to the device (Fig. 2). A needle hub with a single groove was used to deliver the inner phase to the opening of the injection capillary. Needle hubs with two grooves were placed at the intersections of the injection and collection capillary with the square capillary and used to introduce the middle and outer phase, respectively. The device was left 4-5 h to allow the epoxy to fully cure.

3.3 Generation of double emulsions

Three gas-tight glass syringes were filled with the three phases, mounted on syringe pumps (Harvard Apparatus model 11 Elite) and delivered through polyethylene medical tubing with 0.86mm I.D. and 1.52mm O.D. Once the compound jet emerged in the collection tube, double emulsion droplets were started to form (Fig. 2). Droplet formation was monitored using an inverted microscope (XDS-3, GX Microscopes, UK) and Phantom V9.0 high-speed camera interfaced to a PC computer. The pictures were taken at a speed of 3000 frames per second and the videos were processed using ImageJ software.

4. Results and Discussions 4.1 Validation of numerical model

The qualitative validations of the numerical model with conducted experiments are presented in Fig. 3. It is obvious that the CFD model is capable of predicting morphological changes of compound droplets during formation in the collection capillary, including formation of satellite droplets. Due to fast droplet formation, the model does not account for the adsorption of surfactant molecules at the interfaces, which leads to the coalescence of inner droplets in numerical simulations, behaviour not observed in experiments (Fig. 3). Since the droplet formation time was negligible compared to the time needed for the surfactants to reach the equilibrium saturation, the interfacial tension was considered time-independent.

4.2 Experimental results

The effect of phase flow rates on the inner droplet diameter, D1, the outer droplet diameter, D2 and the shell thickness, ts, is presented in Fig. 4. Solid and dashed lines in Figs. 4 a-c are theoretical predictions based on observed droplet generation frequencies and flow rate settings. As can be seen in Fig. 4 c, D1 can be varied by tuning the flow rate ratio Q1 / Q2, where an

increase in Q1 /Q2 leads to a significant increase in D1, while D2 remains almost constant, which results in thinning the shell. In contrast, an increase in Q2 / Q1 increases D2 and reduces D1 and consequently, the shell becomes thicker, as shown in Fig. 4 b. An increase in Q3

considerably reduces both D2 and D1 while ts remains nearly constant (Fig. 4 a). Fig. 4 d gives a comparison between the experimental results and Eq. (10):

! + Q1 + Q2

As can be seen in Fig. 4 d, a very good agreement was achieved between the experimental and predicted values. It should be noted that a relationship between fluid flow rates and drop morphology (D2, D1 and ts) can be affected by the viscosity ratios of the fluids. For example,

Chang and Su (2008) observed an increase in D1 and D2 as Q1 increased, which is in contrast

with the present work and the results of Lee and Weitz (2008). The discrepancy between the experimental results may have arisen from the fact that Chang and Su (2008) used an outer fluid which was about 30 times more viscous that the middle fluid whereas in the present study a viscosity ratio of the outer to middle fluid is 1.3.

4.3 Parametric numerical study of compound droplet formations

In this section parametric numerical investigation of compound droplet formation in glass capillary microfluidic devices is presented. The chosen liquid phase properties are within a range typically found for encapsulation of conventional CO2 solvents such as potassium carbonate and MEA (monoethanolamine) within polymerizable liquid polymers (Aines et al., 2013). The encapsulation of 30 wt% aqueous MEA solution within a UV-curing liquid polymer adhesive (NOA 61) is simulated in Figs. 5, 7, 8, and 9.

4.3.1 Increasing Q3 reduces both D1 and D2

The effect of outer fluid flow rate, Q3, on droplet formation is illustrated in Fig. 5. In the

dripping regime (Q3 / Qsum < 5.17 ), compound droplets are formed very close to the entry of the

collection capillary (L2 / DOrif 1.2, Fig. 5 f) and the size of both inner and outer droplets is

significantly reduced with increasing Q3 (Fig. 5 a-b), which is coupled with an increase in both

droplet generation frequencies (Fig. 5 g). Q3 has a higher impact on D2 than D1, which leads to

the reduction in shell thickness, ts . A short-term increase in the size of inner and outer drops at

Q3 / Qsum 7 in Fig. 5 e was due to transition from dripping to narrowing jetting regime, but both

D1 and D2 resumed its downward trend at Q3 / Qsum > 7.18. A temporary increase in the drop

size can be explained by the variations of Ca3 and We2 presented in Fig. 5 h. At the transition

point, although there is a considerable decrease in D jet , the steep jump in We2 implies a sudden

increase in the velocity of compound jet, while there is a constant increase in Ca3. Therefore, the

velocity gradient at the outer interface decreases and thus, the shearing imposed by the outer phase on the compound jet also decreases, which results in an increase in droplet formation time of both drops, 1/f1 and 1/f2 (Fig. 5 g). Since Qsum= const, an increase in droplet formation

times results in larger D1 and D2 values.

The transition to jetting regime leads to a dramatic increase in break up lengths (Figs. 5 c-d), with L2 /DOrif reaching 8.7 at Q3 /Qsum = 26 (Fig. 5 f). The size of inner and compound droplets continue to decrease in the jetting regime so that the diameter of the compound droplet is less than one half of the orifice diameter at Q3 / Qsum = 26. A reduction in droplet size as a result of

increase in the outer fluid flow rate can be explained by an increase in the viscous stress exerted by the outer fluid to compound jet, which reduces the drop formation time and thus, increases droplet generation frequencies, f1 and f2 to almost 200 Hz at Q3 /Qsum = 26 (Fig. 5 g). Since

f1 = f2 at all flow rates shown in Fig. 5, each compound droplet contains a single internal drop.

A shear stress at the outer interface, t 2,3, is as follows (see the supplemental material S1):

T23 = -a-2 (11)

n Dorif3 (a2 -1)2

where a is given by Eq. (s.12). The variation of t 2 3 with Q3 calculated from Eq. (11) is shown in Fig. 6. An increase in Q3 leads to an increase in t 2 3, but the effect is less pronounced at Q3 / Qsum > 20.

Fig. 5 h shows the variation of the outer flow capillary number, Ca3 with Q3. In dripping regime, where Ca3 < 0.1, the viscous forces are negligible compared to interfacial forces leading to the formation of droplets adjacent to the tip of the injection tube. At Ca3 - 0.1, the viscous and interfacial forces are of the same order of magnitude causing the transition from dripping to narrowing jetting regime. There is an abrupt increase in L1 and L2 with Q3 at the transition

point. The further increase in Q3 leads to additional elongation and thinning of the compound jet. At Ca3 > 0.26, the break-up occurs due to Rayleigh instability and the droplet diameter approaches the jet diameter. The shell thickness in the jetting regime decreases with Q3 and ts / DOrif reaches 0.12 at Q3 / Qsum = 26, which corresponds to ts = 36 ¡¡m . Based on the fluid properties and the simulated geometry, the core-shell droplets are successfully produced at 1.4 <

Q3 / Qsum < 57.9. At Q3 / Qsum < 1.4, the outer flow cannot focus the compound jet through the collection tube and the middle phase wets the collection tube wall. At Q3 / Qsum > 57.9, the outer phase pushes the compound jet back into the injection nozzle and the droplet formation fails.

4.3.2 Increasing Q2 reduces D1 while increasing D2

An increase in Q2 results in higher D2 values, while D1 is reduced and accordingly, the shell becomes thicker (Fig. 7 e). An increase in Q2 results in higher D2 values, while D1 is reduced and accordingly, the shell becomes thicker (Fig. 7 e). An increase in Q2 increases the middle fluid velocity and correspondingly the compound jet inertia. The increase in D2 with Q2 may be attributed to the high inertia of the middle fluid which leads to faster introduction of middle fluid to the droplet before drop detachment. An increase in Q2 results in higher velocity gradient at

the inner interface and thus increases the shear force exerted on the inner jet by the middle fluid. Therefore, the drop generation frequency of inner droplets increases (Fig. 7 g), leading to the formation of smaller inner droplets. According to the mass balance equation: f1 = 6Q1 / (D13). Since Q1 is a constant, f1 varies inversely as the cube of D1. An increase in f2 in Fig. 7 g is a consequence of increase in Q2.

An increase in Q2 results in higher We2 values and longer break-up lengths, L1 and L2, as

shown in Fig.7 h and 7 f respectively. A widening shape of the compound jet is a consequence of higher velocity of the middle fluid compared to that of the outer fluid. Therefore, a drag force at the outer interface tends to slow down and widen the compound jet. As Q2 increases, a

difference in velocity between the middle and outer fluid increases, leading to the formation of increasingly wider and longer compound jet. Since in widening jetting regime, the momentum force of compound jet is responsible for jet breakup, higher Q2 values result in a faster droplet

production and higher f2 values, as shown in Fig. 7 g. Based on the fluid properties used in the simulations, the successful core-shell droplet generation is at 0.58 < Q2 /Q1 < 46.7. For Q2 /Q1 < 0.58, the shell of compound jet is so narrow that the inner fluid ruptures the compound jet and the droplet formation fails. For Q2 /Q1 > 46.7, the compound jet diameter exceeds the orifice

diameter, DOrif and wets the collection tube leading to failure in droplet formation. 4.3.3 Increasing Q1 increases both D1 and D2

Fig. 8 demonstrates the effect of innermost fluid flow rate, Q1, on the drop formation behaviour. The drop formation in the dripping regime occurs over the Q1 /Q2 range from 0.13 to 0.66, corresponding to 0.013< We1< 0.61 (Fig. 8 a-c). An increase in Q1 leads to an increase in both inertial and viscous forces exerted by the inner fluid, which results in an increase in both We1 and Ca1 (Fig. 8 h). Since there is no considerable change in L1 and L2 over the range of 0.13 < Q1 / Q2 < 0.66 (Fig. 8 g), an increase in the momentum of inner fluid results in bigger inner droplet. In the dripping regime Q1 does not have any noticeable effect on the shear stress at the outer interface and thus, D2 increases only slightly with Q1 (Fig 8 e). Since the size of inner drop increases with increasing Q1 and the size of outer drop remains nearly constant, the shell becomes thinner. As Q1 /Q2 exceeds a critical value of 0.66, any further increase in Q1 results in developing a high speed jet of the inner fluid leading to a burst of the middle shell and failure of

drop formation (Fig. 8 d). In addition, for Q1 /Q2 < 0.032, the core-shell formation is periodically interrupted by the formation of shell droplet without any encapsulated inner droplet. A failure in drop formation when Q1 exceeds a certain critical value was also reported earlier

(Chang and Su, 2008). f1 and f2 increased very slightly with Q1, as shown in Fig. 8 f.

The influence of Q1 on drop formation in jetting regime is presented in Fig. 9. As shown in Fig. 9 e, D1 follows the same trend as in dripping regime, but D2 increases more significantly, albeit at a slower pace than D1, so that the shell thickness is reducing. At very low Q1 / Q2 value of 0.0248, where We1 - 0.006, inner droplets are formed in a narrowing jetting regime (due to the high Q2 value of 2.46 ml/h), while the outer droplets are formed in a widening mode, resulting in the formation of large compound drops with very small inner drops (Fig. 9 a). Formation of inner and outer drops in different regimes results in a considerable difference between L1 and

L2, which leads to the generation of multi-cored droplets. As Q1 increases, the inner drop formation switches from narrowing to widening mode, resulting in virtually the same value of L1 and L2 and formation of core/shell droplets. A high speed jet was observed at Q1 /Q2 > 0.62

and We1> 3.83 (Fig. 9 d), similar to that shown in Fig. 8 d. However, in this case the shell rupture occurs after droplet pinch off, further downstream.

4.3.4 Increasing a23 decreases both L1 and L2

Fig. 10 shows the droplet formation behaviour as a function of interfacial tension at the outer interface, a 23 over the range of a 23 from 0.005 to 0.095 N/m. The interfacial force acting on

the outer interface increases with increasing a 23, which results in the reduction of Ca3 from 0.51 to 0.076 and the reduction in We2 from 0.322 to 0.031 (Fig. 10 h). The reduction in outer fluid capillary number and inner fluid Weber number due to increase in interfacial tension was also reported by Lagus and Edd (2013). Higher a 23 values suppress the jetting mode by pulling back the compound jet and switching the droplet formation to dripping regime, as shown in Fig. 10 ad. Consequently, both break-up lengths, L1 and L2, decrease with increasing a 23 (Fig. 10 f), which is consistent with the findings of Herrada et al. (2010). Two different behaviours have been observed on the graphs in Figs. 10 e and g. At a23 =0.005-0.06 N/m, an increase in a23

considerably delays the pinch off events and reduces both f1 and f2 , resulting in an increase in D1 and D2 . An increase in the size of the outer drop is more significant than that of the inner drop causing ts to grow.

A different type of behaviour was observed at high interfacial tensions at the outer interface of 0.060-0.095 N/m (Fig. 10 e and g). Within this range, as a 23 increases both inner and outer

droplets shrink causing the drop generation frequency to increase in order to maintain constant volume flow rates, Q1 and Q2. The reason for this behaviour could be droplet formation in

jetting to dripping transition mode. The same trend in opposite way have been observed in droplet generation during transition from dripping to narrowing jetting regime (section 4.3.1), where further increase in Q3 led to a sudden growth in both D1 and D2 .

4.3.5 Increasing <712 suppresses inner fluid jetting

The effect of interfacial tension at the inner interface, g12 , on the drop formation is shown Fig 11. Starting from g12 = 0.00574 N/m, where both inner and outer droplets are formed in jetting regime, increasing g12 causes a reduction in Ca1 (Fig. 11 h), which implies higher interfacial force compared to shear force of inner jet. Therefore, there is an increasing tendency to pull back the forming inner jet towards the tip of the nozzle. At q12 = 0.01 N/m (Fig. 11 a), the inner and

middle phase jets are of similar length (L1 - L2). As g12 increases to 0.015 N/m, L2 remains

fairly constant while the inner jet shrinks due to increased interfacial tension force at the inner interface which causes development of capillary instability on the inner jet. Consequently the inner jet pinches off faster, leading to smaller D1 values (Fig. 11 e & g). Indeed, for 0.005

< q12 < 0.03 N/m, both L1 and L2 decreases with increasing g12 (Fig. 11 a-c). This effect is more pronounced for the inner fluid, causing the difference between L1 and L2 to increase, which leads to the formation of multi-cored compound droplets (Fig. 11 b).

In the q12 range of 0.015-0.03 N/m, the formed inner droplets, which have the diameter

comparable to that of the compound jet, move downstream inside the compound jet, which results in the appearance of capillary-like waves on the outer interface. Therefore, the outer interface becomes more unstable and consequently the compound jet breaks closer to the nozzle exit, meaning a reduction in L2 (Fig. 11 c). The further increase in g12 over the range from 0.04

to 0.06 N/m leads to reduction in both L1 and D1, where D1 is considerably smaller than the

compound jet diameter. Therefore, the capillary waves previously formed on the outer surface (Fig. 11 c) due to large inner droplet are supressed (Fig. 11 d), and the compound jet pinch off is

delayed leading to an increase in L2. It should be mentioned that the droplet formation at high a12 values is very chaotic leading to formation of multi-cored droplets along with generation of large satellite droplets (Fig. 11 d).

4.3.6 Increasing ju3 increases both D2 and D1

The effect of outer fluid viscosity, ^ 3, on the production of compound drops is shown in Fig. 12.

According to Fig. 12 h, an increase in ^ 3 leads to a linear increase in Ca3 due to higher viscous

force acting on the outer interface, which results in stretching both the inner phase and middle phase jet (Fig. 12. f). The higher level of shearing exerted by outer fluid accelerates the droplet formation, which results in higher droplet formation frequencies, f1 and f2 , and smaller drop

diameters, D1 and D2 (Fig. 12 e & g). As can be seen in Figs. 12 a-b, the variations of D1 and

D2 with ^ 3 are strongly correlated, leading to almost constant shell thickness.

4.3.7 Increasing ju2 increases both L2 and L1

Fig. 13 shows the effect of middle fluid viscosity, ^ 2, on droplet production. As ^ 2 increases,

there is an increasing tendency for transition from dripping to jetting (Fig. 13 a-f). In the current numerical model, the effect of adsorption of surfactants onto the interfaces was not considered, so that droplet coalescence occurs, as shown in Fig. 13 a-c. As shown in Fig. 13 j, the capillary number of middle phase Ca2 increases with increasing ^ 2 , since Ca2 ^ 2 at constant u2 and

a23. The increased ^2 leads to an increase in viscous stress on inner and outer fluids from middle fluid, causing both liquid jets to stretch in the downstream direction. As shown in Fig. 11g, the transition from dripping to jetting caused by an increase in ^ 2 from 0.002 to 0.005 Pas

results in a steep increase in both D2 and D1 and drop formation time. However, there is no

considerable difference in drop formation time and accordingly drop size before and after dripping-to-jetting transition. This fact was pointed out by Suryo et al. (2006) who stated that the droplet size and formation time are insensitive to middle phase viscosity in a certain range of viscosity ratios, ^3 / ^2. The production frequency of compound droplets in the dripping regime

was relatively high, f2 = 120 s-1, compared to the one in jetting mode, f2 = 20 s-1. It can be explained by the fact that at very low ^2 (in dripping regime), due to large viscosity ^3 / ^2 the required force to pinch off the compound jet and make compound droplets is much smaller than that at high ^ 2 values and the compound droplet pinch off accelerates.

4.3.8 Reducing ju1 leads to formation of inner droplets in dripping regime

As illustrated in Fig. 14, D2 is almost unaffected by the innermost fluid viscosity, , due to

negligible effect of on the shear stress at the outer interface. However, there is a steep

increase in D1 as increases from 0.001 to 0.019 Pas , followed by a sharp drop in f1 from

nearly 600 s-1 to 60 s-1. At = 0.001 Pa-s, very small inner droplets are formed in the vicinity

of nozzle exit and there is a large difference in drop generation frequency between inner and outer droplets, resulting in multi-cored droplets. The number of internal droplets can be predicted from the relation: N1 = f1 /f2. E.g., at = 0.001 Pa-s, the number of inner droplets observed is

31 and f1 /f2 = 31.48, whereas at = 0.004 Pa-s, three inner droplets are engulfed by an outer

droplet and f1 /f2 = 3.39. According to (Chu et al., 2007), when the predicted value of N1 is

between two integers, N1 can take either integral value. An increase in to 0.0396 Pas causes an inner jet to develop along with a dramatic growth in D1 , while the outer jet is shortening, and core/shell droplets are formed (Fig. 14 a-d & f). The reduction in L2 can be attributed to the fact that at low values, fast and sequential formation of inner droplets delays the drop pinch off from the compound jet and results in larger L2. By increasing and consequently decreasing the number of inner droplets trapped by shell, the break-up of compound jet facilitates, resulting in shorter L2. The further increase in beyond = 0.0396 Pa-s results in an increase in Ca1,

corresponding to higher viscous forces of the inner fluid and a slight increase in both L1 and L2. However, there is no significant change in f1 and f2 and therefore D1 and D2 .

4.3.9 Fluid densities do not have any impact on flow field

The effect of fluid density on droplet formation is shown in Fig. 15. No difference in flow field was observed as the density of the outer fluid was varied from 600 to 1800 kg / m3. It could be explained by the fact that due to low fluid velocity and low density ratio, the convection term of Navier-Stokes equation which includes the effect of density is negligible. The same behaviour was found when the density of inner and middle fluid was modified.

4.3.10 Decreasing Dorif changes droplet formation regime from dripping to jetting

Fig. 16 shows the effect of orifice size, DOrif, on droplet formation. In Figs. 16 a-b, Q2 and Q3

are relatively high and the drop formation occurs in narrowing jetting regime. Here, an increase

in DOrif from 150 to 300 pm, corresponding to an increase in DOrif /DN from 6 to 12, causes a decrease in both L1 and L2 by 6-7 % as a result of decrease in pulling forces acting on both jets

due to smaller velocities of the middle and outer fluid. Due to decrease of the shear stress acting at the inner and outer interface, the time interval between two consecutive breakup events increase from 13 to 15 ms causing an increase in D1 and D2 by 4-8 % (Table s.2, jetting).

The effect of DOrif on the droplet formation at low flow rates is shown in Figs. 16 c & d. At

DOrif =150 pm (Fig. 16 c), droplets are formed in widening jetting regime since the middle phase

velocity is higher than the outer phase velocity. As a result of the frictional forces at the outer interface, the velocity of the middle phase decreases in the downstream direction causing the tip of the compound jet to grow to a large size. It is now the inertial force rather than shear force that must exceed the interfacial force for the jet to break up into drops. At DOrif =300 pm (Fig. 16 d),

droplets are formed in dripping regime since the outer phase velocity exceeds the middle phase velocity, but shear force at the outer interface is insufficient for narrowing jetting to occur.

4.3.11 Increasing L changes the droplet formation pattern at low flow rates

The effect of distance between the injection and collection capillary, L, on the drop formation is shown in Fig. 17 and Table s.3. As L increases at high flow rates (Fig. 17 a-b), D2 decreases while D1 is almost unaffected, leading to a decrease in the shell thickness (Table s.3). The increased distance L causes the compound jet to slow down before it enters the collection tube.

Therefore, a magnitude of the inertial force of the middle phase decreases compared to interfacial force at the outer interface, which causes a reduction in D2, L1, and L2.

Although L does not have a considerable impact on drop generation behaviour at high flow rates, it dramatically effects the drop generation pattern at low flow rates (Fig. 17 c-d). In Fig. 17 c, the droplets are formed in the widening jetting regime due to very low velocity of the outer phase, which is surpassed by the middle phase velocity. An increase in L considerably reduces the middle fluid velocity at the orifice of the collection tube. Since in this case the middle phase flow rate is relatively low, the viscous force originating from the middle phase is not sufficiently high to elongate the inner fluid and the pinch off happens at a very small distance from the exit orifice. Also, the generation frequency of the outer drops is much higher than that of the inner drops and this non-uniformity in drop generation frequency leads to alternating generation of middle phase droplets and compound droplets. Therefore, in this case the optimum value of L /DOrif lies between 0.3 and 0.83.

5. Conclusion

An in-depth parametric study of compound droplet formation in a three-phase glass capillary device was performed by developing a VOF-CSF numerical model. The model successfully predicted the experimental data and was used to investigate the effect of fluid flow rates, fluid properties and device geometry on the size, morphology and generation frequency of compound droplets. To the best of our knowledge, this is the first CFD model for simulation of three-phase glass capillary device, developed using the exact geometry without any simplifications. The model was capable of reproducing droplet formation in all regimes, as well as complex

phenomena such as satellite and multi-cored droplet formation and alternating generation of single and compound droplets.

The size of both inner and outer droplets was reduced as the outer fluid flow rate, Q3 or the

viscosity of the outer fluid, 3 was increased, which led to the transition from dripping to

narrowing jetting regime at certain critical value of Q3 and 3. An increase in the middle fluid

flow rate, Q2 reduced the size of inner drops and increased the size of outer drops, which led to a

transition to widening jetting regime at the critical Q2 value. As the inner fluid flow rate, Qx

increased in dripping regime, considerably larger inner droplets were observed but the size of outer droplets remained nearly constant. The jetting regime was supressed by increasing the size of the orifice in the collection capillary. When the distance between the injection and collection capillary was too large, at low fluid flow rates the generation of inner and outer droplets was non-synchronised, leading to alternating generation of compound and single droplets. As the interfacial tension at the outer interface increased from 0.005 to 0.06 N/m, the size of compound droplets increased and the jetting regime was highly suppressed. On the other hand, an increase of interfacial tension at the inner interface suppressed jetting of the inner fluid only and led to the formation of multi-cored droplets. The number of inner droplets in multi-cored droplets was consistent with the ratio of generation frequency of inner and outer droplets. As the viscosity of the middle phase increased, jet break up lengths increased until a transition from dripping to jetting occurred. As the viscosity of inner phase increased, the size of inner droplets increased and the number of inner droplets encapsulated within each outer droplet was reduced, but no impact on the size of outer droplets was observed. Due to small droplet sizes and low density ratios, no difference in droplet formation was observed at different fluid densities.

Although the developed model accurately predicted generation of compound droplets in glass capillary devices, due to the presence of two interfaces it was computationally very expensive and time consuming. In the future study, a grid adaption method will be coupled to the current CFD model to refine the meshes along the interface, so as to reduce the size of meshes.

6. Acknowledgement

The authors gratefully acknowledge the financial support for this work by the UK Engineering and Physical Sciences Research Council (EPSRC), project grant: EP/J020184/1and and FP7 Marie Curie iComFluid project grant: 312261. The authors would like to thank Dr. Andrew S. Utada for his helpful discussions.

Nomenclature

Symbol Description Unit

A surface area m2

Ca capillary number -

Co Courant number -

D droplet diameter m

Dn injection nozzle diameter m

DOrif exit orifice diameter m

E' integral constant Pa.m

E/ integral constant m/s

E 2'' integral constant m/s

f' volume fraction -

f droplet generation frequency 1/s

Fb body force N

F. interfacial force N

L break up length m

L ' distance between collection tube and injection nozzle m

n unit normal -

P pressure Pa

Q volume flow rate m3/s

r radial direction m

R radius m

S spreading coefficient N/m

t time s

ts shell thickness m

U mean velocity m/s

We Weber number -

x axial direction m

a radius ratio -

K curvature of interface 1/m

M dynamic viscosity kg/m-s

P density kg/m3

a interfacial tension N/m

0 contact angle -

Subscripts

1 inner phase

2 middle phase

3 outer phase

1,2 interface between inner and middle phases

2,3 interface between middle and outer phases

jet jet

Orif orifice

sum Q1+Q2

w wall

7. References

Aines, R.D., Spaddaccini, C.M., Duoss, E.B., Stolaroff, J.K., Vericella, J., Lewis, J. a., Farthing, G., 2013. Encapsulated solvents for carbon dioxide capture. Energy Procedia 37, 219-224. doi:10.1016/j.egypro.2013.05.105

Brackbill, J.., Kothe, D.., Zemach, C., 1992. A continuum method for modeling surface tension. J. Comput. Phys. 100, 335-354. doi:10.1016/0021-9991(92)90240-Y

Chang, F.-C., Su, Y.-C., 2008. Controlled double emulsification utilizing 3D PDMS

microchannels. J. Micromechanics Microengineering 18, 065018. doi:10.1088/0960-1317/18/6/065018

Chang, Z., Serra, C. a, Bouquey, M., Prat, L., Hadziioannou, G., 2009. Co-axial capillaries microfluidic device for synthesizing size- and morphology-controlled polymer core-polymer shell particles. Lab Chip 9, 3007-11. doi:10.1039/b913703c

Chen, C.-H., Abate, A.R., Lee, D., Terentjev, E.M., Weitz, D. a., 2009. Microfluidic assembly of magnetic hydrogel particles with uniformly anisotropic structure. Adv. Mater. 21, 32013204. doi:10.1002/adma.200900499

Chen, P.W., Erb, R.M., Studart, A.R., 2012. Designer polymer-based microcapsules made using microfluidics. Langmuir 28, 144-52. doi:10.1021/la203088u

Chen, Y., Liu, X., Shi, M., 2013a. Hydrodynamics of double emulsion droplet in shear flow. Appl. Phys. Lett. 102, 051609. doi:10.1063/1.4789865

Chen, Y., Wu, L., Zhang, C., 2013b. Emulsion droplet formation in coflowing liquid streams. Phys. Rev. E 87, 013002. doi:10.1103/PhysRevE.87.013002

Chu, L., Utada, A.S., Shah, R.K., Kim, J., Weitz, D.A., 2007. Controllable monodisperse multiple emulsions. Angew. Chemie 46, 8970-8974. doi:10.1002/anie.200701358

Cramer, C., Fischer, P., Windhab, E.J., 2004. Drop formation in a co-flowing ambient fluid. Chem. Eng. Sci. 59, 3045-3058. doi:10.1016/j.ces.2004.04.006

Cristini, V., Blawzdziewicz, J., Loewenberg, M., 1998. Drop breakup in three-dimensional viscous flows. Phys. Fluids 10, 1781. doi:10.1063/1.869697

Cristini, V., Blawzdziewicz, J., Loewenberg, M., 2001. An adaptive mesh algorithm for evolving surfaces: Simulations of drop breakup and coalescence. J. Comput. Phys. 168, 445-463. doi: 10.1006/jcph.2001.6713

Edris, a, Bergnstähl, B., 2001. Encapsulation of orange oil in a spray dried double emulsion. Nahrung 45, 133-7. doi:10.1002/1521-3803(20010401)45:2<133::AID-F00D133>3.0.C0;2-C

Erb, R.M., Obrist, D., Chen, P.W., Studer, J., Studart, A.R., 2011. Predicting sizes of droplets made by microfluidic flow-induced dripping. Soft Matter 7, 8757-8761. doi :10.1039/c1 sm06231j

Gallarate, M., Carlotti, M.E., Trotta, M., Bovo, S., 1999. On the stability of ascorbic acid in emulsified systems for topical and cosmetic use. Int. J. Pharm. 188, 233-241. doi:10.1016/S0378-5173(99)00228-8

Herrada, M. a., Montanero, J.M., Ferrera, C., Gañán-Calvo, a. M., 2010. Analysis of the dripping-jetting transition in compound capillary jets. J. Fluid Mech. 649, 523-536. doi: 10.1017/S0022112010000443

Herrada, M., Gañán-Calvo, A., Ojeda-Monge, A., Bluth, B., Riesco-Chueca, P., 2008. Liquid flow focused by a gas: Jetting, dripping, and recirculation. Phys. Rev. E 78, 036323. doi: 10.1103/PhysRevE.7 8.036323

Hirt, C.., Nichols, B.., 1981. Volume of fluid (VOF) method for the dynamics of free boundaries. J. Comput. Phys. 39, 201-225. doi:10.1016/0021-9991(81)90145-5

Kim, H.K., Park, T.G., 2004. Comparative study on sustained release of human growth hormone from semi-crystalline poly(L-lactic acid) and amorphous poly(D,L-lactic-co-glycolic acid) microspheres: morphological effect on protein release. J. Control. Release 98, 115-25. doi: 10.1016/j .jconrel.2004.04.020

King, A.G., Keswani, S.T., 1994. Colloid mills: Theory and experiment. J. Am. Ceram. Soc. 77, 769-777. doi:10.1111/j.1151-2916.1994.tb05364.x

Lagus, T.P., Edd, J.F., 2013. A review of the theory, methods and recent applications of high-throughput single-cell droplet microfluidics. J. Phys. D. Appl. Phys. 46, 114005. doi: 10.1088/0022-3727/46/11/114005

Lee, D., Weitz, D. a., 2008. Double emulsion-templated nanoparticle colloidosomes with selective permeability. Adv. Mater. 20, 3498-3503. doi:10.1002/adma.200800918

Maa, Y.-F., Hsu, C., 1996. Liquid-liquid emulsification by rotor/stator homogenization. J. Control. Release 38, 219-228. doi:10.1016/0168-3659(95)00123-9

Martinez, C.J., Kim, J.W., Ye, C., Ortiz, I., Rowat, A.C., Marquez, M., Weitz, D., 2012. A microfluidic approach to encapsulate living cells in uniform alginate hydrogel microparticles. Macromol. Biosci. 12, 946-51. doi:10.1002/mabi.201100351

Muschiolik, G., 2007. Multiple emulsions for food use. Curr. Opin. Colloid Interface Sci. 12, 213-220. doi :10.1016/j .cocis.2007.07.006

Nisisako, T., Okushima, S., Torii, T., 2005. Controlled formulation of monodisperse double emulsions in a multiple-phase microfluidic system. Soft Matter 1, 23-27. doi:10.1039/b501972a

Notz, P.K., Chen, A.U., Basaran, O. a., 2001. Satellite drops: Unexpected dynamics and change of scaling during pinch-off. Phys. Fluids 13, 549. doi:10.1063/1.1343906

Okushima, S., Nisisako, T., Torii, T., Higuchi, T., 2004. Controlled production of monodisperse double emulsions by two-step droplet breakup in microfluidic devices. Langmuir 20, 99058. doi:10.1021/la0480336

Osher, S., Fedkiw, R.P., 2001. Level Set Methods: An overview and some recent results. J. Comput. Phys. 169, 463-502. doi:10.1006/jcph.2000.6636

Pannacci, N., Bruus, H., Bartolo, D., Etchart, I., Lockhart, T., Hennequin, Y., Willaime, H., Tabeling, P., 2008. Equilibrium and nonequilibrium states in microfluidic double emulsions. Phys. Rev. Lett. 101, 164502. doi:10.1103/PhysRevLett.101.164502

Park, J.M., Anderson, P.D., 2012. A ternary model for double-emulsion formation in a capillary microfluidic device. Lab Chip 12, 2672-7. doi:10.1039/c2lc21235h

Perro, A., Nicolet, C., Angly, J., Lecommandoux, S., Le Meins, J.-F., Colin, A., 2011. Mastering a double emulsion in a simple co-flow microfluidic to generate complex polymersomes. Langmuir 27, 9034-42. doi:10.1021/la1037102

Radev, S., Tchavdarov, B., 1988. Linear capillary instability of compound jets. Int. J. Multiph. Flow 14, 67-79. doi:10.1016/0301-9322(88)90034-1

Seo, M., Paquet, C., Nie, Z., Xu, S., Kumacheva, E., 2007. Microfluidic consecutive flow-focusing droplet generators. Soft Matter 3, 986-992. doi:10.1039/b700687j

Shin, S., Juric, D., 2002. Modeling Three-Dimensional Multiphase flow using a level contour reconstruction method for front tracking without connectivity. J. Comput. Phys. 180, 427470. doi:10.1006/jcph.2002.7086

Shirk, K., Steiner, C., Kim, J.W., Marquez, M., Martinez, C.J., 2013. Assembly of colloidal silica crystals inside double emulsion drops. Langmuir 29, 11849-57. doi:10.1021/la4019986

Suryo, R., Doshi, P., Basaran, O.A., 2006. Nonlinear dynamics and breakup of compound jets Nonlinear dynamics and breakup of compound jets. Phys. Fluids 18, 082107. doi: 10.1063/1.2245377

Tryggvason, G., Bunner, B., Esmaeeli, a., Juric, D., Al-Rawahi, N., Tauber, W., Han, J., Nas, S., Jan, Y.-J., 2001. A front-tracking method for the computations of multiphase flow. J. Comput. Phys. 169, 708-759. doi:10.1006/jcph.2001.6726

Utada, A.S., Chu, L., Link, D.R., Holtze, C., Weitz, D.A., 2007. Dripping, jetting, drops, and wetting□: The magic of microfluidics. MRS Bull. 32, 702-708.

Utada, A.S., Lorenceau, E., Link, D.R., Kaplan, P.D., Stone, H.A., Weitz, D.A., 2005.

Monodisperse double emulsions generated from a microcapillary device. Science 308, 537541. doi :10.1126/science .1109164

Vladisavljevic, G.T., Shahmohamadi, H., Das, D.B., Ekanem, E.E., Tauanov, Z., Sharma, L., 2014. Glass capillary microfluidics for production of monodispersed poly (dl-lactic acid) and polycaprolactone microparticles: Experiments and numerical simulations. J. Colloid Interface Sci. 418, 163-70. doi:10.1016/j.jcis.2013.12.002

Vu, T. V., Homma, S., Tryggvason, G., Wells, J.C., Takakura, H., 2013. Computations of

breakup modes in laminar compound liquid jets in a coflowing fluid. Int. J. Multiph. Flow 49, 58-69. doi: 10.1016/j.ij multiphaseflow .2012.10.004

Wilkes, E.D., Phillips, S.D., Basaran, O. a., 1999. Computational and experimental analysis of dynamics of drop formation. Phys. Fluids 11, 3577. doi:10.1063/1.870224

Yabe, T., Xiao, F., Utsumi, T., 2001. The Constrained Interpolation Profile Method for Multiphase Analysis. J. Comput. Phys. 593, 556-593. doi:10.1006/jcph.2000.6625

Yamaguchi, Y., Takenaga, M., Kitagawa, A., Ogawa, Y., Mizushima, Y., Igarashi, R., 2002. Insulin-loaded biodegradable PLGA microcapsules: initial burst release controlled by hydrophilic additives. J. Control. Release 81, 235-249. doi:10.1016/S0168-3659(02)00060-3

Zhang, X., 1999. Dynamics of drop formation in viscous flows. Chem. Eng. Sci. 54, 1759-1774. doi:10.1016/S0009-2509(99)00027-5

Zhou, C., Yue, P., Feng, J.J., 2006. Formation of simple and compound drops in microfluidic devices. Phys. Fluids 18, 092105. doi:10.1063/1.2353116

Table 1. A summary of the boundary conditions applied in this work

Boundary Type

Inlets Mass flow inlet

Outlet Pressure outlet

Symmetry Axis Axis

Wall No-slip wall

Table 1. The properties of the fluids used for experimental double emulsion preparation.

Density Viscosity

Phase (kg/m3) (mPa-s)

Inner phase 1012 1.24

Middle phase 940 10.37

Outer phase 1107 7.91

Inlet (Q,) Symmetry Axis

Fig. 1. Schematic of simulated geometry along with the mesh.

Inner Phase Inlet

Middle Phase Inlet

Outer Phase Inlet

Fig. 2. A fabricated three-phase glass capillary device with DOrif = 255 |im (Up). A region marked by the red dashed line is magnified in the lower Figure. The direction of inner, middle and outer fluid is shown by the blue, green and red arrow, respectively. Supplimentary movie 1.

(t0+ 7.34 ms)

—"TVr^ in

/ —JAy \KJ)

Fig. 3. A qualitative comparison of the numerical (right) and experimental (left) drop generation behaviour. All scale bars are 300 ¡um. Supplementary movie 2.

0 2 4 6 8 02468 10 12 14

Q3/Qsuln Q2/Q1

Q1/Q2 Q2/Q1

- D2 - MC .....D, - MC ♦ D2 - Exp

À Dt - Exp ■ ts

Fig. 4. The size of inner and outer droplets and the shell thickness as a function of phase flow rates: (a) Qi = 2 ml/h, Q2 = 4 ml/h; (b) Qi = 1.5 ml/h, Q3 = 25 ml/h; (c) Q2 = 3 ml/h, Q3 = 12 ml/h; (d) Comparison of experimental data with theoretical predictions. MC in the legend stands for mass conservation.

q°0.8 O

t 6 § 4 2 0

■ D, ♦ D-

250 200 nso '100 50 0

♦ fl - £

5 10 15 20 25 30 Q^Qsum

0 5 10 15 20 25 30 Qî/Qmm

• ♦ L|

f • 1*2 1

0 5 10 15 20 25 30 Q.VQSum

5 10 15 20 25 30

Q3/Qsum

Fig. 5. The variation of flow field with Q3: a) Qs/Qsum = 1.44; b) Qs/Qsum = 5.17; c) Qs/Qsum = 11.48; d) Qs/Qsum = 17.22. Constant parameters: p1 = 1180 kg/m3, p2 = 1170 kg/m3, p3 = 1200 kg/m3, ^ = 0.0396 Pa • s, |2 = 0.0648 Pa • s, |3 = 0.0482 Pa • s, a12 = 0.00574 N/m, a23 = 0.0137 N/m, Q1 = 0.12 ml/h, Q2 = 0.92 ml/h, DOrif = 300 |m.

Fig. 6. The shear stress at the outer interface at the orifice of the collection capillary as a function of Q3. Constant parameters: = 1180 kg/m3, p2 = 1170 kg/m3, p3 = 1200 kg/m3, p1 = 0.0396 Pa • s, p2 = 0.0648 Pa • s, p3 = 0.0482 Pa • s, ct12 = 0.00574 N/m, g23 = 0.0137 N/m, Q1 = 0.12 ml/h, Q2 = 0.92 ml/h, DOrif = 300 pm.

♦ ♦

♦ ♦ ♦

■ ■ ■ ■ ■ ■ ■ ■

- + + + + + + + +

++ + 1

Fig. 7. The variation of flow field with Q2: a) Q2/Qi = 1.34; b) Q2/Qi = 2.69; c) Q2/Qi = 6.72. Constant parameters: p1 = 1180 kg/m3, p2 = 1170 kg/m3, p3 = 1200 kg/m3, ^ = 0.0396 Pa • s, |2 = 0.0648 Pa • s, |3 = 0.0482 Pa • s, a12 = 0.00574 N/m, a23 = 0.0137 N/m, Q1 = 0.46 ml/h, Q3 = 6 ml/h, DOrif = 200 |m. The keys are the same as in Fig. 5.

0 0.2 0.4 0.6 0.8 Q1/Q2

0 0.3 0.6 0.9 1.2 Q1/Q2

0 0.2 0,4 0.6 0,8 Q1/Q2

0 0.3 0.6 0.9 1.2

Fig. 8. The variation of flow field with Q1 in dripping mode: a) Q1/Q2 = 0.0991; b) Q1/Q2 = 0.232; c) Q1/Q2 = 0.661; d) Q1/Q2 = 1.156. Constant parameters: p1 = 1180 kg/m3, p2 = 1170 kg/m3, p3 = 1200 kg/m3, ^ = 0.0396 Pa • s, |2 = 0.0648 Pa • s, |3 = 0.0482 Pa • s, a12 = 0.00574 N/m, a23 = 0.0137 N/m, Q2 = 0.92 ml/h, Q3 = 2.4 ml/h, DOrif = 300 |m. The keys are the same as in Fig. 5.

i 2 Q 1-5 Q 1

0 0.2 0.4 0.6 0.8 0 0.2 0.4 0.6 0.8

Q,/Q2 Qi'Qz

0 0.2 0.4 0.6 0.!

0 0.2 0.4 0.6 O.i Q1/Q2

Fig. 9. The variation of flow field with Q1 in jetting mode: a) Qi/Q2 = 0.0248; b) Qi/Q2 = 0.112; c) Qi/Q2 = 0.372; d) Qi/Q2 = 0.744. Constant parameters: p1 = 1180 kg/m3 , p2 = 1170 kg/m3 , p3 = 1200 kg/m3 , ^ = 0.0396 Pa • s , |2 = 0.0648 Pa • s , |3 = 0.0482 Pa • s , a12 = 0.00574 N/m , a23 = 0.0137 N/m , Q2 = 2.46 ml/h, Q3 = 6 ml/h, DOrif = 200 |m. The keys are the same as in Fig. 5.

0.025 0.05 0.075 0.1 o23( N/m)

0.025 0.05 0.075 0.1

Fig. 10. The variation of flow field with g23 : a) g23 = 0.0057 N/m; b) a23 = 0.02 N/m; c) a23 = 0.06 N/m; d) a23 = 0.095 N/m. Constant parameters: p1 = 1180 kg/m3 , p2 = 1170 kg/m3 , p3 = 1200 kg/m3 , ^ = 0.0396 Pa • s , |2 = 0.0648 Pa • s , |3 = 0.0482 Pa • s, a12 = 0.00574 N/m Q1 = 0.27 ml/h , Q2 = 2.46 ml/h , Q3 = 15 ml/h, DOrif = 200 |m. The symbols on the graphs are the same as in Fig. 5.

0 0.02 0.04 0.06 0.08 *12 (N/m)

0 0.02 0.04 0.06 0.08 <r12 (n/m)

Fig. 11. The variation of flow field with g12 : a) g12 = 0.01033 N/m; b) g12 = 0.015 N/m; c) g12 = 0.03 N/m; d) a12 = 0.06 N/m. Constant parameters: p1 = 1180 kg/m3, p2 = 1170 kg/m3, p3 = 1200 kg/m3, |1 = 0.0396 Pa • s, |2 = 0.0648 Pa • s, |3 = 0.0482 Pa • s, a23 = 0.0137 N/m, Q1 = 0.27 ml/h, Q2 = 2.46 ml/h, Q3 = 15 ml/h, DOrif = 200 |m. The keys are the same as in Fig. 5.

0 0.02 0.04 0.06 0.08 H3 (Pa s)

0 0.02 0.04 0.06 0.08 H3 (Pa s)

0 0.02 0.04 0.06 0.08 lh (Pa-s)

0 0.02 0.04 0.06 0.08 m (Pas)

Fig. 12. The variation of flow field with |3: a) |3 = 0.007 Pa.s; b) |3 = 0.02412 Pa.s; c) |3 = 0.072 Pa.s; d) |3 = 0.08 Pa.s; Constant parameters: p1 = 1180 kg/m3 , p2 = 1170 kg/m3 , p3 = 1200 kg/m3 , = 0.0396 Pa • s, |2 = 0.0648 Pa • s, a12 = 0.00574 N/m, a23 = 0.0137 N/m, Q1 = 0.27 ml/h, Q2 = 2.46 ml/h, Q3 = 2.7 ml/h , DOrif = 200 |m. The keys are the same as in Fig. 6.

0 0.025 0.05 0.075 0.1 H2 (Pa s)

20 16 112

0 0.025 0.05 0.075 0.1 H2 (Pa s)

0.025 0.05 0.075 0.1 (Pa'S)

0 0.025 0.05 0.075 0.1 |i2 (Pas)

Fig. 13. The variation of flow field with |2: a) |2 = 0.002 Pa • s at t0 +0 ms; b) |2 = 0.002 Pa • s at t0 + 1.3 ms; c) |2 = 0.002 Pa • s at t0 +4.3 ms; d) |2 = 0.005 Pa • s; e) |2 = 0.05 Pa • s; f) |2 = 0.09 Pa-s. Constant parameters: p1 = 1180 kg/m3 , p2 = 1170 kg/m3 , p3 = 1200 kg/m3 , = 0.0396 Pa-s , |3 = 0.0482 Pa-s , a12 = 0.00574 N/m, a23 = 0.0137 N/m , Q1 = 0.27 ml/h , Q2 = 2.46 ml/h , Q3 = 2.7 ml/h , DOrif = 200 |m. The keys are the same as in Fig. 5.

0 0.05 0.1 0.15 0.2 Ml (Pa s)

0 0.05 0.1 0.15 0.2 Hi (Pa s)

Fig. 14. The variation of flow field with ^: a) ^ = 0.001 Pa • s; b) ^ = 0.007 Pa • s; c) ^ = 0.0198 Pa • s; d) = 0.06 Pa • s. Constant parameters: p1 = 1180 kg/m3, p2 = 1170 kg/m3, p3 = 1200 kg/m3, |2 = 0.06482 Pa • s, |3 = 0.0482 Pa • s, a12 = 0.00574 N/m, a23 = 0.0137 N/m, Q1 = 0.27 ml/h, Q2 = 2.46 ml/h, Q3 = 2.7 ml/h, DOrif = 200 |m. The keys are the same as in Fig. 5.

1800 kg/m3. Constant parameters: p1 = 1180 kg/m3, p2 = 1170 kg/m3, ^ = 0.0396 Pa • s, |2 = 0.06482 Pa - s, |3 = 0.0482 Pa • s, a12 = 0.00574 N/m, a23 = 0.0137 N/m, Q1 = 0.27 ml/h, Q2 = 2.46 ml/h, Q3 = 2.7 ml/h, DOrif = 200 |m. The keys are the same as in Fig. 5.

Fig. 16. The variation of flow field with DOrif at high flow rates: a) DOrif= 150 |m; b) DOrif= 300 |m. Constant parameters: px = 1180 kg/m3, p2 = 1170 kg/m3, p3 = 1200 kg/m3, |x = 0.0396 Pa • s, |2 = 0.06482 Pa • s, |3 = 0.0482 Pa • s, a12 = 0.00574 N/m, a23 = 0.0137 N/m, Qx = 0.27 ml/h, Q2 = 2.46 ml/h, Q3 = 15 ml/h, DN = 25 |m.. The variation of flow field with DOrif at low flow rates: c) DOrif= 150 |m; d) DOrif = 300 |m. Constant parameters: px = 1180 kg/m3, p2 = 1170 kg/m3, p3 = 1200 kg/m3, |x = 0.0396 Pa • s, |2 = 0.06482 Pa • s, |3 = 0.0482 Pa • s, a12 = 0.00574 N/m, a23 = 0.0137 N/m, Qx = 0.12 ml/h, Q2 = 0.92 ml/h, Q3 = 2.4 ml/h, DN = 25 |m. The keys are the same as in Fig. 5.

Fig. 17. The variation of flow field with L' at high flow rates (Q1 = 0.27 ml/h, Q2 = 2.46 ml/h, Q3 = 15 ml/h): a) L' = 90 |m; b) L' = 250 |m. The variation of flow field with L' at low flow rates (Qx = 0.12 ml/h, Q2 = 0.92 ml/h, Q3 = 2.4 ml/h): c) L' = 90 |m; d) L' = 250 |m. Constant parameters: p1 = 1180 kg/m3, p2 = 1170 kg/m3, p3 = 1200 kg/m3, ^ = 0.0396 Pa • s, |2 = 0.06482 Pa • s, |3 = 0.0482 Pa • s, a12 = 0.00574 N/m, a23 = 0.0137 N/m, DOrif = 300 |m. The keys are the same as in Fig. 5.

Highlights

- Compound droplet formation in glass capillary devices was parametrically studied.

- Good agreement was achieved between experiments and VOF-CSF numerical model.

- Model reproduced well dripping, narrowing jetting, widening jetting and transitions.

- Effects of fluid properties, flow rates and geometry on droplet size were explored.

- Results are applicable for encapsulation of CO2 solvents in core/shell microcapsules.

Graphical Abstract

(for review)