Double emulsion production in glass capillary micro ﬂ uidic device: Parametric investigation of droplet generation behaviour

(cid:1) Compound droplet formation in glass capillary devices was parame-trically studied. (cid:1) Good agreement was achieved between experiments and VOF – CSF numerical model. (cid:1) Model reproduced well dripping, narrowing jetting, widening jetting and transitions. (cid:1) Effects of ﬂ uid properties, ﬂ ow rates and geometry on droplet size were explored. (cid:1) Results are applicable for encapsula- tion of CO 2 solvents in core/shell microcapsules. 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- ﬂ ow and countercurrent ﬂ ow 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 ﬂ ow rates, ﬂ uid properties, and geometry on the size, morphology, and production rate of droplets. As the outer ﬂ uid ﬂ ow rate increased, the size of compound droplets was reduced until a dripping-to-jetting transition occurred. By increasing the middle ﬂ uid ﬂ ow rate, the size of compound droplets increased, which led to a widening jetting regime. The jetting was supressed by increasing the ori ﬁ ce 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 CO 2 solvents within gas-permeable microcapsules.


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 Bergnståhl, 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 CO 2 solvents within a CO 2 permeable polymer are a promising alternative to conventional liquid CO 2 capture materials such as monoethanolamine (Aines et al., 2013). The most challenging step in the production of tailor-made CO 2 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 CO 2 and the CO 2 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. (2005Utada et al. ( , 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 countercurrently 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 drippingto-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.
To the best of our knowledge, this is the first numerical simulation study of double emulsion formation in an axisymmetric threephase 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 CO 2 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.
where U ! is the velocity vector, P is the pressure, and t, μ and ρ are the time, dynamic viscosity and density, respectively. A source term is denoted by F b 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 F b ¼ 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).
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 dynamic viscosity and density in the momentum equation, Eq.
(2) were calculated as follows: 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 σ is the interfacial tension and κ is the local curvature of the interface calculated as follows: wheren is the unit normal defined aŝ

Numerical method
The governing equations were solved by utilizing an unsteady pressure-based segregated algorithm established in finite-volume based commercial software Ansys s Fluent v. 14.0. The discretized 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 pressurevelocity 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, 4, 3, 2 and 1 mm. No difference in behaviour was found for the resolutions of 1 À4 mm. 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 mm, 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 summarized in Table 1.

Dimensionless numbers
Non-dimensional numbers used in this study are Weber numbers of inner phase, We 1 , and that of compound jet, We 2 , and capillary numbers of inner, Ca 1 , middle, Ca 2 and outer phase, Ca 3 where ρ, σ, μ, 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, compound jet, and nozzle, respectively. The values of u jet , u 2 , u 3 , and D jet were directly measured in the vicinity of exit orifice from the numerical data, while u 1 was calculated as u 1 ¼ 4Q 1 = πD 2 N . Since the difference in density between the inner and middle phase is very small, the outer fluid density is used for estimation of We 2 .

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 s 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, M w ¼13,000-23,000 g mol À 1 , 87-89% hydrolyzed, 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 Krüss DSA-100 pendant drop tensiometer.

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 s 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 s 20 G, 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.

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.86 mm I.D. and 1.52 mm 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.

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.

Experimental results
The effect of phase flow rates on the inner droplet diameter, D 1 , the outer droplet diameter, D 2 , and the shell thickness, t s , is presented in Fig. 4. Solid and dashed lines in Fig. 4a-c are theoretical predictions based on observed droplet generation frequencies and flow rate settings. As can be seen in Fig. 4c, D 1 can be varied by tuning the flow rate ratio Q 1 =Q 2 , where an increase in Q 1 =Q 2 leads to a significant increase in D 1 , while D 2 remains almost constant, which results in thinning the shell. In contrast, an increase in Q 2 =Q 1 increases D 2 and reduces D 1 and consequently, the shell becomes thicker, as shown in Fig. 4b. An increase in Q 3 considerably reduces both D 2 and D 1 while t s remains nearly constant (Fig. 4a).  gives a comparison between the experimental results and Eq. (10) As can be seen in Fig. 4d, 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 (D 2 , D 1 and t s ) can be affected by the viscosity ratios of the fluids. For example, Chang and Su (2008) observed an increase in D 1 and D 2 as Q 1 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.

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 CO 2 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-9.
4.2.1. Increasing Q 3 reduces both D 1 and D 2 The effect of outer fluid flow rate, Q 3 , on droplet formation is illustrated in Fig. 5. In the dripping regime (Q 3 = Q sum o 5:17), compound droplets are formed very close to the entry of the collection capillary (L 2 =D Orif 1:2, Fig. 5f) and the size of both inner and outer droplets is significantly reduced with increasing Q 3 (Fig. 5a, b), which is coupled with an increase in both droplet generation frequencies (Fig. 5g). Q 3 has a higher impact on D 2 than D 1 , which leads to the reduction in shell thickness, t s . A short-term increase in the size of inner and outer drops at Q 3 = Q sum 7 in Fig. 5e was due to transition from dripping to narrowing jetting regime, but both D 1 and D 2 resumed its downward trend at Q 3 = Q sum 4 7:18. A temporary increase in the drop size can be explained by the variations of Ca 3 and We 2 presented in Fig. 5h. At the transition point, although there is a considerable decrease in D jet , the steep jump in We 2 implies a sudden increase in the velocity of compound jet, while there is a constant increase in Ca 3 . 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=f 1 and 1=f 2 (Fig. 5g). Since Q sum ¼const, an increase in droplet formation times results in larger D 1 and D 2 values.
The transition to jetting regime leads to a dramatic increase in break up lengths (Fig. 5c, d), with L 2 =D Orif reaching 8.7 at Q 3 =Q sum ¼26 (Fig. 5f). 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 Q 3 =Q sum ¼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, wh ich reduces the drop formation time and thus, increases droplet generation frequencies, f 1 and f 2 to almost 200 Hz at Q 3 =Q sum ¼ 26 (Fig. 5g). Since f 1 ¼ f 2 at all flow rates shown in Fig. 5, each compound droplet contains a single internal drop. A shear stress at the outer interface, τ 2;3 , is as follows (see the Supplemental material S1): where α is given by Eq. (S.12). The variation of τ 2;3 with Q 3 calculated from Eq. (11) is shown in Fig. 6. An increase in Q 3 leads to an increase in τ 2;3 , but the effect is less pronounced at Q 3 =Q sum 420.  There is an abrupt increase in L 1 and L 2 with Q 3 at the transition point. The further increase in Q 3 leads to additional elongation and thinning of the compound jet. At Ca 3 40.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 Q 3 and t s =D Orif reaches 0.12 at Q 3 =Q sum ¼26, which corresponds to t s ¼ 36 mm. Based on the fluid properties and the simulated geometry, the core-shell droplets are successfully produced at 1.4oQ 3 =Q sum o 57.9. At Q 3 =Q sum 1.4, the outer flow cannot focus the compound jet through the collection tube and the middle phase wets the collection tube wall. At Q 3 =Q sum 457.9, the outer phase pushes the compound jet back into the injection nozzle and the droplet formation fails.
4.2.2. Increasing Q 2 reduces D 1 while increasing D 2 An increase in Q 2 results in higher D 2 values, while D 1 is reduced and accordingly, the shell becomes thicker (Fig. 7e). An increase in Q 2 results in higher D 2 values, while D 1 is reduced and accordingly, the shell becomes thicker (Fig. 7e). An increase in Q 2 increases the middle fluid velocity and correspondingly the compound jet inertia. The increase in D 2 with Q 2 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 Q 2 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. 7g), leading to the formation of smaller inner droplets. According to the mass balance equation: f 1 ¼ 6Q 1 =ðD 1 3 Þ. Since Q 1 is a constant, f 1 varies inversely as the cube of D 1 . An increase in f 2 in Fig. 7g is a consequence of increase in Q 2 . An increase in Q 2 results in higher We 2 values and longer break-up lengths, L 1 and L 2 , as shown in Fig. 7h and 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 Q 2 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 Q 2 values result in a faster droplet production and higher f 2 values, as shown in Fig. 7g. Based on the fluid properties used in the simulations, the successful core-shell droplet generation is at 0.58 oQ 2 =Q 1 o46.7. For Q 2 =Q 1 o0.58, the shell of compound jet is so narrow that the inner fluid ruptures the compound jet and the droplet formation fails. For Q 2 =Q 1 446.7, the compound jet diameter exceeds the orifice diameter, D Orif and wets the collection tube leading to failure in droplet formation.

4.2.3.
Increasing Q 1 increases both D 1 and D 2 Fig. 8 demonstrates the effect of innermost fluid flow rate, Q 1 , on the drop formation behaviour. The drop formation in the dripping regime occurs over the Q 1 =Q 2 range from 0.13 to 0.66, corresponding to 0.013o We 1 o 0.61 (Fig. 8a-c). An increase in Q 1 leads to an increase in both inertial and viscous forces exerted by the inner fluid, which results in an increase in both We 1 and Ca 1 (Fig. 8h). Since there is no considerable change in L 1 and L 2 over the range of 0:13 oQ 1 =Q 2 o0:66 (Fig. 8g), an increase in the momentum of inner fluid results in bigger inner droplet. In the dripping regime Q 1 does not have any noticeable effect on the shear stress at the outer interface and thus, D 2 increases only slightly with  Q 1 (Fig. 8e). Since the size of inner drop increases with increasing Q 1 and the size of outer drop remains nearly constant, the shell becomes thinner. As Q 1 =Q 2 exceeds a critical value of 0.66, any further increase in Q 1 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. 8d). In addition, for Q 1 =Q 2 o0.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 Q 1 exceeds a certain critical value was also reported earlier (Chang and Su, 2008). f 1 and f 2 increased very slightly with Q 1 , as shown in Fig. 8f.
The influence of Q 1 on drop formation in jetting regime is presented in Fig. 9. As shown in Fig. 9e, D 1 follows the same trend as in dripping regime, but D 2 increases more significantly, albeit at a slower pace than D 1 , so that the shell thickness is reducing. At very low Q 1 =Q 2 value of 0.0248, where We 1 E0.006, inner droplets are formed in a narrowing jetting regime (due to the high Q 2 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. 9a). Formation of inner and outer drops in different regimes results in a considerable difference between L 1 and L 2 , which leads to the generation of multi-cored droplets. As Q 1 increases, the inner drop formation switches from narrowing to widening mode, resulting in virtually the same value of L 1 and L 2 and formation of core/shell droplets. A high speed jet was observed at Q 1 =Q 2 40.62 and We 1 43.83 (Fig. 9d), similar to that shown in Fig. 8d. However, in this case the shell rupture occurs after droplet pinch off, further downstream.
156. Constant parameters: ρ 1 ¼ 1180 kg=m 3 , ρ 2 ¼ 1170 kg=m 3 ; ρ 3 ¼ 1200 kg=m 3 ; μ 1 ¼ 0:0396 Pa s; μ 2 ¼ 0:0648 Pa s; μ 3 ¼ 0:0482 Pa s; σ 12 ¼ 0:00574 N=m, σ 23 ¼ 0:0137 N=m; Q 2 ¼ 0:92 ml=h; Q 3 ¼ 2:4 ml=h; D Orif ¼ 300 μm. The keys are the same as in Fig. 5. interface increases with increasing σ 23 , which results in the reduction of Ca 3 from 0.51 to 0.076 and the reduction in We 2 from 0.322 to 0.031 (Fig. 10h). 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 σ 23 values suppress the jetting mode by pulling back the compound jet and switching the droplet formation to dripping regime, as shown in Fig. 10a-d. Consequently, both break-up lengths, L 1 and L 2 , decrease with increasing σ 23 (Fig. 10f), which is consistent with the findings of Herrada et al. (2010). Two different behaviours have been observed on the graphs in Fig. 10e and g. At σ 23 ¼0.005-0.06 N/m, an increase in σ 23 considerably delays the pinch off events and reduces both f 1 and f 2 , resulting in an increase in D 1 and D 2 . An increase in the size of the outer drop is more significant than that of the inner drop causing t s 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. 10e and   g). Within this range, as σ 23 increases both inner and outer droplets shrink causing the drop generation frequency to increase in order to maintain constant volume flow rates, Q 1 and Q 2 . 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 Q 3 led to a sudden growth in both D 1 and D 2 .

Increasing σ 12 suppresses inner fluid jetting
The effect of interfacial tension at the inner interface, σ 12 , on the drop formation is shown Fig. 11. Starting from σ 12 ¼0.00574 N=m; where both inner and outer droplets are formed in jetting regime, increasing σ 12 causes a reduction in Ca 1 (Fig. 11h), 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 σ 12 ¼0.01 N/m (Fig. 11a), the inner and middle phase jets are of similar length (L 1 EL 2 ). As σ 12 increases to 0.015 N/m, L 2 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 D 1 values ( Fig. 11e and g). Indeed, for 0.005 o σ 12 o 0:03 N=m, both L 1 and L 2 decreases with increasing σ 12 (Fig. 11a-c). This effect is more pronounced for the inner fluid, causing the difference between L 1 and L 2 to increase, which leads to the formation of multi-cored compound droplets (Fig. 11b).
In the σ 12 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 L 2 (Fig. 11c). The further increase in σ 12 over the range from 0.04 to 0.06 N/m leads to reduction in both L 1 and D 1 , where D 1 is considerably smaller than the compound jet diameter. Therefore, the capillary waves previously formed on the outer surface (Fig. 11c) due to large inner droplet are supressed (Fig. 11d), and the compound jet pinch off is delayed leading to an increase in L 2 . It should be mentioned that the droplet formation at high σ 12 values is very chaotic leading to formation of multi-cored droplets along with generation of large satellite droplets (Fig. 11d).

Increasing μ 3 increases both D 2 and D 1
The effect of outer fluid viscosity,μ 3 , on the production of compound drops is shown in Fig. 12. According to Fig. 12h, an increase in μ 3 leads to a linear increase in Ca 3 due to higher viscous force acting on the outer interface, which results in stretching both the inner phase and middle phase jet (Fig. 12f). The higher level of shearing exerted by outer fluid accelerates the droplet formation, which results in higher droplet formation frequencies, f 1 and f 2 , and smaller drop diameters, D 1 and D 2 ( Fig. 12e and g). As can be seen in Fig. 12a,b, the variations of D 1 and D 2 with μ 3 are strongly correlated, leading to almost constant shell thickness.
4.2.7. Increasing μ 2 increases both L 2 and L 1 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. 13a-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. 13a-c. As shown in Fig. 13j, the capillary number of middle phase Ca 2 increases with increasing μ 2 , since Ca 2 μ 2 at constant u 2 and σ 23 . 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 Pa s results in a steep increase in both D 2 and D 1 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, f 2 ¼120 s À 1 , compared to the one in jetting mode, f 2 ¼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.

Reducing μ 1 leads to formation of inner droplets in dripping regime
As illustrated in Fig. 14, D 2 is almost unaffected by the innermost fluid viscosity, μ 1 , due to negligible effect of μ 1 on the shear stress at the outer interface. However, there is a steep increase in D 1 as μ 1 increases from 0.001 to 0.019 Pa s, followed by a sharp drop in f 1 from nearly 600 s À 1 to 60 s À 1 . At μ 1 ¼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: N 1 ¼ f 1 =f 2 . E.g., at μ 1 ¼0.001 Pa s, the number of inner droplets observed is 31 and f 1 =f 2 ¼ 31:48, whereas at μ 1 ¼0.004 Pa s, three inner droplets are engulfed by an outer droplet and f 1 =f 2 ¼ 3:39. According to , when the predicted value of N 1 is between two integers, N 1 can take either integral value. An increase in μ 1 to 0.0396 Pa s causes an inner jet to develop along with a dramatic growth in D 1 , while the outer jet is shortening, and core/shell droplets are formed ( Fig. 14a-d, f). The reduction in L 2 can be attributed to the fact that at low μ 1 values, fast and sequential formation of inner droplets delays the drop pinch off from the compound jet and results in larger L 2 . By increasing μ 1 and consequently decreasing the number of inner droplets trapped by shell, the break-up of compound jet facilitates, resulting in shorter L 2 . The further increase in μ 1 beyond μ 1 ¼0.0396 Pa s results in an increase in Ca 1 , corresponding to higher viscous forces of the inner fluid and a slight increase in both L 1 and L 2 . However, there is no significant change in f 1 and f 2 and therefore D 1 and D 2 .

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/m 3 . It could be explained by the fact that due to low fluid velocity and low density ratio, the convection term of the 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.2.10. Decreasing D orif changes droplet formation regime from dripping to jetting Fig. 16 shows the effect of orifice size, D Orif , on droplet formation. In Fig. 16a, b, Q 2 and Q 3 are relatively high and the drop formation occurs in narrowing jetting regime. Here, an increase in D Orif from 150 to 300 μm, corresponding to an increase in D Orif =D N from 6 to 12, causes a decrease in both L 1 and L 2 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 D 1 and D 2 by 4-8% (Table S. 2, jetting).
The effect of D Orif on the droplet formation at low flow rates is shown in Fig. 16c, d. At D Orif ¼ 150 μm (Fig. 16c), 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 D Orif ¼300 μm (Fig. 16d), 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.2.11. Increasing L 0 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 0 increases at high flow rates (Fig. 17a,b), D 2 decreases while D 1 is almost unaffected, leading to a decrease in the shell thickness (Table S.3). The increased distance L 0 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 D 2 , L 1 , and L 2 .
Although L 0 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. 17c,d). In Fig. 17c, 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 0 =D Orif lies between 0.3 and 0.83.

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   15. The variation of flow field with ρ 3 . (a) ρ 3 ¼ 600 kg/m 3 ; (b) ρ 3 ¼1000 kg/m 3 and (c) ρ 3 ¼1800 kg/m 3 . Constant parameters: ρ 1 ¼ 1180 kg=m 3 , ρ 2 ¼ 1170 kg=m 3 , μ 1 ¼ 0:0396 Pa s, μ 2 ¼ 0:06482 Pa s, μ 3 ¼ 0:0482 Pa s, σ 12 ¼ 0:00574 N=m, σ 23 ¼ 0:0137 N=m, Q 1 ¼ 0:27 ml=h, Q 2 ¼ 2:46 ml=h, Q 3 ¼ 2:7 ml=h, D Orif ¼ 200 μm. The keys are the same as in Fig. 5. 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, Q 3 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 Q 3 and μ 3 . An increase in the middle fluid flow rate, Q 2 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 Q 2 value. As the inner fluid flow rate, Q 1 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. Fig. 16. The variation of flow field with D Orif at high flow rates. (a)D Orif ¼ 150 μm and (b) D Orif ¼ 300 μm. Constant parameters: ρ 1 ¼ 1180 kg=m 3 , ρ 2 ¼ 1170 kg=m 3 , ρ 3 ¼ 1200 kg=m 3 , μ 1 ¼ 0.0396 Pa s, μ 2 ¼ 0.06482 Pa s, μ 3 ¼0.0482 Pa s, σ 12 ¼ 0:00574 N=m, σ 23 ¼ 0:0137 N=m, Q 1 ¼ 0:27 ml=h, Q 2 ¼ 2:46 ml=h, Q 3 ¼ 15 ml=h, D N ¼ 25 μm. The variation of flow field with D Orif at low flow rates: (c) D Orif ¼ 150 μm; (d)D Orif ¼ 300 μm. Constant parameters: ρ 1 ¼ 1180 kg=m 3 , ρ 2 ¼ 1170 kg=m 3 , ρ 3 ¼ 1200 kg=m 3 , μ 1 ¼ 0:0396 Pa s, μ 2 ¼ 0:06482 Pa s, μ 3 ¼ 0:0482 Pa s, σ 12 ¼ 0:00574 N=m, σ 23 ¼ 0:0137 N=m, Q 1 ¼ 0:12 ml=h, Q 2 ¼ 0:92 ml=h, Q 3 ¼ 2:4 ml=h, D N ¼ 25 μm. The keys are the same as in Fig. 5. Fig. 17. The variation of flow field with L ' at high flow rates (Q 1 ¼ 0:27 ml=h, Q 2 ¼ 2:46 ml=h, Q 3 ¼ 15 ml=h): (a) L ' ¼90 μm; (b) L ' ¼ 250 μm. The variation of flow field with L ' at low flow rates (Q 1 ¼ 0:12 ml=h, Q 2 ¼ 0:92 ml=h, Q 3 ¼ 2:4 ml=h): (c) L ' ¼ 90 μm and (d) L ' ¼250 μm. Constant parameters: ρ 1 ¼ 1180 kg=m 3 , ρ 2 ¼ 1170 kg=m 3 , ρ 3 ¼ 1200 kg=m 3 , μ 1 ¼ 0:0396 Pa s, μ 2 ¼ 0:06482 Pa s, μ 3 ¼ 0:0482 Pa s, σ 12 ¼ 0:00574 N=m, σ 23 ¼ 0:0137 N=m, D Orif ¼ 300 μm. The keys are the same as in Fig. 5