Squeezing Droplet Formation in a Flow-Focusing Micro Cross-Junction

1. Introduction

In the last years, the proliferation of micro-fabrication technologies has provided a broad range of microfluidic applications in engineering. Micro-droplet technology has proven to be a promising and flexible platform for microfluidic functions, such as the production of mono-disperse particles, droplets, bubbles, foams and emulsions with precise control of components and sizes, which can give advancement in chemical, pharmacological, medical and industrial applications. Many microfluidic devices have been designed to generate uniform droplets, including geometry-dominated devices, flow-focusing devices, T-junctions, and co-flowing devices. However, the mechanism of droplet formation in micro-junctions is not yet fully understood. The two-phase flow characteristics are determined by flow conditions, fluids properties and geometry of the micro-device. The first key parameter in the control of the diameter of droplets created within a micro-junction is the flow rate ratio, i.e. the ratio between the flow rate of the dispersed phase to the flow rate of the continue phase [1]. Other key parameters are the geometry of the junction and the properties of the two fluids, such as the capillary number or the viscosity ratio [2]. Several micro-junction geometries are presented in the literature, such as the most studied micro T-junctions [3] and the micro cross-junctions [4]. For these devices, the dependence of the two-phase flow patterns by the flow rate ratio has been shown [5] and correlations have been provided between the droplet diameter to important key parameters, such as the flow rate ratio and the capillary number [6]. The mechanism of droplet formation, called droplet breakup has been largely studied since two decades in the micro-scale. Different mechanisms have been identified, such as squeezing, dripping and jetting [7]. By studying the underlying mechanisms in the droplet breakup, some scaling laws have been established, to predict the size of droplets produced in micro-junctions. However, more experimental data are needed to generalise the results. Based on statistical analysis of a large number of available literature data, Steegmans et al. [8] have shown that none of the scaling models, which are developed to predict droplet formation in a microfluidic T-junction, is general enough to describe the original data and data from other literature sources. Numerical approaches can be complementary to experimental investigation, if validated and integrated with the experimental measurements. The aim of this work is the prediction of the droplet generation in a micro cross-junction under the squeezing regime. The real micro-device considered has a restriction in the junction, i.e. the diameter of the micro-channels in the junction is smaller than the diameters of the micro-channels out of the junction. The dynamics of the droplet creation within the micro-junction is studied numerically by means of a VOF based code in OpenFOAM environment. The simulation results are compared with experimental data obtained on the real cross-junction by means of images obtained by a high-speed camera and micro-PIV measurements. The numerical results are used to build a model for the prediction of the droplet diameter, based on energy balancing during the breakup. A relation between interfacial tension, size of neck cross-section, and the time interval needed for droplet creation is obtained. By this approach, droplet diameters are predicted. The results obtained by the model are in a very good agreement with the the validated numerical simulations. This study can provide useful information for understanding micro-droplet dynamics, providing a basis for optimal design of multi-phase microfluidic devices.

2. Materials and methods

The micro-junction geometry considered for the numerical simulations refers to a real micro-device manufactured by Dolomite Microfluidics, shown in Figure 1 (top). The junction connects four micro-channels with a stadium-shaped cross-section, as shown in Figure 1 (bottom), with a restriction in the junction.

The channel width is $W=390.0\mu$m, while the channel depth is $H=190.0\mu$m. At the junction, the restriction has width $W_j=195.0\mu$m and depth $H_j=190.0\mu$m. The channel aspect ratio is $\beta =H/W=0.4871$, while the aspect ratio in the junction is ${\beta }_j=H_j/W_j=0.974$. The restriction ratio, defined as the ratio between the channel width at the restriction and the channel width is $\gamma =W_j/W=0.5$. For the numerical simulations, a computational domain, shown in Figure 2 has been built, characterised by three inlet branches with length $L_{in}=599.0\mu$m, and one outlet branch with length $L_{out}=1299.0\mu$m. The ratios between the inlet and outlet branches and the hydraulic diameter of the channels, $D_h=266.24\mu$m, are $L_{in}/D_h=2.25$ and $L_{out}/D_h=4.88$ respectively. A polyhedral mesh has been built, Figure 2 (bottom) shows the volume (in m³) of the cells, where the red zones at the three inlets are characterised by larger grid elements, while the outlet region is characterized by smaller elements, with a refinement in the junction where the breakup phenomena occurs, the total number of elements is about 5 milions.

The solver used to perform the numerical simulations is the interFoam solver provided by the opensource code OpenFOAM. This is a solver for two incompressible, isothermal and immiscible fluids, based on a VOF (volume of fluid) phase-fraction-based interface capturing approach.The equations solved are the continuity equation,

$$\nabla ·\mathbf{u}=0$$

(1)

and the Navier-Stokes equation,

$$\frac{\partial \left(\rho \mathbf{u}\right)}{\partial t}+\nabla ·\left(\rho \mathbf{u}\mathbf{u}\right)=-\nabla p+\nabla ·\tau +\rho \mathbf{g}+{\mathbf{f}}_{\sigma }$$

(2)

plus an additional equation to capture the interface between the fluids, the advection of the scalar quantity $\alpha$,

$$\frac{\partial \alpha }{\partial t}+\nabla ·\left(\alpha \mathbf{u}\right)-\nabla ·\left[\alpha (1-\alpha ){\mathbf{u}}_r\right]$$

(3)

The field $\alpha$ is used to distinguish the two fluids and is defined as

$$\alpha =\left\{\begin{array}{cc}1\hfill & \mathrm{in}\mathrm{the}\mathrm{continuous}\mathrm{phase}\hfill \\ 0.5\hfill & \mathrm{at}\mathrm{the}\mathrm{interface}\hfill \\ 0\hfill & \mathrm{in}\mathrm{the}\mathrm{disperse}\mathrm{phase}\hfill \end{array}\right.$$

(4)

To solve these equations, the MULES (multi-dimensional limiter for explicit solution ) algorithm is used in order to guarantee the boundness of the solution and obtain more smeared interfaces. The last term of the right hand side of Eq. 3 is the so-called compression term, which is not negligible only at the interface between the two fluids . In this term is present the relative velocity ${\mathbf{u}}_f$ defined as:

$${\mathbf{u}}_f=min(C_{\alpha }\left|\mathbf{u}\right|,max\left(\right|\mathbf{u}\left|\right))\frac{\nabla \alpha }{|\nabla \alpha |}$$

(5)

For micro-fluidic applications, $C_{\alpha }$ can be taken equal to 1 [10]. In Eq. 2 the source term ${\mathbf{f}}_{\sigma }$ is used to estimate the surface tension forces by means of the Continuum Surface Force (CSF) model, which gives

$${\mathbf{f}}_{\sigma }=\sigma \kappa \nabla \alpha$$

(6)

where $\kappa$ is the curvature, evaluated starting from the volume fraction $\alpha$,

$$\kappa =\frac{\nabla \alpha }{|\nabla \alpha |}$$

(7)

The dispersed phase is water, with density ${\rho }_d$ = 998 kg/m³ and viscosity ${\nu }_d$ = 8.788$·{10}^{-7}$ m²/s. The continuous phase is oil, with density ${\rho }_c$ = 950 kg/m³ and viscosity ${\nu }_c$ = 1.900$·{10}^{-5}$ m²/s.

The dispersed phase enters from the channel along the x-axis and the continuous phase enters from the two channels along the z-axis. The dispersed phase flow rate has been varied in the range $0.6-7.5$ ml/h, while the continuous phase flow rate has been fixed to $7.5$ ml/h. The regime is laminar, as the Reynolds number calculated for the dispersed phase is $R{e}_c=v_c{D}_h/{\nu }_c=0.3$ and $R{e}_d=v_d{D}_h/{\nu }_d$ ranges from 0.5 to 6. The capillary number for the dispersed phase, $C{a}_d={\mu }_D{u}_d/\sigma$ ranges from 0.00005 to 0.0006, where $\sigma =0.04244$ N/m is the surface tension between the two fluids [6]. The capillary number for the continuous phase is $C{a}_c={\mu }_D{u}_d/\sigma =0.0134$. Then, the regime for droplet creation is the squeezing regime, for all the cases in these ranges [6].

Five meshes have been built to check the grid convergence. The convergence has been proven by comparing the velocity and the diameter of the droplets as a function of the number of elements, as shown by Figure 3.

3. Validation of the numerical simulations

A schematic and a picture of the experimental apparatus used in this work are shown in Figure 4. Two syringe pumps (Harvard Instruments PHD 400) were used to control the flow rates of the two working fluids. The pumps were connected to the microfluidic-chip with the micro-junction (Dolomite Microfluidics) which was placed over an inverted microscope (Nikon Eclipse TE2000-U) and observed through a 10× objective lens. High-speed video recordings of the droplets in the microfluidic device were taken using a high-speed camera (Olympus I-Speed I). The illumination was provided by a high-power LED powered by DC current. The shapes and dimensions of the droplets were determined by a custom-made image processing code developed in the Python environment. The two-dimensional velocity field in the mid-plane of the dispersed phase was measured utilizing the micro-PIV technique [11]. For these measurements, polystyrene micro-spheres with a diameter of 1.19 $\mu$m and a density of 1050 kg/m³ were inserted at low concentration in the dispersed phase and used as a passive tracer. The micro-PIV analysis was performed using the open-source library DefocusTracker in the MATLAB environment [12].

In order to validate the simulations, the dimensions of the droplets obtained numerically were compared with the dimensions of the droplets obtained experimentally for the same working conditions. A qualitative comparison is presented in Figure 5 for three different values of the dimensionless flow rates showing an excellent agreement between numerical results and experiments (with a maximum estimated error of 5 %). Moreover, the dynamics of the breakup phenomena in the experiments was compared with the results of the simulations. Figure 6 shows the 3D shapes of the interface colored with the velocity magnitude, to highlight the zone where the acceleration is higher, and compared with the corresponding picture in the experiments. The comparison shows that the simulations can accurately predict the drop geometries observed in the experiment, with an error on the time needed for breakup of about 5.3 %.

Finally, a comparison between the experimental and simulated velocity field in the mid-plane of the dispersed phase was performed. Figure 7 shows the velocity vector map obtained by the numerical simulations (left), together with the velocity vector map obtained by micro-PIV analysis (right), which has been mirrored. The vectors show a vortex in the thread nose, which is well known for droplets moving in two-phase flows [5]. The position and the dimensions of the vortexes obtained by the numerical simulation are in agreement with those obtained by the micro-PIV measurements. Also the velocity magnitude of the analysed fields is in agreement. In summary, we can conclude that the numerical simulations are in good agreement with the experiments both for the prediction of the dynamics of the droplet formation and for the velocity distribution.

4. Interface dynamics during droplet breakup

The dynamics of the droplet during the breakup is analysed in this section, starting from the case with a dimensionless flow rate $Q^{\ast }$, defined as $Q^{\ast }=\frac{Q_d}{Q_c}$, equal to 0.08 and $C{a}_c=6.7·{10}^{-3}$ and increasing the dispersed phase flow rate. During the breakup, three stages can be observed. In a first phase, called filling stage, the dispersed phase evolves as a thread which starts with a hemispherical shape, then the thread increases with the shape of a cylinder with a rounded front end which we will call the thread nose. In the second phase, called necking, the thread starts to be squeezed while the thread nose position moves along the junction. The final phase is the pinch-off, when the neck thickness decreases quickly until the thread surface breaks and the droplet is created. Figure 8 shows the thread interface at four time instants equally spaced from $t=0$ to $t=0.0012$ s, during the first two phases.

The thread at the final instant is also shown in the figure (right) to show the point P1, which represents the position of the thread nose, and the point P2, which represents the position of the neck. Figure 9 shows the evolution in time of P1 (left) and the evolution in time of the thickness of the neck (right) (2P2) as a function of the dimensionless time.

The dimensionless time has been obtained by dividing the time by the time interval between a droplet detachment and the successive droplet detachment. The figure (left) shows that the filling and necking phases last for most of the total time, while P1 increases slowly until $t^{\ast }=0.85$, (black dots), while for $t^{\ast }>0.85$ it increases faster. The figure (right) also shows that the neck thickness decreases slowly for $t^{\ast }<0.9$ and drops very quickly for $t^{\ast }>0.95$. The first time interval corresponds to the filling and necking phases, while the second one characterise the pinch-off stage. Figure 10 shows the evolution of the thread interface with time, from the first instant where the neck starts to be visible until the last time instant before the breakup. The dimensionless time has been obtained in this case dividing by the time duration of the necking stage.

The figure shows that during the necking phase, the squeezing of the thread in the neck is slow and the diameter of the thread can be assumed as a constant for a large time interval before the pinch-off phase. Similar trends can be shown at different flow rate ratios, as shown in Figure 11.

In the Figure 11, the evolution of P1 and P2 are shown for $Q^{\ast }$ = [0.08, 0.15, 0.25, 0.35, 0.65]. The figure shows that P1 increases proportionally to $Q^{\ast }$, in the necking phase, while the neck thickness P2 does not depend on $Q^{\ast }$, i.e. the neck collapses with the same trend in all the cases. This can suggest that the duration of the breakup is an key-parameter of the droplet creation process. During the breakup phenomena, the main forces acting on the drop are shear-stress at the interface and the pressure differences between the two phases. In the squeezing regime the effect of the pressure build-up is the most important [7], because the thread can occlude the micro-channel in the junction, then the pressure in the continuous phase increases, breaking the interface. To investigate this phenomena, the pressure difference between the phases is plotted along the interface during the process, as shown in Figure 12. The plot is made in polar coordinates, where the angle 0 is on the thread nose and the angle 180 is on the inlet section.

In Figure 12, the red dashed lines indicates the position of the restriction in the micro-junction. The figure shows that, when the thread occupies a large portion of the channel, the pressure difference increases, while it decreases when the forming drop has already crossed the restriction in the junction. In a similar way it is possible to describe the shear-stress acting on the thread surface, as shown by Figure 13 for the same time instant as those shown in Figure 12. The plot is made in polar coordinates, where the angle 0 is on the thread nose and the angle 180 is on the inlet section.

Figure 13 shows that the shear is higher when the thread occupies the restriction in the first time instants, as the continuous flow rate reaches the maximum velocity due to the reduction on the cross section. Then, the shear decreases when the cross-section free for continuous phase passage increases. Figure 13 and Figure 13 show that the shear stress is lower than the pressure difference between the phases for all the time instants considered.

5. Model of the droplet formation

The length of the droplets obtained by the simulations in function of the flow rate ratio are reported in Table 1. In the table, ${\dot{Q}}^{\ast }$ indicates the flow rate ratio and $l^{\ast }=l/W_j$ is the dimensionless length of the droplet.

Many papers in the literature show linear correlation between $l^{\ast }$ and ${\dot{Q}}^{\ast }$ for T-junctions [13,14] and cross-junctions [15]. By fitting the dimensionless length obtained by the simulations we obtain the following relation

$$\frac{l}{w_j}=\alpha Q^{\ast }+\beta =0.382Q^{\ast }+0.905$$

(8)

with an average error about 1.3%.

In order to find a model which supports this trend the dynamics of the droplet formation can be considered, following the approach shown in [16] for predicting the time interval which involves the necking stage. The approach is applied to the entire necking stage, as it affects the the final droplet size more than the filling stage, as shown by [17]. The control volume considered in this approach is shown in Figure 14.

The control volume is a cylinder containing the volume of the thread at the beginning of the necking stage, as shown by the figure. This region is inside the restriction in the junction. The final droplet volume is assumed to be given by the product of the dispersed phase volume flow rate multiplied by the time lapse between the beginning of the necking phase and the pinch-off. The droplets created in the junction show a circular shape by the high-speed camera images which record the the droplets on a plane $xy$ containing the direction of the flow, as shown in Figure 5. This means that the volume of the created droplet is a sphere in the case a droplet diameter (that we will call length of the droplet $l_d$) smaller than the height of the channel H, and is a cylinder if $l_d>H$. Then, the volume of the droplet is $V_d=4/3\pi r_d^3$ if $l_d<H$, or $V_d=\pi r_d^{2}H$ if $l_d>H$. Assuming that the time interval between the beginning of the droplet creation and the time instant when the droplet detaches from the thread is $\delta t$, the volume of the droplet is

$$V_d=V_0+{\dot{V}}_d\delta t=\frac{\pi l_d^{2}H}{4}$$

(9)

where ${\dot{V}}_d$ is the dispersed volume flow rate and $V_0$ is the volume of the thread at the end of the filling stage, i.e. at the beginning of the necking stage, that is the stage described by the modelling in this section. Then, the length of the droplet $l_d$ after the detachment can be obtained by the dispersed phase volume flow rate

$$l_d={\left(\frac{4(V_0+{\dot{V}}_d\delta t)}{\pi H}\right)}^{1/2}$$

(10)

The time interval for the droplet creation during the necking stage can be obtained by the energy balance approach introduced by [16]. The surface tension energy stored in the thread before the droplet detachment is balanced by the energy difference between the inlet and the outlet sections in the control volume after the droplet detachment,

$$\sigma S_d=M_d\left(\frac{pc}{{\rho }_d}+\frac{u_d^2}{2}-\frac{p_d}{{\rho }_d}\right)$$

(11)

where $\sigma$ is the surface tension between the two phases, $S_d=2\pi rl$ is the surface of the droplet in the region of the neck approximated by a cylinder with radius r and length l, and $M_d=\pi r^{2}l{\rho }_D$ is its mass. If $p_c$ and $p_d$ are the pressure in the continuous and dispersed phase respectively, and $u_d$ is the velocity of the droplet after the detachment, on the right of the equation we have the energy of the droplet after creation. Equation 11 can be rewritten as

$$\frac{u_d^2}{2}=\frac{2\sigma }{r{\rho }_d}+\frac{p_d-p_c}{{\rho }_d}$$

(12)

and considering that $\frac{p_d-p_c}{{\rho }_d}=\frac{\sigma }{{\rho }_{d}r}$ one obtains the velocity of the dispersed phase:

$$u_d={\left(\frac{6\sigma }{{\rho }_{d}r}\right)}^{1/2}$$

(13)

Assuming that $u_d=l/\delta t$ then one obtains

$$\delta t=l{\left(\frac{{\rho }_{d}r}{6\sigma }\right)}^{1/2}$$

(14)

Substituting equation 14 in equation 10 one obtains

$$l_d={\left(\frac{4}{\pi H}\left(V0+V_{d}l{\left(\frac{{\rho }_{d}r}{6\sigma }\right)}^{1/2}\right)\right)}^{1/2}$$

(15)

Dividing the droplet length by the width of the junction $W_j$ and rewriting the equation in a dimensionless form one obtains

$$l^{\ast }=l_0^{\ast }+a\dot{(}{Q}^{\ast }{)}^{1/2}W{e}_c^{1/4}$$

(16)

where $W{e}_c=\frac{{\rho }_c{v}_c^2{W}_j}{\sigma }$ is the Weber number referred to the continuous phase in the junction and $v_c=\frac{4{\dot{V}}_c}{\pi W_j^2}$ is the superficial velocity of the continuous phase in the junction. In equation $a={\left(\frac{l^2{\rho }_d}{6H^2{\rho }_c}\right)}^{1/4}$ is a constant with $l=4W_j$ and $l_0^{\ast }=0.8$ is obtained by the initial volume of the thread. Equation 5 is in a very good agreement with the measurements, as shown by Figure 15.

In the figure, the dimensionless length of the droplets obtained by the numerical simulations are shown by the symbols, while the continuous line refers to the results obtained by Eq. 8 and the dashed line represents the linear fit 5. The error between the numerical results and the Eq. 5 is less than 1%, showing a better fit with respect to the linear correlation. This result is in agreement with what has been shown by [16,18]. The dependence on the Weber number evaluated in the junction emphasises the influence of the restriction of the junction which forces the continuous phase to flow at a higher velocity in the region of the neck. Similar results have been shown by Liu et al. [19] for cross-junctions where the two channels carrying the continuous phase have a smaller width with respect to the central channel for the dispersed phase inlet. This results can be extended to micro cross-junctions with restriction and flow-focusing device, showing to be potentially used for optimal design of micro-machines for particles and droplets manipulation.

6. Conclusions

The dynamics of the droplet formation in a micro cross-junction with a restriction has been studied by an integrated approach between CFD simulations in the OpenFOAM environment, experimental measurements, and an up-scaled model of the droplet breakup. The results of the simulations have been validated by a comparison with experimental measurements by means of high-speed camera images of the droplets and velocity measurements in the dispersed phase by micro-PIV. By the numerical simulations, the forces acting on the dispersed phase during the droplet formation and the droplets diameters have been obtained for different dimensionless flow rates. These results have been used for the development of a model of the droplet breakup under the squeezing regime, and a novel correlation between the dimensionless length of the droplet and the dimensionless flow rate has been obtained. These results can be generalised to similar flow-focusing cross-junctions, showing that the methodology introduced in this paper can be used for the optimization and control of droplets production in microfluidics applications.