## Abstract

A model was developed for the analysis of steady and unsteady bubble generation frequencies in a microfluidic flow-focusing device. In both cases, the generation frequency depends on the downstream influence of bubbles created by the flow-focusing device, which provides a source of feedback. For steady-state generation of bubbles, this feedback explains the relationships among frequency, gas pressure and liquid flow rate in our experiments as well as gives a physical explanation for previously observed frequency relations. The role of feedback is also exploited to explain unsteady behaviour; we develop a numerical model and analytical expressions that agree with experimental measurements.

## 1. Introduction

A wide variety of microfluidic devices are available to produce well-controlled and monodisperse foams, emulsions and even multiple emulsions, with volume fractions of the dispersed phase varying from 0 to 1 (see Umbanhowar *et al*. 2000; Ganan-Calvo & Gordillo 2001; Anna *et al*. 2003; Garstecki *et al*. 2004; Okushima *et al*. 2004; Utada *et al*. 2005). The increased control of these flows comes at a cost—system volume throughputs are significantly smaller than with conventional devices. For many applications of microfluidics, such as those in biology and medicine, the decreased volume is beneficial, since it reduces the required amount of expensive protein or reagent (see the review by Beebe *et al*. 2002). To create industrially useful amounts of foam or emulsion, however, the volumetric flow rates produced by microfluidic devices must be increased.

One route to produce larger quantities from microfluidic devices is to run them in parallel. Designing massively parallel microfluidic channels is easy using current fabrication techniques, but the behaviour of microfluidic devices under highly parallel operation is poorly understood. Generating monodisperse foams and emulsions in parallel systems has proved challenging, in part due to the significant crosstalk and feedback in these systems (see, for instance, Barbier *et al*. 2006).

The primary focus of studies on the generation of multiphase flows is on the dynamics at the fluid–fluid interface. This viewpoint is understandable, as the mechanics of thread break up is an old and fundamentally interesting problem (see the review by Eggers 1997). However, understanding how these devices will operate in parallel requires a complete description of the device, from fluid inlets to multiphase outlet. This approach is especially true for parallel systems, where an increased resistance in one channel can affect the flow rates in all of the other channels.

We use the specific example of a microfluidic flow-focusing bubble generator to demonstrate the interrelation between the behaviour at the fluid–fluid interface and the downstream flow behaviour, including steady-state and time-dependent behaviours. We start by giving a brief review of the principles of multiphase flow generation in microfluidic devices, in particular the behaviour at the fluid–fluid interface. Next, we describe experiments performed to quantify the relationship between the generation frequency of individual drops and the input flow parameters in a flow-focusing device. Remarkably, the generation frequency even in this simple device depends both on the dynamics at the fluid–fluid interface and on the resulting multiphase flow generated from the device. A quantitative theory is developed for the bubbling frequency as a function of the flow and geometrical parameters.

## 2. Multiphase flow generation

The underlying principle in all microfluidic flow-focusing generators is the injection of one phase, the discrete phase, into a second, immiscible continuous phase. The two most common microfluidic implementations of this idea are the T-junction (see Thorsen *et al*. 2001), where the discrete phase is injected perpendicular to the continuous phase (figure 1*a*), and the flow-focusing junction (see Ganan-Calvo & Gordillo 2001; Anna *et al*. 2003), where the discrete phase is injected into a coflowing continuous phase (figure 1*b*).

The dynamics underlying both the T-junction and flow-focusing devices are controlled by competing stresses—surface tension, pressure, viscous and inertial—at the fluid–fluid interface. As the flow-focusing geometry lends itself to more straightforward experimental, numerical and theoretical description, we focus on the behaviour of these devices, although significant progress has also been made in understanding the operation of T-junction devices (for instance Garstecki *et al*. 2006).

Anna & Mayer (2006) described three distinct regimes in the formation of multiphase flows in planar flow-focusing devices, in the order of increasing capillary number: geometrically controlled break up, dripping, and jetting (figure 2). Garstecki *et al*. (2005) showed that controlled break up occurs when the thread of the discrete phase is comparable in diameter to the flow-focusing constriction. The thread grows until it blocks flow by the continuous phase, at which point the break up of the thread is mediated by volumetric contraction of the thread until pinch-off. The thread then retracts to its original position and the process repeats. Droplets formed in this regime are similar in size to the constriction and are highly monodisperse. Dripping is similar to geometrically controlled break up, except that the break up is not directly controlled by the geometry of the constriction, but instead driven by surface tension (see Eggers 1993).

Jetting is characterized by steady droplet formation and pinch-off from a thread that is significantly smaller and extends considerably further than the flow-focusing junction. Here, break up is controlled by a Rayleigh–Plateau-like capillary instability of the long, thin thread (see Plateau 1849; Lord Rayleigh 1879). Utada *et al*. (2007) found two different regimes of jetting, ‘narrowing jetting’ and ‘widening jetting’, characterized by the downstream width of the narrow thread. The narrowing jetting features a jet that narrows downstream of the constriction and the subsequent droplet sizes are comparable to the thread size. The widening jetting is characterized by an increasing width of the thread downstream from the flow-focusing junction, which produces relatively large, polydisperse droplets.

## 3. Steady bubble generation

### (a) Experimental set-up and measurements

Flow-focusing devices (figure 1*c*) are fabricated in polydimethylsiloxane (PDMS) using standard soft lithography techniques following Duffy *et al*. (1998). We use solutions of glycerol and water combined with Tween-20 surfactant (0.1% wt) as the continuous phase. The solution viscosity is adjusted by changing the relative concentrations of glycerol and water (table 1). The liquid flow rate is controlled using a syringe pump. We use compressed air that is controlled by a pressure regulator as the discrete phase.

Measurements are taken on an inverted microscope using a fast camera, with a typical frame rate of 4000 frames s^{−1}. Analysis of these movies is performed using image analysis software to first find the bubble size and position in each movie frame and then stitch positions together to produce complete tracks for each bubble. From the evolution of bubble position and size with time, all of the relevant parameters for the bubble are calculated. For example, the velocity is measured by using the numerical derivative of the bubble position with respect to time and the frequency is calculated from the time delay between successive bubbles to pass a fixed reference point in the channel. Finally, the volume is estimated as the product of the visible area of the bubble and the height of the channel, since the discrete phase is shaped similar to a disc.

### (b) Experimental results

The bubble velocity (figure 3) is observed to be approximately equal to the average fluid velocity. For large bubbles in a circular channel, Bretherton (1961) found that the discrepancy between the average fluid velocity and bubble velocity is expected to be proportional to *Ca*^{2/3}, where *Ca*=*ηu*/*γ* is the ratio of viscosity *η* and velocity *u* to surface tension *γ*. While the channel is not circular in cross section in our system, the bubbles are large enough to be geometrically confined along the vertical dimension and produce a thin wetting layer, so this analysis is not entirely inappropriate. A generalization to account for the rectangular cross section and the motion of large bubbles has been given by Wong *et al*. (1995), and confirms this idea.

The steady-state production of bubbles in such a device is characterized by two parameters—the volume of each bubble and the frequency of generation. The bubble volume has been studied experimentally by Garstecki *et al*. (2004) in a planar device and numerically by Jensen *et al*. (2006) for the case of a low Reynolds number axisymmetric flow. In both of these studies, volume was found to be proportional to gas pressure *P* and inversely proportional to liquid flow rate *q*_{ℓ} and viscosity *η*, *V*_{b}∝*P*/*ηq*_{ℓ}. This relation is explained in terms of geometrically controlled pinch-off. The pinch-off time is inversely proportional to liquid flow rate, and the gas flow rate *q*_{g} is set by viscous drag on the gas thread to be *q*_{g}∝*P*/*η* due to pushing forward the downstream liquid column of viscosity *η*. Our measurements of bubble volume confirm the scaling of this relationship, as shown in figure 3*b*.

Garstecki *et al*. (2004) found that frequency of bubble formation *f* was directly proportional to the flow rate and gas pressure *f*∝*Pq*_{ℓ}, but our measurements (varying the flow rate with a fixed gas pressure), reported in figure 4*a*, clearly deviate from this linear relationship at high flow rates. This discrepancy can be understood in terms of the extra pressure drop produced in the channel by the flow of liquid through the channel *P*_{c} as well as the pressure drop induced by the bubbles in the channel (see figure 5 for a schematic of the notation used). For a rectangular channel with cross-sectional dimensions *h* and *w* and length *L*, the pressure drop for a fluid of viscosity *η* and flow rate *q*_{ℓ} is(3.1)

For simplicity, we rewrite this equation as *P*_{c}=*αq*_{ℓ}, where *α* is the hydrodynamic resistance. Raven *et al*. (2006) extended the frequency relationship to include the channel pressure drop *P*_{c} with *f*∝(*P*−*P*_{c})*q*_{ℓ} for low pressures and *f*∝1/*q*_{ℓ} for high pressures. This modified relation arises from the two time scales present in bubble formation: the thread pinch-off time, which is proportional to 1/*q*_{ℓ}; and the thread growth time, which is observed experimentally to be proportional to 1/(*q*_{ℓ}(*P*−*P*_{c})). Direct measurements of the evolution of gas thread volume (figure 6) confirm the idea that there are two distinct time scales in the thread pinch-off dynamics in our system as well (labelled as a, b in figure 6). There is initially slow volume growth, corresponding to the thread growth up to the flow-focusing constriction followed by a rapid volume increase corresponding to bubble formation and pinch-off in the constriction. In our experiments, the thread growth time is significantly longer than the bubble formation and pinch-off and dominates in determining the frequency of bubble generation. The frequency relation *f*∝*q*_{ℓ}(*P*−*P*_{c}) explains the observed variation for different viscosities and flow rates, but does not capture the dependence on gas inlet pressure *P* (figure 4*b*).

### (c) A model for bubbling frequency

To move beyond an empirical model for the bubble generation frequency, we must understand what controls the time scale of thread growth. We can model the initial thread growth by a column of gas at fixed pressure *P* growing into a fluid-filled tube. Neglecting the effects of surface tension and interface shape, a simple force balance sets the flow rate of gas *q*_{g} to be *q*_{g}=*P*/*α*, where *α* is the hydrodynamic resistance of the liquid through the channel. If, in addition to the gas flow, the liquid flows with a volumetric rate *q*_{ℓ}, the gas flow rate is reduced to *q*_{g}=*P*/*α*−*q*_{ℓ}=(*P*−*P*_{c})/*α*, giving a total flow rate *q*=*q*_{g}+*q*_{ℓ}. For the flow-focusing device, the gas can only fill a finite volume, *V*_{f}, before the effects of the flow-focusing junction dominate, which sets a time scale for slow growth of *τ*_{g}=*V*_{f}/*q*_{g}. When thread growth dominates bubble generation, this estimate gives a frequency generation of *f*=(*P*−*P*_{c})/*αV*_{f}.

This simple model is in disagreement with the observed scaling of frequency, which depends on both pressure drop and flow rate, and not simply pressure drop (figure 4*b*; Garstecki *et al*. 2004; Raven *et al*. 2006). This discrepancy can be explained by realizing that the observed experimental frequencies correspond to steady bubble trains, not simply fluid-filled channels. As has been observed by Raven & Marmottant (2006) in a foam generator, the bubbles downstream also contribute to the pressure drop in the channel and can dramatically alter the flow behaviour. The total pressure drop Δ*P* for a steady flow in a channel of length *L* may then be written as , where Δ*P*_{b,i} is the additional pressure drop of the *i*th bubble in the channel. This extension to the model allows the generation frequency to be calculated as(3.2)

If the stream is assumed to be in a single period, steady-state mode, the sum in equation (3.2) can be replaced by *N*Δ*P*_{b}, where each of the *N* bubbles contributes the same pressure drop Δ*P*_{b}. The bubble velocity is proportional to the average fluid velocity, so for a channel of height *h*, width *w* and length *L*, the total number of bubbles, *N*, in the channel is *N*=*fLhw*/*q*_{ℓ}. Combining the modified pressure drop from a stream of steady bubbles with the simple model for frequency yields a relation for the steady-state generation frequency(3.3)where *V*_{c}=*Lhw* is the volume of the channel. In the limit that *q*_{ℓ}≪*V*_{c}Δ*P*_{b}/*αV*_{t}, equation (3.3) can be well approximated by *f*_{0}=*q*_{ℓ}(*P*−*P*_{c})/*V*_{c}Δ*P*_{b}.

A general relation for the pressure drop across a bubble Δ*P*_{b}, for a bubble comparable in size to the channel width in a rectangular channel is not known to us. Feurstman *et al*. (2007) determined a relation for Δ*P*_{b} for bubbles long compared to the channel width in rectangular channels that depends on the bubble length, liquid flow rate and surface tension. Adzima & Velankar (2006) likewise found a relation for liquid droplets in a liquid–liquid flow. The pressure drop of the bubble should depend on size, velocity, viscosity and surface tension of the bubble and viscosity and flow rate of the fluid. As the viscosity of air is small enough to neglect and the surface tension is similar for all of our experiments, we expect the observed Δ*P*_{b} behaviour to be dominated by bubble velocity and size effects. We do not have a direct measurement of the additional pressure drop due to a bubble, but we can invert equation (3.3) to arrive at an estimate of Δ*P*_{b}. This step yields, in the limit of small *q*_{ℓ}, Δ*P*_{b}=*q*_{ℓ}(*P*−*P*_{c})/*f*_{0}*V*_{c}, giving Δ*P*_{b} in terms of properties that have been experimentally measured. We expect the pressure drop to be proportional to the viscosity and liquid flow rate, which is itself proportional to *P*_{c}, so a plot of our experimental Δ*P*_{b} estimate normalized by *P*_{c} should yield the bubble volume dependence of Δ*P*_{b} (figure 7). The data are approximately linear with *V*_{b}, which allows us to make the empirical approximation Δ*P*_{b}=*ψP*_{c}*V*_{b}/*V*_{c}, where the dimensionless proportionality constant *ψ* depends on surface tension and channel geometry.

The final piece to clarifying the scaling relation for frequency is the relation for volume scaling found by Garstecki *et al*. (2004), *V*_{b}=*V*_{t}*P*/*P*_{c}, where *V*_{t} is the proportionality constant in the volume relation (figure 3*b*). Thus, the pressure drop can be rewritten as Δ*P*_{b}=*ψPV*_{t}/*V*_{c}. The final approximate relation for frequency becomes(3.4)

This relation agrees well with our measurements over a range of pressures, flow rates and viscosities (figure 7).

## 4. Unsteady bubble generation

One experimental test of the role of feedback in the bubble generation frequency is to measure this frequency with varying numbers of bubbles in the outlet channel. This test can be practically accomplished by measuring the start-up behaviour of an initially bubble-free channel. To this end, measurements of a start-up period were taken by raising the liquid velocity above the rate at which the gas thread retreats from the outlet and, once the device is cleared of bubbles, lowering the flow rate to the desired value.

The start-up behaviour of a flow-focusing device from such an experiment is presented in figure 8*a*. The frequency of bubbles during start up oscillates from an initially high value and then slowly settles towards a stable generation frequency. One unexpected feature of the measurement is an initial increase in the bubble generation frequency, as one might have expected bubbles to increase the flow-generated pressure and thus decrease the frequency. This feature can be understood by looking at the average bubble velocity during generation (figure 8*a*), where it is seen that the liquid flow rate decreases slowly from an initially large value. This fluctuation is due to compressibility in the syringe—most likely owing to trapped small air bubbles.

A numerical model can be developed to predict the time behaviour of a flow-focusing junction attached to a single outlet channel. The model is based on generating bubbles at a frequency dictated by equation (3.2). These bubbles are then carried by the flow at the mean flow velocity *u* until they reach the end of the channel. With the positions of all of the bubbles in the channel, the total pressure drop can be calculated to yield the new generation frequency. This process is then iterated to determine the temporal behaviour of the system. When taking into account the variable flow rate caused by compressibility of the system, the model does a reasonable job of replicating the observed oscillations (figure 8*a*).

The numerical model lends itself to an analytical approximation. When the number of bubbles in the channel is large, it is reasonable to approximate the discrete bubbles by a continuous density of bubbles *ρ*(*x*), i.e. the number per length. For bubbles moving at a constant velocity *u*, there is a simple relation between a continuous frequency *f*(*t*) and density, namely *ρ*(*x*)=*f*(*t*−*x*/*u*)/*u*. Similarly, the total pressure drop due to bubbles in the channel of length *L* can be calculated as . Combining these elements, we arrive at a continuous equation for the generation frequency(4.1)

In addition to the steady-state solution, one analytical solution to this integral equation is an offset decaying exponential *f*=*f*_{0}+*A*e^{(i−Δ)ωt}. Substituting this function in equation (4.1), one finds that *f*_{0} is the steady-state solution for frequency, given in equation (3.4). There is also a particular set of solutions for *ω* and *Δ* given by(4.2)(4.3)where *t*_{0}=*L*/*u* is the bubble travel time through the channel and *β*=Δ*P*_{b}/*αV*_{t}. The amplitude *A* of frequency oscillations is a free parameter and must be set by the initial conditions.

The frequencies that satisfy equation (4.3) form a family of solutions to the bubble oscillation problem. The smallest exponential decay—corresponding to the longest lived solution—occurs for the smallest non-zero value of *ω*. Comparing this solution with the numerical solution for bubble start up shows a strong agreement after the initial decay (figure 8*b*). The initial discrepancy is most likely due to the initially empty starting conditions in the numerical simulation. The agreement between the numerical and the analytical models exists until the number of bubbles in the channel becomes small, at which point the underlying assumption that the bubbles can be treated continuously is no longer valid.

## 5. Conclusions

We have furthered the understanding of the generation frequency behaviour in a microfluidic flow-focusing bubble generator. Frequency is controlled by two time scales—the time taken for a thread to fill the flow-focusing constriction and the time required for the thread to pinch-off to form a bubble. The growth time is controlled by the pressure difference between the gas thread and the surrounding liquid, which in turn is controlled by the liquid flow rate and the density of bubbles downstream from the constriction. For steady-state bubble production, the frequency of bubble production can be approximated by *f*∝*q*_{ℓ}(*P*−*P*_{c})/*P*, which is similar to experimental results obtained by others.

For single flow-focusing devices, downstream effects are most important when the flow is pressure driven, as the channel flow primarily affects the absolute upstream pressure. In geometrically controlled flow-focusing devices, this feedback stabilizes the flow—increased bubble production increases the channel resistance, which in turn lowers the bubble production rate. The negative feedback in the channels damps frequency transients. For massively parallel devices, however, these feedback effects become important even for flow rate driven devices, as these would typically be supplied from a single syringe pump. This manner of forcing would cause the flow rate at a particular junction to be set by the resistance through the entire parallel arm, so that device feedback can have a pronounced effect, as seen by Barbier *et al*. (2006), even when all components of the system are incompressible. In designing parallel devices for monodisperse foam and emulsion formation, care must be taken to mitigate the downstream effects of bubbles or droplets on flow, especially at channel junctions (see Engl *et al*. 2005; Jousse *et al*. 2006). Feedback, even in parallel systems, need not be simply a liability in the operation of microfluidic devices. For example, Prakash & Gershenfeld (2007) have taken advantage of the nonlinear interactions at channel junctions to produce a variety of microfluidic logic operations and oscillators. The nonlinear operation of multiphase flow generators might be similarly leveraged.

## Acknowledgments

We would like to acknowledge the support from the Harvard MRSEC (DMR-0213805) and the Schlumberger Future of Research programme.

## Footnotes

One contribution of 11 to a Theme Issue ‘New perspectives on dispersed multiphase flows’.

- © 2008 The Royal Society