Elucidating the mechanism of step-emulsification
Andrea Montessori, Marco Lauricella, Sauro Succi, Elad Stolovicki, David A. Weitz
EElucidating the mechanism of step-emulsification
Andrea Montessori ∗1 , Marco Lauricella , Sauro Succi , Elad Stolovicki , and DavidWeitz Istituto per le Applicazioni del Calcolo CNR, via dei Taurini 19, Rome, Italy Center for Life Nano Science@La Sapienza, Istituto Italiano di Tecnologia, 00161 Roma, Italy Institute for Applied Computational Science, Harvard John A. Paulson School ofEngineering And Applied Sciences, Cambridge, MA 02138, United States School of Engineering and Applied Sciences, Harvard University, McKay 517 Cambridge, MA02138, USA Department of Physics, and School of Engineering and Applied Sciences, Harvard University,Pierce 231, 29 Oxford St. Cambridge, MA 02138, USAJune 22, 2020
Abstract
Three-dimensional, time-dependent direct simulations of step emulsification micro-devices highlight twoessential mechanisms for droplet formation: first, the onset of an adverse pressure gradient driving a back-flow of the continuous phase from the external reservoir to the micro-channel. Second, the striction ofthe flowing jet which leads to its subsequent rupture. It is also shown that such a rupture is delayed andeventually suppressed by increasing the flow speed of the dispersed phase within the channel, due to thestabilising effect of dynamic pressure. This suggests a new criterion for dripping-jetting transition, based onlocal values of the Capillary and Weber numbers.
Step emulsification (SE) has captured significant interest in the recent years as a viable microfluidic techniquefor the controlled production of liquid droplets [1, 2]. Among others, one of the main appeals of the SE techniqueis the prospect of producing large volume rates of the dispersed phase, which are out of reach for previoustechniques, such as flow-focusers [3, 4, 5, 6].The basic idea behind SE is to exploit the pressure drop due to a sudden channel expansion (step) to inducethe pinch-off of the dispersed phase, leading to droplet formation [7]. Albeit conceptually straightforward,the details of the process depend on a number of physical and geometrical parameters, primary the capillarynumber Ca and the aspect ratio h/w of the height versus width of the micro channel cross-section (see Figure1). Such parameters dictate the shape of the droplet and the transition between the dripping and the jettingregimes[7, 8, 9]. Although of primary importance, the Capillary number (viscous dissipation/surface tension)does not capture the full picture and needs to be complemented by other dimensionless groups, namely theWeber number (inertia/surface tension) and/or Reynolds numbers (inertia/viscous dissipation).Despite major technological advances, the theoretical description and the numerical simulation of micro-channel emulsification is still under development.In this Letter, we present direct numerical simulations of the fully three-dimensional, time-dependent Navier-Stokes equations for a specific step-emulsification micro-device, in order to elucidate the basic fluid phenomenaunderpinning the step-emulsification process. The simulations highlight two essential mechanisms: i) the back-flow of the continuous phase from the external reservoir to the confined micro channel, driven by an adversepressure gradient, ii) the resulting striction of the flowing jet within the channel and its subsequent rupture. Itis also shown that such a rupture is delayed or even suppressed upon increasing the flow speed of the dispersedphase within the channel, due to the stabilising effect of dynamic pressure.In order to simulate the droplet breakup in a recently proposed class of step emulsification devices [10, 2, 8]we solve the multi-component Navier-Stokes equations, using an extensions of the Lattice Boltzmann method[11, 12, 13, 14, 15, 16, 17]. See Supplemental Material at 2 for a detailed description of the numerical modelemployed [18, 19, 20]. The device, made of polydimethylsiloxane (PDMS), is used for producing water in oilemulsions and a sketch is reported in figure 1. The water flows through the device inlet, and splits into hundredsof step-emulsifier nozzles with rectangular cross sections. The PDMS device is submerged in quiescent oil, thecontinuous phase, with the nozzles pointing upwards. The dispersed phase (water), is then pumped through ∗ Electronic address: [email protected] ; Corresponding author a r X i v : . [ phy s i c s . c o m p - ph ] J un igure 1: Sketch of the nozzle geometry in the simulation box, along with the imposed boundary conditions(top panel a). The adopted conditions reproduce a periodic array of independent nozzles, which is consistentwith the geometry of the volcano device (bottom panel b). Here, the dispersed phase (red) is pumped throughthe device, forming mono-disperse drops in a reservoir containing a continuous immiscible phase (cyan).the device and forms monodisperse drops, whose sizes are proportional to the nozzle height ( h ). We wish topoint out that, being the device submerged in quiescent oil, there is no net flow of the continuous phase, thisin contrast with the emulsification system employed in [21]. In this work, we simulate a single nozzle out of2Figure 2: a) sequence of the break process in the dripping regime in a step emulsification nozzle. The dispersedphase (water), is pumped through the device and forms monodisperse drops whose sizes are proportional to thenozzle height. b) Sequence of the step emulsifier in the jetting regime Ca d / h Dripping J e tt i n g Figure 3: Dimensionless droplet diameter ( d/h ) versus the Capillary number for a nozzle aspect ratio h/w =1 /
5. The dots are the normalized droplet diameters predicted by the numerical simulations, while the trianglesstand for the experimental values of the normalized diameters. In the dripping region, for Ca between 0 .
002 and O (10 − ) the average value of the droplet diameter is ∼ µm ( d/h ∼ ∼ .
2, as in [8], [22, 23, 24, 25].First, we perform a comparison with the experiments reported in [8], by investigating the effect of the mainnon-dimensional parameters i.e., the capillary number ( Ca = ρU ν/σ ) and the Weber number (We= ρU L/σ ),where U is the velocity of the dispersed phase at the inlet, ρ is the density of the dispersed phase, L is acharacteristic length defined as √ wh (see fig. 1), σ is the oil-water surface tension and ν is the kinematicviscosity of the dispersed phase.In our simulations the dispersed phase flows through nozzles with a rectangular cross section of 135 × µm and length L = 810 µm , with a characteristic nozzle aspect ratio h/w = 1 /
5, being h and w the height andthe width of the microchannel, respectively. Experiments show that the droplet sizes are nearly independent ofthe flow rate over an extended range (between 12 and 70 mL/min , 0 . < Ca < . d = 567 ± µm . In this range of flow rates, thus for Capillary numbers in the range O (10 − ) ÷O (10 − ), the step emulsifier has been shown to operate in the dripping regime, while for larger flow rates,i.e. for Capillary numbers larger than a critical value, a transition from dripping to jetting regime occurs,which is characterized by the production of much larger and polydisperse droplets [2, 10]. The dripping-jettingtransition occurs whenever the droplet does not break up anymore and starts ’ballooning’ (i.e., above a critical We performed additional simulations in which the same range of capillary numbers was spanned by changing the viscosity ν and or the surface tension σ . The simulation results show that the transition from dripping to jetting occurs almost at the samecritical Capillary number, corresponding to lower flow rates as the viscosity of the inner phase is increased. p ∗ ), (b) velocity ( u ∗ ), (c) vorticity ( ω ∗ ) and (d) stress fields ( S ∗ yz ), in a y-z midplanetaken between the two walls separated by h , from the focusing stage (1,2) to the pinching (3) and finally breakup(4). (e) Zoom of the droplet structure and associated flow field at breakup time. In this simulation Ca = 0 . h/w = 1 / . < Ca < . d ∼ µm ( d/h ∼ Ca ∼ . U/h . This sequence displays atypical elongational flow structure, especially in figure 4, panel c3, which stretches the jet until rupture.Finally, in panel (d) we show the component of the stress field in the yz mid-plane, in units of U/h . Thesequence shows that the stress field is highly localised around the oil-water interface, with a null-point right atthe pinch position [26]. The present analysis is in line with previous studies [27], in which the breakup is notinterpreted as due to a Plateau-Rayleigh instability [28], but rather to the back flow of the continuum phase,triggered by the adverse pressure gradient which arises in correspondence with the focusing of the water jet.We wish to emphasize that our analysis is fully dynamic and three-dimensional, i.e it does not rely on anyquasi-static assumption [7], nor on any axial-symmetry of the flow [29].
Dripping to jetting transition.
Most experimental studies of step-emulsification report a dripping to jetting transition above a critical cap-illary number Ca crit ∼ O (10 − ) [1, 8, 10]. However, the underlying mechanisms behind such transition are stillunder investigation. Here, we wish to point out that the transition to the jetting regime is facilitated by thecontribution of the dynamic pressure ρ in u in / ρ in and u in the local density and velocity of the dispersedphase in the pinching region. Such dynamic pressure withstands the effects of the negative-curvature in thepinch region. Due to the pinching effect, the local flow speed within the dispersed phase significantly exceedsthe inlet velocity. Hence, the local capillary number attains larger values, of the order of 0 . −
1, whenever thenominal capillary number of the dispersed phase reaches its critical value around 0 .
01. Note that the nominal4 * Figure 5: Normalized velocity magnitude and vector field in the dripping nozzle. The counterflow in thecontinuous phase within the nozzle, is clearly evidenced by the quiver plot. The solid white line identifies theinterface between the continuous and the dispersed phase while the red line denotes the walls of the nozzle. Theinset in the leftmost panel, highlights the backflow occurring inside the nozzle. P in + r u in /2 + s / d P out P out d Figure 6: Sketch of the static and dynamic pressure components in and out of the neck region.capillary number is computed with the imposed velocity at the inlet.The local acceleration of the flow field inside the pinch region is clearly visible in Figure 5, which reportsthe flow field inside the neck region of the dispersed phase.As pinching progresses in time, the flow speed inside the pinching region grows accordingly, so that, at somepoint, inertial effects can no longer be neglected.To clarify the point, let us write the dynamic force balance under flow conditions, namely: P in + ρ in u in / P out + ρ out u out / − σδ (1)where subscripts "in" and "out" refer to the neck region (see fig.6), inside and outside the water jet, respectively.Note that the minus sign in front of the surface tension reflects the negative curvature (see fig. 6).In equation 1, δ is the characteristic length scale of the neck region, which is found to be smaller butcomparable with the channel height h . This is plausible, because δ > h is not feasible since the neck diametercannot be larger than the height of the nozzle, while δ (cid:28) h signals the imminent breakup.This expression shows that the inner dynamic pressure adds to the surface tension in withstanding the outerpressure. As a result, it is natural to extend the definition of Capillary number so as to include the contributionof the dynamic pressure, namely: K = µu in σ + ρu in δ/ ≡ Ca in W e in (2)where we have neglected the outer velocity since u out /u in ∼ δ/w < Ca in = 0 .
015 and
W e in ∼ .
22, showing that the dynamic pressure is stillsub-dominant with respect to the capillary pressure σ/δ . On the other hand, in the case of jetting, (see panel(c) in Fig. 2)
W e in ∼
1, indicating that the jetting regime is entered whenever the dynamic pressure becomescomparable or higher than the capillary pressure.We wish to point out that since δ ∼ h , Ca in ∼ wh Ca , which is precisely the quantity controlling the dripping-jetting criterion discussed in [21]. In this Letter, we noted that such criterion should also take into account thecontribution of dynamic pressure, which becomes dominant in the pinch region.5 Acknowledgments
The research leading to these results has received funding from the European Research Council under theEuropean Union’s Horizon 2020 Framework Programme (No. FP/2014- 2020)/ERC Grant Agreement No.739964 ("COPMAT").
The LB immiscible multicomponent model is based on the following lattice Bhatnagar-Gross-Krook (BGK)equation: f ki ( ~x + ~c i ∆ t, t + ∆ t ) = f ki ( ~x, t ) + Ω ki ( f ki ( ~x, t )) , (3)where f ki is the discrete distribution function, representing the probability of finding a particle of the k − th component at position ~x and time t with discrete velocity ~c i . The lattice time step is taken equal to 1, and i the index spans the lattice discrete directions i = 0 , ..., b , where b = 26 for a two dimensional nine speed lattice(D3Q27). The density ρ k of the k − th fluid component is given by the zeroth order moment of the distributionfunctions ρ k ( ~x, t ) = X i f ki ( ~x, t ) , (4)while the total momentum ρ~u is defined by the first order moment: ρ~u = X i X k f ki ( ~x, t ) ~c i . (5)The collision operator Ω ki results from the combination of three sub-operators, namely [16]Ω ki = (cid:0) Ω ki (cid:1) (3) h(cid:0) Ω ki (cid:1) (1) + (cid:0) Ω ki (cid:1) (2) i . (6)Here, (cid:0) Ω ki (cid:1) (1) is the standard BGK operator for the k − th component, accounting for relaxation towards a localequilibrium (cid:0) Ω ki (cid:1) (1) f ki ( ~x, t ) = f ki ( ~x, t ) − ω k (cid:16) f ki ( ~x, t ) − f k,eqi ( ~x, t ) (cid:17) , (7)where ω k is the relaxation rate, and f k,eqi ( ~x, t ) denotes local equilibria, defined by f k,eqi ( ~x, t ) = ρ k w i (cid:18) ~c i · ~uc s + ( ~c i · ~u ) c s − ( ~u ) c s (cid:19) (8)Here, w i are weights of the discrete equilibrium distribution functions, c s = 1 / p (3) is the lattice sound speed[ ? ]. In this model, (cid:0) Ω ki (cid:1) (2) is a perturbation operator, modeling the surface tension of the k − th component.Denoting by ~F the color gradient in terms of the color difference (see below), this term reads (cid:0) Ω ki (cid:1) (2) f ki ( ~x, t ) = f ki ( ~x, t ) + A k | ~F | " w i ( ~F · ~c i ) | ~F | − B i , (9)with the free parameters A k modeling the surface tension, and B k a parameter depending on the chosen lattice[17, 18]. The above operator models the surface tension, but it does not guarantee the immiscibility betweendifferent components. In order to minimize the mixing of the fluids, a recoloring operator (cid:0) Ω ki (cid:1) (3) is introduced.Following the approach in Ref. [18], being ζ and ξ two immiscible fluids, the recoloring operators for the twofluids read as follows (cid:16) Ω ζi (cid:17) (3) = ρ ζ ρ f i + β ρ ζ ρ ξ ρ cos( φ i ) X k = ζ,ξ f k,eqi ( ρ k , (cid:16) Ω ξi (cid:17) (3) = ρ ξ ρ f i − β ρ ζ ρ ξ ρ cos( φ i ) X k = ζ,ξ f k,eqi ( ρ k ,
0) (10)where β is a free parameter and cos( φ i ) is the cosine of the angle between the color gradient ~F and the latticedirection ~c i . Note that f k,eqi ( ρ k ,
0) stands for the set of equilibrium distributions of k − th fluid evaluated settingthe macroscopic velocity to zero. In the above equation, f i = P k f ki . The LB color gradient model has beenenriched with the so called regularization procedure [19, 20, ? ], namely a discrete Hermite projection of thepost-collisional set of distribution functions onto a proper set of Hermite basis. The main idea is to introduce6 set of pre-collision distribution functions which are defined only in terms of the macroscopic hydrodynamicmoments.All the higher-order non-equilibrium information, often referred to as ghosts [12], is discarded.Inequations, the regularized LB reads as follows: f ki ( x i + c i ∆ t, t + ∆ t ) = R f k,neqi ( x, t ) ≡ f k,eqi − ∆ tω k ( f k,regi − f k,eqi ) (11)where f k,regi is the hydrodynamic component of the full distribution f ki (see [20]) for the k − th color, and R is the regularization operator. The above equation shows that the post-collision distribution, of a 4 th -orderisotropic lattice, is defined only in terms of the conserved and the transport hydrodynamic modes, namelydensity ρ , current ρ~u and momentum-flux tensor Π = P i f i ~c i ~c i . References [1] S. Sugiura, M. Nakajima, N. Kumazawa, S. Iwamoto, M. Seki, Characterization of spontaneoustransformation-based droplet formation during microchannel emulsification, The Journal of Physical Chem-istry B 106 (36) (2002) 9405–9409. doi:10.1021/jp0259871 .URL https://doi.org/10.1021/jp0259871 [2] C. Priest, S. Herminghaus, R. Seemann, Generation of monodisperse gel emulsions in a microfluidic device,Applied Physics Letters 88 (2) (2006) 024106. doi:10.1063/1.2164393 .URL https://doi.org/10.1063/1.2164393 [3] P. Garstecki, H. A. Stone, G. M. Whitesides, Mechanism for flow-rate controlled breakup in confinedgeometries: A route to monodisperse emulsions, Physical review letters 94 (16) (2005) 164501.[4] P. Garstecki, M. J. Fuerstman, H. A. Stone, G. M. Whitesides, Formation of droplets and bubbles in amicrofluidic t-junction-scaling and mechanism of break-up, Lab on a Chip 6 (3) (2006) 437–446.[5] M. Costantini, C. Colosi, P. Mozetic, J. Jaroszewicz, A. Tosato, A. Rainer, M. Trombetta, W. Święszkowski,M. Dentini, A. Barbetta, Correlation between porous texture and cell seeding efficiency of gas foaming andmicrofluidic foaming scaffolds, Materials Science and Engineering: C 62 (2016) 668–677.[6] S. L. Anna, N. Bontoux, H. A. Stone, Formation of dispersions using "flow focusing" in microchannels,Applied physics letters 82 (3) (2003) 364–366.[7] R. Dangla, E. Fradet, Y. Lopez, C. N. Baroud, The physical mechanisms of step emulsification, Journal ofPhysics D: Applied Physics 46 (11) (2013) 114003.[8] E. Stolovicki, R. Ziblat, D. A. Weitz, Throughput enhancement of parallel step emulsifier devices by shear-free and efficient nozzle clearance, Lab on a Chip (2018). doi:10.1039/c7lc01037k .URL http://dx.doi.org/10.1039/C7LC01037K [9] A. Ofner, D. G. Moore, P. A. RÃŒhs, P. Schwendimann, M. Eggersdorfer, E. Amstad, D. A. Weitz,A. R. Studart, High-throughput step emulsification for the production of functional materials using aglass microfluidic device, Macromolecular Chemistry and Physics 218 (2) (2017) 1600472–n/a, 1600472. doi:10.1002/macp.201600472 .URL http://dx.doi.org/10.1002/macp.201600472 [10] N. Mittal, C. Cohen, J. Bibette, N. Bremond, Dynamics of step-emulsification: From a single to a collectionof emulsion droplet generators, Physics of Fluids 26 (8) (2014) 082109. doi:10.1063/1.4892949 .URL https://doi.org/10.1063/1.4892949 [11] A. Montessori, M. Lauricella, M. La Rocca, S. Succi, E. Stolovicki, R. Ziblat, D. Weitz, Regularized latticeboltzmann multicomponent models for low capillary and reynolds microfluidics flows, Computers & Fluids167 (2018) 33 – 39. doi:https://doi.org/10.1016/j.compfluid.2018.02.029 .URL [12] F. Higuera, S. Succi, R. Benzi, Lattice gas dynamics with enhanced collisions, EPL (Europhysics Letters)9 (4) (1989) 345.[13] A. Montessori, P. Prestininzi, M. La Rocca, S. Succi, Lattice boltzmann approach for complex nonequilib-rium flows, Physical Review E 92 (4) (2015) 043308.[14] R. Benzi, M. Sbragaglia, S. Succi, M. Bernaschi, S. Chibbaro, Mesoscopic lattice boltzmann modeling ofsoft-glassy systems: theory and simulations, The Journal of Chemical Physics 131 (10) (2009) 104903.715] S. Succi, The Lattice Boltzmann Equation: For Complex States of Flowing Matter, Oxford UniversityPress, 2018.[16] S. Leclaire, M. Reggio, J.-Y. Trépanier, Numerical evaluation of two recoloring operators for an immiscibletwo-phase flow lattice boltzmann model, Applied Mathematical Modelling 36 (5) (2012) 2237–2252.[17] T. Reis, T. Phillips, Lattice boltzmann model for simulating immiscible two-phase flows, Journal of PhysicsA: Mathematical and Theoretical 40 (14) (2007) 4033.[18] S. Leclaire, M. Reggio, J.-Y. Trépanier, Isotropic color gradient for simulating very high-density ratios witha two-phase flow lattice boltzmann model, Computers & Fluids 48 (1) (2011) 98–112.[19] A. Montessori, G. Falcucci, P. Prestininzi, M. La Rocca, S. Succi, Regularized lattice bhatnagar-gross-krookmodel for two-and three-dimensional cavity flow simulations, Physical Review E 89 (5) (2014) 053317.[20] R. Zhang, X. Shan, H. Chen, Efficient kinetic method for fluid simulation beyond the navier-stokes equation,Physical Review E 74 (4) (2006) 046703.[21] Z. Li, A. Leshansky, L. Pismen, P. Tabeling, Step-emulsification in a microfluidic device, Lab on a Chip15 (4) (2015) 1023–1031.[22] I. Kobayashi, S. Mukataka, M. Nakajima, Effects of type and physical properties of oil phase on oil-in-wateremulsion droplet formation in straight-through microchannel emulsification, experimental and cfd studies,Langmuir 21 (13) (2005) 5722–5730.[23] M. Stoffel, S. Wahl, E. Lorenceau, R. Höhler, B. Mercier, D. E. Angelescu, Bubble production mechanismin a microfluidic foam generator, Physical review letters 108 (19) (2012) 198302.[24] K. van Dijke, I. Kobayashi, K. Schroën, K. Uemura, M. Nakajima, R. Boom, Effect of viscosities ofdispersed and continuous phases in microchannel oil-in-water emulsification, Microfluidics and nanofluidics9 (1) (2010) 77–85.[25] G. T. Vladisavljević, I. Kobayashi, M. Nakajima, Effect of dispersed phase viscosity on maximum dropletgeneration frequency in microchannel emulsification using asymmetric straight-through channels, Microflu-idics and Nanofluidics 10 (6) (2011) 1199–1209.[26] J. Eggers, Universal pinching of 3d axisymmetric free-surface flow, Phys. Rev. Lett. 71 (1993) 3458–3460. doi:10.1103/PhysRevLett.71.3458 .URL https://link.aps.org/doi/10.1103/PhysRevLett.71.3458 [27] V. van Steijn, C. R. Kleijn, M. T. Kreutzer, Flows around confined bubbles and their importance intriggering pinch-off, Physical review letters 103 (21) (2009) 214501.[28] J. Eggers, E. Villermaux, Physics of liquid jets, Reports on Progress in Physics 71 (3) (2008) 036601.URL http://stacks.iop.org/0034-4885/71/i=3/a=036601http://stacks.iop.org/0034-4885/71/i=3/a=036601