Volume-of-fluid simulations in microfluidic T-junction devices: Influence of viscosity ratio on droplet size
11 Volume-of-fluid simulations in microfluidic T-junction devices: Influence of viscosity ratio on droplet size
Mehdi Nekouei and Siva A. Vanapalli
Department of Chemical Engineering, Texas Tech University, Lubbock, TX, 79409, USA
Abstract
We used volume-of-fluid (VOF) method to perform three-dimensional numerical simulations of droplet formation of Newtonian fluids in microfluidic T-junction devices. To evaluate the performance of the VOF method we examined the regimes of drop formation and determined droplet size as a function of system parameters. Comparison of the simulation results with four sets of experimental data from the literature showed good agreement, validating the VOF method. Motivated by the lack of adequate studies investigating the influence of viscosity ratio ( ) on the generated droplet size, we mapped the dependence of drop volume on capillary number (0.001 < Ca < 0.5) and viscosity ratio (0.01 < < 15). We find that for all viscosity ratios investigated, droplet size decreases with increase in capillary number. However, the reduction in droplet size with capillary number is stronger for < 1 than for > 1. In addition, we find that at a given capillary number, the size of droplets does not vary appreciably when < 1, while it increases when > 1. We develop an analytical model for predicting droplet size that includes a viscosity-dependent breakup time for the dispersed phase. This improved model successfully predicts the effects of viscosity ratio observed in simulations. Results from this study are useful for the design of lab-on-chip technologies and manufacture of microfluidic emulsions, where there is a need to know how system parameters influence droplet size.
1. Introduction
Droplet-based microfluidics where fluid volumes down to picoliters are manipulated has witnessed a remarkable growth due to applications in biochemical analysis and material synthesis . In these applications, it is necessary to produce droplets of controlled size. Droplet size is an extremely important parameter as it controls the efficiency of encapsulation of individual cells and biomolecules and rates of reaction
1, 9, 10 . From a fundamental point of view, droplet size dictates the mixing dynamics
1, 9-13 , flow resistance , breakup , coalescence and collective
24, 25 behavior. Since experimental conditions including flow rates, fluid properties and channel dimensions, may vary for different applications it is important to understand drop formation and develop predictive models of how system parameters influence droplet size. A widely used microfluidic geometry for producing droplets is the T-junction device, where the continuous and dispersed phase flowing orthogonally meet at a junction producing droplets
27, 28 . Several parameters have been shown to influence droplet size in T-junction devices . These include dimensionless parameters such as flow rate ratio (
𝑄 = 𝑄 𝐷 𝑄 𝐶 , where Q D and Q C are volumetric flow rates of dispersed and continuous phases, respectively), capillary number ( 𝐶𝑎 = 𝐶 𝑈 , where C and U are the viscosity and velocity of the continuous phase, and is the interfacial tension), capillary number of dispersed phase ( 𝐶𝑎 𝐷 = 𝐷 𝑈 𝐷 , where D is the viscosity of dispersed phase, U D is the inlet velocity of dispersed phase), Reynolds number ( 𝑅𝑒 = 𝐶 𝑈𝑤 𝐶 𝐶 where C is the density of continuous phase and w C is the width of main channel), viscosity ratio ( = 𝐷 𝐶 ) and density ratio ( = 𝐷 𝐶 , where D is the density of dispersed phase). In addition, geometrical parameters such as width ratio ( 𝑊 = 𝑤 𝐷 𝑤 𝐶 , where w D is the width of side channel) and height ratio ( 𝐻 = ℎ𝑤 𝐶 , where h is the height of the channel) can also influence droplet size. Several experimental studies have investigated the dependence of drop size on system parameters in T-junction devices . Among these, here we discuss those that are relevant to our work (see Table 1). Garstecki et al . focused on the effect of flow rate ratio on droplet size and found a linear relationship between drop length (L) and flow rate ratio . Van Steijn et al . investigated how the dimensions of the main and side channel in the T-junction influence droplet size . Similar to Garstecki et al. , for a given T-junction geometry, they also found linear dependence of drop volume (V) on flow rate ratio. However, they showed that at a given flow rate ratio, drop size also depends on width ratio (W) and height ratio (H). In both these studies, experiments were conducted at low capillary numbers (Ca < 0.01). Christopher et al . measured droplet size across a wider range of capillary numbers (0.001 < Ca < 0.5) while maintaining fixed flow rate ratios . In addition to linear dependency of droplet size on flow rate ratio, they found that the droplet size decreases by increasing capillary number. The parameter range covered by these experimental studies is shown in Table 1. Investigating experimentally, the effect of system parameters on droplet size has limitations. For example, experimentally it is difficult to produce droplets for viscosity ratios greater than unity, making it unclear how viscosity ratio influences droplet size. This feature is evident from Table 1, where studies were limited to viscosity ratio less than 1. Thus, fixing one control parameter while changing the other parameters is not always guaranteed in experiments, lending access to a narrow system parameter space. Moreover, it is difficult to experimentally measure the three- Table 1. Experimental studies of droplet formation in microfluidic T-junction devices, whose results on droplet size have been compared to VOF simulations. Study Ca Q W Focus of study Garstecki et al. , 2006 [30] Ca<0.01 0.01 1, the viscous stress of dispersed phase can impact drop production, studies covering a broader range of viscosity ratios need to be pursued. In this study, we use numerical simulations based on VOF method to study the influence of system parameters on drop size (see Table 2). In the first part of our study, we carried out numerical simulations to predict regimes of drop formation and the generated droplet size. We compared the simulation results with data from experimental studies shown in Table 1, allowing us to assess the capability of VOF method and validate it. In the second part of the study, we investigated the influence of viscosity ratio on droplet size and find it to play an important role. Subsequently, we develop a model that predicts droplet size considering the effect of viscosity ratio.
2. Numerical Simulation 2.1. Volume-of-fluid method and its implementation
Three-dimensional simulation of droplet formation in T-junction geometries was performed using volume-of-fluid (VOF) method. VOF is a Eulerian method of multiphase flow simulations where fluid properties such as viscosity and density are smoothed and the surface tension force is distributed over a thin layer near the interface as a body force. In VOF, a phase fraction parameter, , is used to indicate the presence of each phase at every location of the simulation domain. In our simulation, =1 for phase 1 (i.e. continuous phase), = 0 for phase 2 (i.e. dispersed phase) and 0< <1 in the interface region. In VOF, the governing equations including continuity (Eqn.1), momentum balance (Eqn. 2) and phase fraction equations (Eqn. 3), are solved simultaneously. 𝛁. 𝑼 = 0 (1) 𝜕𝜌 𝑏 𝑼𝜕𝑡̅ + 𝛁. (𝜌 𝑏 𝑼𝑼) = −𝛁𝑝 + 𝛁. 𝑻 + 𝜌 𝑏 𝒇 + 𝑭 𝑠 (2) 𝜕𝛼𝜕𝑡̅ + 𝛁. (𝛼𝑼) = 0 (3) In Eqns. 1-3, U is the velocity vector field, p is the pressure field, T is the deviatoric stress tensor ( 𝑻 = 2𝜇𝑺 − 2𝜇(𝛁. 𝑼)𝑰/3 , where
𝑺 = 0.5(𝛁𝑼 + 𝛁𝑼 𝑇 ) and I is identity matrix), f is gravitational force. Parameters b and b are bulk viscosity and density are based on the weighted average of the distribution of phase fraction: 𝜇 𝑏 = 𝛼𝜇 𝐶 + (1 − 𝛼)𝜇 𝐷 (4) 𝜌 𝑏 = 𝛼𝜌 𝐶 + (1 − 𝛼)𝜌 𝐷 (5) The last term on the right-hand side of Eqn. 2, F S , represents the continuum surface tension force (CSF) and is nonzero only on the interface. This force term is defined as 𝑭 𝑠 = (∇𝛼) where is curvature ( = ∇. ( ∇𝛼 ∇𝛼 ) ). We solved the momentum balance equation in conjunction with the continuity equation using the Pressure Implicit with Splitting of Operators (PISO) method . In this method, the velocity field is predicted and then corrected to advance the pressure and velocity fields in time. In this work we used three PISO iterations. Once the velocity field was found, Eqn. 3 was solved to find the phase fraction. Even though Eqns. (1-3) yield velocity and phase fraction at every cell in the domain, the location of the interface needs to be identified with high resolution. To achieve this, we used a two fluid formulation where the contribution of each phase to the velocity of interface is considered, i.e. 𝜕𝛼𝜕𝑡̅ + 𝛁. (𝛼𝑼 𝐶 ) = 0 (6) 𝜕(1−𝛼)𝜕𝑡̅ + 𝛁. ((1 − 𝛼)𝑼 𝐷 ) = 0 (7) where U C and U D are velocity vector fields of continuous and dispersed phase respectively. Here, we assumed that velocity of each phase has a contribution to the convection of interface based on their phase fraction, i.e. 𝑼 = 𝛼𝑼 𝐶 + (1 − 𝛼)𝑼 𝐷 (8) Equation (6) can be rearranged and used as phase fraction equation : 𝜕𝛼𝜕𝑡̅ + 𝛁. (𝛼𝑼) − 𝛁. (𝛼(1 − 𝛼)𝑼 𝑟 ) = 0 (9) where 𝑼 𝑟 = 𝑼 𝐷 − 𝑼 𝐶 , is called as the compression velocity. Eqn. 9 has a new convective term, compared to Eqn. 3. This term is only present at the interface and vanishes in the pure phases. The compression velocity is given by: 𝑼 𝑟 = 𝒏 𝒇 min [𝐶 | 𝜑𝑺 𝒇 | , 𝑚𝑎𝑥 ( |𝜑||𝑺 𝒇 | )] (10) 𝒏 𝒇 is the cell normal flux, the ratio ( |𝜑||𝑺 𝒇 | )] is the magnitude of velocity where and S f are the cell face volume flux and surface area, respectively. C is a user-specified compression factor that can vary from zero to four. Using larger compression factor results in thinner interface. However, Hoang et al. who studied microfluidic droplet break-up showed that compression factor greater than one causes parasitic currents at interface , therefore we chose C = 1 in this study. In order to assure stability and convergence of the simulation we used an adaptive time step method. At the beginning of each iteration, a new time step was calculated based on the Courant number (Co), which is defined as: 𝐶𝑜 = 𝑼 𝑓 .𝑺 𝑓 𝑑 𝑓 .𝑺 𝑓 ∆𝑡̅ (11) where d f is the distance between two neighbor cells and U f is the velocity at the surface of the cell and ∆𝑡̅ is time step. Then based on the identified courant number, a new time step size was calculated in order to keep Co less than a predefined limit. The value of Courant number indicates how much each fluid element is displaced over one time step. For example, when Co=1 it means that, each element of the fluid moves in a distance of one grid size in one time step. Courant number was evaluated at each computational cell including bulk fluid and the interface in each iteration. All the numerical simulations were performed in the open-source code OpenFOAM. VOF was implemented in the interFoam solver as reported by Deshpande et al . . Several works have used interFoam for solving incompressible two phase flows
19, 44, 61-63 . Among those papers Hoang et al . and Nieves-Remacha et al . have used this solver for simulation of two phase flows at micro-scale
19, 44, 61 . In OpenFOAM, governing equations were discretized by finite volume center-based method. Boundedness of phase fraction parameter was controlled by a Total Variation Diminishing (TVD) method implemented in OpenFOAM called MULES (Multidimensional Universal Limiter with Explicit Solution) . Figure 1. Drop formation at a microfluidic T-junction. (a) Schematic of the microfluidic T-junction highlighting the relevant system parameters. The definition of these parameters is provided in the main text. Representative images showing the three modes of droplet breakup (b) Squeezing, Ca=0.01 (c) Dripping, Ca=0.03 (d) Jetting, Ca=0.07. The different breakup behaviors were obtained by varying the capillary number, while keeping other system parameters fixed (Q=0.2, =0.1, W=1 and H =1). a) b) c) d) The schematic of T-junction’s geometry used in the simulation study is shown in Figure 1. The system has two inlets for continuous and dispersed phases and one outlet. We set = 0 at the inlet of the dispersed phase and = 1 at the inlet of the continuous phase. A constant velocity was set for both the inlets in the system. No slip boundary condition was applied at the walls. All the simulations were carried out for a density ratio, = 1 and Reynolds number is Re < 0.1, therefore inertia is negligible. In this study, time is scaled by 𝑤 𝐷 𝑈 𝐷 . In order to predict the values of phase fraction near the walls, i.e. the interface, the normal of interface is related to the wall’s normal and tangential unit vectors: 𝑛 𝑖𝑛𝑡𝑒𝑟𝑓𝑎𝑐𝑒 = cos(𝜃) 𝑛 𝑤𝑎𝑙𝑙 + sin (𝜃)𝑡 𝑤𝑎𝑙𝑙 (12) where n interface is the normal of interface, is the wall’s static contact angle and n wall and t wall are the unit normal and unit tangential vectors to the wall, respectively. The walls of the channels were considered as fully wetted by the continuous phase; therefore we set the contact angle as zero at the walls in this study. The boundary condition at the outlet was set as atmospheric pressure. The simulation domain was meshed using hexahedral cells. Adjacent cells to the wall were refined five times to capture lubrication film around the droplet. The domain meshing was carried out using ANSYS Gambit version 2.4.6. Since determining the size of the droplets is the main goal of this study, grid size of the mesh becomes important. We carried out numerical simulations on a representative T-junction geometry to evaluate the sensitivity of grid size (d) on Figure 2. Effect of mesh size on the generated droplet volume in a T-junction geometry, evaluated at two different capillary numbers. The conditions of the simulation are Q=0.5, =1, W=0.5 and H=0.33. the droplet volume. The width and height ratio were set to be W=0.5 and H=0.33 respectively. We performed simulations for a fixed flow rate ratio (Q=0.5) and two capillary numbers (Ca=0.01 and 0.07). We measured the droplet volume (V) by counting the number of grid points in which ≤ 0.5. The results are shown in Fig. 2. We find that for w C /d > the variation in droplet size is smaller than 4%, while for w C /d > it’s smaller than 1%. Therefore, a grid size of w C /d=
100 was chosen in this study. To clarify whether the upstream and downstream channel lengths of T-junction are long enough for the flow to become fully developed, we computed the entrance length in our system. The entrance length for channel flows under laminar conditions (Re<2000) is given by 𝐿 𝑒 = 𝑑 ℎ ( 0.61 + 0.035𝑅𝑒 + 0.056𝑅𝑒) where d h is the hydraulic diameter of the channel. We find that the maximum entrance length, corresponding to highest Re, needed to achieve a fully developed flow is just a small fraction of the upstream (L e =0.11 L Upstream ) and downstream lengths (L e =0.02 L Downstream ) used in the simulation. Therefore, the choice of entrance lengths used in the simulation is sufficient to obtain fully developed flow. All the simulations were performed in parallel on a Linux cluster (Hrothgar cluster at High Performance Computing Center, Texas Tech University and Stampede cluster, Texas Advanced Computing Center) by employing 48 to 96 processors. Based on initial test simulations, we determined that in order to prevent spurious currents and nonphysical behavior of interface, courant number and interface courant number need to be kept below 0.1. Typical clock time for the formation of one droplet ranges from 22 – 40 hrs depending on the number of processors used.
3. Results and discussion 3.1. Regimes of Drop Formation
Previous studies have shown that three distinct regimes – squeezing, dripping and jetting - exist for the dispersed phase behavior in a T-junction geometry
11, 38, 66, 67 . In this section, we perform simulations to examine whether VOF can capture the different regimes of drop production. We also compare our simulation results to the experimental data of Tice et al . who identified these regimes as a function of flow rate ratio and capillary number. Similar to experiments, we observe the three behaviors in our VOF simulations as shown in Fig. 1. We find that in the squeezing regime, the dispersed phase blocks the main channel significantly and breakup occurs in the vicinity of the T-junction. In the dripping regime, the dispersed phase only blocks the main channel partially, and penetrates further from the T-junction and the droplets are produced at a fixed spatial location in the microchannel. In contrast to the dripping regime, in the jetting regime, we find that after the generation of the first droplet the thread of the dispersed phase continues to move forward and a second droplet is produced further downstream. At high capillary numbers, Ca > 0.1, we observed that dispersed phase co-flows with the continuous phase forming parallel stream without any droplet being formed. These different behaviors, we observed with VOF simulations are consistent with those observed in experiments
11, 27, 67 , phase field and lattice-Boltzmann simulations . To test whether our VOF simulation can quantitatively predict the transitions between these regimes as a function of system parameters, we compared the results from our simulation with the experimental data of Tice et al. . In their study, experiments were conducted for a geometry with W=H=1, for two different viscosity ratios =0.1 and 10, and for 0.1 < Q < 10 and 10 -3 < Ca < 10 -1 . Simulations were performed at the same parameter values as the experiments. Figure 3 shows the map of the regimes as a function of the two control parameters - Ca and Q. A good agreement is observed between the simulation results and the experimental data. As shown in Figure 3, for relatively small capillary number and flow rate ratio, droplets are produced in the squeezing regime. However, by increasing Ca and/or Q, the regime of drop formation changes to dripping and jetting. The dripping region is much narrower compared to squeezing and jetting. By increasing the flow rate ratio, we find that the transition from squeezing to dripping and dripping to jetting occurs at a lower capillary number. This observation indicates that at higher flow rate ratios dispersed phase penetrates much more into the main channel such that the continuous phase can break it up only in the downstream region. By comparing results for two different viscosity ratios in Fig. 3, we observe that the boundaries defining the transition between regimes are steeper for =10 than =0.1. This means that the transition between regimes occurs at significantly smaller capillary numbers for higher viscosity ratio. When the viscous force of dispersed phase increases, it is difficult for the continuous phase Figure 3. Regimes of drop formation in a microfluidic T-junction device with W=1, H=1. VOF simulations (closed symbols) were conducted for (a) =0.1 and (b) =10. Experimental data (open symbols) are from Tice et al . . The dashed lines are drawn to guide the eye. to fragment it. The viscous force of dispersed phase depends on both the viscosity and velocity of dispersed phase. Therefore, for high viscosity ratios or high velocities of dispersed phase, drops are generated in jetting regime. In the previous section, we have quantified the regimes of drop formation in the T-junction device with VOF simulations. We note it is difficult to define sharply the transition between the regimes. Therefore, in order to benchmark the capabilities of VOF simulations more precisely we compared VOF predictions of drop size against experimental data. As mentioned earlier, we choose the three experimental data sets reported in Table 2 to validate the VOF method. Figure 4 compares the data reported by Garstecki et al . and the numerical simulations in a fixed T-junction geometry at two different viscosity ratios, =0.01 and 0.1. The experiments were conducted at three different continuous phase flow rates, Q c = 0.0028, 0.028 and 0.28 L/s. Garstecki et al . plotted the droplet length normalized with the main channel width, L/w c , as a function of flow rate ratio, Q. They found that the droplet length is mostly constant at low flow rate ratios but increases almost linearly at high flow rate ratios. Our simulation results also show similar trends and are in good agreement with their experimental data. The experimental data obtained by Garstecki et al. in Fig. 4 pertains to a single T-junction geometry. To explore the capability of VOF simulations to predict droplet size generated in other T-junction geometries, we relied on the experimental study by Van Steijn et al . . They used three different geometries with different width and height ratios, (W, H) = (0.33, 0.33), (0.67, 0.17) and Figure 4. Effect of flow rate ratio on the size of droplets for (a) =0.01 and (b) 0.1. Experimental data (open symbols) are from reference [30] and the VOF data are represented by closed symbols. The geometrical parameters are W=0.5, H=0.33 and the flow rates of continuous phase are 0.0028, 0.028 and 0.28 µL/S. The prediction based on Eqn. 14 (see main text) is shown by ( ). (1.33, 0.11). We used the same geometries in our study and the same viscosity ratio of =0.1. We varied the flow rate ratio (0.2 < Q < 6) by fixing the continuous phase flow rate and Figure 5. Effect of T-junction geometry on the generated droplet size for =0.1. Experimental data (open symbols) are from reference [31]. VOF data (closed symbols) were obtained for Ca ~ O(10 -3 ). The dashed lines are drawn using equations in Figure 2 of reference [31]. V / h W C Q ⚫ / ◯ W=0.33, H=0.33 (Sim/Exp) ▲ / △ W=0.67, H=0.11(Sim/Exp) ◆ / ◇ W=1.33, H=0.17 (Sim/Exp)
Figure 6. Effect of capillary number on the size of droplets for W=1, H=0.5 and =0.01. Experimental data (open symbols) are from reference [29] and VOF data is shown by closed symbols. changing the flow rate of dispersed phase. Van Steijn et al. plotted normalized droplet volume, V/hw C2 , as a function of flow rate ratio for the three mentioned geometries (see Fig. 5). They found that the width of the side channel influences the volume of droplet significantly. At a fixed flow rate ratio, the volume of produced droplet in a T-junction with a wider side channel is larger. For example, at Q = 5 the droplet volume in T-junction with W=1.33 is four times larger than in W=0.33. As shown in Fig. 5, our simulation results predict the same trend of data as the experimental findings. The experimental data of Garstecki et al. and van Stein et al. pertain mostly to low capillary numbers (Ca ~ O(10 -3 )). We therefore conducted VOF simulations at high capillary numbers and compared the results with that of the experimental data reported by Christopher et al. . They employed one set of fluids, =0.01 and one geometry, W=1 and H=0.33 in their experimental study. They fixed the flow rate ratios at, Q=0.05, 0.25 and 0.5, and changed the capillary number by varying inlet velocity of continuous phase. We performed numerical simulations at the same parameter values as experiments. As shown in Figure 6, our numerical findings are consistent with the experimental data. By increasing the capillary number for any fixed value of flow rate ratio the droplet volume decreases. On the other hand, by increasing flow rate ratio at fixed capillary number the volume of droplet increases. Overall, comparison of simulation results with the three sets of experimental data reveals that our VOF methodology and choice of simulation parameters (e.g. grid size, compression factor, Courant number) are suitable for making measurements of droplet size in microfluidic T-junction devices. Figure 7. Images from the VOF simulation (for Q=0.8, Ca=0.05, =0.1, W=1 and H=0.33) showing the droplet formation process, highlighting the filling and squeezing stages. During the filling stage, the dispersed phase occupies the main channel. Because of the blockage of the main channel by the dispersed phase, the upstream pressure builds, squeezing the neck, which connects the dispersed phase to droplet. In this section, we discuss a popular model for predicting droplet size at low capillary number in T-junction devices, as originally proposed by Garstecki et al . . Given that this simple model ignores the effect of fluid viscosity ratio, we subsequently study the influence of viscosity ratio using VOF simulations and highlight the need to incorporate this parameter into analytical predictions of droplet size. At low capillary numbers, Garstecki et al. assume that a drop is produced in a two-stage process. In the first stage, the dispersed phase starts to enter the main channel and occupies a portion of the channel. This stage is called the filling stage. In the second stage, called the squeezing stage, the neck which connects the drop to the dispersed phase is squeezed by the continuous phase, ultimately pinching-off the dispersed phase. Figure 7 shows snapshots from the VOF simulation depicting these two stages of droplet generation. Figure 8. Dependence of droplet volume on the capillary number for different viscosity ratios. Other parameters were fixed at Q = 0.3, W = 1 and H = 0.33. Using this conceptual picture, Garstecki et al . proposed scaling arguments to obtain a prediction for droplet size . They argued that during the filling stage, the extent to which the dispersed phase fills the main channel is approximately equal to the width of the main channel, w c . In the squeezing step, the droplet is being filled by an amount equal to U D t sq where U D is the dispersed phase velocity and t sq is the time needed to squeeze and rupture the dispersed phase finger. By adding these two contributions to the overall droplet length, we have 𝐿 = 𝑤 𝐶 + 𝑈 𝐷 𝑡 𝑠𝑞 . Assuming that the dispersed phase is being squeezed by a rate proportional to the continuous phase inlet Figure 9. The effect of viscosity ratio on the thickness of dispersed phase and the timescales during droplet formation. The simulations were performed for Ca = 0.01, Q=0.3, W=1 and H=0.33 (a) Snapshots of mid-plane view of the dispersed phase as a function of time. Yellow arrows represent thickness of the neck, d neck , over time. Regions corresponding to filling and squeezing are highlighted (b) Time-evolution of dispersed phase thickness for two viscosity ratios of =0.1 and 5. The dashed vertical line demarcates the duration of filling and squeezing stages. (c) Plot of squeezing and filling times as a function of viscosity ratios for Ca = 0.003 (circles) and Ca=0.01 (triangles). Filling Squeezing a) c)b) d n ec k / w C tt fill t sq □ =0.1 =5 velocity (U C ) gives an estimate of 𝑡 𝑠𝑞 = 𝑑 𝐶 𝑈 𝐶 , where d C is an undetermined length scale. Using this estimate of t sq, the equation for droplet length can be rewritten as 𝐿𝑤 𝐶 = 1 + 𝛽𝑄 (14) where = d C /w C ~ O(1) is a fitting parameter that depends on the T-junction geometry. For example, for the data in Fig. 4, = 1.13 ± 0.16 and 1.67 ± 0.41 for = 0.01 and 0.1 respectively. The effect of the T-junction geometry was explicitly considered by Van Steijn et al. , as reflected by the data of Fig. 5. The low capillary number model of Garstecki et al. relates the droplet size to simply the flow rate ratio and does not consider fluid viscosity ratio. Likewise, the compilation of studies in Tables 1 and 2 did not cover a broad range of viscosity ratios and were conducted with fluids for which ≤ 1. We therefore investigated the importance of viscosity ratio by varying from 0.01 to 15. In the simulations, was varied by changing the viscosity of dispersed phase and keeping the continuous phase viscosity constant. Additionally, simulations were conducted at a fixed flow rate ratio (Q=0.3) in a device geometry with W=1 and H=0.33. We chose this flow rate ratio as it allows a wider access into the squeezing and dripping regimes. However, at higher viscosity ratios ( =10, 15), we obtained limited data, because the operating window for squeezing/dripping regime is small as discussed earlier in section 3.1. The capillary number corresponding to these simulations was 0.001 < Ca < 0.02. Fig. 8 shows the droplet size as a function of Ca for different viscosity ratios. For all viscosity ratios, by increasing the Ca, the size of droplet decreases suggesting that increasing the viscous stress of continuous phase produces smaller droplets. However, the reduction in droplet size with Ca is stronger for < 1 than for > 1. In addition, we find that at a given Ca, the size of droplets does not vary appreciably when < 1, while it increases when > 1. To understand the observations on the effect of viscosity ratio, we followed the time evolution of the width of the dispersed phase neck (d neck ) during the breakup process, as shown in Fig. 9a. Here d neck is measured at the channel location where ultimately the pinch-off occurs. Fig.9b shows how the non-dimensional neck thickness, d neck /w D, varies with time for = 0.1 and 5. As expected, we observe that the neck thickness increases and then decreases with a maximum. We define the duration in which the neck thickness continues to increase and then decrease as t fill and t sq respectively, since these two time scales effectively correspond to the filling and squeezing stages in the droplet generation process. From the data in Fig. 9b, we observe that at a fixed capillary number, the duration of filling stage for the two viscosity ratios is the same while the duration of the squeezing stage is different. In fact, as shown in Fig. 9c, when explored across the entire viscosity range (0.01 ≤ ≤ 15) for varying capillary number, we find t fill is independent of viscosity ratio and is larger when the capillary number is low. In contrast, t sq increases with viscosity ratio, but decreases with increase in capillary number. The observed dependence of t fill and t sq on Ca and help to explain our simulations results shown in Fig. 8. Since both t fill and t sq decrease with Ca, the drop volume decreases with increase in Ca. At a fixed Ca, since t sq increases with viscosity ratio, while t fill is constant we observe larger droplet size for higher viscosity ratio. Our results thus far highlight that the viscosity ratio affects the dispersed phase breakup behavior and that the squeezing duration increases with viscosity ratio, leading to larger droplet sizes. The drop breakup model of Garstecki et al. discussed earlier (see Sec. 3.3) estimates the squeezing time as 𝑡 𝑠𝑞 = 𝑑 𝑐 𝑈 𝐶 . Clearly, the effect of viscosity ratio is missing in this estimate of squeezing time. Here we address this gap and develop an improved model for predicting the generated droplet size. In our model, similar to previous work
30, 31 , we assume the overall drop volume is resulting from volumetric contributions during the filling and squeezing stages, i.e.
𝑽 = 𝑽 𝒇𝒊𝒍𝒍 + 𝑽 𝒔𝒒 (15) To estimate the volume of the droplet at the end of the filling period, we use Eqn. 16, which was derived by Christopher et al . . Eqn.16 was obtained by conducting a force balance that includes continuous phase shear stress, upstream fluid pressure and Laplace pressure jumps across the front and rear of the droplet. (1 − √ 𝑉 𝑓𝑖𝑙𝑙 ℎ𝑤 𝐶2 ) = 𝐶𝑎√ 𝑉 𝑓𝑖𝑙𝑙 ℎ𝑤 𝐶2 (16) Eqn.16 is consistent with our simulation result that t fill is independent of viscosity ratio but dependent on capillary number (c.f. Fig. 9). In addition, as
𝐶𝑎 → 0 , Eqn.16 predicts that the non-dimensional fill volume 𝑉 𝑓𝑖𝑙𝑙 ℎ𝑤 𝐶2 → 1 , which is consistent with the low capillary number model of Garstecki et al . To estimate V sq , we determine the overall squeezing time (T sq ) by summing the break-up time of a viscous thread (t b ) and the squeezing time considered by Garstecki et al. , i.e. 𝑇 𝑠𝑞 = 𝑡 𝑠𝑞 + 𝑡 𝑏 (17) To determine t b , we consider the break-up dynamics of Newtonian filaments that has been investigated in several studies . These studies report the viscous filament break-up time to be 𝑡 𝑏 = 𝑘 µ 𝐷 𝑑 𝐶 , with 6 < k < 33. This time scale results from the competition between surface tension trying to shrink the liquid filament and the viscous stresses in the filament opposing it. Incorporating the expressions for t sq and t b into Eqn. 17 we obtain 𝑇 𝑠𝑞 = 𝑑 𝐶 𝑈 𝐶 [1 + 𝑘𝜆𝐶𝑎] (18) Given the volumetric contribution to the drop size during the squeezing period is 𝑉 𝑠𝑞 = 𝑇 𝑠𝑞 𝑄 𝐷 , we have 𝑉 𝑠𝑞 ℎ𝑤 𝐶2 = 𝛽𝑄 + 𝐶𝑎 𝐷 (19) where = 𝑘 𝑤 𝐷 𝑑 𝐶 𝑤 𝐶2 . Solving Eqns. (15), (16) and (19) together gives the final prediction for the generated droplet volume in T-junction devices taking into account the viscosity ratio of the fluids. In Figure 10, we compare the predictions of our improved model with the simulation data across a wide range of viscosity ratio and for two geometries with different aspect ratios H=0.33 and 1. For the first geometry, H=0.33, the viscosity ratio varies between 0.01 and 15 and for the second geometry, H=1, we have =0.1 and 10. Two sets of fit parameters and were used for these two geometries. We find excellent agreement between our model predictions and the simulation data capturing different viscosity ratios. For each set of data with the same Ca, we found the best-fit values of . For the first geometry, H=0.33, this value varies between 7 to 36 for the entire set of data with an average of 21.96±8.6. For the second geometry, H=1, the best-fit value of varies between 8 and 19 with an average of 14.8±3.69. These finding are consistent with the Figure 10. Parity plot showing the prediction of the improved model incorporating viscosity ratio, versus the simulation data for a) W=1, H=0.33 and 0.01≤ ≤15 and b) W=1, H=1 and =0.1 and 10. In this plot, V is the volume of the droplet normalized by ℎ𝑤 𝐶2 . Models predictions reported by Garstecki et al. and Van Steijn et al. are shown by horizontal solid and dotted lines, respectively [30,31]. The grey lines indicate 15% deviation from the expected true value. The insets show the model predictions (dotted lines) as a function of Ca D for different Ca values: ▲ , ◆ ○ □ ■ a) b) range of k values reported in the literature (5 < k < 33) since we have 𝑤 𝐷 𝑑 𝐶 𝑤 𝑐2 ~ O(1). According to McKinley and Tripathi, the k-value depends on the longitudinal stress in the viscous thread, which is a function of the axial curvature of the filament . In our study, the rate at which the dispersed phase is injected into the droplet changes when the capillary number is varied. This can cause the axial curvature of the dispersed phase finger to vary with capillary number and therefore it can result in variability of k-values, as observed when we generated the parity plot (Fig. 10). We also found the best fit parameter β = 1.18±0.4 for the first geometry, and β =1.4±0.77 for the second geometry. These values are reasonable estimates since we get =1.8 and 1.15 for the first and second geometry respectively, when evaluated with the theoretical model proposed by Van Steijn et al. that includes the effect of T-junction geometry . An important outcome of our analysis is that it explains why at a given (low) capillary number, the generated drop volumes are the same for < 1 but varies for > 1 (see Fig. 8). As shown in Fig. 11, when we plot the contributions of t sq and t b to the overall squeezing time as a function of viscosity ratio, we find that for < 1, the contribution of viscous thread breakup time is negligible compared to the squeezing time. This implies that the produced droplets will have the same size for < 1. Alternatively, we observe that for > 1, the breakup time contributes more making the drop volumes dependent on viscosity ratio. Figure 11. Dependence of normalized contributions of breakup and squeezing time-scales with viscosity ratio. The dashed vertical line indicates that the contribution of the breakup time to the overall squeezing time is small for ≤ 1. t b / T s q t s q / T s q Conclusions
We implemented the volume-of-fluid (VOF) method for investigating droplet generation in microfluidic T-junction devices. We examined the different regimes of dispersed phase behavior and found the VOF simulations to be in good agreement with experimental maps of the regimes. We also performed simulations at low and high capillary numbers, as well as in different T-junction geometries. The predicted droplet size was in good agreement with experimental reports in the literature. Taken together, these results validate our VOF methodology and choice of simulation parameters. Our simulation results show that for > 1, the droplet size depends on viscosity ratio. Previous theoretical models have not considered this dependence of viscosity ratio on the generated droplet size. We developed an improved model and incorporated this dependence by including the break-up time needed to fragment the high-viscosity dispersed phase. This improved model predicts successfully the influence of drop volume on capillary number as well as viscosity ratio. More broadly, we find that VOF is a useful tool for droplet-based microfluidics and has the potential to reveal new insights into the fluid physics of drop formation and behavior, which can have important technological applications. Acknowledgments
We thank the National Science Foundation (CAREER Grant No. 1150836 and AIR-TT Grant No.1445070) for supporting this work. We are grateful to Jerzy Blawzdziewicz and Michael Loewenberg for useful discussions.
References
1. H. Song, D. L. Chen and R. F. Ismagilov,"Reactions in droplets in microfluidic channels", Angew Chem Int Ed Engl (44), 7336-7356 (2006). 2. S. Y. Teh, R. Lin, L. H. Hung and A. P. Lee,"Droplet microfluidics", Lab Chip (2), 198-220 (2008). 3. X. C. I. Solvas and A. deMello,"Droplet microfluidics: recent developments and future applications", Chem Commun (7), 1936-1942 (2011). 4. H. N. Joensson and H. A. Svahn,"Droplet Microfluidics-A Tool for Single-Cell Analysis", Angew Chem Int Edit (49), 12176-12192 (2012). 5. R. Seemann, M. Brinkmann, T. Pfohl and S. Herminghaus,"Droplet based microfluidics", Rep Prog Phys (1) (2012). 6. A. Dewan, J. Kim, R. H. McLean, S. A. Vanapalli and M. N. Karim,"Growth kinetics of microalgae in microfluidic static droplet arrays", Biotechnology and bioengineering (12), 2987-2996 (2012). 7. J. Avesar, T. B. Arye and S. Levenberg,"Frontier microfluidic techniques for short and long-term single cell analysis", Lab on a Chip (13), 2161-2167 (2014). 8. M. A. Khorshidi, P. K. P. Rajeswari, C. Wählby, H. N. Joensson and H. A. Svahn,"Automated analysis of dynamic behavior of single cells in picoliter droplets", Lab on a Chip (5), 931-937 (2014). 9. H. Song, J. D. Tice and R. F. Ismagilov,"A microfluidic system for controlling reaction networks in time", Angewandte Chemie (7), 792-796 (2003).
10. A. M. Huebner, C. Abell, W. T. Huck, C. N. Baroud and F. Hollfelder,"Monitoring a reaction at submillisecond resolution in picoliter volumes", Analytical chemistry (4), 1462-1468 (2011). 11. J. D. Tice, A. D. Lyon and R. F. Ismagilov,"Effects of viscosity on droplet formation and mixing in microfluidic channels", Analytica chimica acta (1), 73-77 (2004). 12. P. Paik, V. K. Pamula and R. B. Fair,"Rapid droplet mixers for digital microfluidic systems", Lab on a Chip (4), 253-259 (2003). 13. W. S. Wang and S. A. Vanapalli,"Mechanisms of mass transport during coalescence-induced microfluidic drop dilution", Physical Review Fluids (6), 064001 (2016). 14. M. J. Fuerstman, A. Lai, M. E. Thurlow, S. S. Shevkoplyas, H. A. Stone and G. M. Whitesides,"The pressure drop along rectangular microchannels containing bubbles", Lab on a Chip (11), 1479-1489 (2007). 15. S. S. Bithi and S. A. Vanapalli,"Behavior of a train of droplets in a fluidic network with hydrodynamic traps", Biomicrofluidics (4), 044110 (2010). 16. C. N. Baroud, F. Gallaire and R. Dangla,"Dynamics of microfluidic droplets", Lab on a Chip (16), 2032-2045 (2010). 17. S. A. Vanapalli, A. G. Banpurkar, D. van den Ende, M. H. Duits and F. Mugele,"Hydrodynamic resistance of single confined moving drops in rectangular microchannels", Lab on a Chip (7), 982-990 (2009). 18. A. Leshansky, S. Afkhami, M.-C. Jullien and P. Tabeling,"Obstructed breakup of slender drops in a microfluidic T junction", Physical review letters (26), 264502 (2012). 19. D. A. Hoang, L. M. Portela, C. R. Kleijn, M. T. Kreutzer and V. van Steijn,"Dynamics of droplet breakup in a T-junction", Journal of Fluid Mechanics (2013). 20. D. Link, S. L. Anna, D. Weitz and H. Stone,"Geometrically mediated breakup of drops in microfluidic devices", Physical review letters (5), 054503 (2004). 21. M. Sun and S. A. Vanapalli,"Generation of chemical concentration gradients in mobile droplet arrays via fragmentation of long immiscible diluting plugs", Analytical chemistry (4), 2044-2048 (2013). 22. S. S. Bithi, W. S. Wang, M. Sun, J. Blawzdziewicz and S. A. Vanapalli,"Coalescing drops in microfluidic parking networks: A multifunctional platform for drop-based microfluidics", Biomicrofluidics (3), 034118 (2014). 23. B. Bhattacharjee and S. A. Vanapalli,"Electrocoalescence based serial dilution of microfluidic droplets", Biomicrofluidics (4), 044111 (2014). 24. P. Janssen, M. Baron, P. Anderson, J. Blawzdziewicz, M. Loewenberg and E. Wajnryb,"Collective dynamics of confined rigid spheres and deformable drops", Soft Matter (28), 7495-7506 (2012). 25. S. S. Bithi and S. A. Vanapalli,"Collective dynamics of non-coalescing and coalescing droplets in microfluidic parking networks", Soft matter (25), 5122-5132 (2015). 26. T. Ward, M. Faivre, M. Abkarian and H. A. Stone,"Microfluidic flow focusing: Drop size and scaling in pressure versus flow‐rate‐driven pumping", Electrophoresis (19), 3716-3724 (2005). 27. J. Nunes, S. Tsai, J. Wan and H. Stone,"Dripping and jetting in microfluidic multiphase flows applied to particle and fibre synthesis", Journal of physics D: Applied physics (11), 114002 (2013). 28. G. Christopher and S. Anna,"Microfluidic methods for generating continuous droplet streams", Journal of Physics D: Applied Physics (19), R319 (2007). 29. G. F. Christopher, N. N. Noharuddin, J. A. Taylor and S. L. Anna,"Experimental observations of the squeezing-to-dripping transition in T-shaped microfluidic junctions", Phys Rev E Stat Nonlin Soft Matter Phys (3 Pt 2), 036317 (2008). 30. P. Garstecki, M. J. Fuerstman, H. A. Stone and G. M. Whitesides,"Formation of droplets and bubbles in a microfluidic T-junction-scaling and mechanism of break-up", Lab Chip (3), 437-446 (2006). 31. V. van Steijn, C. R. Kleijn and M. T. Kreutzer,"Predictive model for the size of bubbles and droplets created in microfluidic T-junctions", Lab Chip (19), 2513-2518 (2010).
32. V. van Steijn, M. T. Kreutzer and C. R. Kleijn,"μ-PIV study of the formation of segmented flow in microfluidic T-junctions", Chemical Engineering Science (24), 7505-7514 (2007). 33. J. Xu, S. Li, J. Tan and G. Luo,"Correlations of droplet formation in T-junction microfluidic devices: from squeezing to dripping", Microfluidics and Nanofluidics (6), 711-717 (2008). 34. T. Glawdel, C. Elbuken and C. L. Ren,"Droplet formation in microfluidic T-junction generators operating in the transitional regime. I. Experimental observations", Physical Review E (1) (2012). 35. T. Glawdel, C. Elbuken and C. L. Ren,"Droplet formation in microfluidic T-junction generators operating in the transitional regime. II. Modeling", Physical Review E (1), 016323 (2012). 36. T. Fu, Y. Ma, D. Funfschilling, C. Zhu and H. Z. Li,"Squeezing-to-dripping transition for bubble formation in a microfluidic T-junction", Chemical engineering science (12), 3739-3748 (2010). 37. J. Sivasamy, T.-N. Wong, N.-T. Nguyen and L. T.-H. Kao,"An investigation on the mechanism of droplet formation in a microfluidic T-junction", Microfluidics and nanofluidics (1), 1-10 (2011). 38. M. De Menech, P. Garstecki, F. Jousse and H. Stone,"Transition from squeezing to dripping in a microfluidic T-shaped junction", journal of fluid mechanics , 141-161 (2008). 39. H. Yang, Q. Zhou and L.-S. Fan,"Three-dimensional numerical study on droplet formation and cell encapsulation process in a micro T-junction", Chemical Engineering Science , 100-110 (2013). 40. Y. Yan, D. Guo and S. Wen,"Numerical simulation of junction point pressure during droplet formation in a microfluidic T-junction", Chemical Engineering Science , 591-601 (2012). 41. S. Van der Graaf, T. Nisisako, C. Schroen, R. Van Der Sman and R. Boom,"Lattice Boltzmann simulations of droplet formation in a T-shaped microchannel", Langmuir (9), 4144-4152 (2006). 42. M. N. Kashid, A. Renken and L. Kiwi-Minsker,"CFD modelling of liquid–liquid multiphase microstructured reactor: Slug flow generation", Chemical Engineering Research and Design (3), 362-368 (2010). 43. H. H. Liu and Y. H. Zhang,"Droplet formation in a T-shaped microfluidic junction", J Appl Phys (3) (2009). 44. D. A. Hoang, V. van Steijn, L. M. Portela, M. T. Kreutzer and C. R. Kleijn,"Benchmark numerical simulations of segmented two-phase flows in microchannels using the Volume of Fluid method", Computers & Fluids , 28-36 (2013). 45. L. Sang, Y. Hong and F. Wang,"Investigation of viscosity effect on droplet formation in T-shaped microchannels by numerical and analytical methods", Microfluidics and nanofluidics (5), 621-635 (2009). 46. C. W. Hirt and B. D. Nichols,"Volume of fluid (VOF) method for the dynamics of free boundaries", Journal of computational physics (1), 201-225 (1981). 47. D. Gueyffier, J. Li, A. Nadim, R. Scardovelli and S. Zaleski,"Volume-of-fluid interface tracking with smoothed surface stress methods for three-dimensional flows", Journal of Computational Physics (2), 423-456 (1999). 48. D. J. Benson,"Volume of fluid interface reconstruction methods for multi-material problems", Applied Mechanics Reviews (2), 151-165 (2002). 49. B. Van Wachem and A.-E. Almstedt,"Methods for multiphase computational fluid dynamics", Chemical Engineering Journal (1), 81-98 (2003). 50. V. Badalassi, H. Ceniceros and S. Banerjee,"Computation of multiphase systems with phase field models", Journal of Computational Physics (2), 371-397 (2003). 51. R. Folch, J. Casademunt, A. Hernández-Machado and L. Ramirez-Piscina,"Phase-field model for Hele-Shaw flows with arbitrary viscosity contrast. I. Theoretical approach", Physical Review E (2), 1724 (1999). 52. M. De Menech,"Modeling of droplet breakup in a microfluidic T-shaped junction with a phase-field model", Physical Review E (3), 031505 (2006). 53. S. Chen and G. D. Doolen,"Lattice Boltzmann method for fluid flows", Annual review of fluid mechanics (1), 329-364 (1998). 54. D. Yu, R. Mei, L.-S. Luo and W. Shyy,"Viscous flow computations with the method of lattice Boltzmann equation", Progress in Aerospace Sciences (5), 329-367 (2003).
55. C. K. Aidun and J. R. Clausen,"Lattice-Boltzmann method for complex flows", Annual review of fluid mechanics , 439-472 (2010). 56. R. R. Nourgaliev, T.-N. Dinh, T. Theofanous and D. Joseph,"The lattice Boltzmann equation method: theoretical interpretation, numerics and implications", International Journal of Multiphase Flow (1), 117-169 (2003). 57. J. Brackbill, D. B. Kothe and C. Zemach,"A continuum method for modeling surface tension", Journal of computational physics (2), 335-354 (1992). 58. H. Rusche, Imperial College London (University of London), 2003. 59. E. Berberović, N. P. van Hinsberg, S. Jakirlić, I. V. Roisman and C. Tropea,"Drop impact onto a liquid layer of finite thickness: Dynamics of the cavity evolution", Physical Review E (3), 036306 (2009). 60. S. S. Deshpande, L. Anumolu and M. F. Trujillo,"Evaluating the performance of the two-phase flow solver interFoam", Computational science & discovery (1), 014016 (2012). 61. M. J. Nieves-Remacha, L. Yang and K. F. Jensen,"OpenFOAM Computational Fluid Dynamic Simulations of Two-Phase Flow and Mass Transfer in an Advanced-Flow Reactor", Ind Eng Chem Res (26), 6649-6659 (2015). 62. A. Q. Raeini, M. J. Blunt and B. Bijeljic,"Modelling two-phase flow in porous media at the pore scale using the volume-of-fluid method", Journal of Computational Physics (17), 5653-5668 (2012). 63. F. Raees, D. Van der Heul and C. Vuik, Report No. 1389-6520, 2011. 64. O. OpenCFD,"The Open Source CFD Toolbox", User Guide, OpenCFD Ltd (2009). 65. C. J. Pipe, T. S. Majmudar and G. H. McKinley,"High shear rate viscometry", Rheologica Acta (5-6), 621-642 (2008). 66. A. Gupta and R. Kumar,"Effect of geometry on droplet formation in the squeezing regime in a microfluidic T-junction", Microfluidics and Nanofluidics (6), 799-812 (2010). 67. T. Thorsen, R. W. Roberts, F. H. Arnold and S. R. Quake,"Dynamic pattern formation in a vesicle-generating microfluidic device", Physical Review Letters (18), 4163-4166 (2001). 68. A. Gupta, S. S. Murshed and R. Kumar,"Droplet formation and stability of flows in a microfluidic T-junction", Appl Phys Lett (16), 164107 (2009). 69. G. H. McKinley and A. Tripathi,"How to extract the Newtonian viscosity from capillary breakup measurements in a filament rheometer", Journal of Rheology (1978-present) (3), 653-670 (2000). 70. D. T. Papageorgiou,"On the breakup of viscous liquid threads", Physics of Fluids (1994-present) (7), 1529-1544 (1995). 71. J. Eggers,"Universal pinching of 3D axisymmetric free-surface flow", Physical Review Letters (21), 3458 (1993). 72. M. P. Brenner, J. R. Lister and H. A. Stone,"Pinching threads, singularities and the number 0.0304", Physics of Fluids (1994-present) (11), 2827-2836 (1996). 73. L. E. Rodd, T. P. Scott, J. J. Cooper-White and G. H. McKinley,"Capillary Break-up Rheometry of Low-Viscosity Elastic Fluids", Applied Rheology (1) (2005).(1) (2005).