## 1 Introduction

Emulsions, foams, dispersions are example of multiphase fluids. These materials are ubiquitous in nature and in industry: e.g., concentrated emulsions, colloidal dispersion, granular flows, in the food and in the cosmetic industry [1]. In many of these applications the dispersion is polydisperse and poorly controlled due to the difficulties encountered during the emulsification process. Renewed interest in emulsification field has been brought by the demonstration that microfluidic devices may produce highly monodisperse droplet at the micrometer size [2]. This control of the droplet size distribution has been demonstrated to be very useful in new applications and for developing new technologies. Monodisperse capsules, multiple emulsions [3,4] or janus particles [5] have been produced using microfluidic devices. The perfect control of the size and velocity of droplet allows one to provide measurements of interest to the biological and chemical communities. Droplets in microfluidic devices has been used to study chemical reaction kinetics on millisecond time scales [6] by injecting reactive chemicals together to form drops and in developing platform for crystallisation and phase diagram [7–9]. Droplets are particularly advantageous for this sort of study because reagents are not dispersed beyond the boundary of the drop. In this article, we deal with the mechanisms involved in the droplet formation. We first present an overview of the study of droplets formation in microfluidic devices. We then report some experiments in various geometries and we present a simple analytic model allowing us to predict the nature of the flow. Section 3 is devoted to cylindrical geometry. Section 4 deals with rectangular geometries which are commonly encountered in microfluidic devices.

## 2 Biphasic flow in microfluidic devices

Biphasic flow has been widely studied in microfluidic devices [10]. The wetting properties [11] of the microchannel crucially control the patterns obtained [12]. Direct emulsions are produced in glass devices whereas, an inversed one is produced in PDMS devices. The nature of the flow does also depend upon the geometry of the microchannel. Using a flow focusing geometry, Anna et al. [2] have studied drop formation in liquid–liquid systems as a function of flow rates and flow rate ratios of the two liquids. They present a phase diagram including one regime where drop size is comparable to orifice width and a second regime where drop size is dictated by the diameter of a thin “focused” thread, so drops much smaller than the orifice are formed. They point out that both monodisperse and polydisperse emulsions can be produced. Drops are not the single pattern that may be obtained in microfluidic devices. Jets or truncated jets [13] are commonly encountered for high values of the non-wetting phase flow rate. Surface tension and the viscosity ratio rule the drop size and the transition between the parallel flow regime and the drop regime [14]. The understanding of the droplet formation requires the use of flow stability analysis or of numerical computation [15]. Flow stability analysis has been widely used to apprehend break-up scenario in unbounded flowing systems. Determining the conditions required to get an absolute or a convective instabilty allows one to predict the zone of droplets and jets production. Absolute instability corresponds to disturbances growing and propagating both in the downstream and upstream directions. A continuous jet cannot thus exist and typically, drops are released intermittently either right at the injection nozzle or at a finite distance from it to form a dripping jet. At the opposite, convective instability corresponds to perturbations that propagate downstream while they grow, allowing for a long continuous fluid thread to persist. Several experimental studies support this picture and the link between the absolute/convective transition of the instability to the dripping/jetting transition in the observed spatio-temporal behaviour of the biphasic flow [16–24]. In the following, we will apply these considerations to flows in cylindrical microfluidic devices and in rectangular cross-section microfluidic devices.

## 3 Drops and jets in microfluidic devices with a circular cross-section

### 3.1 Experimental results in cylindrical geometries

We generate a jet in a cylindrical glass capillary of inner radius R_{c} [25], using as a nozzle a glass capillary of square cross-section with a tapered end (see Fig. 1).

The outer dimension of this square capillary is very close to the inner diameter of the cylindrical tube which ensures good alignment and centering. R_{c} is in the 200–500 μm range, whereas the radius of the tapered orifice of the square tube is set between 20 and 50 μm using a pipette–puller set up. Syringe pumps are used to inject an inner fluid of viscosity η_{i} at a rate Q_{i} in the square capillary and the outer fluid of viscosity η_{e} at a rate Q_{e} through the cylindrical capillary. This leads to coaxial injection at the tapered orifice. We observe flow patterns which vary significantly with operational (Q_{e}, Q_{i}), geometrical (R_{c}), and system parameters (η_{i}, η_{e}, surface tension). Fig. 2 displays the typical outcome of an experiment where the flow rates are varied for a given system (here the inner solution is 50 in weight glycerine in water solution with η_{i} 55 mPa s and the outer one a silicone oil for which η_{e} 235 mPa s). A droplet regime is found for low Q_{i}, with either droplets emitted periodically right at the nozzle symbol (open circle) or non spherical plug like droplets resulting from the instability of an emerging oscillating jet (filled gray circle). Jets are found in the bottom right corner of Fig. 2 with different visual aspects: wavy jets with features that are convected downstream (open square), and for larger values of Q_{i}, straight jets (filled square) that persist throughout the cylindrical capillary. For large values of the external flow rate Q_{e}, we observe what we call jetting: thin and rather straight jets (open diamond) that extend over some distance in the capillary tube before breaking into droplets at a well-defined and reproducible location. This jet length increases with Q_{i} for a fixed value of Q_{e}.

### 3.2 Theoretical description of the droplets jet transition in cylindrical geometries

We now attempt to model these phenomena analytically. To reach this aim, we will perform a linear analysis of the stability of the flow and determine the conditions required to get an absolute or a convective instability. Beyond the importance of viscous forces (the Reynolds numbers are typically small to moderate), the essential contrast with most of the previous studies which focused on unbounded flows is the major role of the microchannel walls which induce parabolic flow profiles and strongly affect the development of perturbations.

In order to get an analytical description, we proceed with approximations. We neglect inertial effects (i.e., as the Reynolds number is small in most of our experiments), and we use lubrication theory (i.e., we formally assume that the wavelengths of the perturbations are long compared to the capillary radius) which has been shown to be remarkably insightful in somewhat related situations [26].

### 3.3 Lubrication analysis for the cylindrical geometry

We first consider a cylindrically symmetric geometry (see Fig. 1). In a capillary tube of radius R_{c} are flown two immiscible and incompressible liquids, an “inner” fluid of viscosity η_{i} is injected at rate Q_{i} in the stream of an “external” liquid of viscosity η_{e} flowing at rate Q_{e}.

In the unperturbed state, the flow is unidirectional with an inner fluid jet of radius ${r}_{\text{i}}^{\text{o}}$ and pressure gradients ∂_{z}P_{e} and ∂_{z}P_{i} in the two fluids that are constants. Together with the boundary conditions at the surface of the jet (continuity of the velocity field and tangential shear stress) and at the walls of the geometry (no slip condition) and local force balance leads to:

$${Q}_{\text{e}}^{\text{o}}=\frac{{R}_{\text{c}}^{2}{r}_{\text{i}}^{\text{o}2}}{4{\eta}_{\text{e}}}\left(\frac{{r}_{\text{i}}^{\text{o}2}}{{R}_{\text{c}}^{2}}-1-2\frac{{r}_{\text{i}}^{\text{o}2}}{{R}_{\text{c}}^{2}}\phantom{\rule{0.25em}{0ex}}\mathrm{ln}\left(\frac{{r}_{\text{i}}^{\text{o}}}{{R}_{\text{c}}}\right)\right)\left({\partial}_{z}{P}_{\text{i}}^{\text{o}}-{\partial}_{z}{P}_{\text{e}}^{\text{o}}\right)-\frac{\mathrm{\pi}{\left({R}_{\text{c}}^{2}-{r}_{\text{i}}^{\text{o}2}\right)}^{2}}{8{\eta}_{\text{e}}}{\partial}_{z}{P}_{\text{e}}^{\text{o}}$$ | (1) |

$${Q}_{\text{i}}^{\text{o}}=\frac{\mathrm{\pi}{r}_{\text{i}}^{\text{o}4}}{2{\eta}_{\text{e}}}\phantom{\rule{0.25em}{0ex}}\mathrm{ln}\left(\frac{{r}_{\text{i}}^{\text{o}}}{{R}_{\text{c}}}\right)\left({\partial}_{z}{P}_{\text{i}}^{\text{o}}-{\partial}_{z}{P}_{\text{e}}^{\text{o}}\right)-\frac{\mathrm{\pi}{r}_{\text{i}}^{\text{o}4}}{8{\eta}_{\text{i}}}{\partial}_{z}{P}_{\text{i}}^{\text{o}}+\frac{\mathrm{\pi}{r}_{\text{i}}^{\text{o}2}{\left({r}_{\text{i}}^{\text{o}2}-{R}_{\text{c}}^{2}\right)}^{2}}{4{\eta}_{\text{e}}}{\partial}_{z}{P}_{\text{e}}^{\text{o}}\text{.}$$ | (2) |

$${P}_{\text{i}}^{\text{o}}-{P}_{\text{e}}^{\text{o}}=\frac{\Gamma}{{r}_{\text{i}}^{\text{o}}}\text{,}$$ | (3) |

We perform a linear stability analysis and consider the spatial-temporal response of the system to small z-dependent cylindrically symmetric perturbations δQ_{e}, δQ_{i}, ∂_{z}δP_{e}, ∂_{z}δP_{i} and δr_{i}. We make the perturbations proportional to e^{(ikz+ωt)} with k and ω complex numbers. As indicated in Section 1, we restrict our analysis to the lubrication approximation, assuming formally that the perturbation wavelength is larger than the capillary radius R_{c}. In this framework, the expressions obtained for the unperturbed flow can still be used locally, so that the local perturbations in the flow rates read:

$$\delta {Q}_{\text{e}}=\frac{\partial {Q}_{\text{e}}^{\text{o}}}{\partial \left({\partial}_{z}{P}_{\text{e}}^{\text{o}}\right)}{\partial}_{z}\delta {P}_{\text{e}}+\frac{\partial {Q}_{\text{e}}^{\text{o}}}{\partial \left({\partial}_{z}{P}_{\text{i}}^{\text{o}}\right)}{\partial}_{z}\delta {P}_{\text{i}}+\frac{\partial {Q}_{\text{e}}^{\text{o}}}{\partial {r}_{\text{i}}^{\text{o}}}{\partial}_{z}\delta {r}_{\text{i}}\text{,}$$ | (4) |

$$\delta {Q}_{i}=\frac{\partial {Q}_{\text{i}}^{\text{o}}}{\partial \left({\partial}_{z}{P}_{\text{e}}^{\text{o}}\right)}{\partial}_{z}\delta {P}_{\text{e}}+\frac{\partial {Q}_{\text{i}}^{\text{o}}}{\partial \left({\partial}_{z}{P}_{\text{i}}^{\text{o}}\right)}{\partial}_{z}\delta {P}_{\text{i}}+\frac{\partial {Q}_{\text{i}}^{\text{o}}}{\partial {r}_{\text{i}}^{\text{o}}}{\partial}_{z}\delta {r}_{\text{i}}\text{,}$$ | (5) |

However, an important difference is that, as the radius of the jet varies, the Laplace law now requires that the two pressure gradients are different:

$${\partial}_{z}\left(\delta {P}_{\text{i}}-\delta {P}_{\text{e}}\right)=-\Gamma \left({\partial}_{z}\left(\frac{\delta {r}_{\text{i}}}{{R}_{\text{c}}^{2}}+{\partial}_{z}^{2}\delta {r}_{\text{i}}\right)\right)\text{,}$$ | (6) |

_{z}

^{2}δr

_{i}to the curvature in the flow direction.

This set of equations is equivalent to computation of the linear response in δr_{i}, through a second order expansion in powers of $\epsilon ={R}_{\text{c}}/{L}_{z}$ of the velocity field, where ${L}_{z}=2\mathrm{\pi}/k$ is the characteristic length scale involved in the z direction. We recall that in the lubrication approximation L_{z} is larger than the capillary radius R_{c} which implies that ε ≪ 1. Note however, that an additional approximation is required to get these expressions. Specially, in the calculation of the O(ε^{2}) term, where it is equivalent to neglect the pressure gradient created by the zeroth order in ε velocity compared to the one generated by the interface curvature along z. We will quantify more precisely this approximation at the end of this section.

Mass conservation of the incompressible fluids allows to close the system of equations:

$$\delta {Q}_{\text{e}}+\delta {Q}_{\text{i}}=0\text{,}$$ | (7) |

$${\partial}_{t}\left(\mathrm{\pi}{\left({r}_{\text{i}}^{\text{o}}+\delta {r}_{\text{i}}\right)}^{2}\right)={\partial}_{z}\delta {Q}_{\text{e}}\text{.}$$ | (8) |

$$x=\frac{{r}_{\text{i}}^{\text{o}}}{{R}_{\text{c}}},\phantom{\rule{1em}{0ex}}\lambda =\frac{{\eta}_{\text{i}}}{{\eta}_{\text{e}}},\phantom{\rule{1em}{0ex}}\tilde{k}={r}_{\text{i}}^{\text{o}}k,\phantom{\rule{1em}{0ex}}\tilde{\omega}=\frac{\omega 16{\eta}_{\text{e}}{R}_{\text{c}}}{\Gamma}\text{,}$$ |

$$\tilde{\delta {Q}_{\text{i}}}=\frac{8{\eta}_{\text{e}}}{\mathrm{\pi}{R}_{\text{c}}^{4}{\partial}_{z}{P}_{\text{i}}^{\text{o}}}\delta {Q}_{\text{i}},\phantom{\rule{1em}{0ex}}{\partial}_{z}\tilde{\delta {P}_{\text{i}}}=\frac{{\partial}_{z}\delta {P}_{\text{i}}}{{\partial}_{z}{P}_{\text{i}}^{\text{o}}}\text{,}$$ |

$$\tilde{\delta {Q}_{\text{e}}}=\frac{8{\eta}_{\text{e}}}{\mathrm{\pi}{R}_{\text{c}}^{4}{\partial}_{z}{P}_{\text{e}}^{\text{o}}}\delta {Q}_{\text{e}}\phantom{\rule{1em}{0ex}}\text{and}\phantom{\rule{1em}{0ex}}{\partial}_{z}\tilde{\delta {P}_{\text{e}}}=\frac{{\partial}_{z}\delta {P}_{\text{e}}}{{\partial}_{z}{P}_{\text{e}}^{\text{o}}}\text{.}$$ | (9) |

$$\delta \tilde{{Q}_{\text{e}}}=a\left(x\right){\partial}_{z}\delta \tilde{{P}_{\text{e}}}+b\left(x\right)\delta x+c\left(x\right)\left({\partial}_{z}\delta \tilde{{P}_{\text{i}}}-{\partial}_{z}\delta \tilde{{P}_{\text{e}}}\right)\text{,}$$ | (10) |

$$\delta \tilde{{Q}_{\text{i}}}=d\left(x\right){\partial}_{z}\delta \tilde{{P}_{\text{e}}}+e\left(x\right)\delta x+f\left(x\right)\left({\partial}_{z}\delta \tilde{{P}_{\text{i}}}-{\partial}_{z}\delta \tilde{{P}_{\text{e}}}\right)\text{,}$$ | (11) |

$$a\left(x\right)=-{\left(1-{x}^{2}\right)}^{2}\phantom{\rule{1em}{0ex}}d\left(x\right)=2{x}^{2}\left({x}^{2}-1\right)-{\lambda}^{-1}{x}^{4}\text{,}$$ |

$$b\left(x\right)=4x\left(1-{x}^{2}\right)\phantom{\rule{1em}{0ex}}e\left(x\right)=8{x}^{3}-4{\lambda}^{-1}{x}^{3}-4x\text{,}$$ |

$$c\left(x\right)=2{x}^{4}-2{x}^{2}-4{x}^{4}\phantom{\rule{0.25em}{0ex}}\mathrm{ln}\left(x\right)\phantom{\rule{1em}{0ex}}f\left(x\right)=4{x}^{4}\phantom{\rule{0.25em}{0ex}}\mathrm{ln}\left(x\right)-{\lambda}^{-1}{x}^{4}\text{.}$$ |

$$\tilde{\omega}=\frac{-i\tilde{k}\xb7{x}^{3}\mathrm{Ka}E\left(x\right)+F\left(x\right)({\tilde{k}}^{2}-{\tilde{k}}^{4})}{D\left(x\right)}$$ | (12) |

$$\mathrm{Ka}=\frac{-{\partial}_{z}{P}^{\text{o}}{R}_{\text{c}}^{2}}{\Gamma}$$ |

$$E\left(x,\lambda \right)=-4x+\left(8-4{\lambda}^{-1}\right){x}^{3}+4\left({\lambda}^{-1}-1\right){x}^{5}\text{,}$$ |

$$F\left(x,\lambda \right)={x}^{4}\left(4-{\lambda}^{-1}+4\phantom{\rule{0.25em}{0ex}}\mathrm{ln}\left(x\right)\right)+{x}^{6}\left(-8+4{\lambda}^{-1}\right)+{x}^{8}\left(4-3{\lambda}^{-1}-\left(4-4{\lambda}^{-1}\right)\mathrm{ln}\left(x\right)\right)\text{,}$$ |

$$D\left(x,\lambda \right)={x}^{9}\left(1-{\lambda}^{-1}\right)-{x}^{5}\text{.}$$ |

The first term of the dispersion equation $-i\xb7{x}^{3}\mathrm{Ka}\left(E\left(x\right)/D\left(x\right)\right)\tilde{k}$ summarizes the convection kinematics and the mass conservation; the second term $\left(F\left(x\right)/D\left(x\right)\right)\left({\tilde{k}}^{2}-{\tilde{k}}^{4}\right)$ describes the convection effect. E and F are positive functions in cylindrical geometry. Ka is a genuine capillary number as it is the ratio between viscous forces ${\eta}_{\text{e}}{V}_{\text{e}}{R}_{\text{c}}=-{\partial}_{z}{P}^{\text{o}}{R}_{\text{c}}^{3}$ and capillary forces ΓR_{c}. We have however used a notation different from the usual Ca to call the reader's attention to the fact that Ka is a capillary number at the scale of the capillary R_{c} rather than at the scale of the jet. Obviously in contrast with our approach, studies of unconfined jets focus on capillary numbers defined using either the average jet velocity or the velocity at the surface of the unperturbed jet.

The system is stable if all ω_{r}s are negative. This is not the case. Indeed, low real k values, corresponding to large wavelengths, are unstable since F(x) is a positive value. Note that the instability comes from the k^{2} term which is related to the curvature in the cross-section of the jet. Decreasing the jet radius size, decreases the interfacial area and energy cost and promotes instability. The term k^{4}, related to the curvature in the flow direction, is a stabilizing term since undulations in the flow direction increase the interfacial energy cost. With some ω_{r} positive, an initially localized perturbation generates a growing distortion.

As the perturbation grows, a leading edge profile is selected by the flow, such that the long-time profile is dominated by the mode K_{r} corresponding to the largest growth rate w_{r} and an external velocity for the envelope. This leads to the following characteristics for the selected perturbations [27–29]:

$${v}^{\ast}=\frac{{\omega}_{\text{r}}^{\ast}}{{k}_{\text{i}}^{\ast}},\phantom{\rule{0.50em}{0ex}}\frac{\partial {\omega}_{\text{r}}^{\ast}}{\partial {k}_{\text{r}}}=0\phantom{\rule{1em}{0ex}}\text{and}\phantom{\rule{1em}{0ex}}{v}^{\ast}=\frac{\partial {\omega}_{\text{r}}^{\ast}}{\partial {k}_{\text{i}}}\text{.}$$ | (13) |

^{∗}is negative and convectively unstable if all the solutions v

^{∗}are positive. In the former case, a continuous jet can no longer be formed with perturbations propagating backwards. This leads typically to droplets being released intermittently to form a dripping jet, either right at the injection nozzle or at any further point of the entrance zone where the inner stream adjusts its size to the imposed environment. At the opposite, if all v

^{∗}s are positive, the disturbances which grow in time are simultaneously convected downstream and a continuous jet can persist in the system. Fig. 3 gives a schematic representation of the two cases.

Solutions of the dispersion equation that satisfy Eq. (13), lead to four wave vectors that are independent of the imposed flow (or capillary number):

$${\tilde{k}}_{1}^{\ast}=+\sqrt{\frac{9+3\sqrt{7}}{24}}+\mathrm{i}\sqrt{\frac{-1+\sqrt{7}}{24}}\text{,}$$ |

$${\tilde{k}}_{2}^{\ast}=-\sqrt{\frac{9+3\sqrt{7}}{24}}+\mathrm{i}\sqrt{\frac{-1+\sqrt{7}}{24}}\text{,}$$ |

$${\tilde{k}}_{3}^{\ast}=+\sqrt{\frac{9+3\sqrt{7}}{24}}-\mathrm{i}\sqrt{\frac{-1+\sqrt{7}}{24}}\text{,}$$ |

$${\tilde{k}}_{4}^{\ast}=-\sqrt{\frac{9+3\sqrt{7}}{24}}-\mathrm{i}\sqrt{\frac{-1+\sqrt{7}}{24}}\text{.}$$ |

$${\omega}_{\text{r}}^{\ast}=\frac{\mathrm{Ka}E\left(x\right){x}^{3}\tilde{{k}_{\text{i}}^{\ast}}+F\left(x\right)\left(5+\sqrt{7}/18\right)}{D\left(x\right)}\text{.}$$ | (14) |

^{∗}always positive. We thus focus on the disturbances associated to ${k}_{3}^{\ast}$ and ${k}_{4}^{\ast}$, which correspond to the same physics as they differ only by a phase factor. The velocity of the associated perturbations is given by:

$$\tilde{{v}^{\ast}}=\frac{\mathrm{Ka}E\left(x,\lambda \right){x}^{3}-{C}_{1}F\left(x,\lambda \right)}{D\left(x\right)}\text{,}$$ | (15) |

For sufficiently low capillary numbers Ka, the corresponding v^{∗} is negative (and ω_{r}^{∗} > 0), whereas in the opposite limit of high flow speeds and large Ka, v^{∗} becomes positive. This suggests an absolute to convective instability transition as the flow rate is increased, with the associated transition from dripping to a continuous jet. We thus reach a rather simple analytical prediction for the transition, $\tilde{{v}^{\ast}}(\mathrm{Ka},x,\lambda )=0$, plotted on Fig. 4 in the (x, Ka) plane describing operational conditions for various values of the viscosity ratio system λ = η_{i}/η_{e}. This plot can be envisioned as a dynamic behaviour diagram with a dripping and a jet regime. For a given λ, increasing the capillary number Ka (i.e., the normalized pressure drop) or the confinement x always eventually lead to a continuous jet. This is physically sound as increasing Ka corresponds to convecting away the perturbations faster, while increasing the confinement x results in slowing down the development rate of the perturbations due to the proximity of the walls. Decreasing λ = η_{i}/η_{e} increases the “droplets” regime at the expense of the “jet” regime.

At this stage, we may comment the validity of our approximations. Two different approximations have been made. First, we used a dispersion equation derived in the lubrication approximation framework. This implies that the parameter $\varepsilon ={r}_{\text{i}}^{\text{o}}/{L}_{z}$ is smaller than 1. The selected wave vector corresponds to a characteristic length in the z direction, L_{z}, equals to ${L}_{z}=2\mathrm{\pi}\sqrt{24/9+3\sqrt{7}}{r}_{\text{i}}^{\text{o}}$. This leads to a value of $\varepsilon ={r}_{\text{i}}^{0}/{L}_{z}=0.0659/x$. Assuming that the use of the dispersion equation obtained in the lubrication approximation framework is valid for ε < 0.2, we can conclude that this hypothesis is valid for x > 0.33. Second, we neglected the recirculations inside the perturbated jets. Indeed, in the calculation of the O(ε^{2}) term, we neglected the pressure gradient created by the zeroth order in ε velocity compared to the one generated by the interface curvature along z. This is valid for x > 0.33 if ${\eta}_{\text{i}}/{\eta}_{\text{e}}<1$ and for x > 0.6 if ${\eta}_{\text{i}}/{\eta}_{\text{e}}>1$. These ranges of parameters correspond thus to the range where our diagrams are fully valid. Note that this point will prevent us to compare quantitatively our results with the ones obtained in unbounded geometries [22,30].

### 3.4 Comparison with experiments

We now quantitatively compare our predictions to experimental data obtained for two surface tensions (see Fig. 5) and two capillary radii R_{c}s (see Fig. 6).

Clearly, given the approximations involved, our simple model describes very well the experimental data with no adjustable parameters, and appears as a powerful predictive tool. A certain level of disagreement is expected for weak confinement (small Q_{i} large Q_{e}) as the lubrication approximation is formally invalid in this case. Our model indeed overestimates Ka at the transition for vanishing x, as demonstrated by comparison to the exact result of Gañán-Calvo [30] for unbounded creeping flows. Inertial effects may also slightly alter the picture for the largest outer flow rates. We also report in Figs. 7 and 8 all the experimental data obtained for a given viscosity ratio, but for various surface tensions, various capillary radii, and various flow rates. This shows the relevance of our description in terms of the dimensional variables Ka and x. An additional interest of this mapping onto the (Ka, x) plane of this large set of data is that it collapses relatively well the different types of flows observed within the droplet and jet regimes as can be seen by the grouping of the symbols. In short, droplets correspond to small capillary number Ka and small confinement ratio x, while plugs require larger values of x. At higher values of Ka, increasing the confinement ratio x shifts the behaviour from jetting (with emission of droplets at a large finite distance from the nozzle) to stable jets.

## 4 Drops and jets in microfluidic devices with a rectangular cross-section

### 4.1 Experiments in rectangular geometries

In this section, we turn to the most commonly encountered geometry in microfluidics, namely that of a microchannel with a rectangular cross-section.

Our microfluidic devices are fabricated using soft lithography technology [31]. Polydimethylsiloxane (PDMS) channels or PDMS–glass channels are used. The microdevices have two inlet arms which meet at a T junction with a funnel design (see Fig. 9). The flow patterns are observed in the outlet channel after the T junction. The inlet channels are connected via tubing to syringes loaded with the fluids. Syringe pumps allow us to control the flow rates of the liquids. In this study the immiscible fluids used are aqueous and oil solutions. Oils are silicone oils (Rhodorsil) of different viscosities and hexadecane (Prolabo). Aqueous solutions are dilute solutions of Sodium dodecyl sulphate (SDS) (Merck) above the critical micellar concentration (cmc) in water (cmc = 0.24 wt%) or in a mixture of glycerine (Prolabo) and water. Adding glycerine allows us to match the optical index of the aqueous solution with silicone oil and thus, enables us to perform fluorescent confocal microscopy.

Fig. 10 shows the flow pattern diagram obtained with a mixture of 50 wt% of water with SDS at 2 cmc and 50 wt% of glycerine as aqueous phase and hexadecane in a 100 μm × 100 μm microchannel. The viscosity of the aqueous phase is of 7 cP and the oil phase is of 3 cP. During the experiment, three typical flow patterns have been identified: droplets formed at the T junction (DTJ), parallel flows or jets (PF) and parallel flows which break into droplets inside the channel (DC). Droplets are obtained in the upper left corner of Fig. 10, i.e., at low water flow rates and substantial oil flow rates. Two kinds of regimes must be distinguished. In the DTJ regime represented by (○), droplets are not always produced with reproducibilty in size whereas in the DC regime (●), monodisperse droplets are produced. In the latter case, a parallel flow (i.e., a jet) takes place at the entrance of the channel and becomes unstable after a short distance, breaking up into droplets. In the lower right, corner jets are observed. To investigate the three dimensional structure of the flow, confocal fluorescence microscopy experiments are performed. In this case rhodamine (Sigma) is added to the aqueous phase and hexadecane is used as the oil phase. Fig. 11 displays cross-section images obtained for PF in a 100 μm × 100 μm microchannel. The channel walls represented here have been drawn and added to the picture. The bottom wall is in glass and the three others are in PDMS. At the entrance of the funnel, hexadecane starts wrapping around the water. This is due to the difference of wettability of the two fluids with the PDMS walls. The water flow cross-section continues evolving until the flow mainly evolves along the propagation axis. Depending on the flow rates, jets or DC are then observed. PF is a truncated cylinder jet located on the glass at the PDMS glass corner. The pressure in a cross-section is uniform inside each fluid and differs from one fluid to the other. The Laplace law imposes that the free fluid surface has a single curvature radius.

## 5 Extension of our model to rectangular geometry

In this section, we try to extend the previous analysis to a microchannel with a rectangular cross-section. Two situations have to be analyzed. The inner jet may or may not be squeezed by the geometric confinement. We first focus on situations where the jet touches the upper and lower walls.

### 5.1 Squeezed jet in rectangular geometry

We consequently consider as the reference geometry a jet squeezed by the geometry. We assume that the contact angle θ between the internal and external phases on the glass is constant. In the unperturbated state, the free fluid–fluid interface has a shape involving a single radius curvature (see Fig. 12).

We wish to analyze again the stability of such a jet to axial perturbations without having to move to 3D numerical simulations. We therefore proceed with an uncontrolled approximation and consider only axially symmetric perturbations of the radius of the jet independent of the polar angle. In the small perturbations case and for a squeezed jet, the perturbated jet remains squeezed by the geometric confinement. This point is very important, since it sets the shape of the perturbations. Indeed, it ensures that the curvature radius remains constant in the geometry cross-section. The first term in the Laplace equation (Γ∂_{z}δr_{i}) is thus nil. As a consequence the dispersion equation in the approximation of the lubrication reads:

$$\tilde{\omega}=\frac{-i\tilde{k}\xb7{x}^{3}\mathrm{Ka}E\left(x\right)-F\left(x\right)\left({\tilde{k}}^{4}\right)}{D\left(x\right)}\text{.}$$ | (16) |

_{r}s are negative and the 2D jet is thus stable. Confining the jet and squeezing it thus ensures its stability. Perturbations leading to the interfacial area reduction are unable to develop. The same results have been previously derived by Miguel et al. [32]. Note also that this scenario of stable “2D” jet where the area reduction is killed by the confinement had been experimentally shown by Dollet et al. [33] in a flow focusing design. This results explain why coflow situations are often encountered in microfluidic devices [13,34].

### 5.2 Non-squeezed jet in rectangular geometry

We now deal with situations where the jet is not squeezed by the wall. The flow studied is a wetting jet on the glass plate as observed in confocal experiments. In our calculation we will only consider axisymmetric perturbation of the radius of the jet. We assume that the contact angle θ between the internal and external phases on the glass is constant. With this assumption the center of the jet depends on the radius r and θ as depicted in Fig. 13.

To calculate the front velocity of the propagation we need to calculate the variation of the internal flow rate δQ_{i} and the external one δQ_{e} as a function of a small variation of the radius (δr) and of the pressure in both fluids (δP_{e} and δP_{i}):

$$\delta {Q}_{\text{e}}=\frac{\partial {Q}_{\text{e}}}{\partial \nabla {P}_{\text{e}}}\partial \delta {P}_{\text{e}}+\frac{\partial {Q}_{\text{e}}}{\partial \nabla {P}_{\text{i}}}\partial \delta {P}_{\text{i}}+\frac{\partial {Q}_{\text{e}}}{\partial r}\partial \delta r$$ | (17) |

$$\delta {Q}_{\text{i}}=\frac{\partial {Q}_{\text{i}}}{\partial \nabla {P}_{\text{e}}}\partial \delta {P}_{\text{e}}+\frac{\partial {Q}_{\text{i}}}{\partial \nabla {P}_{\text{i}}}\partial \delta {P}_{\text{i}}+\frac{\partial {Q}_{\text{i}}}{\partial r}\partial \delta r$$ | (18) |

_{e}and ∇P

_{i}are the pressure gradients in the internal and external phases, respectively.

The flow rates are numerically calculated using the Stokes equation in the geometry of Fig. 13. The boundary condition used is no-slip at the wall and the continuity of the velocity and of the shear stress at the fluid–fluid interface [35,36]. Calculations are made using the scientific software Scilab developed by INRIA. The numerical procedure used to calculate the flow profile is described in precision elsewhere [37]. We used a Cartesian regular mesh whose size is 1 μm × 1 μm. To calculate the different terms of Eqs. (17) and (18) we calculate the flow rates for both fluids for the different conditions described in Fig. 14. We calculate then the flow rates at the transition between the absolute and convective instabilities for different pressure gradients.

## 6 Comparison with experimental data

The absolute/convective transition is calculated taking into account the geometry of the channel, the viscosity of both phases and the contact angle θ measured either by confocal microscopy pictures or by the transmission picture (to precision on this measure see Ref. [37]). Fig. 15 shows results with silicone oil sytems. Lines are theoretical transitions between absolute and convective and symbols (● and ○) are experimental transitions obtained for two viscosities (100 mPa s and 50 mPa s) and two surface tensions (Γ = 15 mN/m and Γ = 30 mN/m).

Results are in good agreement without any adjustable parameter and the theoretical transitions move in the same way as the experimental ones. The same conclusions are reached when the external phase is hexadecane with (◊) and without () Span 80 (see Fig. 16). When the internal phase viscosity is changed to 7 mPa s, good agreement is obtained with hexadecane (), silicone oil at 20 mPa s (▾) and 100 mPa s (●).

## 7 Conclusion

In this work, using the lubrication approximation and focusing on low Reynolds numbers, we have studied the stability of a pressure-driven jet confined in a capillary. Analyzing the transition from a continuous jet to a dripping system in terms of a convective/absolute transition, we find the jet remains continuous at high capillary numbers and breaks down into droplets at lower flow rates, i.e., pressure gradients. Analytical formulas are provided in cylindrical capillaries for a given system where viscosities and surface tension are known. The influence of the various system and operation parameters has been highlighted. At the expense of further approximations, we have explored other microchannel geometries, that are more relevant in the microfluidics context given fabrication processes. Microchannels of rectangular cross-section promote droplets at strong confinement when contrasted with their cylindrical counter-parts, most likely because corners allow easier fluid transfer. The outlook deals with the modelization of the volume of the droplets. Numerical simulations using level set methods are under way in order to address this problem.