Parameters influencing the droplet formation in a focusing microfluidic channel

. In the present work a detailed numerical study of the parameters influencing the droplet formation in a flow-focusing microfluidic device are made. First, an extensive verification of the simulations with data from the literature is carried out. Influence of parameters like viscosity and inflow velocity are compared with the results from literature showing a good agreement. Some differences are attributed to the different numerical techniques used: in the present work a pure volume-of-fluid method is used, while in the reference study this method is combined with the level-set method. As a second step of the verification of the present model, a comparison with experimental data from the literature was carried out which shows a very good agreement. After the verification was completed, eight new simulations are carried out covering a wide range of velocities of the continuous phase u c . In these simulations the velocity of the discrete phase u d remains unchanged. The variation of the continuous phase velocity reveals that with increasing the value of u c , respectively the value of the capillary number Ca , the droplet length reaches a point of saturation, i.e. a point where the droplet length does not decrease any more. For the present setup this saturation occurs for Ca > 0,03. On the other hand, when the velocity of the continuous phase goes towards very low values ( Ca < 0,01 for the present setup), the droplet size increases significantly. Further, it was found that for increasing capillary numbers Ca above a value around 0,015 for water/oil and above 0,025 for water + 40% glycerol / oil systems, a transmission from the dripping towards the jetting regimes of droplet formation occurs. It was shown that when the viscosity of the continuous phase increases, higher total pressure jumps in the droplet occur, also leading to the formation of smaller droplets.


Introduction
In recent years, microfluidic devices have emerged as a novel tool for the realization of different biological, chemical and medical processes. The decreased reagents volumes, usually in the nano/picolitres range, can significantly reduce reaction times, costs and energy consumption for certain processes [1]. Droplet-based microfluidics is a novel technique, which allows the encapsulation of different chemical or biological compounds into single picolitre-droplets. This makes it possible for the isolation and control of individual reactions, protecting them from their surroundings and unwanted interactions with other components. These small microreactors have allowed an easily implementable and relatively cheap approach for a broad range of processes including cell lysis [2], DNA purification [3], polymerase chain reaction (PCR) [4] and other. Therefore, understanding the mechanisms behind droplet generation is crucial for the control and efficiency of a droplet-based microfluidic system.
Droplets formation usually occurs when two immiscible liquids intersect each other. The process is achieved either by active (using electric fields) or by passive (using pressuredriven flow and channel geometry) methods. The specific design of the microfluidic channels makes it possible for a dispersed phase, usually water, to be sheared by another continuous phase (hydrocarbon oils, fluorocarbon oils etc. [5]) and to produce uniformsized drops. Different geometry designs are possible, the most frequently used in the praxis are shown in Figure 1: T-junction [6][7] (a), flow-focusing [8][9] (b), and co-flowing [10] (c). The present work studies devices working on the flow-focusing principle. Prior knowledge of droplet size, shape, formation frequency or pressure drop are usually required for the work with devices based on droplet-based microfluidic. Several experimental studies have already investigated the mechanisms behind droplet generation in microfluidic channels [11][12][13][14][15]. Those measurements give an adequate information about the needs of a certain experimental setup; it is not possible however, to include the influence of all effective parameters on the system. A feasible way to collect in advanced the needed information is by utilizing a predictive CFD model. A couple of numerical studies on droplet-based microfluidics were carried out in the past few years making use of different numerical methods (level set (LS) [16], Lattice Boltzmann method (LBM) [17][18] and other). In this article the droplet formation in a flow-focusing microfluidic channel is investigated numerically by utilizing the volume of fluid (VOF) method. First, we validate our results with data from numerical and experimental studies from the literature. In a second part, using the same geometry as in [19], a parameter study is carried out. It consists of eight  new velocities of the continuous phase while at the same time the velocity of the dispersed  phase was kept unchanged. Thus, the critical capillary number Ca, where the transition  between the dripping and the jetting regimes occurs, can be determined. The pressure  variation along the center line of the channel for the two flow regimes, is presented and discussed.

Mathematical model, geometrical setup, boundary conditions and numerical details
In the present study three-dimensional simulations of droplet formation in a flow-focusing geometry as shown in Fig. 2 In the equations above, U is the velocity vector field and p is the pressure field. FS represents the surface tension force and is defined as shown in Eqn. 4 and is nonzero only at the interface between the two phases: Here N represents the curvature (N = ∇ ∇. ( ∇ߙ /~∇ߙ~)). Only one such transport equation (Eqn. 3) is solved since the volume fraction of the other phase can be inferred from the limitation: where the index 'c' stands for continuous and 'd' for dispersed phase. The viscosity P and the density U are based on the weighted average of the phase fraction: In order to verify the results from the VOF method, a geometry of a flow-focusing microchannel, already investigated by Sontti and Atta [19] is considered. The lengths and dimensions of the square cross-section, the inlet and outlet channels are presented in Figure  2.
The continuous phase (water/water + 40% glycerol) is introduced through the two side channels and the dispersed phase (oil: octane +2,5% SPAN 80) is entered from the main (central) channel. Table 1 summarizes the physical properties of the two oil-water systems considered for the verification of the model and the later calculations in our work. As in [19], the information is obtained from the experimental results of Yao et al. [20]. A finite volume method based CFD solver from ANSYS Fluent 16 with the implemented Volume of Fluid technique is used to solve the system of time-dependent partial differential equations. For the boundary conditions, constant velocity block profile was utilized for both continuous and dispersed phase inlets. The phase fraction D is set to D =0 at the two inlets of the continuous phase and D = 1 at the inlet of the dispersed phase.
The walls of the channels were considered as fully wetted by the continuous phase, therefore the contact angle was set as zero at the walls. No slip boundary conditions are applied at the walls. Pressure boundary was specified at the outlet of the main channel. Table 2 summarizes the varied boundary conditions for each case investigated. The first two cases serve as a verification of our model. All other cases are completely different from [19] and cover a large range of new velocities for the continuous phase. The purpose is to investigate in detail how these velocities influence the flow regimes (dripping (D) or jetting (J)), which will be introduced in the following sections. The inlet velocity of the dispersed phase (oil) was kept constant at ud = 0,0185 m/s for all 10 cases.  Making use of the mesh sensitivity analysis made by Sontti and Atta [19], we considered a near wall mesh refinement with a first row of thickness of 6 μm and continuing with a 13% size increase towards the middle of the channel. The middle elements reach 35 μm length size, and an overall of 390 000 elements are required for the hexahedral equidistant mesh. Figure 3 shows a zoomed view of the hexahedral elements in the cross region where the two fluids meet.

Verification and experimental comparison
As mentioned before, in order to examine the efficiency of the VOF model, first a verification of the results with the work of Sontti and Atta [19] is made. The authors of [19] utilized a combination of the VOF and LV (Level Set) methods, which compared to the pure showed more accurate methodology in capturing the interface between the two fluids. A detailed comparison between the VOF and the CLSVOF models is made by Keshavarzi et al. [21]. As a next step, a visual comparison with the experimental work of Wu et al. [22] was carried out, in order to further confirm the accuracy of the model. For our verification, cases 1 and 2 from Table 3 are considered and compared with identical simulations from [19]. Figure 4 shows the comparison of the droplet formation in the middle plane of the domains through time in the two studies for case 1. Both results show very close time scales regarding the generation of oil (blue) in water (red) droplets. The whole process can be divided into three stages. In the first, so-called filling stage (0,020s -0,040s), the dispersed phase is injected into the main channel. At some point the growing oil-front blocks the flow from the side channels causing the upstream pressure to increase until it reaches a critical value where the continuous phase begins to squeeze the interface [23]. In the second stage (0,040s -0,055s) the dispersed phase (oil) is still being injected into the droplet at a constant flow rate, while the neck starts to collapse. This accelerates the last pinch-off stage, in which droplet detaching occurs (around 0,060s). The three stages described above form, the so-called dripping regime in a microfluidic device [24]. Fig. 4. Comparison of the droplet formation in our work (left) and Sontti and Atta [19] (right) for case 1. Results from [19] are reproduced with permission Similar observations are made in the second case (case 2 of Table 2), in which water +40% glycerol is utilized as continuous phase. The comparison is shown in Fig. 5. For this case the generated droplets have smaller volumes and higher formation frequencies compared to pure water (case 1). At the same operating condition for oil-water+40% glycerol system, an earlier droplet breakup (0,051 s) is also observed. This is due to the higher shear force, 'cutting' the droplet, connected with the higher viscosity of the red colored fluid in the second case. A larger difference in viscosity makes it possible to increase the shear at the boundary between the fluids thus increasing the generation frequencies. [19] (right) for case 2. Results from [19] are reproduced with permission.

Fig. 5. Comparison of the droplet formation in our work (left) and Sontti and Atta
A quantitative comparison of the two cases considered is also carried out and summarized in Table 3. The nondimensional droplet length with respect to channel width (= LD/Wc) and the droplet volume VD are utilized for this purpose. For the calculation of VD Equation 8 is used, with f standing for the frequency of droplet formation in Hertz and Qd being the volume flow rate of the dispersed phase: A relative difference of around 13-14% is observed for the volume of the generated droplets and around 10% for their length with respect to the channel width. This is in a reasonably good agreement with the CLSVOF method presented in [19], as the VOF technique often tends to smear the step profile of the interface over several mesh cells because of numerical diffusion, as it has been demonstrated in [22].
To further confirm the accuracy of the utilized VOF method, an additional comparison with the experimental work of Wu et al [22] was carried out. Numerical simulations were performed with the same fluid properties and channel dimensions, considered by the mentioned authors. Figure 6 demonstrates the visual comparison of droplet development.
Here, the velocity of the continuous phase is uc = 0,00252 m/s, the velocity of the discrete phase -ud = 0,00084 m/s and the surface tension is σ = 30 mN/m.

Influence of continuous phase velocity and significance of the capillary number for the transition of regimes
In order to understand better how droplet size for the flow-focusing microdroplet generator could be controlled more precisely, eight additional cases with varied continuous phase velocities (at constant dispersed flow rate of ud = 0,0185 m/s) were carried out, extending the range of simulated values by Sontti and Atta [19]. First, the value of uc was decreased two times, compared to cases 1 and 2, to further evaluate the characteristics of the drops in the dripping regime (cases 3 and 4 from Table 2). Figure 7(a) shows the pressure variation along the center line of the channel for these two cases. Due to the increase in the internal friction of the continuous phase in case 4, a bigger total pressure in the domain (around 40%) is observed when compared to case 3. In any pressure profile, each peak and its width indicate a droplet position and shape, respectively. Due to their smaller size and hence, larger surface curvatures, the droplets from case 4 have on average approx. 15% higher pressure differences with their surrounding compared to case 3. Table 2: (a) the total pressure profiles along the channel centerline ; (b) a three dimensional iso-surface evolution of droplet shapes for case 3; (c) a three dimensional iso-surface evolution of droplet shapes for case 4 Similar to cases 1 and 2, the droplet length in cases 3 and 4 decreases with increasing continuous phase viscosity, when comparing the three-dimensional iso-surface evolutions. As seen in Figure 7(b) and Figure 7(c) both cases are again in the dripping flow regime, described in section 3.1.

Fig. 7 Comparison of case 3 and 4 from
In the next step the velocity of the continuous phase is increased four times, compared to cases 1 and 2 (again for ud = 0,0185 m/s). This way we were able to further investigate the effect of uc on the droplets and answer the question if the dripping regime take up even at higher velocity values. Compared to Figure 7, even larger differences in the total pressure can be seen between cases 5 and 6; these are shown in Figure 8(a). The bigger curvature radii lead to bigger resulting capillarity forces and hence, to three times higher pressure differences between the inside and the outside of the droplets as compared to cases 3 and 4. A noticeable pressure gradient across the droplet itself can also be seen. This is due to the difference in radius of curvature from back to nose of a droplet, shown also in the theoretical analysis made by Abiev et al. [25]. From Figures 8(b) and 8(c) a clear change in the shape and size of droplets is visible. Due to the increased "cutting'' shear force on the dispersed phase, droplet length and volume for case 6 (when compared to case 5) decrease and higher formation frequencies occur. In this, so called jetting regime, a longsuspended column of the disperse fluid flows in the central channel as shown in Figures  8(b) and 8(c). The breakup of drops does not occur due to shear stress but due to the pressure drop across the emerging drop [26]. The question where the transition from the dripping to the jetting regimes occurs has been discussed in the literature [27]. One of the critical parameters thoroughly investigated in droplet producing microfluidic devices is the capillary number, described in Equation 9: where μc is the dynamic viscosity of the continuous phase. The number identifies the ratio of viscous to interfacial forces. For low capillary numbers, capillary forces dominate and therefore big droplets will be formed, whereas for high capillary numbers viscous forces generate small drops, increasing their surface energy. Similar to the Reynolds (Re) number, the Ca number can be used to determine a transition point between two flow regimes. One way to do this is to show the geometrical changes in the generated droplets as a function of Ca. Figure 9 shows the non-dimensional length of the droplets for the 10 cases as a function of the Capillary number. In the Figure it can be seen that the largest changes of the drop size occur up to Ca ≈ 0,02. The Figure shows that for even larger Ca numbers (Ca > 0,03) the droplet size stagnates, i.e. it does not change significantly. Such a qualitative behavior is observed also in other studies with other geometries (e.g. Figure 6 from [28]).
At the left-hand side of the figure, the size of the droplets increases significantly when the Ca value decreases under the value 0,01. W. For water (blue) this value is around Ca ≈ 0,015 and for water + 40 % glycerol (red) it is Ca ≈ 0,025. The blue and red symbols show the same qualitative trend: larger slope for smaller Ca numbers and an almost constant droplet size for larger Ca numbers. The difference in the transition points of the two fluids is attributed to the different fluid properties. Although the Ca number of the continuous phase is the most important parameter that indicates the transition, there are other properties (viscosity of the dispersed phase, surface tension) which also influence the flow regimes. Fig. 9 The nondimensional droplet length with respect to channel width as a function of Ca It is important to note that for the definition of the Ca number, the average domain velocity -and not the inlet velocity -of the continuous phase was utilized. Table 4 summarizes the averaged domain values for the continuous phase which are used in Equation 9 and Fig. 9.

Conclusion
In this work the droplet generation in a flow-focusing microfluidic device has been investigated. The continuous phase (water, or water + 40 % glycerol) was introduced through the two side channels and the dispersed phase (oil: octane +2,5% SPAN 80) was entered from the main channel. For all simulations the VOF method was utilized. In the first part, a verification of the technique with the results from Sontti and Atta [19] has been made. Here, two cases were simulated, in the first, pure water was used as a continuous phase, and in the second -water+40 % glycerol which leads to increased viscosity. A good agreement between the two works, with maximum differences of 13% is found which are attributed mainly to the different numerical techniques. To further confirm the accuracy of the utilized VOF method, an additional comparison with the experimental work of Wu et al [22] was carried out which shows a very good agreement.
In the second part of the presented study, the influence of velocities of the continuous phase on the droplet formation is investigated. Two main flow regimes are observeddripping and jetting, which are distinguished by the size of the droplet (for jetting, the droplets are smaller than the width of the channel). The transition between these regimes is investigated. This is done for a wide range of velocities uc (0,1400 m/s > uc > 0,00925 m/s). The large number of studied regimes enabled the following main observations for the investigated setup: x At very large velocities uc (Ca > 0,03) droplet size reaches a point of saturation, hence droplets do not decrease their size with increasing velocity uc any more.
x At very small velocities uc, (Ca < 0,01) the observation is that droplet size still increases strongly with decreasing velocity. x There is a critical Ca number, which serves as a transition point between the two described flow regimes (dripping and jetting). For the cases with water as a continuous phase, a critical value around Ca ≈ 0,015 is observed, whereas for the cases with water + 40 % glycerol it is Ca ≈ 0,025. This difference is due to the influence of the dispersed phase properties. x It was found that at higher flow rates (velocities) higher total pressures and pressure jumps in the droplets occur, also leading to the formation of smaller drops. The same effect was observed with increasing continuous phase viscosity.
As a whole, the present study shows that the VOF method is a reliable technique for the simulation and prediction of droplet generation in flow-focusing channels. It will allow the future study of other diverse setups and of various fluid combinations.