Liquid flow-focused by a gas: jetting, dripping and recirculation
Miguel A. Herrada, Alfonso M. Ganan-Calvo, Antonio Ojeda-Monge, Benjamin Bluth, Pascual Riesco-Chueca
LLiquid flow-focused by a gas: jetting, dripping and recirculation
Miguel A. Herrada, Alfonso M. Ga˜n´an-Calvo, AntonioOjeda-Monge, Benjamin Bluth, Pascual Riesco-Chueca
E.S.I, Universidad de Sevilla.Camino de los Descubrimientos s/n 41092 Spain. (Dated: January 5, 2014)
Abstract
The liquid cone-jet mode can be produced upon stimulation by a co-flowing gas sheath. Mostapplications deal with the jet breakup, leading to either of two droplet generation regimes: jet-ting and dripping. The cone-jet flow pattern is explored by direct axisymmetric VOF numericalsimulation; its evolution is studied as the liquid flow-rate is increased around the jetting-drippingtransition. As observed in other focused flows such as electrospraying cones upon steady threademission, the flow displays a strong recirculating pattern within the conical meniscus; it is shownto play a role on the stability of the system, being a precursor to the onset of dripping. Closeto the minimum liquid flow rate for steady jetting, the recirculation cell penetrates into the feedtube. Both the jet diameter and the size of the cell are accurately estimated by a simple theoreticalmodel. In addition, the transition from jetting to dripping is numerically analyzed in detail in someillustrative cases, and compared, to good agreement, with a set of experiments. a r X i v : . [ phy s i c s . f l u - dyn ] M a y . INTRODUCTION The controllable production of small flowing geometries is a crucial challenge for chemicalengineering [1] and bio-industry [2, 3]. Drops, bubbles, jets and recirculation cells in themicro-scale provide a useful platform for diverse technical applications. Here we concentrateon cone-jet flow patterns, looking at the streamline geometry in the jetting to drippingtransition. Recirculating flow is shown to take place just before the transition. This mayhave been disregarded in some experimental setups, whose prime concern was the studyof the cone and jet geometry, or the analysis of drop generation. Its study requires eitherspecific flow visualization techniques or numerical simulation methods. The conditions forrecirculation are extremely interesting, both as an indicator phenomenon associated to thejetting-dripping threshold, and as an attractive technological feature.Small droplet generation by means of co-flowing immiscible fluid streams has becomewidespread. The intrinsic smallness of the output droplets generally leads to small Reynoldsnumber flows. Thus, a number of classic studies back to Taylor [4] including the recentlyblooming field of co-flowing microfluidics take the low Reynolds number assumption forgranted. For example, a simple scheme (a straight tube surrounds a coaxial, more slen-der tube, two immiscible fluids being fed through each tube) has been thoroughly exploredby Suryo and Basaran [5] using a computer simulation: a locally extensional flow sponta-neously develops at the tip of the elongated drop drawn by the co-flowing liquid, causingthe seemingly continuous ejection of extremely small droplets by tip streaming under cer-tain parametrical combinations. This is obtained without the burden of complex geometry.It is hoped, therefore, that low Reynolds number co-flowing small-droplet generation maybecome a hydrodynamic standard in a near future. Different setups have been proposedwhere the flow is driven by an external straining flow. Among them, the elegant analyticalsolution by Zhang [6] points to parametrical combinations where extremely thin fluid jets,even down to the molecular scale, could be continuously reached. Those jets, if confirmed,would yield unimaginably small droplets upon breakup.An interesting research field is concerned with the behavior of electrified cones and dropsof leaky dielectric fields. The problem of freely-suspended liquid droplets deforming dueto an applied electrostatic field was examined by Haywood, Renksizbulut and Raithby [7].Collins et al. [8] reported simulations and experiments supplying a comprehensive picture2f the mechanisms of cone formation, jet emission and break-up that occur during EHD tipstreaming from a liquid film of finite conductivity. Lac and Homsy used a boundary integralmethod to describe the axisymmetric deformation and stability of a viscous drop in a steadyelectric field [9]. In the present paper, however, only non-electrified fluids will be considered.In spite of their generality and tractability, low Reynolds number flows are constrainedby the requirement that the overall flow velocity does not grow above a certain threshold toensure that inertial forces remain negligible. This constraint limits the overall productivityof low-Re systems. Co-flowing with inertia was successfully explored, aiming at a reductionof the issued bubble diameter, by Oguz and Prosperetti[10]. Subsequently, new perspec-tives were open by the emergence of moderate-high Reynolds number flow-focusing [11] asa high-productivity alternative to low Reynolds number co-flowing systems. Compared toother co-flowing techniques, flow focusing (FF) stands today as a mature microfluidic stan-dard yielding steady capillary jets or droplets whose size is well below the scale of the flowboundaries. As originally conceived, FF aimed at the generation [11] of continuous steadymicro-jets upon focusing by a co-flowing gas stream at moderate-high Reynolds number.Furthermore, FF was shown [12] to produce perfectly monodisperse microbubble streamswhen the co-flowing current is a liquid. A slight variation of the concept was subsequentlyintroduced by Takeuchi et al. [13] to produce microbubbles. When the axisymmetric ge-ometry originally proposed was reduced to a planar topology [14], particularly suitable formicrofluidics, the scientific literature production on flow focusing underwent an enormousboost[15, 16, 17]. In addition, axisymmetric multiple-phase FF leading to compound co-axialmicrojets[11] has been developed by other authors to produce microcapsules in a microfluidicsetup at relatively low [18] and moderate Reynolds number[19].The technological applications of FF were evident from its very inception. A crucial ad-vance resulted from the combination of massive production (high production rates of micro-scopic fluidic entities) and accurate tailoring. Depending on the geometry and arrangementof the involved fluids (a decision determining which interfaces are to be created), the natureof the fluids involved, gas or liquid, and the system geometry, an output including nearlymonodisperse micro-droplets, bubbles or complex capsules can be obtained at an unprece-dentedly controllable rate. Surface tension becomes a paramount ally in the conformationof discrete (generally spherical) fluidic units. Capillary jets have long ago been observedto give rise to continuous drop streams at fast emission rate upon Rayleigh axisymmetric3reakup. Here, although surface tension is negligible compared to other driving forces in theglobal scale, it becomes the main driving agent for jet instability and breakup. Obtaininga jet is therefore the precondition for the creation of a fluid domain with higher velocitiesand smaller dimensions at no cost in terms of control; and surface tension is free to performits conformation task in this new scale. Thus, as first proposed in FF, the steady capillarythin jet conformed by pressure forcing by an immiscible co-flowing fluid provides favorablelocal conditions, a suitable environment for the generation of bubbles, capsules or droplets.A FF capillary jet is driven by three main agents: fluid inertia, viscosity and surfacetension. Owing to the simplicity of the slender jet geometry, which asymptotically ren-ders all forces strictly additive in one dimension, FF can be scaled with the help of twodimensionless parameters: (i) the inertia to surface tension forces ratio (Weber number)and (ii) the viscous to surface tension ratio (Capillary number). Other classic numberssuch as Reynolds are combinations of the former. Nevertheless, as early noticed [11], anintrinsic feature of FF, namely the presence of a focusing fluid, gives rise to supplementaryinfluences issuing from the correlation between the properties of the focusing and focusedfluids. In particular, when a liquid is being focused by a gas, the gas sheath flows muchfaster than the liquid jet at the exit orifice (Fig. 1). Thus, in addition to the extensionalviscous forces at the neck of the meniscus, transversal viscous diffusion of momentum causesa non-trivial axial velocity profile. Some simplifying assumptions have been adopted, yield-ing accurate first order solutions [11, 20, 21]. However, they are not applicable to predictcritical phenomena like the onset of steady jetting, or the jetting-dripping transition, as afunction of the working parameters. These problems have been made analytically tractableat the expense of a drastic geometry simplification, i.e. assuming infinite jet slenderness[22, 23, 24]. Such simplified models are predictive in a variety of situations, but FF systemsexhibit an intrinsically three-dimensional meniscus from which the jet or the small dropletsissue. Simultaneous modelling of the meniscus and the jet goes beyond the scope of presenttheoretical frameworks. Thus, numerical simulation or experiment are the only avenue todiscern the physics of the fluid emission and its parametric conditions. Some further insightcan be gained by general scaling laws. This is the approach chosen for the research presentedhere.Many authors have applied numerical simulation [25] to this class of problems, whereit has supplied welcome information on the droplet dynamics of complex flows [26]. A4ignificant number of studies have been proposed on microfluidic FF devices; occasionalcomparison with experimental data is provided to validate the numerical model. A liquid-liquid configuration for the production of microemulsions has been simulated [27] to goodagreement with experiments [14]. Other authors have considered the microbubbling setup[17, 28], where good experimental fit is also obtained[15].In this work, we make use of numerical simulation, with some experimental support,to study the generation of a liquid jet focused by gas in an axisymmetric FF device, atmoderate-high Reynolds numbers. The jet diameters obtained in the simulation are in goodagreement with our experiments and scaling laws [11], a fact that fully validates our hypothe-ses. Among other findings, we determine the flow rate at the jetting-dripping transition fortwo combinations of Reynolds and Weber numbers. We observe that either the jet or thecusp-like meniscus are responsible for the global instability of the system, which drives it intowell defined dynamical cycles (global dripping). A detailed description of the flow patternsheds light on the physics of the jetting-dripping transition and the peculiar appearance ofthese two regimes in co-flow problems, as opposed to faucet jetting and dripping.One of the key findings of the simulation is the occurrence, under favorable drivingconditions, of a recirculation cell in the meniscus. This is in perfect analogy to recircu-lating meridian fluid flows observed inside Taylor cones when electrospraying liquids withsufficiently large values of both the viscosity and the electrical conductivity [29, 30]. Addi-tionally, in experiments aiming at the production of tip streaming patterns in liquid-liquidtwo-dimensional FF (surfactant treated interface), the streamline image of fluorescent par-ticles seeding the flow of the internal, aqueous liquid during thread formation, was shownto consist of symmetric recirculation vortices [31]. In both types of motion, either driven bythe electrical stresses acting at the cone surface or by the external focusing flow, the liquidflows towards the meniscus (cone) tip, along the generatrix, and away from it along the axis.The problem under consideration here is comparable to these other instances of recirculatingcell, because the driving action of the gas sheath, which causes a strong tangential forcingat the interface, plays a similar role to either tangential electric stress[29, 30], surfactant-aided liquid-liquid interaction[31], or purely surfactant-driven tip-streaming.[32, 33] In thiswork, scaling arguments are developed to describe the size and occurrence of purely FFrecirculation cells. 5
IG. 1: Simulated boundary geometry and fluid flow domains.
II. GOVERNING EQUATIONS AND BOUNDARY CONDITIONS
The axisymmetric flow-focusing device and the computational domain used are sketchedin Fig. 1. A constant liquid flow rate Q l , flowing through a capillary tube (outer diameter D = 2 R , inner diameter D = 2 R ), is forced through a coaxial round orifice of diameter D = 2 R (nozzle) located at an downstream distance H from the tube outlet. The liquidstream is drawn by a constant flow rate Q g of focusing gas stream discharging throughthe nozzle into a infinitely large chamber. The gas flow is assumed incompressible, inasymptotic consistency with the low pressure drop at the exit orifice, a condition prevailingwhen maximum droplet size monodispersity is required. Therefore, the incompressible,axisymmetric and unsteady Navier-Stokes equations in cylindrical ( z, r, φ ) coordinates areused to describe the time evolution of both fluids.Fig. 1 also shows the boundary conditions: a) at the liquid inlet, z = − z l , a Hagen-6oiseuille profile, U l ( r ) = V [1 − ( r/R ) ] , is specified; b) at the gas inlet, z = 0 , R < r < R ,a uniform axial flow, U g ( r ) = V , is assumed. This assumption is based on the following:the gas Reynolds number under the conditions explored is relatively high (of the order of,or above 100), while most flow-focusing devices have a gas inlet length L g which is notmuch bigger than the width (cid:52) R = R − R , so that the relevant dimensionless number forboundary layer development, ρ g U g ( (cid:52) R ) / ( µ g L g ) is sufficiently above unity; c) on all solidwalls we assume no slip and no penetration u = ; d) at the axis r = 0 a symmetry conditionis applied; e) the outlet discharge chamber has been modelled as a rectangular box, z = z out and r = r out being two open surfaces where the pressure is set to zero. This assumption isdiscussed later on.Note that the corresponding gas and liquid flow rates can be derived from the inletvelocity field: Q l = (cid:90) R πrU l ( r ) dr ; Q g = (cid:90) R R πrU g ( r ) dr. (1)Parametric studies of the dimensionless variables involved are carried out next. Thevelocity field u = ( u, v ) is scaled with the mean gas velocity at the nozzle V = Q g / ( πR ),while length is scaled with the nozzle radius R , time t with R/V , and pressure p with ρ g V , ρ g being the density of the focusing gas. A single geometrical configuration is considered in thiswork, characterized by the following aspect ratios: R /R = 0 . R /R = 1 . R /R = 3 . H/R = 1 and
L/R = 0 .
75. We have chosen a liquid-gas combination where ρ l ρ g = 833 . , µ l µ g = 55 .
55 (2) ρ and µ being the density and viscosity of the liquid (subindex l ) and the gas ( g ). Thischoice is representative of the experimental jetting of air-focused water. The problem isgoverned by the following dimensionless parameters: Reynolds and Weber numbers Re = ρ g V Rµ g , (3) W e = ρ g V Rσ ; (4) σ being the surface tension between the two phases. Q is defined as the flow rate quotient: Q = Q l /Q g . (5)7or a given value of Re and W e we wish to analyze the formation of a steady liquid jet andthe dependence of the flow on the quotient Q . In particular, we identify the minimum valueof Q , Q ∗ ( Re, W e ), below which the liquid jet ceases to be steady and a dripping regime isobserved in the simulation. The regime is considered to be steady (and the jet convectivelyunstable) if the liquid meniscus remains steady for a sufficiently large period of time.We should point out that in order to focus a jet of liquid by gas, moderate-high Reynoldsnumbers are needed. We consider in detail two different conditions for the focusing gas:Case 1: Re = 465 . W e = 8 . Re = 931 . W e = 32 . III. NUMERICAL PROCEDURE
In order to predict the interface geometry during the time-solution, several techniqueshave been used, falling into one of three categories. These are: (i) interface tracking meth-ods, including a moving mesh, [34] (ii) front tracking and particle tracking schemes, [35]and (iii) interface capturing methods, including volume of fluid (VoF)[36, 37] and level settechniques.[38] We chose a VoF method consisting of two parts: an interface reconstructionalgorithm to approximate the interface from the set of volume fractions and a VoF transportalgorithm to determine the volume fraction at the new time level from the velocity field andthe reconstructed front. The basic method is robust and flexible, and is based on widelyused VoF schemes [40, 41, 42, 43].For convenience and with the aim of making our results readily reproducible for others,we have used the well tested commercial solver FLUENT v 6.3 (laminar unsteady) to resolvethe discretized mass continuity, momentum conservation, and the liquid volumen fractionequations in the mesh depicted in Fig. 2, generated by commercial code GAMBIT inFLUENT v 6.2. Observe that the smallest quadrilateral elements lie between the needleedge and the nozzle, where the liquid meniscus is located, and in the near-axis region, wherewe expect the development of the liquid jet. The basic mesh should be sufficiently refinedto capture, in the absence of the liquid, the strong velocity gradients experienced by thegas flow at the orifice region. In the grid shown in Fig. 2 the minimum cell radial andaxial lengths are (∆ z ) min = (∆ r ) min = 0 .
02. Several numerical tests with smaller size mesh8
IG. 2: Grid of the domain under study. A denser mesh is provided in areas where the interfaceis expected to lie. To avoid numerical diffusion of the interface, the interface region is defined witha higher density of nodes. cell have shown that this level of accuracy is comfortably sufficient to describe the gas flowpattern for the two cases considered ( Re = 465 .
83 and 931 . z ) min = (∆ r ) min = 0 .
01. Finally, only for the more difficult cases (case 2 with Q small),the results where recomputed in a finer grid with (∆ z ) min = (∆ r ) min = 0 . z out = 10 and r out = 3 .
5, to minimizeartificial boundary effects in the results obtained. On the other hand, the z position wherethe inlet boundary conditions for the liquid are imposed, has been located sufficiently faraway from the needle edge, at z = z l = −
3. This choice has been made, as will be shownlater, because a liquid recirculation cell intrudes upon the capillary tube when Q decreases.Therefore, in order to impose a Hagen-Poiseuille profile for the liquid velocity as a wellposed inlet boundary condition, this boundary should be set sufficiently far upstream fromthe recirculating region.Tracking the interface between the phases is accomplished by solving a continuity equationfor the volume fraction of one of the phases using an explicit time-marching scheme. The restof the equations are solved implicitly. The time steps selected were fixed and sufficientlysmall to ensure that the global Courant number based on the mesh cell size, the meanvelocity in the cell and the time step was always smaller than one. Regarding the spatialdiscretization of the equations, the third-order modified MUSCL scheme[44] is used to obtainthe face fluxes whenever a cell is completely immersed in a single phase. When the cell is nearthe interface, the CICSAM algorithm is used[45]. The pressure corrections are computedwith the body forces weighted scheme and the pressure-velocity coupling in segregated solveris treated with the PISO method [46]. All under-relaxation factors are set to one to avoidany numerical masking of fade-out effects in our physical problem. IV. NUMERICAL RESULTS. DISCUSSION
A fruitful interpretation of the results obtained needs to be situated in the frame of theliterature on the dripping faucet. The book by Shaw [47] gave rise to a rich and insightful se-ries of studies, among them major contributions by Fuchikami et al. [48], Ambravaneswaranand co-workers [49, 50, 51], and Coullet, Mahadevan and Riera.[52] To discriminate betweenthe jetting and dripping modes, it is helpful to make use of the categories introduced byAmbravaneswaran and co-workers: [50]1. the dimensionless limiting length L d /R from the capillary edge to the extremity ofthe first drop at detachment. 10. the ratio of the distance L s between the centers of mass of the drop that is about toform and the previously formed drop, and L d
3. the ratio of the volume of the drop that is about to form, V d , to that of the drop thatis attached from the capillary, V p When undergoing the transition from dripping to jetting, the first parameter undergoesa sudden increase, while the two other ones experience an opposite trend. In general, thedripping mode is characterized by bulky drops, relatively distant from each other, whosediameter is considerably larger than the jet diameter from which they detach.The usual categories applicable to faucet dripping need some adaptation before beingused in a co-flow problem, where there is considerable stretching of the cone-jet and droplettrain by the coaxial current: the drops are deformed, their radial extent is limited, andthey may undergo secondary breakup (particularly so in the dripping regime, whose bulkydrops are more vulnerable to shear) after detachment from the filament. Therefore, under co-flowing, the classical aspect of the jetting and dripping regimes is modified, and the transitionbetween them is not sharp: such features are confirmed by experiment, as explained below.This is the reason why the behavior of the meniscus can be used as a further indicator ofthe jetting mode. Dripping leads to a pulsating meniscus, each detached drop giving riseto recoil and oscillating; while, in jetting, the detachment of the drops does not cause anyfluctuation of the meniscus and jet (see Fig. 19). In a full dripping regime, these pulses areperfectly regular (see Fig. 14); in most cases a slender unsteady liquid ligament detachesfrom the meniscus and breaks up into droplets of heterogeneous size [53]. However, atthe onset of dripping (a situation which will be labelled “incipient dripping”), completelyirregular fluctuations of the meniscus are observed [50].Accordingly we begin by studying the formation of a steady (convectively unstable) liquidjet in the FF device. Initially, the capillary needle is filled with liquid up to z = 0 while therest of the domain is filled with gas. We start the simulation from rest ( u = v = 0) in thewhole domain except at the inlet sections, where velocity profiles are prescribed. Fig. 3(a)-3(h) shows the formation of a steady liquid jet for case 1 and Q = 0 . α = 0 . θ ( t ), has reached a constant value in time, θ ( t ) = θ o .2. Both the jet diameter at the nozzle inlet, d in ( t ), and at the nozzle exit, d out ( t ), shouldreach a steady regime or a stable oscillating regime around a mean value; these quan-tities are of course to stay above zero. This amounts to excluding jet breakup inthe nozzle region, a feature associated with a non-slender jet and possible drippingbehavior (unsteady meniscus-jet).Here, d in ( t ) and d out ( t ) are computed at each time step, by integrating radially the liquidvolume fraction, α , at the nozzle inlet, z = 2, and at the exit, z = 2 . d in ( t ) = 2 (cid:115) (cid:90) α ( t, z = 2 , r ) rdr ; d out ( t ) = 2 (cid:115) (cid:90) α ( t, z = 2 . , r ) rdr. (6)For sufficiently large Q , as illustrated in Fig. 3, both d in and d out evolve towards a steadyvalue.However, oscillations of these two quantities are observed when Q is reduced. For exam-ple, Fig. 4 shows the time oscillation of d in and d out for case 1 and Q = 0 . θ , and the flow structure inside the meniscus as afunction of Q . Since the flow is unsteady, we will use a mean value of the jet diameter atthe nozzle inlet and outlet as defined by:¯ d in = 1 T (cid:90) T + t i t i d in ( t ) dt, ¯ d out = 1 T (cid:90) T + t i t i d out ( t ) , dt (7)where t i is a time position once a steady jet has developed and T is a time period longenough to ensure a significative mean value. For example, selecting t i = 0 and T = 500leads to ¯ d in = 0 . d out = 0 . r r r (b)(a)(e) (d)(c) (h)(g) (f) t=147 t=168t=231t=252t=210t=357 t=315t=1000 FIG. 3: A sequence of snapshots from the simulated growth of an eventually steady jet (case 1, Q = 0 . The procedure is the same for the two cases under consideration. The simulation isstarted from rest with a value of Q sufficiently high to obtain a steady jet. Then, Q isreduced and the solution is monitored in time until a new steady jet is obtained. Fig. 5shows the stabilized liquid-gas interface for case 1 and different Q . It should be pointed outthat Q = 0 . Q ∗ for steady jetting: Q ∗ = 0 . Q values once a steady regime isreached. The smallest jetting flow rate here is Q ∗ = 0 . d in and ¯ d in are plotted as a function of Q for (a) case 1 and (b) case 2. To complete the picture, Fig. 8 shows the dependence of13
100 200 300 400 5000.2450.24550.2460.2465 d i n t d ou t FIG. 4: Time dependence of the jet diameter, d in (dashed line, left ordinates) and d out (continuousline, right ordinates) for case 1 and Q = 0 . d out as compared to d in . the meniscus angle θ with Q for the two cases. In both examples, θ becomes smaller as Q decreases (smaller flow rate quotient implies stronger focusing action). Just before dripping,as the liquid flow rate is reduced, the angle appears to become independent from Q : theinterface geometry becomes invariant (local hydrostatic balance). The smallest angles areobtained in case 2. This is to be expected since the normal pressure forces produced by thegas stream, which cause the focusing flow, are larger for that case.Analyzing in more detail the structure of the flow inside the liquid meniscus in jettingmode, in the lower- Q range, a meniscus recirculation cell is observed, in analogy with otherco-flowing systems [5, 31] and Taylor cones [29]. Fig. 9 shows instantaneous streamlinesfor case 1 and different values of Q . The recirculation increases when Q decreases, the cellpenetrating into the capillary needle. Fig. 10 depicts instantaneous streamlines for case 2and four different values of Q . Again, a recirculation region appears before the meniscus jetsystem ceases to be steady. The size of the recirculation can be calculated by finding thetwo z -positions where the velocity at the axis becomes zero.14 −1.5−1−0.500.511.5 −1.5−1−0.500.511.5z r (a) (b)(c) (d) Q =0.004 Q = 0.0008Q = 0.0006 Q = 0.0004
FIG. 5: Shape of the liquid-gas interface as a function of Q (case 1) Fig. 11 shows the velocity at the axis, v axis , as a function of z for five different valuesof Q : (a) case 1 and (b) case 2. It is worth observing that v axis is roughly uniform insidethe capillary needle, its value being given by the Hagen-Poiseuille expression; in the nozzleregion it increases owing to the focusing effect of the gas stream, which creates the issuingjet. A region where u axis decreases is located in the meniscus region, between the capillaryand the nozzle; note that when Q decreases, a local minimum of u axis is observed at a givenposition z = z min in the meniscus region. If Q is sufficiently small, v axis becomes negativenear the local minimum in a region delimited by the two z positions, z and z , where v axis = 0. Therefore, the size of the recirculation region, S r , observed in Figs. 10 and 11,can be computed as S r /R = s r = z − z . There is a threshold value of Q , Q r , below whicha recirculation pattern is observed. At the threshold flow rate v axis = 0 at z = z min = z and v axis > s r as function of Q for (a) case 1 and (b) case 2. Looking back at Fig. 9,note that the size of the recirculation cell increases as Q decreases. In situations of incipient15 r (a) (b)(c) (d) Q = 0.0012 Q = 0.0004Q = 0.0002 Q = 0.0001
FIG. 6: Shape of the liquid-gas interface as a function of Q (case 2) −4 −3 −2 d −4 −3 −2 (a) (b) FIG. 7: Mean jet diameter at the inlet (dashed line) and at the exit (solid line) of the nozzle orificeversus the flow rate quotient Q . (a) Case 1, (b) Case 2. −4 −3 −2 θ −5 −4 −3 −2 (a) (b) FIG. 8: Angle θ (degrees) versus flow rate quotient Q . (a) Case 1, (b) Case 2.FIG. 9: Instantaneous streamlines at four different flow rates quotients Q for case 1: (a) Q = 0 . Q = 0 . Q = 0 . Q = 0 . recirculation ( Q smaller than but similar to Q r ) this growth appears to be linear, as derivedlater from dimensional arguments. In Fig. 12, the discrete points ’o’ have been obtaineddirectly from the simulations. The dashed lines are linear interpolations computed in therecirculating regime, s r >
0. The linear interpolation is not only in good agreement with17
IG. 10: Instantaneous streamlines with four different flow rate quotients Q for case 2: (a) Q =0 . Q = 0 . Q = 0 . Q = 0 . −1 0 1 2−0.0200.020.04 u a x i s z −1 0 1 2−0.0100.010.020.03 zQ = 0.004Q = 0.0016Q = 0.0008Q = 0.0006Q = 0.0004 Q = 0.0012Q = 0.0008Q = 0.0004Q = 0.0002Q = 0.0001 (a) (b) FIG. 11: Velocity at the axis versus z for several flow rate quotients Q : (a) Case 1 and (b) Case 2. s r
0 0.0005 0.001 0.00150123 Q (a) (b)
FIG. 12: Size of the recirculation s r as a function of Q : (a) Case 1 and (b) Case 2. the data but also provides a reliable approximation to compute Q r . The estimations are Q r = 0 . Q r = 0 . Q is: s r ∼ A ( Q r − Q ) . (8)This means that s r is proportional to a back flow rate Q B /Q g = Q b = ( Q r − Q ) for agiven set of fluid properties, geometry, and gas flow Reynolds number. The relative locationof the jetting threshold Q ∗ and the recirculation threshold Q r is, in both cases, Q ∗ < Q r .Therefore, recirculation can be taken as a dripping-precursor: as far as can be gathered fromour simulation, it always precedes global instability of the meniscus-jet system.Finally, some results are presented with flow rate values below the jetting threshold. For Q < Q ∗ , with Q close to Q ∗ , our simulations show the flow to exhibit different behaviors ina sequence: a period where a thin jet breaks up in the nozzle region alternates with otherperiods where a thin jet breaks up downstream of the nozzle. The irregular time behaviorof the flow for Q < Q ∗ , but Q close to Q ∗ ( incipient dripping ) can be observed in Fig. 13,where d i (a) and d out (b) are shown as a function of time for case 1 and Q = 0 . Q < Q ∗ but Q sufficiently different from Q ∗ , the flow behavior becomes more regular andperiodic with a unique dripping frequency.As anticipated at the first part of this section, our results, while belonging to the dripping-regime, are untypical in that the co-flowing current gives rise to axial stretching of the19
500 100000.050.10.150.20.250.3 t d i n d ou t (a) (b) FIG. 13: Time evolution of (a) d in and (b) d out in a irregular dripping regime for case 1 and Q = 0 . r r r r (b)(a)(e) (d)(c) (h)(g) (f) t=0 t=49t=105t=119t=77t= 161 t= 133t= 189 FIG. 14: A sequence of the dripping regime for case 1 and Q = 0 .
200 400 600 80000.10.20.30.4 t d i n d ou t (a) (b) FIG. 15: Time evolution of (a) d in and (b) d in in a dripping regime for case 1 and Q = 0 . jet and drops, so that unusual breakup geometries result. The radial extension of thedrops is limited, and the axial stretching gives rise to secondary breakup, immediately afterdetachment. The pattern observed in Fig. 14 points to a dripping regime: it is perfectlyperiodic (transient jetting can therefore be excluded) and each period is associated with thefilling up of a drop, its breakup from a thinning filament, and the recoil of this filament. Fig.14 shows a complete time sequence of a dripping process (case 1, Q = 0 . T ∼
210 for each cycle. Fig. 15 plots d i (a) and d out (b) as a function of time for this case. Initially, a liquid meniscus is growingwith no jet production, and d in = d out = 0. Subsequently, a liquid jet issues and d in and d out become positive. Both diameters reach a maximum at a certain time and then start todecrease. Finally, the jet breaks into droplets, d in and d out are set to zero and the processbegins anew. In spite of the observed differences, the dripping process in this case is quitesimilar to regular faucet dripping, the time period being mainly imposed by the filling ofthe meniscus until reaching a critical volume.21 y y y (b)(a)(e) (d)(c) (h)(g) (f) t=0 t=139.7t=244.6t=261.3t=209.6t=314.4 t=279.5t=419.2 FIG. 16: A sequence of the dripping regime for case 2 and Q = 0 . A similar situation is observed in case 2. Fig. 16 shows a complete time sequence of aregular dripping process with Q = 0 . T ∼
500 for each cycle (see d i (a) and d out (b) inFig. 17 as a function of time). A. Influence of the BCs and the spatial and temporal resolution on the numericalresults
The numerical problem addressed is quite complex: it involves a high speed stream ofgas discharging through a nozzle into a infinity large chamber plus a meniscus-liquid jetwhich may break into droplets within the finite numerical domain. This complexity leads todifferent time and spacial scales associated to a plurality of interacting physical phenomena(jet breakup due to capillary and Kelvin-Helmholtz instabilities, mixing layer instabilities22
500 1000 150000.050.10.150.2 t d i n d ou t (a) (b) FIG. 17: Time evolution of (a) d in and (b) d in in a dripping regime for case 2 and Q = 0 . in the main gas stream). Therefore, an accurate analysis of the jet breakup is difficult toachieve. On account of it, though our VOF method is fully reliable in qualitative terms asa predictor of jet breakup and drop formation, we have devoted this subsection to checkthat our numerical results (meniscus-jet shape as a function of Q for different setups) areindependent of the selected BCs and numerical meshes.As indicated above, the most problematic simulation choice is setting p = 0 at theoutlet boundaries, since any jet or a drop crossing the boundary is influenced by the strongand artificial restriction that the pressure remains fixed. Our choice is a simplification( p = 0) which takes advantage of the essentially parabolic character of the equations. Asecond option has been explored, the so-called outflow conditions (assuming uniformity, i.e.Neumann type), but they give rise to a false constraint on the flow pattern, because theyimply that the gas flow is coaxial. A minisymposium held in 1994 on the open boundarycondition problem in incompressible flow, by Sani and Gresho [54] led to the concludingremark: “We have made some attempts at shedding more light on the difficult and unresolvedarea of seeking good OBCs for incompressible flow simulations. It has been an exercise infrustration and we are not thrilled with the results obtained”.23 r (b)(a) FIG. 18: The effect of the p = 0 boundary condition on the liquid-gas interface comparing todownstream boundary locations: z = 6 (thicker line) and z = 10 (thinner line):(a) data as in case1 and Q = 0 . Q = 0 . There is an evident inaccuracy involved in our p = 0 choice: the pressure jump associatedto an interface will lead to high local pressure inside the liquid jet or droplets. However,this assumption can be reconciled with our aim, which is not a study of the breakup processand its transient geometries. We are addressing a wider scale: the cone-jet flow pattern, andthe general drop generation regime. To show that the distortion caused by this artificial BCis local and does not modify the global behavior at the cone-jet region, some explorationas been carried out. It can be shown that setting the external boundary sufficiently fardownstream from the nozzle region, at z out = 10, the meniscus-jet is not affected by the24 IG. 19: The effect on the liquid-gas interface of the spatial resolution: (a) a general view of case 2, Q = 0 . z ) min = (∆ r ) min = 0 . z ) min = (∆ r ) min = 0 . boundary condition. To show this, we have considered the worst scenario: we choose largeliquid flow rates Q and weak gas flow (case 1). Fig. 18 shows the stabilized liquid-gasinterface for case 1 and two different values of Q , computed in the original domain and ina shorter one. In Fig. 18(a), the jet does not break up within any of the two numericaldomains and the jet and meniscus interface in the nozzle region is evidently not affected bythe artificial p = 0 boundary condition. The influence of the artificial BC is confined to afew diameters upstream of the downstream boundary. In the case considered in Fig. 18(b),the jet is breaking up into drops within the large domain. Even in this case, the meniscusand jet interface in the nozzle region is not affected by the artificial boundary condition.Let us now show the consistency of the model by comparing the results in two differentmeshes. In this case again, we have considered the worst scenario, by selecting smaller valuesof the liquid flow rate ( Q small) and a large gas flow (case 2), since thinner jets are obtained25
20 40 60 80 100 12002468 v o r t i c i t y ω F r equen cy c on t en t (a)(b) FIG. 20: a) Comparison of the time evolution of the vorticity magnitude at the point ( z = 9 . r = 0 . Q = 0 . t = 0 .
014 (solid line) and with ∆ t = 0 . in these cases. Fig. 19 (a) shows a instantaneous picture of a steady (convectively unstable)liquid jet breaking up into drops (case 1, Q = 0 . ω ∼ .
46 is related to the passage of the vortices generated in the gas mixing layer. Thereis a secondary characteristic frequency, ω ∼ .
56, associated to the interaction between themixing layer and the vorticity wake of the drops. This interpretation was strengthened byrecomputing the flow without the liquid jet: the frequency content of the vorticity signal atthe same observation point only showed the main frequency peak ω ∼ .
46 .
V. COMPARISON WITH ANALYTICAL MODELS AND SCALING LAWS
The first predictive model for the jet diameter d j at the orifice exit [11] assumes thatviscous and capillary effects are small enough compared to liquid inertia. This demandslarge enough Reynolds and Weber numbers of the liquid jet, in reasonable agreement withmost experimental conditions (common solvents including water, down to the micron scale).In this limit, the overall pressure difference ∆ P = ∆ P ( Q g ) (pressure difference between thegas inlet and the gas outlet) imposed in the downstream direction (i.e. through the orifice),transmitted to the liquid stream by normal surface stresses, is converted into kinetic energy,so that ∆ P (cid:39) ρ l v (cid:39) Q l π d j , (9)27
10 100 10000.1110100 Q l /Q o d j / d o FIG. 21: Jet size measured at the entrance of the nozzle using ¯ d in ( (cid:50) correspond to case 1 and ∗ to case 2) and at the exit using ¯ d out ( o correspond to case 1 and × to case 2) compared to thetheoretical prediction (continuous line). which readily gives d j = (cid:18) ρ l π ∆ P (cid:19) / Q / l . (10)Furthermore, the jet is assumed sufficiently small compared to the orifice diameter D suchthat no only it does not touch the orifice borders, but also the boundary layers of thefocusing fluid (gas) at the orifice and at the jet’s surface are sufficiently small compared tothe corona defined between the jet and the orifice. This is why D does not enter expression(10). Neither does D have any direct influence on the jet diameter; only as a parameterdetermining the liquid flow rate Q l .Interestingly, if viscous effects and surface tension are neglected, and we assume d j <
1. Correction for surface tension effects
The liquid surface tension reduces the effective pressure drop ∆ P l in the liquid stream as∆ P l = ∆ P − σ/d j . (14)Consequently, the jet velocity decreases and its diameter increases accordingly. The re-sulting expression for the non-dimensional jet diameter d j /d o , neglecting third order termsproportional to O ( d o /d j ) <<
1, reads: d j /d o = (8 /π ) / ( Q l /Q o ) / + 1 / . (15)In other words, the second order correction of the jet diameter d j to account for surfacetension effects is asymptotically equal to d o /
2. Correction for liquid viscosity effects (extensional stresses)
Assuming that the extensional viscous forces in the liquid are smaller than inertia, thebalance of the different terms of the momentum equation, including the second order termsof the expansion, leads to the following order of magnitude for the correction to the firstorder diameter (10): d e = O (cid:34) d µ (cid:18) Q max Q l (cid:19) / (cid:35) (16)30 . Correction for tangential stresses owing to the gas stream In the same way, the diameter correction (decrease) owing to the momentum injected bythe much faster gas stream through the jet surface is of the order of: d g = O (cid:18) µ g U g D ∆ P (cid:19) / (17)where µ g and U g are the gas viscosity and velocity. The latter is of the order of U g ∼ O (∆ P/ρ g ) / , where ρ g is the gas density.The relative weight of these three corrections provides information on the importance ofthe surface tension and the viscosity of the liquid and gas phases. Interestingly, for mostcommon solvents, these relative weights are of order unity. This happens to be the casewhen measuring the relative importance of the surface tension and the gas tangential stresseffects for water focused by any gas at standard conditions. Therefore, since both correctionsare opposite, the best agreement with experimentally measured jet diameters and numericalsimulations is obtained, interestingly enough, using the first order expression (10), or itsalternative forms (11-13).
4. Correction owing to the nozzle flow pattern
The jet diameter as measured at the nozzle may also differ from the simplest theoreticalprediction given by Eq. (10) because of local flow effects. A complex but symmetric structuredevelops owing to the coexistence of (i) a core potential flow and (ii) the detachment of aradially convergent boundary layer at the inner lip of the nozzle. In any real situationwhere the gas viscosity is non-zero and the continuum hypothesis holds, this flow patternis not aptly described by the pure potential flow through a round orifice given by Morseand Feshback[55] (page 1294) for a stationary discharge. The potential flow solution ischaracterized by an axial velocity distribution with a minimum value at the axis, v ( r = 0) =2 Q g / ( πD ) (half the average velocity through the orifice), and an infinite value at r = D/ Q g being the theoretical gas flow rate discharged. The actual flow geometry is characterizedby the well known vena contracta effect, a consequence of the radial momentum carriedby the collapsing potential flow, which slips at the nozzle border owing to the boundarylayer. The vena contracta flow exhibits an axial velocity distribution which echoes the31otential flow solution, showing a local minimum velocity at the axis, and a maximumvalue at the streamlines coming just from the outside of the boundary layer detached at theorifice (see Fig. 9). The immediate consequence of this particular flow structure is that thetransversal pressure gradients are negligible only sufficiently far downstream of the innerlip of the nozzle: in fact, they become negligible at the axial downstream station wherethe vena contracta effect ends, i.e. where the streamlines become almost parallel. It occursrelatively close to the inner orifice plane, at an approximate D/ A. Scaling of the recirculation zone
For a given gas flow rate Q g and orifice diameter D = 2 R , the typical gas velocity close tothe meniscus surface can be estimated as V = Q g / ( πR ). Given the small ρ = ρ g /ρ l valuesin liquid jets focused by gas, liquid velocities are much smaller than V everywhere. As theliquid approaches the neck, the boundary layer will collapse (Fig. 22). This implies that atleast a liquid flow rate Q R ∼ U s δ l (18)would be drawn into the jet in the absence of recirculation ( U s is the velocity of the interface,that can be obtained from V , and the densities and the viscosities ratios[56]). On thecontrary, whenever Q l < Q R , part of Q R must have been recirculated back into the meniscus(Fig. 22). Therefore Q R can be interpreted as the minimum flow rate for no recirculation(scaled as Q r = Q R /Q g ).The boundary layer in the liquid meniscus is confined. It grows along the cone duringlengths comparable to R (the orifice radius) till the apex of the meniscus is reached. In thisarea, the gas speed gradients are steep: δ l ∼ ( µ l R/ρ l U s ) / . (19)Whenever there is recirculation, the peripheral boundary layers merge at the meniscus apexand give rise to a jet, whose initial radius at the neck will accordingly be of the same order.32 δ l δ r S D RD ≡ l Q B Q s U l Q FIG. 22: Sketch of the recirculation zone, showing boundary layers, cell size ( S r ), and typicalvelocities. In the absence of liquid emission, maximum recirculation will be observed. Experimentally,however, a dripping instability will occur before reaching this limit. In the opposite case (norecirculation), the boundary layers do not merge, and an inviscid core should be observed atthe neck. The threshold flow rate for recirculation can therefore be estimated as Q R ∼ U s δ l ,a result which happens to be independent of the gas velocity. In effect, by definition ofthe meniscus boundary layer, the viscous stress µ l U s /δ l must be of the same order as themomentum convection ρ l U s /R , so that, interestingly: Q R ∼ R µ l ρ l = ⇒ Q r ∼ ρµRe (20)This scaling is fully confirmed by the numerical simulations: the values of Q r · Re are 0.6768for case 1 and 0.6596 for case 2, deviating by less than 2.6% from the scaling predictions.Assume now the recirculation cell to be S r in axial length. The backflow Q B = Q R − Q l will come to rest within a length of the order S r . In this length, viscous momentum diffusionshould slow down the flow and deflect both the incoming flow injected by the feeding tube33nd the recirculated flow at the axis (Fig. 22). Thus, viscous and inertia forces shouldbalance within that length S r : in other words, the liquid Reynolds number associated toaxial lengths of order S r should be of order unity so as to deflect the unidirectional flowissuing from the feed tube (Hagen-Poiseuille). This is in analogy to the entry length or exitlength in laminar pipe flow. Two cases need to be considered, depending of the relative sizeof the cell compared to the feed tube radius R : • When S r < R , viscous stress, of the order O ( µ l Q B S − r ), balances inertia, O ( ρ l Q B S − r ),which leads to S r ∼ ρ l Q B µ − l . • When S r > R , viscous stress, O ( µ l Q B R − ), balances inertia, O ( ρ l Q B R − S − r ), leadingagain to S r ∼ ρ l Q B µ − l .Interestingly enough, again, the length of the recirculation flow is independent of the gasflow for any given geometry. The latter scaling can be expressed in non-dimensional termsas: s r ≡ S r /R ∼ µ l ( Q R − Q l ) ρ − l R − . (21)Using equation (18), one may write: s r = C − C Re R , (22)where Re R = ρ l Q l / ( µ l R ) is a Reynolds number of the liquid flow, and C and C areconstants which depend on the geometry only (i.e., R /R , H/R , etc.). In our case, we haverepresented all our measured s r values from numerical simulations versus Re R in Fig. (23).Linear fitting to all points leads to C = 2 .
636 and C = 0 . . Q r − Q and Re as well, as: s r = kµρ − ( Q r − Q ) , (23)where k is again a constant depending on the geometry only, in full agreement with expression(8), as anticipated by experiments. VI. EXPERIMENTS
In the following, we provide experiments corresponding to the same local geometricalparameters as in cases 1 and 2 in the vicinity of the exit orifice. The basic flow focusing34
10 20 30 4000.511.522.53 Re R s r S r (Re R )=C −C *Re R Re=931.66Re=698.75Re=465.83
FIG. 23: Recirculation cell size s r as a function of Re R : dots, squares and stars are obtained bynumerical simulation; the line is a theoretical prediction resulting from dimensional arguments.An additional series of simulations have been performed for an intermediate gas flow condition( Re = 698 .
75 and
W e = 18 .
31) to assess the validity of the scaling proposed: note the gooddegree of collapse obtained. The small deviations can be attributable to the small differences inthe geometry of the cone for different gas flow conditions. chamber is a box consisting of five aluminum faces and one clear methacrylate face. It is 5cm by 5 cm by 5.65 cm, with its longest side along the capillary/orifice axis. The chamberis situated with the methacrylate face horizontal and pointing upwards, the capillary beinglocated parallel to this face. The orifice is made in a stainless steel orifice disk attached tothe box side, perpendicular and opposite to the capillary. The disk is 4.0 mm in diameterwith a thickness of 75 µ m and an orifice of diameter 0.200 mm. Both the air tube andthe capillary enter through the face opposite the orifice. After the capillary tube has beenaligned with the orifice, the distance H from the tube to the orifice can be simply adjustedby carefully sliding the capillary in its housing on the opposite face to the orifice disk. H is measured with a microscope through the methacrylate face. Fig. 24 shows some views35 μ m200 μ m150 μ m FIG. 24: (Left) Experimental tube-orifice set up as numerically simulated in this work ( D =200 µ m, D = 150 µ m, H = 100 µ m; here, ∆ P = 10 KP a , Q l = 3 mL/h). (Right) Photographsof experimental conditions with twice the distance from the feeding tube to the exit orifice, usinga different tube material (fused silica): (A) jetting ( D = 200 µ m, D = 150 µ m, H = 200 µ m,∆ P = 30 KP a , Q l = 6 . Q l = 2 . of the feeding tube-orifice setup as seen through the thick methacrylate window (inevitableliquid spills leave behind some debris on the inner face of the window causing a blurredimage). In particular, Fig. 24 (left) shows the geometry numerically simulated in this work.After setting H and ensuring that the capillary is perfectly coaxial with the nozzle orifice,the pressure is set using a pressure gauge and a pressure meter. A water flow rate is thensupplied using a syringe pump (Cole-Palmer 74900 Series) with a 20 ml syringe. The systemis given sufficient time to relax until either a characteristic steady or unsteady flow is present.This can be checked by illuminating the jet that exits the orifice or by looking at the meniscuswhen the distance H is 0.100 mm or greater. Unsteady jet flow appears very faint to thenaked eye and contains thin streaks of water along with large scattered spray. This is insignificant contrast to steady jet flow, which has bright illumination as a result of a finer,concentrated stream with uniform characteristics. In experiments where the meniscus wasvisible, it was also possible to discriminate steady versus unsteady flow (see Fig.24 (right)A -jetting- and B -dripping-), in perfect correlation with the spray observations: a steadymeniscus had sharp edges and a clear, unwavering glass-like appearance (Fig. 24 A, seesteady jet reflected in the metal plate), while an unsteady meniscus had blurred edges andflickered (Fig. 24 B, no jet is visible at all). Both the jet test and meniscus test displayed36 Re l We l H=0,100mm UpH=0,100mm DownNumerical, jettingNumerical, drippingCase 2Case 1 (gas pressure increasing)
FIG. 25: Jetting -dripping transition in the { Re l , We l } plane. Diamonds: experimentally deter-mined conditions (filled symbols: liquid flow rate decreasing -“Down”; open symbols: liquid flowrate increasing-“Up”. In most cases, both “Up” and “Down” points coincide). Circles: numericallytested conditions. Filled circles: jetting conditions. Open circles: dripping conditions. clear and abrupt transitions between the two states. Once unsteady flow is establishedfor a given pressure, the rate determined by the syringe pump is increased in steps of 0.1ml/hr. After each flow rate increase, a 30 seconds waiting period was established, so as toensure that the system had relaxed and all the readings were accurate. This period has beenchosen after it was found that 15 s was not enough to observe fluctuations in the system:occasionally a steady regime would revert back to an unsteady one after the 15 s period. The30 s delay has proven long enough to accurately characterize the flow. Accordingly, the ratewas increased until the unsteady jet sharply transitioned to a steady one; at this point, theflow rate was read from the syringe pump and recorded as the minimal flow rate (increasing,or “up”; Fig. 25). Keeping this same steady flow rate the process is then reversed to findthe minimum decreasing (or “down”) flow rate (steps of 0.1 ml/hr and intervals of 30 s until an unsteady regime developed). When the flow became unsteady, a rate 0.1 ml/hrabove the reading on the syringe pump was recorded, since the rate which produced the last37teady flow (i.e. minimum flow rate) was one step (0.1 ml/hr) higher. The resulting valuewas recorded as the minimal (decreasing, or “down”) flow rate. This process is repeated forvarying pressures and distances of H to get an accurate mapping of minimal jetting flowrates as a function of varying geometry and flow conditions. Following this procedure, wecollected the experimental data plotted in Fig. 25 for H/R = 1. The gas (air) pressure ∆ P increases as indicated by the arrow.Six conditions numerically tested for cases 1 and 2 are plotted in Fig.25. In order tomake our results readily translatable in most of the capillary jet stability literature (whichuses the jet radius as a characteristic length), we may introduce liquid Reynolds and Webernumbers consistent with previous definitions and using scaling law (10): Re l = (cid:18) π (cid:19) / (cid:18) ρ l Q l ∆ Pµ l (cid:19) / , W e l = (cid:18) π (cid:19) / (cid:18) ρ l Q l ∆ P σ (cid:19) / (24)As it follows from the plot, using these definitions, jetting or dripping conditions are accu-rately predicted by the numerical model. This lends additional support to the use of fullVOF simulation analysis in flow focusing systems. VII. CONCLUSIONS
The cone-jet geometry associated with flow focusing has been handled by a diversityof tools, numerical, experimental and theoretical. Order-of-magnitude estimations followfrom dimensional arguments: such procedures contribute a valuable theoretical framing andprovide the scaling criteria for data representation. Analytical approaches are generallybased on the consideration of a perfectly cylindrical infinite jet, a simplification that ignoresthe influence of the meniscus (a source of instability) and the role of streamline convergenceor divergence in the jet. Experiments are burdened by the diversity of influencing parametersand visualization difficulties associated with the small scale of the meniscus and jet.In this paper, experimental results are backed up by a numerical simulation based on VOFelements. Numerical schemes allow a more systematic exploration of the parametric influ-ence. In addition, the shortcomings of theoretical models (unavoidable in a situation wherethe geometry of the fluid domain is complex, as in a cone-jet flow pattern) are overcome,and a detailed description of the streamlines can be readily obtained.The key results of the above exploration are the following:38
The theoretical scaling leading to jet diameter estimates is confirmed by the simulation.The expressions for flow focusing scales, notwithstanding their simplicity, are thereforeto be considered a reliable shortcut for the prediction of jet dimensions. • The complete sequence from meniscus growth to jet emission (jetting regime) and tothe sequential filling of drops (dripping regime) is portrayed in detail. • The jetting-dripping transition is documented in detail, both by experiment and sim-ulation. A two-branch structure is observed in the plot, showing the simultaneousinfluence of the jet and the meniscus as instability sources. Incipient dripping (Fig.13) is shown to give rise to highly irregular fluctuations; while fully developed dripping(Fig. 15) produces perfect cycles of drop detachment. • A recirculation cell is identified in the jetting regime at the meniscus tip. This oc-currence appears to be linked to intensive forcing by the gas sheath, leading to highinterface velocity along the meridians; the issuing jet is unable to convey all of themobilized flow, so that a return flow around the axis is observed. The recircula-tion cell grows as the liquid flow rate is reduced: eventually, dripping conditions arereached. Similar recirculation cells have been observed in electrospray cones, underthread emission, and in liquid-liquid two-dimensional flow focusing, assisted or not bya surfactant [5, 29, 31]. All of the recirculation instances reported thus far appearto share a common attribute: strong interfacial forcing, either electric, capillary orhydrodynamic. • A reliable scaling is provided, identifying the parametric conditions where recirculationis to be expected and estimating the size and flow rates of the cell.A key feature in the flow pattern explored is the recirculation cell, and its conceptuallink to the merging of the boundary layers which grow from the meniscus edge and fusetogether at the neck of the jet. Controllable recirculation is an extremely attractive feature,providing adjustable residence times within a very simple flow setup. The cell can be viewedas a flow trap or reactor, where biosynthesis or chemical operations take place in a protectedenvironment; the liquid flow rate can be increased to flush the recirculation products.An additional focus deals with the peculiarities of the jetting and dripping regimes underthe influence of a co-flowing sheath current. The aspect of the jet and droplet train and the39ynamics of the meniscus (an indicator of dripping) are a contribution to a problem whosecomplexity forbids a global theoretical approach.
VIII. ACKNOWLEDGMENTS
This work has been supported by the Spanish Ministry of Science and Education, projectnumber DPI2004-07197, and partially supported by European Commission through grantCOOP-CT-2005-017725. Thorough discussions of one of the authors (AMGC) with Dr.Joan Rosell are warmly acknowledged. 40
1] O. A. Basaran. Small-scale free surface flows with breakup: Drop formation and emergingapplications.
AIChE J. , 48:1842–1848, 2002.[2] H. Song, D. L. Chen, and R. F. Ismagilov. Reactions in droplets in microfluidic channels.
Angew. Chem. Int. Ed. , 45(44):7336 – 7356, 2006.[3] A. D. Griffiths and D. S. Tawfik. Miniaturising the laboratory in emulsion droplets.
TrendsBiotechnol. , 24(9):395–402, 2006.[4] G. I. Taylor. The formation of emulsions in definable fields of flow.
Proc. R. Soc. London,Ser. A , 146:501, 1934.[5] R. Suryo and O. A. Basaran. Tip streaming from a liquid drop forming from a tube in aco-flowing outer fluid.
Phys. Fluids , 18:082102, 2006.[6] W. W. Zhang. Viscous entrainment from a nozzle: Singular liquid spouts.
Phys. Rev. Lett. ,93:184502, 2004.[7] R. J. Haywood, M. Renksizbulut, and G. D. Raithby. Transient deformation of freely-suspended liquid droplets in electrostatic fields.
AIChE J. , 37(9):1305–1317, 1991.[8] R. T. Collins, J. J. Jones, M. T. Harris, and O. A. Basaran. Electrohydrodynamic tip streamingand emission of charged drops from liquid cones.
Nature Phys. , 4:149–154, 2007.[9] E. Lac and G. M. Homsy. Axisymmetric deformation and stability of a viscous drop in asteady electric field.
J. Fluid Mech. , 590:239, 2007.[10] H. N. Oguz and A. Prosperetti. Dynamics of bubble-growth and detachment from a needle.
J. Fluid Mech. , 257:111–145, 1993.[11] A. M. Ga˜n´an-Calvo. Generation of steady liquid microthreads and micron-sized monodispersesprays in gas streams.
Phys. Rev. Lett. , 80:285–288, 1998.[12] A. M. Ga˜n´an-Calvo and J. M. Gordillo. Perfectly monodisperse microbubbling by capillaryflow focusing.
Phys. Rev. Lett. , 87:274501, 2001.[13] S. Takeuchi, P. Garstecki, D. B. Weibel, and G. M. Whitesides. An axisymmetric flow-focusingmicrofluidic device.
Adv. Mater. , 17:1067–1072, 2005.[14] S. L. Anna, N. Bontoux, and H. Stone. Formation of dispersion using flow-focusing in mi-crochannels.
Appl. Phys. Lett. , 87:364, 2003.[15] P. Garstecki, I. Gitlin, W. DiLuzio, E. Kumacheva, H. A. Stone, and G. M. Whitesides. ormation of monodisperse bubbles in microfluidic flow-focusing device. Appl. Phys. Lett. , 85:2649, 2004.[16] M. Seo, Z. H. Nie, S. Q. Xu, M. Mok, P. C. Lewis, R. Graham, and E. Kumacheva. Continuousmicrofluidic reactors for polymer particles.
Langmuir , 21:11614 –11622, 2005.[17] M. J. Jensen, H. A. Stone, and H. Bruus. A numerical study of two-phase Stokes flow in anaxisymetric flow-focusing device.
Physics of Fluids , 18, 2006.[18] A. S. Utada, E. Lorenceau, D. R. Link, P. D. Kaplan, H. A. Stone, and D. A. Weitz. Monodis-perse double emulsions generated from a microcapillary device.
Science , 308:537–541, 2005.[19] C. Berkland, E. Pollauf, D. W. Packa, and K. Kim. Uniform double-walled polymer micro-spheres of controllable shell thickness.
J. Control. Release , 96:101–111, 2004.[20] Luc´ıa Mart´ın-Banderas, A. Rodr´ıguez-Gil, A. Cebolla, S. Ch´avez, T. Berd´un-Alvarez, J. M.Fernandez-Garcia, M. Flores-Mosquera, and A. M. Ga˜n´an-Calvo. Flow focusing: a versatiletechnology to produce size-controlled and specific-morphology microparticles.
Small , 1:688–692, 2005.[21] L. Mart´ın-Banderas, M. Flores-Mosquera, P. Riesco-Chueca, A. Rodr´ıguez-Gil, A. Cebolla,S. Ch´avez, and A. M. Ga˜n´an-Calvo. Towards high-throughput production of uniformly en-coded microparticles.
Adv. Mater. , 18:559–564, 2006.[22] A. M. Ga˜n´an-Calvo and P. Riesco-Chueca. Jetting-dripping transition of a liquid jet in alower viscosity co-flowing immiscible liquid: the minimum flow rate in flow focusing.
J. FluidMech. , 553:75–84, 2006.[23] A. M. Ga˜n´an-Calvo, M. A. Herrada, and P. Garstecki. Bubbling in unbounded coflowingliquids.
Phys. Rev. Lett. , 96:124504, 2006.[24] A. M. Ga˜n´an-Calvo. Absolute instability of a viscous hollow jet.
Phys. Rev. E , 75:027301,2007.[25] D. Erickson. Towards numerical prototyping of labs-on.chip: modeling for integrated microflu-idics devices.
Microfluid Nanofluid , 1:301–318, 2005.[26] V. Cristini and Y.-Ch. Tan. Theory and numerical simulation of droplet dynamics in complexflows-a review.
Lab Chip , 4:257–264, 2004.[27] M. M. Dupin, I. Halliday, and C. M. Care. Simulation of a microfluidic flow-focusing device.
Physical Review E , 73, 2006.[28] M. W. Weber and R. Shandas. Computational fluid dynamics analysis of microbubble forma- ion in microfluidic flow focusing devices. Microfluid Nanofluid , 3:195–206, 2007.[29] A. Barrero, A. M. Ga˜n´an-Calvo, J. D´avila, A. Palacio, and E. G´omez-Gonz´alez. Low andhigh Reynolds number flows inside Taylor cones.
Phys. Rev. E , 58:7309–7314, 1998.[30] M. A. Herrada and A. Barrero. Self-rotation in electrocapillary flows.
Physical Review E(Statistical, Nonlinear, and Soft Matter Physics) , 66(3):036311, 2002. URL http://link.aps.org/abstract/PRE/v66/e036311 .[31] S. L. Anna and H. C. Mayer. Microscale tipstreaming in a microfluidic flow focusing device.
Phys. Fluids , 18:121512, 2006.[32] J. M. Fernandez and G. M. Homsy. Chemical reaction-driven tip-streaming phenomena in apendant drop.
Phys. Fluids , 16:2548–2555, 2004.[33] R. Krechetnikov and G. M. Homsy. On physical mechanisms in chemical reaction-driventip-streaming.
Phys. Fluids , 16:2556–2566, 2004.[34] D. B. Kothe and W. J. Rider. A comparison of interface tracking methods.
Fluid Dynamics ,pages 19–22, March 1995.[35] Salih Ozen Unverdi and Gr´etar Tryggvason. A front-tracking method for viscous, incom-pressible, multi-fluid flows.
J. Comput. Phys. , 100(1):25–37, 1992. ISSN 0021-9991. doi:http://dx.doi.org/10.1016/0021-9991(92)90307-K.[36] Dalton J. E. Harvie and David F. Fletcher. A new volume of fluid advection algorithm: thestream scheme.
J. Comput. Phys. , 162(1):1–32, 2000. ISSN 0021-9991. doi: http://dx.doi.org/10.1006/jcph.2000.6510.[37] Eugenio Aulisa, Sandro Manservisi, and Ruben Scardovelli. A mixed markers and volume-of-fluid method for the reconstruction and advection of interfaces in two-phase and free-boundaryflows.
J. Comput. Phys. , 188(2):611–639, 2003. ISSN 0021-9991. doi: http://dx.doi.org/10.1016/S0021-9991(03)00196-7.[38] S. Tanguy and A. Berlemont. Application of a level set method for simulation of dropletcollisions.
International Journal of Multiphase Flow , 31:1015–1035, 2005.[39] For convenience and with the aim of making our results readily reproducible from others, inthis work we have used the well tested commercial solver FLUENT v 6.3. Mesh generation isaided by the software GAMBIT in FLUENT v 6.2.[40] C.W. Hirt and B.D. Nichols. Volume of fluid (vof) method for the dynamics of free boundaries.
J. Comput. Phys. , 39:201, 1984.
41] P. Heinrich. Nonlinear numerical model of landslide-generated water waves.
Int. J. Engrg.Fluid Mech. , 4:403, 1991.[42] A. Tomiyama, A. Sou, H. Minagawa, and T. Sakaguchi. Numerical analysis of a single bubbleby vof method.
JSME Int. J , 36:51, 1993.[43] B. Lafaurie, C. Nardone, R. Scardovelli, S. Zaleski, and G. Zanetti. Modelling merging andfragmentation in multiphase flows with surfer.
J. Comput. Phys. , 113:134, 1994.[44] B. Van Leer. Toward the ultimate conservative difference scheme. iv. a second order sequel togodunov’s method.
J. Comput. Phys. , 32:101–136, 1979.[45] O. Ubbink.
Numerical prediction of two fluid systems with sharp interfaces . PhD thesis,Imperial College of Science, Technology and Medicine, London, 1997.[46] R. I. Issa. Solution of implicitly discretized fluid flow equations by operator splitting.
J.Comput. Phys. , 62:40–65, 1986.[47] R. Shaw.
The Dripping Faucet as a Model Chaotic System . Aerial, Santa Cruz, 1984.[48] D. R. Chen, D. Y. H. Pui, and S. L. Kaufman. Dripping faucet dynamics by an improvedmass-spring model.
J. Phys. Soc. Japan , 68:3259–ˆu3270, 1999.[49] B. Ambravaneswaran, S. D. Phillips, and O. A. Basaran. Theoretical analysis of a drippingfaucet.
Phys. Rev. Lett. , 85:5332, 2000.[50] B. Ambravaneswaran, H. J. Subramani, S. D. Phillips, and O. A. Basaran. Dripping-jettingtransitions in a dripping faucet.
Phys. Rev. Lett. , 93(3):034501, 2004.[51] H. J. Subramani, H. K. Yeoh, R. Suryo, Q. Xu, B. Ambravaneswaran, and O. A. Basaran.Simplicity and complexity in a dripping faucet.
Phys. Fluids , 18:032106, 2006.[52] P. Coullet, L. Mahadevan, and C. S. Riera. Hydrodynamical models for the chaotic drippingfaucet.
J. Fluid Mech. , 526:1–17, 2005.[53] E. Villermaux, P. Marmottant, and J. Duplat. Ligament-mediated spray formation.
Phys.Rev. Lett. , 92(7):074501, 2004.[54] R. L. Sani and P. M. Gresho. Resume and remarks on the open boundary condition minisym-posium.
International Journal for Numerical Methods in Fluids , 18(10):983–1008, 1994.[55] P. M. Morse and H. Feshback.
Methods of theoretical physics . McGraw-Hill Pub. Co., 1953.[56] A. M. Ga˜n´an-Calvo. Electro-flow focusing: The high conductivity, low viscosity limit.
Phys.Rev. Lett. , 98:134503, 2007., 98:134503, 2007.