Drag-induced breakup mechanism for droplet generation in dripping within flow focusing microfluidics☆

Ping Wu ,Zhaofeng Luo ,Zhifeng Liu ,Zida Li,Chi Chen ,Lili Feng ,Liqun He ,*

1 Department of Thermal Science and Energy Engineering,University of Science and Technology of China,Hefei 230027,China

2 School of Life Sciences,University of Science and Technology of China,Hefei 230027,China

3 Department of Mechanical Engineering,Hong Kong University of Science and Technology,Hong Kong,China

Keywords:Microfluidics Droplet Dripping Flow focusing Thread Force balance

ABSTRACT Based on viscous drag-induced breakup mechanism,a simple model was proposed to predict the dripping droplet size as a function of controllable parameters in flow focusing micro devices.The size of thread before breakup was also investigated through laminar flow theory.Experiments and numerical simulations by VOF are carried out simultaneously to validate the theoretical analysis,showing that droplet size decreases rapidly with the increase of the flow rate ratio and capillary number.

1.Introduction

Micro droplets provide versatile ways for chemical and biological analysis such as multiple emulsion[1],polymerase chain reaction(PCR)[2-5],drug delivery[6-8],microreactors[9,10]and so on.Most of these applications demand highly uniform droplets in size,and flow-focusing microfluidic device is one of the most widely utilized method for such droplet formation.In these micro channels there are usually two immiscible fluids such as water and oil,and the inner liquid is squeezed to form a thread along the centerline by wrapping viscous oil.The thread head will become unstable when flowing through an orifice and releases monodisperse droplets by spontaneous generation.A lot of work have been carried out[11-14]since this approach was first implemented by Anna et al.[15],and experiments have shown that fluid flow undergoes several distinct flow regimes including jetting and dripping as the capillary number alters[16,17].In the jetting pattern,the water thread will be extended to a long distance out of the orifice,and the thread is subjected to the well-known Rayleigh-Plateau capillary instability that is dominatingly responsible for droplet formation.As for the dripping pattern,the thread will experience a sudden expansion near the exit of orifice,causing velocity redistribution and pressure drop in the flow field,and the detailed breakup mechanism appears more complicated and remains controversial to date.Several researches believed that the capillary instability also plays important role in droplet breakup in dripping,and they speculated that droplet radius should be proportional to the thread size[16,18].However,other investigations suggested that the droplet/bubble pinch-off period in dripping is weakly dependent on inter facial tension and thus it is difficult to solely associate the droplet/bubble rupture process with capillary instability as expected[11,19-24].In general,the dripping regime produces more homogeneous droplets than that of jetting regime,and no simple theoretical model was yet reported to explain the collapse mechanism and predict the droplet size distribution[25].

This paper attempts to present an intuitional force-balance model to elucidate the droplet breakup mechanism in dripping mode and forecast the droplet size as a function of flow conditions.Numerical investigation suggests that viscous drag contributes appreciably in droplet detachment at extremely low Reynolds and Weber numbers[26].From this standpoint,we describe and calculate holding and detaching forces acting on a growing droplet,and show that the dynamical competition between them determines the droplet deformation and breakup,as well as the final droplet size.Experiments and numerical simulations are carried out simultaneously to validate our theoretical analysis.

2.Model Development

2.1.Drag-induced droplet breakup mechanism

In our flow focusing microfluidic device shown in Fig.1(a),water and oil act as inner and outer phases.With the net action of oil on water,the liquid-liquid interface is squeezed strongly towards the center line and water forms a thread in the entrance of the orifice.In dripping regime,the thread expands abruptly at the very exit of the nozzle,and interfacial tension forces the thread head into spherical shape.However,the spherical droplet is easily confined into disk-like ones when the droplet size grows larger than channel depth,as shown in Fig.1(b).The sudden expansion leads to deceleration and stagnation of the spherical tip,and viscous drag from outer fluid acting analogously as those co- flowing capillary dropformation devices is a main detaching force.At low Reynolds number,viscous drag on the emerging droplet stretches and thins the neck,and becomes more important than capillary instability in the role of droplet breakup.

We now propose a relatively simple hydrodynamic model to explain the droplet detachment when only one droplet is in birth at the dripping thread head.As illustrated in Fig.1(a),Rdropis droplet radius,and the shape of water jet flow in orifice region is close to cylindrical with a diameter dj.In the growth of a droplet,forces acting on the growing droplet volume τ can be divided into holding and detaching ones.During the expansion,the interfacial tension Fσacting around thread surface is a resistance preventing the droplet expanding and detaching.In the presence of the neck connecting the droplet,the droplet shape of it is not entirely spherical,leading to the fact that the net effect interfacial tension acting on the growing droplet is opposite to the flow direction and can be perceived as a body force principally preventing the drop neck breakup.In contrast,the Stokes-like force FSconsisting of viscous and pressure stress performs as the most principal drag impelling the thread breakup.The linear momentum force FLcaused by inner fluid momentum is also considered as detaching forces accelerating the breakup.Generally,the inertial force is always negligible since it is about 10?4-10?7order of magnitude less than others in these low Reynolds number microfluidics(R e≈0.001-0.5).In the beginning,interfacial tension plays a dominate role and holds the droplet as a slender tube[27].With the growth of droplet,detaching forces become more important.When these detaching forces achieve balance with the holding one,the forming droplet begins to break away from the thread and tends to develop as spherical shape by the interfacial tension.This drag-induced breakup mechanism performs analogously to the breakup of a pendant drop under gravitational force.

2.2.Force balance model

Based on the above statement,one can write the force balance equation to determine the critical condition of droplet detachment as

Interfacial tension is generally described as

Here the cylindrical thread holds the drop,and σ is the interfacial tension coefficient.

Before detachment the droplet bears viscous drag increasingly from the outer fluid and it can be approximately calculated according to Stokes'Law at low Reynolds number.We must realize that the growing droplet is easily confined to disk-like ones due to a high width-to-depth ratio of the channel,and the fluid flow in microfluidics is always simplified as a planar flow.Nevertheless,the disk-like droplets would transit into spherical ones as the flow rate ratio(q=Qo/Qw)increases,resulting in those droplets whose size is even smaller than the channel depth.Our experiments confirm that spherical droplets are generated when q is higher than 30.So we describe the viscous drag according to a disk-like(2D)or spherical(3D)droplet,which is approximately determined by flow rate ratio.In addition,the internal circulation flow within micro droplets can be always neglected owing to the use of surfactant,and the Stokes flow assumption is acceptable in both cases.

In 2D case,the disk-like droplets with a depth h is considered as a cylinder and we obtain the viscous drag of per unit length on the grounds of the unbounded O seen flow around a two-dimensional cylinder as[28]:

With the consideration of microchannel wall effect,the net Stokes drag acting on a growing droplet can be described according to the Brenner's study[29]as

Here,Vo=2Qo/(wdownh),l=wdown/2,and β is a constant referring to the shape factors,β=1.004 is always appropriate for the case of parallel plate flow with a obstacle in the center line[29].

Fig.1.(a)Droplet formation in flow focusing micro geometry and forces acting on a growing droplet τ;(b)illustration of a disk-like droplet confined between the top and bottom wall.

In the 3D case,the classical Stokes drag acting on a spherical droplet is

Likewise,the net Stokes drag is calculated with wall effect considered[29],

Thus the viscous drag is determined by flow rate ratio as

The effect of the inner fluid momentum crossing the boundaries can be regarded as linear momentum force approximately described as[30]

Similar to co- flowing,we define a capillary number characterizing the relative importance of viscous stresses and capillary pressure across the interface between the two phases as

Ultimately,we obtain the prediction of final droplet size through Eqs.(1)-(7):

2.3.Thread size

Eqs.(2)and(10)indicate that djmust be calculated out first when estimating the interfacial tension force.Unlike the calculation of viscous drag as planar flow,the interfacial tension force is calculated according to three-dimensional(3D) flow case because the thread size is always smaller than the width and depth of orifice,especially in dripping regime(shown in Fig.4(a)).It is meaningful that the thread shape would always keep cylindrical in the rectangular cross-section of orifice region and 3D fluid flow is more reasonable in the calculation of the thread size.

It is observed that the radius of water thread Rjet(Rjet=dj/2)thins slightly just before achieving the force balance in the orifice region,and is determined mainly by flow rate ratio[16,18].The shape of inner fluid in the rectangular orifice region is close to cylindrical,and the assumption of fully developed inner flow is adopted.However,the forward movement of water-oil interface makes the velocity boundary condition little different from the conventional pipe flow.We prove in the 2D flow case(seen in Appendix)that the parabolic velocity profile of inner fluid should be modified by adding a constant,which equals the interface velocity between two liquid.The present velocity profile is coincident with the prediction of Utada et al.[16],and matches the viscous shear stress across the interface.Extrapolated to 3D case,the water velocity is written as

where the integral curve L:y2+z2=, uois the oil velocity distribution obtained from the traditional Poiseuille flow in rectangular pipe[31],

We reach the conclusion that the first term of the in finite series in the expression above was far more dominant than the others and the relative tolerance between Eq.(13)and Eq.(14)below is O(10?5).For n=0,then

where w and h denote respectively to the width and depth of the orifice channel shown in Fig.1(b).Since the cylindrical thread size is constant through the orifice,the pressure difference caused by interfacial tension between two phases is also nearly invariable,suggesting that the pressure gradient of two phases equals each other.After the simplification,a lot of work can be relieved in the integration for calculating the thread size with sufficient accuracy.Thereupon,Rjetis also calculated through the given inlet flow rate ratio as

In Section 5.2 we will give a reasonable simplified expression for the calculation of thread size through numerical approximation from Eq.(15).

Now we can summarize that the flow rate ratio determines the thread size,and the final droplet size can be estimated from the force balance breakup mechanism described as Eq.(1).The prediction is obtained entirely from the inlet parameter, fluid material properties and channel size but not scaling relation.

3.Numerical Simulation

3.1.The VOF method

To confirm the model above and acquire more detail of the breakup mechanism,numerical simulation is carried out.Volume of Fluid(VOF)method is one of the most powerful models to track the motion of interface in multiphase[32-36].

The simulation is implemented in FLUENT 6.3.26.Here we define oil and water as primary and second phases,and the variableαrefers to the volume fraction of second phase added to global equation.Volume fraction is calculated by

Generally,the source term Sαdoes not exist in our case.The advection of the volume equation relies on the flow field,and the revised momentum equation with surface tension as a body force is expressed as

Density and viscosity in each cell are calculated as the following averages:

Due to Brackbill's contribution[37,38],surface tension force can be added in momentum equation as a body force:

where the curvature of the interface κ can be defined as

Discretization of flow equation is by a second order upwind format with mesh element generated as 1 μm.To stabilize the iteration,we set a small fixed time-step as0.1μs initially.After hundreds of iterations,we increase the time-step slightly.The Courant number is defined in VOF model as

We set the global maximum Courant number to 0.5 to accelerate the convergence rate while maintaining the accuracy of calculation.

Near the interface,geometric reconstruction scheme is adopted to obtain more accurate interface shape.FLUENT provides standard geometric reconstruction scheme to implement the interpolation near the fluid- fluid interface using a piecewise-linear approach.The first step in this reconstruction scheme is calculating the position of the linear interface relative to the center of each partially- filled cell,based on the information about the volume fraction and its derivatives in both directions.The second step is calculating the advecting amount of fluid through each face using the computed linear interface representation and information about the normal and tangential velocity distribution on the face.The third step is calculating the volume fraction in each cell using the balance of fluxes calculated during the previous step.

3.2.The importance of meshes

The computational results,especially the accuracy of the interface shape and curvature,are directly affected by the grids.To find out the most appropriate size of the elements,four cases of 2D simulation with different mesh sizes(Δx=10 μm,4 μm,2 μm and 1 μm)are conducted as shown in Fig.2.It suggests that the coarsest grid produces laminar stratified flow of oil-water-oil rather than droplets.Though droplet generation can be observed in the case of Δx=4 μm,the interface shape obtained in the upstream of the orifice is chaotic and it is illogical with polydisperse droplets instead of monodisperse ones.When the element size is 1 μm the perfect results as well as the independence of grid are achieved.Therefore all the subsequent simulations are conducted on the grids with Δx=1 μm.

Results of the simulations make the process of focusing and breakup plainly visualized,and the size of thread and droplet can be calculated from volume fraction distribution of second phase.Flow parameters,material property and channel geometry is conveniently adjustable,with various droplet size acquired to validate the theoretical analysis.

4.Experimental Approach

Flow focusing microfluidic chips are fabricated via soft-lithography and replica molding[39]to demonstrate droplet formation.Polydimethylsiloxane(PDMS)and its curing agent(SYLGARD?184,Dow Coring,Michigan USA)is fully mixed with a mass ratio of 10:1,and then poured on a monolayer mold made by SU-8 Resists(MicroChem Corp.,USA)over a silicon wafer.The polymer is peeled off after baked for 1 h under 90°C,exposed to UV-light for 3 h,and then immediately bonded to a smooth PDMS face.All the depth of microchannel is h=37.6 μm and the width of inlet channel is win=184 μm.Two widths of the channel orifice(w=80,100 μm)with different outlets(wdown=200,400 μm)are used in the experiments.

Fig.2.Test of the independence of simulation on grids(ρo=980 kg·m?3,ρw=1000 kg·m?3,μo=0.048 Pa·s,μw=0.001 Pa·s,σ =0.023 N·m?1,Q w=2 μl·h-1,Q o=20 μl·h-1).

Fig.3.(a)Microfluidic chip in experiments.(b)Detailed size of flow-focusing microchannel in vertical view.(c)Image of the channel cross-section in scanning electron microscopy(SEM).

Fig.4.Dispersed phase distribution in the cross-section perpendicular to the direction of flow in the orifice region.(ρo=980 kg·m?3,ρw=1000 kg·m?3,μo=0.048 Pa·s,μw=0.001 Pa·s,σ=0.023 N·m?1,Q w=2μl·h-1).(b)Horizontal cross-section view of simulation results of the detailed droplet breakup process.(q=5,Q w=2μl·h?1,μo=0.048 Pa·s,μw=0.001 Pa·s,ρo=980 kg·m?3,ρw=1000 kg·m?3,σ=0.023 N·m?1).

Fig.5.(a)Theoretical and experimental thread size versus flow rate ratio.(ρo=980 kg·m?3,ρw=1000 kg·m?3,μo=0.048 Pa·s,μw=0.001 Pa·s,σ =0.023 N·m?1,Q w=10 μl·h-1).(b)Theoretical and experimental thread size versus viscosity ratio and orifice size.(ρo=980 kg·m?3,ρw=1000 kg·m?3,μw=0.001 Pa·s, flow rate ratio q=5,Q w=2 μl·h-1,w and h denote the width and depth of the orifice).

As illustrated in Fig.3,mineral or silicon oil is used as the continuous liquid and injected into the channel through lateral inlets by a syringe pump(New Era Pump Systems Inc.USA),while deionized water acting as the disperse phase,is transfused through the central inlet.Qoand Qware the inlet flow rate of oil and water,respectively.Dripping regime is observed through an optic microscope by a CCD camera throughout the experiments by controlling the flow rate ratio q ranging from 1 to 40.Different droplet sizes are reported and compared with those of theoretical and simulation results.

5.Results and Discussion

5.1.Simulated breakup process of droplet formation

In the VOF simulation,we firstly monitor the phase distribution in the cross-section of orifice to explore the shape of the thread shown in Fig.4(a).It is strongly validated that cylindrical-shape thread is acceptable when the flow rate ratio is higher than 1 in dripping regime.Furthermore,we show the detailed viscous drag-induced breakup process in numerical simulation as Fig.4(b),where we can see that the drop begins to develop near the orifice exit,and the thread size is almost invariable before the droplet detaches.The result is in good agreement with the assumption in theoretical calculation of thread size.

Fig.6.(a)Droplet diameter versus flow rate ratio.(ρo=980 kg·m?3,ρw=1000 kg·m?3,μo=0.048 Pa·s,μw=0.001 Pa·s,σ =0.023 N·m?1,Q w=2 μ l·h-1).(b)Theoretical droplet size versus capillary number.For silicon oil 1,q=1,μo=0.01 Pa·s,σ =0.0115 N·m?1,Ca=0.01;for mineral oil 2,q=5,μo=0.048 Pa·s,σ =0.023 N·m?1,Ca=0.02;for salad oil 3,q=5,μo=0.098 Pa·s,σ =0.0269 N·m?1,Ca=0.03;for silicon oil 4,q=10,μo=0.5 Pa·s,σ =0.0461 N·m?1,Ca=0.052.

5.2.Thread size

Fig.5(a)presents the results of prediction and experiments of thread size.High oil-water flow rate ratio means inner fluid would occupy less volume proportion in orifice resulting in more slender thread.It indicates as in the theoretical analysis that the larger the thread size,the stronger the holding interfacial tension force,and the harder the detachment resulting in larger droplets.Plugging droplets are generated when the oil flow rate is lower than that of water,and dripping occurs more obviously when q＞1.The careful measurement data validated the theoretical result from Eq.(15)about thread size.Fig.5(b)shows the effect of viscosity ratio and orifice size,and it is found that the thread size decreases rapidly with the increase of viscosity ratio and the narrowing of orifice.The viscosity difference between two phases causes the inner fluid to develop a different velocity profile from that of the outer fluid,and it exerts influence on the flow rate ratio determining the final thread size.

Since the flow rate ratio q,viscosity ratio λ and orifice size contribute appreciably on thread size according to the results of Figs.5(a)and(b),we can give a rational expression about thread size as

Numerical calculation gives the constant A=1.301,and Eq.(24)has sufficient precision for estimating the value of dj.

5.3.Prediction of droplet size

The force balance model above explained that droplets would pinch off when detaching effect surmount the holding force around the spherical interface.The prediction of droplet size versus flow ratio from Eq.(10)is also in good agreement with both numerical simulation and experiments as shown in Fig.6(a).As mentioned above,larger thread size leads to stronger interfacial tension and longer detachment time with bigger droplets achieved.Furthermore,one can observe that the dripping will perform tempestuously as flow ratio increases over 30 with tinier spherical droplets generated,yet the homogeneity of the size would not be guaranteed.We believe that drag-induced breakup mechanism would give way to capillary instability in higher flow rate,and jetting regime occurs when droplet size is so close to the thread.

Fig.6(b)plots the theoretical droplet size versus capillary number,and the curve tendency is coincident with those previous researches[24].Larger capillary number means the droplet would bear more viscous drag,resulting in shorter time to achieve the force balance for collapse and smaller size.In the corresponding experiments we use four different oils as continuous phase(silicon oil with lower and higher viscosity,mineral oil and salad oil)and adjust the flow rate to obtain different capillary numbers to explore its effect.Results of experimental images in the inset of Fig.6(b)confirm that the droplet size decreases with the increase of capillary number.

6.Conclusions

By virtue of viscous drag-induced breakup mechanism,a simple model is proposed to predict the droplet size as a function of controllable parameters.The droplet breakup is explained in light of hydrodynamics at low Reynolds number and the integral momentum equation.Limitations exists in the model such as it fails in the instability-dominate breakup jetting regime at the highest flow rate,but results of the final droplet size prediction is still in remarkable agreement with the experimental measurements and numerical simulations,suggesting that higher flow rate ratio and capillary number will generate finer emulsion droplets in dripping regimes.This theoretical analysis shows that droplet size in flow focusing microfluidics is predictable and controllable,which is meaningful for the control of droplet formation,operating condition selection and the microfluidic geometry excogitation.

Nomenclature

Ca capillary number(=μouo/σ)

Co Courant number(=vfluidΔt/Δx)

djthread diameter,μm

FLlinear momentum force,N

FSStokes-like force,N

Fσinterfacial tension,N

h channel depth,μm

l half width of the downstream channel,μm

Q flow rate,μl·h-1

q flow rate ratio of oil and water(=Qo/Qw)

Rdropdroplet radius,μm

Re Reynolds number(=2ρoVoRdrop/μo)

Rjetthread radius,μm

Δt time step size in VOF simulation,s

uooil velocity distribution in orifice region,m·s-1

Voaverage velocity of oil in orifice region,m·s-1

V average slip velocity of the oil-water interface,m·s-1

vfluidmaximum fluid velocity in VOF simulation,m·s-1

w width of orifice,μm

wdownwidth of the downstream channel,μm

Δx grid size in VOF simulation,μm

α volume fraction of water phase

β constant referring to the shape factors in the modification of viscous drag

λ viscosity ratio of oil and water(=μo/μw)

μ viscosity,Pa·s

ρ density,kg·m?3

σ interfacial tension coefficient,N·m-1

Subscripts

o oil

w water

Appendix

In the simplest 2D case,the microchannel is taken as a planar structure when the dimension in the thickness direction is far less than the other two,and the typical parallel flow holds.In the following,the velocity distribution of oil and water is derived,and the fact that the penetration of water will lead to the increase of oil pressure gradient is validated.

Since the Reynolds number is extremely low in this microfluidics and the stratified flow is nearly steady in the orifice,the Navier-Stokes equation can be simplified as the Stokes equation.As illustrated in Fig.A1,both the velocity profile is parabolic and only about y due to the assumptions of steady,laminar and fully developed flow.Non-slip boundary condition is considered near the channel wall,and zero shear stress must satisfied within water phase at y=0.Across the water-oil interface(y=±Rjet),the velocity u and shear stress τ of two liquid must be continuous,so the equation is written as

Fig.A1.Parabolic velocity profile of oil and water flow in orifice region.

from above we obtain the velocity profile as

Two conclusions can be drawn from the above and extended to 3D flow case.First,we note that the expression of oil velocity distribution is the same in form as the traditional inner flow of single fluid between parallel plates.It suggests that the central penetration of water does not influence the form of oil velocity profile.However,it will cause an increase of oil pressure gradient to reach mass conservation of oil when the inlet flow rate is unaltered.Second,the velocity of water phase consists of the traditional parallel flow velocity and the slip velocity of water-oil interface.