Numerical simulation of droplet formation in a microchannel device

The formation of droplets is a phenomenon of particular importance in the development of industrial emulsions. The quality of these compounds is associated with droplet size and stability over time. Anna et al. (2003) developed a methodology named  ̈flow focusing ̈ to improve droplet formation processes for engineering applications. In this work, Computational Fluid Dynamics (CFD) based techniques are used to assess the capacity of a pseudo-2D numerical model to reproduce the formation of water droplets within silicon oil, as obtained during experiments. Average time of droplet onset obtained via numerical analysis was 1.5 times larger than observed experimentally, whereas droplets convection velocity and diameter predictions differed by 40-45% and 60%, respectively. Nevertheless, the calculated velocity profiles downstream of the discharge slot reproduced the expected free-jet shear layer according to outer/inner flow ratio.


INTRODUCTION
Microfluidics applications have spread out since two decades ago.It is a promising area of mechanics that extends to several disciplines such as physics, chemistry, medicine, cosmetics, biotechnology, and food industry, among others.The main appeal of microfluidics is the possibility of creating small components manufactured to produce moderate volume of complex fluids with highly controlled properties.
In recent years, the phenomenon of droplet formation has been researched for the purpose of studying its stability, specifically in micro-emulsions, related to the conservation of their sizes.Many monodisperse emulsions require a composition of uniform size drops from 0.1 microns to millimeters, with a coefficient of variation CV (defined as (average standard deviation/average droplet size)*100) usually less than 5% (Kobayashi & Nakajima, 2009).
Droplet formation at microscale levels may be accomplished through: a) flow focusing, characterized by confluence channels, i.e., three or more lines parallel to the entrance channel are subject to a contraction; and b) "T" junctions which are based on microscale cross-flow techniques.
These technologies have been studied both empirically and simulated numerically (Zhou et al. 2006, Anna et al. 2003).The application of computational tools reduces the typically high costs in the experimental setting, associated with the manufacture and precise manipulation of microscale geometries and associated fluids at this scale.The numerical simulation provides also qualitative and quantitative information to develop models to characterize the relation between dependent and independent variables, by working as a "virtual laboratory" (computer + model).Anna et al. (2003) developed a model experiment considering the "flow focusing" in a microfluidic device for monodisperse and polydisperse emulsions.
They illustrated the qualitative aspects needed to control the size and distribution of the drop as the inner and outer flow rates of the two liquids vary.
Their experimental microchannel is flat and manufactured through soft lithography methods.The experiment consists in forming water droplets (dispersed phase) in oil (continuous phase).A stream of water enters the inner channel and a stream of oil enters through the two outer channels.The two phases are forced to flow through a small orifice which is located downstream of the three input channels.The pressure and shear viscous forces exerted by the fluid form a thin thread, which breaks down into the orifice or just downstream of it as the thread becomes critically unstable.This paper presents a numerical model aimed to reproduce the phenomenon of "flow focusing" in a pseudo-2D geometry mimicking the experimental facility of Anna et al. (2003), with the main purpose of evaluating the prediction capacity of current numerical tools using inexpensive hardware.

NUMERICAL ANALYSIS
This study was carried out using ANSYS-CFXv12.The computer equipment used was a Hewlett-Packard model with an Intel ® Core (TM) 2 Quad 2.83GHz CPU and RAM 3.5 GB.Storage memory needed to store the generated files is 1TB.The approximate time for a simulation resolution is 30 days.

Geometry Creation
The geometry used was developed with similar lengths to those of the experiment by Anna et al. (2003) as shown in Figure 1.
It consists of three parts: • The co-flow channels, consisting of three channels through which fluids enter to be analyzed, silicone oil enters into the external channels and water flows in the inner channel.

•
Contraction, with a diameter of 43.5 microns.At this point, the drop forms.

•
The discharge channel, is located downstream of the contraction and is where the fluids coming through the co-flow channel exit.

Mesh Generation
CAD geometry is built out of blocks, as shown in Figure 2, to facilitate the structured meshing, shown in Figure 3.

Domain verification analysis
As is customary, domain and mesh verification analyses were performed to assess the appropriateness of the numerical set up.The domain length, downstream of the discharge slot, was set at 500, 600 and 750 µm, to verify the independence of results with this parameter.
Considering that the contraction width ¨D¨ is 43.5 µm, for every length, plots of velocity profiles were made at 1D, 2D and 4D downstream from the contraction and were compared to each other, as shown in Figures 4, 5 and 6, respectively.
The minimum difference was found between centerline velocities between the 600 µm and 750 µm domain lengths.This difference was around 5%, for which the former domain length was adopted for the final simulations.

Mesh verification analysis
The Grid Convergence Index (GCI) methodology (Celik, 2008) was used for mesh verification purposes.
The GCI was applied at three levels of mesh density (coarse, medium and fine) in order to estimate the order of convergence (Table 1).The rate of convergence is computed on shear stress, pressure drop and maximum speed at the contraction where the droplet formation occurs.
The GCI for both grids (medium-to-coarse and fine-to-medium) is less than 3%, however the fine-medium mesh is applied for two-phase model simulations, since relative error is less than 1% (Table 2).

Simulation Model 2.1 BOUNDARY CONDITIONS
The geometry has three inlet channels and one exit.The input is divided into: external channels, through which silicone oil is injected at 4.3E-2 m/s, and an internal channel, through which the water is injected at 1.52E-3 m/s.
An average pressure of 0 Pa is set at the outlet (Figure 7).Symmetry is established on the bottom and top panels, to guarantee a pseudo-2D model, using one volume element in the third dimension.
Non-slip boundary conditions are set for the rest of the boundaries.

INITIAL CONDITIONS
A fully converged steady state single-phase model, using oil as the working fluid, is used as the initial condition for the simulation of the transient two-phase model.The transient analysis is simulated with a time step of 1•10 -6 s.

Numerical model
A second-order high resolution scheme is applied to discretize the convective term in the Navier-Stokes equation.
A homogeneous segregated two-phase model with a sharpening interphase algorithm is used to capture the liquid-liquid interphase.
The Young-Laplace equation is modeled through the Brackbill (1992) approximation of a continuum surface force, using a surface tension of 0.001 N/m.

RESULTS AND ANALYSIS
The comparison between numerical results and experiments is based on the droplet diameter (D) and droplet-to-droplet distance (L) at their inception just after the contraction.

DATA EXTRACTED FROM THE EXPERIMENTAL WORK ANNA ET AL. (2003)
Experimental data was previously processed from Anna et al. (Figures 3-j, 2003), to estimate experimental L and D values (see Figure 8).
The measurement is undertaken using Pixel Profile® ( 2007), an open source digitalization program.First, a line is drawn between two points (droplets) and the program estimates the distance as a number of pixels between the start and end points and convert this into continuous length with reference on the contraction width.Therefore, the distance between the first droplets is 2.09E-4 m (see Figure 8).
The speed of the emerging droplet is estimated as the relation of the total flow and the contraction cross-section area (Eqn 1), obtaining 0.82 m/s. (1) The time-lag between consecutive droplets is approximately 2.53E-4 s, and is calculated by Eqn 2, as: (2) With the help of Pixel Profile® the level of brightness-shade (named Hue) of the experiment's digital picture is extracted along a line joining two consecutives droplets, as shown in Figure 9.The average diameter obtained is 8.1E-5 m (Figures 10).

NUMERICAL RESULTS
In order to obtain the time-lag between consecutive droplets, several monitor points are used to track the transient response of the simulation.These monitoring points are located upstream and downstream of the contraction, as shown in Figure 11.Large print-shots are detailed in the plot of the volume fraction of water to show the shade contour around the given monitor point when the peak occurs.
Large images are detailed at the point where the volume fraction of water is higher in the peaks of the graph.In addition, it shows explicitly the extended area in which the image was taken and the color scale represents the change of volume fraction in the domain.
It is clear that even at Point 1, which is upstream of the contraction, shear layer instability sets and triggers the process of droplet formation downstream.
At Point 2, the plot (Figure 14) depicts a modulation of the instability, by reducing its amplitude and wavelength, as part of the focusing phenomenon.Between points 2 and 3, the shear layer break-up occurs, giving rise to the droplet formation, as shown in Fig. 15.
The section between points 2 and 3, allows capturing the elapsed time between separated droplets.
Section at point 4 shows the high diffusion experienced by the droplet as it travels downstream (see Figure 16).This is a known and expected behavior since the free jet maximum velocity should attenuate as the mixing layer broadens (see Figure 17).

ANALYSIS OF RESULTS
The shear layer between both fluids passing through the contraction sets the framework to start and perturb the instability that gives rise to the droplet formation.A more detailed view of the instability is shown in Figure 18.Furthermore, Figure 19 shows the formation of the water droplet.It is observed that the drop is a mixture of two fluids, while the force of cohesion of the water acts to create a concentration gradient from the center of the drop to the edge.
The droplets obtained from the simulation are not pure water but within them there is a significant concentration of oil up to 60%, as show in Figures 13-16  More defined droplet bounds are found by setting the droplet limits at water fraction of 10%.In such a case, droplet contours are depicted in Figure 21.
The average diameter obtained is 2.9E-5 m which differs about 64% from experiments.

CONCLUSIONS
• Amethodology was developed to simulate the droplet formation process in water-oil emulsions, through a pseudo-2D multiphase segregated numerical model.

•
The numerical model captured the droplet break up, which occurs within the contraction of the microchannel discharge, as a consequence of the flow focusing mechanism, apparently tuned-up by the two-fluid slip velocity.• Droplets form every 5.0E-4 s, corresponding to an increase of 49% over the experiments.

•
The numerically simulated droplets depict a water core with a large oil concentration gradient surrounding it, as a consequence of moderate numerical diffusion associated to the spatial discretization of the domain.

•
The computed droplets average diameter was 64% smaller than obtained in experiments.

RECOMMENDATIONS
The following recommendations are proposed for future projects related to this research: • Extend simulation scenarios to remaining experiments of Anna et al. (2003).

•
Perform a sensibility analysis on the surface tension to achieve a closer fit with the experiments of Anna et al. (2003), since surfactant was used in the experiments.

•
Extend numerical analysis to 3D geometry to simulate the phenomenon as close to reality as possible.Nevertheless, special care must be taken, since this condition will require more memory and simulation time.

Figure 9 .
Figure 9. Line to measure the diameter

Figure 13 .
Figure 13.Volume fraction of water as a function of time at point 1.

Figure 14 .
Figure 14.Volume fraction of water as a function of time at point 2.

Figure 15 .
Figure 15.Volume fraction of water as a function of time between points 2 and 3.

Figure 16 .
Figure 16.Volume fraction of water as a function of time at point 4.

Figure 17 .
Figure 17.Velocity profiles in the discharge channel

Figure 18 .Figure 19 .
Figure 18.Instability of the water layer in the contraction (I) . The average volume fraction at Point 1 is 0.43 due to it being closer to the internal input channel.It´s acknowledged that the strong diffusion occurs downstream of Point 2, by the growth of the multiphase shear layer, which makes the droplets reach monitor Point 4 with an average volume fraction of water of 0.255.The Figure20depicts a superposition of the volume fraction as a shade contour and as xy plot, which allows estimating the droplet diameter obtained numerically and see the presence of both fluids within the droplet bounds.