Numerical Investigation of Double Emulsion Formation in non-Newtonian Fluids using Double Co-Flow Geometry

The paper presents a 2D axisymmetric numerical simulation in three phases to investigate a double Co-Flow microfluidic device's ability to produce double emulsions in both Newtonian and non-Newtonian ambient fluids. The Volume of fluid (VOF) method was utilized to perform an investigation of the creation of a double emulsion in a double Co-Flow geometry. The study utilized a model to examine how the size and generation frequency of double emulsions are impacted by various factors such as the velocity of the phases, viscosities, interfacial tension, and rheological properties of non-Newtonian fluids. The model predicted the process of emulsification successfully in dripping and jetting


Introduction
Double emulsions are a type of dispersions that involve three phases, where a droplet of one liquid is surrounded by a shell of another liquid, and then this combination is dispersed into another continuous phase.Double emulsions have a middle phase that acts as a protective barrier or semi-permeable shell, separating the inner droplet from the continuous phase.This unique structure makes double emulsions useful in various applications such as controlled drug-delivery systems (Lao et al. 2009;Chong et al. 2015;Dluska et al. 2019) food science (Matos et al. 2014;Aditya et al. 2015;Nelis et al. 2019), pharmaceutical (Li and Jiang 2018), ultrasonic contrast agents (Klibanov 2005;Levenstein et al. 2016), singlecell screening (Brouzes et al.), biotechnology (Lu et al. 2015;Yildirim et al. 2016) and cosmetics industries (Yoshida et al. 1999;Matos et al. 2013;Stasse et al. 2020).
In the past two decades, a variety of microfluidic methods have emerged for producing droplet-in-droplet structures, particularly double or multiple emulsions, which have gained widespread adoption.This section provides a concise overview of the different techniques employed to fabricate multiple emulsions.
Various approaches have been devised for generating double or multiple emulsions in microfluidic devices, which can be classified into two types: the "one-step" technique and the "two-step" technique.The choice of technique depends on the flow properties and speeds of the different phases, and is determined by where the inner phase is broken up.
Applying a shear force Leading to the creation of double emulsions in a single step.This is achieved by injecting the inner and middle fluids through a capillary tube in one path while introducing the outer phase via the outer coaxial area in the reverse path.
Abate et al. (Abate et al. 2011) presented a double flow-focusing microfluidic device that was able to make multiple emulsions with controllable shell thickness in one-step.Their one-step configuration could generate very thin-shelled thickness that is useful for capsule synthesis applications.Nabavi et al. (Nabavi et al. 2017) presented another one-step configuration with the potential to create multi-core double emulsion droplets.They discovered that the outer phase's capillary number by lowering it or by raising the middle phase's capillary number can result in compound droplets with multiple cores.Although a one-step configuration has some advantages, including the creation of high monodisperse compound droplets in a high frequency of generation and do not require partial wettability like two-step methods, they are expensive and have an arduous procedure of fabrication in comparison with soft lithography methods (Morita et al. 2016).Akamatsu et al. (Akamatsu et al. 2015) demonstrated a one-step approach for manufacturing W/O/W forms utilizing membrane integrated glass capillary.The key aspect of the integrated device is its capacity to reduce the dimension of droplets when the emulsions created at the capillary move through the membrane.During this process, the majority of the oil shell layer is removed as well, resulting in double emulsion droplets with thin shells.The process of making double emulsions can be segmented into three stages: first, creating the middle oil phase; second, allowing the pinch-off process to take place and the inner phase to penetrate the middle phase; and third, breaking up the resulting double emulsion structure.Nam et al. (Nam et al. 2016) utilized a similar structure in their work, in which they described an effective microfluidic technique for encapsulating Pt nanoparticle stabilized by polyvinylpyrrolidone (PVP) in photocurable double-emulsion droplets with semipermeable thin shells.Arriaga et al. (Arriaga et al. 2015) devised a revolutionary one-step configuration for manufacturing double emulsions, using two successive flow-focusing devices to create double emulsions in a single step.The proposed structure can generate double emulsions with a slender shell thickness of approximately 10% of the diameter.
As mentioned, the main applications of multiphase and microfluidic flows are medicine, pharmacy, drug delivery, and imaging.Given that blood or non-Newtonian fluids in general are used in the majority of the aforementioned applications, it is necessary to study microfluidic and multi-phase flows by considering the non-Newtonian properties of fluids.
Blood, as the most important non-Newtonian fluid for medical engineering researchers, has a shear thinning behavior especially in small conduit (Sattari and Hanafizadeh 2019), which is why researchers are looking for fluids that follow this behavior.The generation of compound droplets in non-Newtonian fluids is crucial for their potential applications in biological contexts.Although numerous studies have concentrated on the creation or motion of non-Newtonian compound droplets or jets in Newtonian fluids (Toose et al. 1999;Mohsin et al. 2012); the investigation of this constructions in non-Newtonian fluids is critical given their promising applications.However, to our best knowledge, thus far, no research has conducted a comprehensive examination on the production of double emulsions in a double co-flow geometry in non-Newtonian fluids.Therefore, with the aim of bridging this research gap, this study analyzes the production of double emulsions in both non-Newtonian and Newtonian ambient fluids utilizing a Double Co-Flow Geometry.
Our research has specific goals of examining how the continuous phase's rheological parameters impact the under-controlled manufacturing and manipulation of double emulsions.Based on the information we have at present, this study is the first to simulate the rheological impacts of the continuous phase on double emulsion production using a double co-flow geometry, which could be a fundamental step in the use of compound droplets in fields such as biotechnology and drug delivery.

Governing equations
The governing equations are Eq.(1) for mass conservation and Eq.(2) for momentum conservation: where  , U ⃗ ⃗ , ,  and  represent density, velocity vector, time, pressure, and viscosity, respectively.The symbol F  stands for surface tension, which operates at the border of two immiscible fluids.Given the small length scale of a few microns, the gravitational force has been ignored.
The Carreau viscosity model can be displayed by the following equation (Byron and Carreau 1968;Carreau et al. 1968) (3): where  is the quantified viscosity, Υ ̇ is the applied shear rate, η 0 is the viscosity limit when Υ ̇ → 0, and η ∞ is the viscosity at unlimited shear rate, which was defined as the solvent viscosity.
The parameter  is expressed in units of time, and the constant value parameter n is acquired for each non-Newtonian fluid using curve fitting.
The VOF approach solves a supplementary continuity equation to monitor the interface for the volume fraction of phase .The phase advection is represented by α p , as follows: The volume fraction of a particular phase is utilized to determine the fraction of each computational cell occupied by that phase.
The model being presented here involves the use of subscripts 1 to represent the internal phase, 2 to represent the middle phase, and 3 to represent the outer phase.In this model, the equivalent dynamic viscosity and density can be calculated using the following relationships: The diffuse interface representation is used to calculate the interfacial surface force, F  , by utilizing the CSF (Continuum Surface Force) method (Brackbill et al. 1992): In the given equation, σ represents the surface tension, and κ represents the interface's local curvature, the calculation for this is as follows: where  ̂ is the unit normal specified as,

Non-dimensional numbers
In this study, the impact of phase flow rates was assessed by analyzing three parameters: the Weber number (We), which compares inertia forces to surface tension, the capillary number (Ca), which measures the relative influence of viscous forces versus interfacial forces, and the Reynolds number (Re), which determines the ratio of inertial to viscous forces.Reynolds number is estimated to describe non-Newtonian fluid flows using the following relation (Ohta et al. 2020): where  denotes the equivalent diameter of the channel and  12 and  23 are respectively surface tension between inner and middle and middle and outer phases, and u is superficial velocity, which is considered equal to the inlet velocity of phases.The subscribes "in," "mid," and "out" are used in reference to the inner, middle, and outer phases, respectively.

Geometric design, fluid properties, and modeling procedure
The microfluidic device's geometry, which includes three nested capillary tubes, mesh   with other solution properties, are tabulated in Table 1.The study examined different physical properties, such as viscosity and surface tension between the inner and middle and middle and outer phases, along with the velocities of the three-phase inlet, which were considered as input parameters.Table 2 represents a list of the variety of parameters that were examined in this study.Notice that phase densities were disregarded since it was determined that they do not significantly impact the creation process and properties of double emulsions (Nabavi et al. 2015).axisymmetric laminar model.To achieve this, we utilized a time-dependent pressure-based segregated algorithm, which is based on the finite volume technique.We utilized the model to conduct a study that involved analyzing various parameters related to the formation of compound droplets.The QUICK technique was used to estimate the discretized momentum equation to give the most accurate result.The PRESTO approach, which directly calculated the pressure term on the interfaces to eliminate interpolation mistakes, was used to achieve the pressure term interpolation (Nabavi et al. 2015).The Semi-implicit method for pressure linked equations (SIMPLE) scheme is used to accomplish the coupling of the pressure and velocity while the pressure interpolation is achieved by applying the Geo-Reconstruction algorithm.To lower the computational cost, an adjustable time step technique was used employing the max Courant number,   , with the specified quantity of 0.35.Constant inlet velocities were provided for all three inlets.The pressure outlet was used as the outlet boundary.As the wettability of the external wall had no impact on the encapsulated structure in the intended geometry, the non-wetted wall circumstances were employed to establish all other boundaries.

Grid convergence test
To study the effects of lattice resolution on the interfaces between immiscible fluids and mechanism of droplet breakup, a mesh investigation with elements of four different sizes was carried out.The diameter of the inner droplet used as the benchmark for the mesh analysis.Fig. 3 shows the four mesh types mentioned above along with their corresponding volume fraction contours, for the corresponding process circumstance of   = 0.003,   = 0.6, and   = 30.As noted by previous researchers (Bashir et al. 2011;Liu et al. 2016), the size of the elements used does not have a notable influence on the interface or the mechanism of breakup.Table 3 provides details about various meshes used in the study, including their statistics and the inner droplet diameter as shown in Fig 3 .This table confirms the claim made in the study that mesh independence is achieved, as the inner droplet diameter only shows a 1% variation between the smallest and greatest mesh resolutions.The following numerical studies were performed on Mesh 3, which has 33,747 elements in total because of its low dependence on topological behavior of the encapsulation process and low computational cost.

Model validation
Although numerous studies have confirmed the VOF model's accuracy, we compare our model with Shao et al. (Shao et al. 2013) in both qualitative and quantitative ways.The two main components of their double emulsion generation were three co-axial capillaries and an assisting body.In this arrangement, the internal fluid is infused into the middle capillary through the innermost capillary tube, and the two capillaries are then combined to inject the inner fluid into the outer capillary.Dripping instability then causes the double emulsion to form downstream of the outer capillary.Our numerical simulation successfully predicted the process of forming compound droplets before and after breakup, as shown in

Dripping regime
The axisymmetric geometry is used to simulate the creation of double emulsions through two co-axial arrangements of the same level.During the dripping regime, the inner fluid is infused into the developing droplet of the middle phase through the internal channel, resulting in the formation of a liquid filament of the inner phase as seen in

Jetting regime
As the dispersed phase's flow rate enhances in the microchannel, there is a transition in droplet formation from a dripping mode to a jetting mode.The fluid filament split away downstream of the channel during the jetting regime due to jet instability.It is commonly understood that the dripping regime occurs when the continuous phase's capillary number is small Because the force caused by viscosity cannot dominate surface tension .The jetting regime develops as the capillary number for continuous phase increases, making the force of viscous shear the prevailing force.The process used to make double emulsions is similar to the dripping regime.In other words, the inner phase of the jet stream infuses the    As shown in Fig. 8, when the value of "  " is enhanced, there is a corresponding increase in the size of the inner droplet (Fig 8a).However, the effect of Weber number on the size of the outer droplet is not as remarkable as it is for the inner droplet, and the size of the outer droplet remains approximately constant (Fig 8b).In addition, as "  " increases, the gap between the inner and outer droplet sizes reduces, resulting in a decline in the thickness of the droplet (Fig 8c), which is defined as (  −   )/2.Therefore, structures with a thick or thin shell can be achieved in a significant range only by changing the flow rate of the internal phase.Figs 8 & 9 show that the inner and outer diameters of the droplets generated in the Newtonian fluid (water) are greater than those of the selected non-Newtonian fluids, and that these diameters decrease as the concentration of CMC solutions rises.As is well known, the generated emulsions have the lowest diameter in the 0.5% Wt CMC solution.As a result, the generated double emulsions will have the largest thickness in Newtonian fluid, the lowest thickness in 0.5% Wt CMC solution, and the thickness of emulsions falls as the CMC solution's concentration enhances.As it is clear from Fig 8d , with the enhance of internal phase flow rate, the frequency of double emulsion structure formation increases with a very small slope.This means that more double emulsions are formed by rising the flow rate of the internal phase over a certain period of time.Also, the frequency of the formation of double emulsions in non-Newtonian fluids is higher than in Newtonian fluid, because in non-Newtonian fluids, the diameter of the droplets is smaller and the separation of the droplets occurs earlier.In many applications, such as modern drug delivery systems, where numerous drops are needed quickly, this condition is highly desirable.

Effect of inner phase flow rate
Fig 9 shows a series of images captured just before and after the point at which the liquid breaks up.As the Weber number (  ) increases, the liquid takes on a shape more resembling a jet and the location at which the liquid breaks up and separates into smaller droplets moves downstream.As the Weber number enhances, it is evident that the role of inertia becomes relatively more significant in contrast to the surface tension force.This shift in balance leads to an enhance in the length of the liquid thread, mainly because of the enhance in inertial forces.Consequently, this elongates the liquid cone that remains on the inner capillary tube and leads to a decline in the size of the resulting inner droplet.The

Effect of middle phase flow rate
A rise in middle phase velocity outcomes in lower D in and greater D out (Fig 10a and Fig 10b,respectively).Obviously, the shell thickness increases as well (Fig 10c).The relative effects of the force of surface tension, which tends to reduce the droplet surface, and the drag force caused by viscosity, which tends to deform the inner droplets using middle phase friction, increase as the capillary number of the middle phase increases.Therefore, the shear stress that is applied on the inner droplet increases in a similar way, so that the droplet separation process occurs faster and the dimensions of the inner droplets are reduced.An enhance in middle phase velocity also outcomes in a greater velocity gradient, and as a consequence, higher shear rate exerted on the inner liquid by the middle phase.These small drops easily penetrate into the middle phase and form a double emulsion.Rising the speed of the fluid jet leads to a relative enhance in the dimensions of the outer droplets due to the increase in the length of the fluid string.Thus, from the provided explanations, it is evident that the thickness of the encapsulated structure rises with an enhance in the velocity of the

Effect of outer phase flow rate
The effect of outer phase velocity on double emulsion creation is displayed in Figs 12 & 13.The jetting regime for inner and outer droplets for CMC 0.4% Wt and CMC 0.5% Wt fluids can be observed for   > 3 and   ≥ 3.35, respectively (Figs 12a and b).By enhancing the flow rate of the continuous phase, the size of both inner and outer droplets decreases dramatically (Figs 12a and b), which results in a higher frequency of compound droplets generation (Fig 12d).The diameter of the outer droplet is more affected by the velocity of the outer phase, resulting in a decrease in shell thickness (Fig 12c).The reduction in droplet size can be explained by the fact that an enhance in the outer phase's capillary number results with increased shear rates at the outer interface, owing to the rising velocity gradient.Additionally, the decrease in the diameter of the middle phase jet generates a higher velocity gradient in the vicinity of the inner interface (the interface between the inner and middle liquids), resulting in higher shear rates and a lower diameter of the inner droplet.The decrease in the dimensions of the droplets is due to the fact that with the enhance in the speed of the continuous phase and the subsequent enhance in capillary number and increase in velocity gradient, the continuous phase applies more shear stress to the middle phase and, as a consequence, from the middle phase to the inner phase, which results in a reduce in the thickness of the forming jet, facilitating the separation of droplets, growing the frequency of the formation of droplets, and finally reducing the final dimensions of the droplets being formed.And as is well established, the dimensions of the outer and inner droplets decrease as the concentration of the non-Newtonian fluid increases.

Conclusion
The aim of this research paper was to use the VOF numerical model to study the creation of double emulsion in a double Co-Flow geometry.By changing the physical properties and phase velocities of the system, the study explored the formation of these droplets in Newtonian vs. non-Newtonian ambient fluids.The model convincingly validated by the experimental data and was utilized to study phase flow rates' effects on the features and behavior of compound droplets.The study investigated how different input parameters influence the properties of double emulsions, such as size of the outer and inner droplets, thickness, and frequency of creation.Additionally, the study demonstrated that it's possible to alter the features of the inner phase droplet by adjusting the velocity of the internal phase, irrespective of whether the ambient fluid is Newtonian or non-Newtonian.However, in both Newtonian and non-Newtonian ambient fluids, altering the size of the outer droplet alone is a more difficult process.The study employed three dimensionless groups (namely, we in , Ca mid , and Re out ) to examine how effective forces influence the size of compound droplets.The findings revealed that increasing we in resulted in larger inner droplets, while greater Re out led to smaller inner droplets.However, as we in increases, the outer droplet size remains almost constant, whereas as Re out increases, the outer droplet size decreases, while increasing Ca mid results in bigger outer droplets and smaller inside droplets in both Newtonian and non-Newtonian ambient fluid.Also, the results showed that in non-Newtonian fluids, smaller droplets are formed than in Newtonian fluid, and with the enhance in the concentration of the non-Newtonian fluid, or, in other words, with the increase of  0 and λ and the decrease of n, the diameter of the double emulsions formed

Declaration of Competing Interest
domain, and the boundary conditions, is depicted in Fig 1.The capillary tubes located inside, in the middle, and outside have varying internal radiuses of 10, 30, and 70 , respectively.To ensure precise formation of the double emulsions, the inner and middle capillary tips are positioned at the same location.Moreover, to reduce the impact of the outflow boundary, the outer capillary tube is kept 40 times as long as the outside capillary's radius.

Fig 1 .
Fig 1.This diagram depicts the double Co-Flow geometry for the fabrication of double emulsions, including the meshed area and boundary conditions.

Fig. 2
Fig. 2 illustrates the properties of non-Newtonian fluids' rheology, which were gathered

Fig 2 .
Fig 2. Schematic of the viscosity of shear-thinning CMC solutions as a function of shear rate.

Fig 3 .
Fig 3. (Left) Four different resolution meshes and (Right) the contour of volume fraction for each mesh.The equivalent process circumstances are   = 0.003,   = 0.6, and   = 30

Fig 4 .
Fig 4. The experimental measurements resulted in droplet diameters of 1433  and 2003 , while the numerical simulation predicted inner and outer droplet diameters of 1403  and 1983 , respectively.As a result, the errors for the inner and outer droplet

Fig 5 .
Fig 5. steps for making a double emulsion (a) First, the jet of the inner phase moves within the Fig 6a.The middle phase's developing droplet is subsequently penetrated by the inner filament (Fig 6b).Because of dripping instability, a double emulsion eventually forms at the downstream end of the channel (Fig 6c).

Fig 6 .
Fig 6.Steps in the process of creating a double emulsion under the dripping regime: (a) The filament of emerging middle phase jet (Fig 7a), the necking pattern causes the inner jet stream to break up into one or more droplets (Fig 7b), and the inner droplets are encapsulated by a number of splits in the outer jet stream (Fig 7c).

Fig 7 .
Fig 7. Steps in the process of creating a double emulsion under the jetting regime: (a) The inner phase's

Fig. 8
Fig. 8 demonstrates how the presence of "  " affects various parameters including inner and outer double emulsion sizes (  and   ), thickness (t = (  −   )), and the frequency

Fig 8 .
Fig 8. (a) Inner phases droplet size, (b) outer phases droplet size, (c) thickness, and (d) formation frequency as a function of We in .
diameter of the double emulsion decreases with increasing non-Newtonian fluid concentration.It should be noted that there are two ways to interpret this condition.One is the change of base viscosity of fluids by adding water soluble polymer concentration and the other is non-Newtonian properties of continuous phase fluids.In other words, the diameter of double emulsions decreases as viscosity increases in non-Newtonian fluids of the Carreau viscosity model, simultaneously with increases in  0 , λ and decreases in n.Another conclusion is that the successful creation of the double emulsion structure depends on selecting the right non-Newtonian fluid concentration.
middle phase.It is a well-established fact that the frequency of double emulsion formation exhibits an almost linear increment with the rise in the flow rate of the middle phase.This outcome is due to the rise in the separation velocity of the middle phase droplets, resulting from the reduction in the size of the inner droplets (Fig 10d).Additionally, as was mentioned in the section before, the diameter of the double emulsion decreases with increasing non-Newtonian fluid concentration.According to Fig 10 and Fig 11, the formation of double emulsions in non-Newtonian fluid CMC 0.3% Wt and Ca mid > 0.16 occurs in the jetting regime.The jetting regime also occurs in CMC 0.4% Wt and CMC 0.5% Wt fluids when Ca mid ≥ 0.14 and Ca mid ≥ 0.11, respectively.

Fig 10 .
Fig 10.(a) Inner phases droplet size, (b) outer phases droplet size, (c) thickness, and (d) formation frequency as a function of   .

Fig 12 .
Fig 12. (a) Inner phases droplet size, (b) outer phases droplet size, (c) thickness, and (d) formation frequency as a function of Re out .

Fig 14 .
Fig 14.(a) Inner phases droplet size, (b) outer phases droplet size as a function of   , (c) Inner phases droplet size, and (d) outer phases droplet size as a function of   .
decreases.In general, this research demonstrates how phase characteristics and flow rates impact the creation of double emulsion in a double Co-Flow microfluidic configuration.The findings of this study can be utilized to produce various types of liquid-liquid encapsulated structures in both Newtonian and non-Newtonian fluids.The findings indicate the development of a microfluidic configuration capable of producing double emulsions with highly controllable adjustments.This device has the potential to be applied in a variety of fields, including drug delivery, fluid transportation, and synthesis.To expand on these findings, future research could focus on a more comprehensive analysis of the effects of phase properties, non-Newtonian fluid rheology, and capillary tube diameter through a parametric study.CRediT authorship contribution statement Vahid Mollania Malakshah: Writing -original draft, Software, Validation, Data curation, Visualization, Investigation.Pedram Hanafizadeh: Supervision, Conceptualization, Methodology, Writing -review & editing.

Table 2 .
The assortment of input variables used in the design of the numerical investigation.

Table 3 .
Details about different mesh types and their inner droplet sizes.