Microfluidic devices for generating droplets have found a huge range of applications in recent years including monitoring of chemical reaction using ultra-fast mixing [1, 2], the formation of nanoparticles [3, 4], sample delivery at X-ray Free Electron Lasers (XFELs) [5–7] and single-molecule detection [8], and for drug delivery [9]. Understanding and modelling droplet generation in these devices is key to controlling droplet volume, frequency, and monodispersity [10]. For example, the application of droplet-based microfluidics for XFEL sample delivery necessitates a highly monodisperse distribution of droplets with a well-characterised and tuneable frequency [11, 12]. It is also often desirable to minimise the size of individual droplets and to be able to reliably deliver them to precise locations at high frequencies. In this paper we present a systematic simulation study exploring the upper limit of droplet generation in terms of frequency/flow rate whilst maintaining a small droplet volume—particularly relevant to XFEL sample delivery [13, 14].

Microfluidic droplet generators can be categorised as either passive or active. Passive droplet generators exploit the flow characteristics and geometry to induce breakup of the liquid stream whilst active droplet-based microfluidics external sources of energy such as acoustic waves or energy and magnetic or electric fields [15]. The simplicity and comparatively low-cost of passive microfluidic droplet generators make them highly attractive in cases where significant modifications to the droplet characteristics are not required ‘on-the-fly’. Of the wide range of passive droplet generators that have been developed, these can be further classified as cross-flow (or “T-junction”), co-flow, or flow-focusing devices [16]. To maximise control over the droplet size and monodispersity we have chosen to simulate droplet generators of the flow-focusing category [17]. Compared to T-junction and co-flow droplet generators, flow focusing devices provide a wider range of parameters to control droplet characteristics and thus create opportunities to tailor the output without, for example, modifying the wall wettability [18].

The basic elements of the microfluidic T-junction droplet generator are that a buffer/continuous phase flows within a channel and receives the droplet forming sample/continuous phase injection from a side channel intersecting at 90°. These devices have a tendency to form ‘plug-like’ droplets whose dimensions scale in proportion to the flow rate and size of the channels. The advantages of T-junction droplet generators are ease of fabrication and reliability; however, it is difficult to achieve very small droplet sizes using these devices due to the practical limits on channel size. The high resistance associated with having small channels also creates issues in terms of blockage of those channels. To address this issue it has been suggested the channel walls can be chemically modified to reduce their hydrophobicity and facilitate high-flow rates [17].

Coflowing microfluidic droplet generators meanwhile consist of two coaxial channels. The coaxial geometry within these devices is often formed by a capillary tube being placed inside of a rectangular channel. The capillary then supplies the dispersed phase/sample whilst the surrounding rectangular tube contains the continuous phase fluid. The flow rate of the surrounding sheath fluid can be adjusted to control the breakup of the dispersed phase fluid [17]. Related to this, a recent study [19] found that the capillary number of the continuous phase in coflowing droplet generators plays a key role in controlling the breakup regime for droplet formation. Whilst coflowing devices offer greater flexibility, for example, in controlling breakup versus surface tension and Rayleigh instability, the added complexity and restrictions in terms of device geometry can potentially limit their applicability in the context of sample delivery.

Flow-focusing droplet generators comprise a junction in which the dispersed phase is hydrodynamically focused by the continuous phase [17]. The geometry of these devices can be either 3-D axisymmetric, or quasi-2D planar [15]. Axisymmetric devices provide a dispersed phase through a cylindrical inlet, which is hydrodynamically focussed as both phases pass through a circular baffle. Alternatively, planar geometries consist of one dispersed phase inlet surrounded by two continuous phase inlets, contained within a maze-like structure of constant height that focuses the dispersed phase using the continuous phase fluid [20]. Flow focusing devices can be manufactured to high spatial resolution via lithography, allowing them to achieve small feature sizes. Flow-focusing droplet generation can be primarily inertial, for example, when the outer channels are directed backwards towards the fluid intersection, or geometric, as in the case of ‘geometry-controlled’ droplet generation. The combination of inertial and geometric effects are a key feature of flow-focusing devices and the mechanism for droplet generation critically depends on parameters such as the capillary number, Weber number, relative viscosity of the dispersed and continuous phase as well as the flow rate.

Flow focusing devices vary widely in design, from typical planar flow focusing with a straight neck [21–25], through variations on a simple rectangular cross motif [26, 27], to coaxial flow focusing [28], and the closely related coflowing droplet generators [29]. The channel diameters tend to vary from 10’s of µm [26] to the upper limit of microfluidics of hundreds of µm [30]. In the literature, flow rate is regularly the independent variable when studying flow regimes and mechanics. Flow rates can range from single digit µL/min [23], into the 10’s like most studies [20], or even many orders higher [29].

Previous studies have examined both the potential applications and mechanics of microfluidic flow focusing droplet generators [26, 28, 31, 32]. The high flow rates and range of surface tensions explored in the present study are rare, and even less commonly achieve comparable droplet frequencies [33]. This study shows that these conditions are conducive to highly frequent and monodisperse droplets. It is envisaged that the results from this study will translate into the design optimisation of microfluidic devices for droplet-on-demand sample delivery and rapid mixing [34, 35].

Passive microfluidic flow-focusing droplet generators can be classified in a number of different ways; most commonly, however, they are categorised according to the capillary number, Ca, of the continuous phase. This dimensionless number describes the ratio of viscous forces to interfacial tension forces that contribute to the fluid flow features C a = µ V σ (1) where µ is dynamic viscosity, V is velocity, and σ is surface tension.

The Weber number, We, is also used to categorise flow-focusing microfluidic devices in the literature: W e = ρ L V 2 σ (2) where ρ is density and L is the characteristic length, or hydraulic diameter of the channel.

However, this is less commonly used then the capillary number since We describes the ratio of inertial forces to interfacial tension forces in a flow and inertia has typically been assumed to be of secondary importance in flow-focusing microfluidics [9]. The common neglect of the effects of inertia in flow-focusing devices can be attributed to the fact that the Reynolds number, Re, R e = ρ L V µ = W e C a (3) which represents the ratio of inertial to viscous forces, is much less than unity in these devices. However, we will observe that under the specific conditions simulated here the effects of inertia are, in fact critical leading to higher values of Re than are typically reported for flow-focusing devices. Hence, in the present case, We was selected as the primary variable since the key characteristics of the devices e.g. droplet diameter and droplet dispersity show the most important variations when plotted as a function of We.

The features of droplet generating microfluidic devices which are typically of most interest are the jet length, droplet breakup position, droplet size, droplet frequency and droplet size uniformity (monodispersity). Experimentally, jet length, breakup position, and droplet size are typically determined from images of the device whilst in operation. Whilst jet length and breakup position are typically defined as simple distance measurements from the point of intersection of the continuous and dispersed phase, droplet size can be characterised in multiple ways, for example, droplet diameter, volume, and diameter coefficient of variation (CV). Because these parameters are normally determined from 2D images of the devices, it can often be hard in experiment to measure the precise droplet volume, however in simulation this quantity can be accurately determined to the mesh resolution of the simulation. CV is defined as: CV(%) = S D D ¯ × 100 , where S D is the standard deviation and D ¯ the mean droplet diameter.

Flow-focusing microfluidic devices have historically been classified as corresponding to one of three regimes: squeezing, dripping, or jetting, depending on the corresponding Ca for the device. More recently however, two additional regimes have been identified, classified as ‘tip-streaming’ and ‘tip-multi-breaking’ [36, 37]. These two additional regimes have been observed when moving to higher values of Ca and are distinguished by the relative influence of viscosity or inertia on droplet generation associated with these devices.

Within the squeezing regime droplets are generated periodically. The process is initiated by the orifice being filled by the dispersed phase. Next, the continuous phase pressure builds up and “squeezes” the dispersed phase into a droplet. The droplet is pushed out of the channel by the continuous phase, building up the pressure of the dispersed phase. Finally, the dispersed phase blocks the orifice again, and the cycle repeats. The upper Ca limit for the squeezing regime has been said to vary between 3 × 10⁻³ and one depending on the channel size, flow rate ratio, and viscosity [38], but is typically on the order of Ca ∼10^–2 [15, 26, 31].

Once the viscosity is sufficiently large that it can overcome surface tension, breakup of the dispersed phase transitions to the dripping regime [15, 19, 28, 38]. In dripping, Ca is typically larger than it is for the squeezing regime, starting at 10^–2 [15, 19, 38]. Within the dripping regime drag forces are often the dominant cause of droplet formation and droplet diameters are often reported as being highly monodisperse [9, 15, 19, 38].

As the capillary number is further increased, the jetting regime is reached. In this regime, droplet size grows and with it the polydispersity of the droplet distribution. The main cause of droplet generation in the jetting regime is typically described as Rayleigh-Plateau instabilities [15, 17]. The capillary number associated with the continuous phase, Ca_c, of jetting is limited to less than ∼0.5. The corresponding onset length of this regime is purported to be equivalent to at least three times the orifice diameter [19].

Whilst jetting was traditionally viewed as the final regime for flow-focusing droplet generation, more recently tip-streaming has been classified as a separate regime [36] occurring when a thin jet under shear or elongational flow, breaks up into extremely small (sub-micron) droplets mediated by Rayleigh-Plateau instabilities. To achieve this the capillary number should be at least 0.5 with a corresponding viscosity ratio, λ = µ_d/µ_c, of at least 0.1, where µ_d and µ_c are the viscosities of the dispersed and continuous phases respectively. To reach this regime, surfactant is typically employed, with the amount of surfactant to be added being inversely proportional to the interfacial surface tension and proportional to droplet diameter [15]. Additionally, the Re of both the dispersed and continuous phase in this regime is << 1. In this regime viscous forces dominate droplet formation, and the corresponding distribution of droplets has been reported as being highly monodisperse.

Even more recently [37] a fifth regime labelled the “tip-multi-breaking” regime has been identified. Unlike the tip-streaming regime, here Ca is normally less than for either tip-streaming or jetting but is greater than for the squeezing and dripping regimes. In the tip-multi-breaking regime, the viscosity ratio is still small, but the breakup is dominated by inertia. Cyclical oscillation of the dispersed liquid causes sequentially sized droplets to be released in this regime.

The systematic variation of surface tension variation is routinely investigated by starting with a pair of immiscible liquids, such as water and oil, which have a high interfacial surface tension, and then decreasing the tension via the addition of surfactant. The combination of the presence of surfactants combined with microfluidics can introduce dynamic variations in the local surfactant kinetics on the boundary of the droplets. The kinetics depend on the local shear forces, surfactant concentration gradients, and the difference between equilibrium surfactant distribution time and droplet breakup time [39–41]. These effects may be approximated by utilising velocity dependant and spatially dependant surface tension values for the droplets during simulation [19, 42, 43].

Largely due to the limits of computational resources most of the early simulation studies of droplet forming microfluidic devices approximated the three-dimensional fluid flow as a two-dimensional process. Continuous advancement in computing capabilities and the sophistication of simulations now provide the opportunity to explore, in a 3D simulated environment, the effect of varying parameters such as surfactant concentration and flow rate. The results of these studies can in turn provide general information allowing for the optimised fabrication and operation of droplet forming microfluidic devices. In simulation, detailed quantitative analysis can also be performed and more extreme situations explored to better understand the practical limits of empirically derived laws governing droplet formation.

Here we carry out a series of three-dimensional simulations of flow-focusing microfluidics with a continuous phase composed of oil and a dispersed phase composed of water. We explore the high-frequency droplet generation regime in the context of varying interfacial surface tension and geometry. In comparison to earlier simulation studies which have explored the effect of surface tension on droplet formation—the dimensions of our devices are small—matching current capabilities for microfluidic fabrication, the flow rates are higher, and the range of interfacial surface tensions are larger than has been looked at previously [33]. In particular, the flow rates were chosen to facilitate Gas Dynamic Virtual Nozzle (GDVN) jetting conditions [14]. It is envisaged that the results will translate into design optimisation of microfluidic devices for sample delivery and rapid mixing.

Simulations were performed using COMSOL Multiphysics software, versions 5.6 and 6.0. The behaviour of the two-phase fluids was studied using the Two-Phase Flow Multiphysics module, which couples the Laminar Flow module responsible for the fluid flow conditions and Level Set module describing the moving interface. The Laminar Flow module was selected to model the constraints present in thin microfluidic systems, which often have low Reynolds numbers due to their small characteristic length scales [17].

The COMSOL Multiphysics FEM software package provides for both adaptive and static mesh simulation approaches. The former yields infinitesimal liquid-liquid interface thickness results by placing mesh element points at the interface. This interface calculation is then repeated for each time point. To allow for dispersed phase separation, a static mesh model, such as ‘Level Set’ or ‘Phase Field’, is required. Considering that the “Level Set” model is the more computationally efficient of the two, it is the more commonly employed in order to minimise calculation time.

The governing equations used by the Level Set model to represent a two fluid water and oil system are the incompressible Navier-Stokes equations. Eq. 4 is called the continuity equation and it ensures fluid incompressibility by requiring that fluid cannot converge to, or diverge from, any single point according to conservation of mass, ∇ ∙ u = 0 , (4) where u is the velocity vector. Eq. 5 is the momentum equation, which fulfills the requirement for changes in fluid flow to result from pressure and other forces whilst maintaining conservation of momentum, ρ δ u δ t + u ∙ ∇ u = ∇ ∙ − p I + μ ∇ u + ∇ u T + F g + F s t + F e x t + F , (5) where ρ, t, p, I, μ, F_g , F_st , F_ext , and F are density, time, pressure, the identity matrix, dynamic viscosity, the force of gravity, the surface tension force, the external force, and any applicable additional forces, respectively. Note that only the Phase Field model uses F_ext . F_g is equal to ρ g where g is the gravity vector. Since there are no additional forces, F is not applicable here. F _st is defined by Eq. 6: F s t = σ δ κ n + δ ∇ s σ , (6) where σ is surface tension, δ is the Dirac-delta function at the interface, ∇ _s is the surface gradient operator, 𝜅 is the curvature of the interface, and n is the vector orthonormal to the interface. It can be seen from Eq. 6 that F_st comprises contributions from the continuum method (57), derived from the Young-Laplace equation, and the single-phase Marangoni force. Eq. 7 is added into the self-consistent equations by the Level Set model for the purpose of tracking the fluid interface: δ Φ δ t + u ∙ ∇ Φ = γ ∇ ∙ ε ∇ ϕ − ϕ 1 − ϕ ∇ ϕ ∇ ϕ , (7) where γ is the reinitialization parameter and ε is the “parameter controlling interface thickness”. The interface controlling parameter is not definitively equal to the thickness of the interface, it only controls a function that determines the thickness. The default value for the parameter is half of the maximum element size within the region intersecting the interface.

Simulations were performed using COMSOL Multiphysics software, versions 5.6 and 6.0. A square channel was simulated with symmetry constraints applied along the y and z axes. This reduced the simulation time by 75% whilst maintaining all essential features of the device. Typical simulation times were around 30 h for some of the devices simulated within this study which comprised around 420,000 elements. Some simulations, e.g., those with especially small neck widths below 25 µm could take up to 5 days to complete due to the increase in the density of mesh elements. The computer used for simulations is equipped with two intel Xeon e5-2687w v2, 8 core, 3.40 GHz CPUs, 128 GB RAM and the simulations each occupied approximately 20 GB of storage space.

All of the results presented here are from simulations which converged to a stable solution for droplet generation. The mesh resolution was optimised and validated to confirm that using a finer mesh did not result in variations in the output of the simulation model. Simulations were carried out assuming a constant slip length. Some studies reported in the literature apply a ‘no slip’ condition, whilst this results in a computationally more efficient model, it comes at the cost of accuracy and therefore was not adopted in the present case. A schematic illustrating the representative geometry used for this study is shown in Figure 1. The channels used throughout the device were square in order to mimic the two most common fabrication geometries encountered for flow-focusing devices, axisymmetric, and planar. A key feature of the flow-focusing devices simulated here are that they include a flow baffle in which the thin neck restricts the flow of the water/oil mixture. The incorporation of a comparatively long outlet channel allowed for the study of a wide range of breakup lengths.

Figure 1A shows the water and oil inlets; water was injected at a volumetric flow rate of Q_d = 20 μL/min whilst oil was injected at a volumetric flow rate of Q_c = 100 μL/min. The resulting flow rate ratio was thus: φ = Q_d/Q_c = 1/5. All simulations were carried out under standard temperature and pressure (STP). For water the viscosity and density were µ_d = 1.0 mPa•s and ρ_d = 1.0 g/cm³ respectively, for oil they were µ_c = 4.1 mPa•s and 1.78 g/cm³, yielding a constant viscosity ratio of λ = 0.24. The main parameters using in the simulation are shown in Table 1.

Solid-liquid interfaces are subject to flow restrictions from friction, approximated by the slip condition, and capillary forces which can be determined from the contact angle. The contact angle used here is 135, corresponding to most common materials used in microfluidics (e.g., SU-8) and the slip condition is Navier Slip, or partial slip, with a slip length factor of 1. It is also worth noting that the distance from the channel wall to the point of zero velocity, defined as the “slip distance,” in COMSOL is set as a constant equal to the constant “slip length factor” multiplied by the smallest mesh element in the simulation. Time-dependant simulations under these conditions were carried out with a time-increment of 10 µs. The total simulation time was 10 ms, during which the droplet generation was observed to have stabilised. Here, a constant viscosity is used whilst the interfacial surface tension, σ, is varied and evaluated dynamically during the simulation as a function of position in the x-direction: σ = σ 0 1 − 0.2 ∙ cos 2 π x L n (8) where L n is the length of the neck and σ₀ is the initial surface tension.

The rationale behind this selection is three-fold. First, the surface tension is a variable, and the systems being simulated match the criteria for surfactant laden droplet generation, in which surfactant, and hence surface tension, vary both spatially and temporally along the length of the neck. Secondly, the specific function used was selected for its computational simplicity, it can be readily interpreted as a spatial approximation of normal dynamic surface tension [43]. Finally, because the function describing varying surface tension doesn’t involve any feedback from the fluid flow results, it can be evaluated and updated at each and every cycle without significantly impacting calculation speed.

Figure 2 shows representative isosurface renderings of droplet generation for a range of different initial surface tensions and neck sizes. These results were used to quantify the droplet onset length and volume.

Using finite element simulation, a range of interfacial surface tensions and geometries have been explored for flow-focusing droplet forming microfluidic devices. In general, the breakup onset lengths observed here are typically longer than previously observed for flow-focusing droplet generators [21, 26, 27] which normally use wider channels. In addition, particularly when using low surface tension values and small neck widths, the droplets formed had a range of non-spherical shapes. The confinement of droplets in in the long, very narrow, outlet channel used here acts to inhibit capillary instabilities [46] such as Rayleigh-Plateau (RP) instabilities which have their growth suppressed with respect to distance but typically cause breakup in open channels. Thus, compared to previous studies using open, or at least, much wider outlet channels than are used here, the breakup length may be expected to be longer. Confinement within a narrow channel also influences the droplet shape. When the surface of the droplets is near a wall, the droplet surface velocity is lowered compared to the bulk of the droplet due to friction. This in turn can lead to the droplet shape becoming non-spherical and would be expected to be particularly apparent in droplets with low surface tension and for very narrow channels—consistent with our observations from simulation.

A further general result from the analysis was that the breakup onset lengths increase, and the maximum pressure decreases as the neck width increases. This is accompanied by an increase in the monodispersity and droplet frequency. Observations of droplets from the simulation suggest that the sudden increase in onset length which occurs for a neck width of 31.25 µm (corresponding to a large CV) appears to occur because of a transition between a shear-dominated jetting mode and a surface tension-dominated jetting mode. This sudden jump may also be influenced by the Re since it seems probable that Re would decrease at the neck section as neck widths increase, and the overall fluid behaviour of the device becomes more like that expected for tip streaming—associated with a larger onset length. Future work will examine this interesting transition between the different regimes in more detail via a finer sampling of channel widths.

Whilst there have been previous studies of droplet generation in flow-focusing devices [9, 21, 42, 47, 48], these have generally used a much lower flow rate to channel cross-section ratio. This results in much lower droplet velocities, which in turn, would make them unsuitable for many sample delivery applications where material needs to be replenished quickly. As a result, the Weber numbers and Reynolds numbers reported here are substantially higher than are typically reported in the literature for flow-focusing microfluidic devices.

The interfacial surface tension was varied in a simulated flow-focusing geometry with a square cross-section and analysis of the key features of the simulation was performed. Comparing the model characteristics with those reported in the literature for the different fluid flow regimes, it was determined that jetting and tip streaming, best describe the model behaviour. The focus of the present study was on high-frequency formation of small droplets with low polydispersity. This resulted in an unusually high ratio for the total flow rate to channel width. By probing the effects of interfacial surface tension over a much larger range then previously reported we found that, surprisingly, a standard water/oil mixture is the least suited to high frequency droplet generation. The highest frequency of droplets was found to be 56 kHz observed for a surface tension of 118.75 mN/m. This had low polydispersity (CV was just 3.1%) and a comparatively small droplet diameter of 22 µm. The geometrical effect of varying the neck width was also examined. By optimising the neck width, a value of 37.5 µm for the neck was observed to result in a droplet generation frequency of 61 kHz and a surface tension of 52 mN/m.

Having completed a comprehensive suite of simulations, the results of this study will now be used to guide future work confirming these observations in experiment. Previous studies investigating the influence of surfactants in inkjet printing have revealed that the distribution of surfactant at the droplet interface may be inhomogeneous [39]. Although the range of surface tensions investigated in this earlier study was different to the present study, changes in surfactant concentration before, during, and after droplet formation could potentially have a significant influence on droplet behaviour in the present case. An extrapolation of the results between the 37.6 and 25 µm necks for surface tension values between 52 and 118.75 mN/m, indicates that droplet frequencies approaching 90 kHz whilst maintaining small droplet size and low polydispersity, should be readily achievable, although it remains to be tested whether this could be replicated in experiment.