Internal flow in droplets within a concentrated emulsion flowing in a microchannel

Droplet microfluidics, in which droplets act as individual nanoliter- to picoliter-sized reactors, has enabled a wide range of high-throughput screening applications including digital polymerase chain reactions and directed evolution of enzymes.^1–4 Important to the efficient mixing of reagents inside the droplets, the internal flow and recirculation of fluids inside single droplets flowing in a channel have been studied extensively and are relatively well understood.^5–9 The recirculation or rotational patterns in the droplets are driven by differences in the velocity of the droplet and the continuous phase.¹⁰ It has been shown that such patterns depend on multiple parameters including the capillary number, the size of the droplet relative to the channel, the choice and concentration of surfactants, as well as the viscosity ratio between the droplet and the continuous phase.^5,6,11

Although the fluid dynamics of single droplets is well understood,^{5,6,8,12–18} to the best of our knowledge, no work has probed the flow field inside droplets of concentrated emulsions at high volume fractions (>80%) in confined microchannels. Understanding how these emulsions flow is important for droplet microfluidics especially during incubation and interrogation steps where the emulsion is often pre-concentrated prior to manipulation.^19,20 The ability to predict and control the flow of concentrated emulsions will also be useful for further increasing the throughput of droplet microfluidic systems, as it removes the need to handle large volumes of continuous phase required in dilute emulsions for processing the same number of droplets. Prior work on concentrated emulsions has focused primarily on their bulk rheological properties.^21–23 Nevertheless, understanding how individual droplets behave within the emulsion is important as each droplet acts as an individual reactor and can carry a different reaction in droplet microfluidics applications. It is also critical for the fundamental understanding of how the microstructure of the emulsion and the local flow field at the single-drop level relate to its bulk rheological properties.

In this paper, we describe the use of microscopic particle image velocimetry (μPIV) to quantify the flow inside individual droplets within a concentrated emulsion having volume fraction φ ∼ 85% flowing as a monolayer in a straight microfluidic channel. Specifically, we investigate the effect of confinement (namely, channel width relative to the droplet size or the number of rows of droplets across the width of the channel) and viscosity ratio on the internal flow patterns. We identify two features in the flow that were never described before: (1) The flow patterns inside the droplets are dependent on the confinement, as well as the location of the droplets in the channel. (2) The flow patterns inside the droplets are independent of the viscosity ratio for the range of values tested. Based on previous studies of the flow of a single droplet in a microchannel, the flow in droplets within a concentrated emulsion in a microchannel is likely to be three-dimensional.^5,6,14 In this work, we focus on two-dimensional measurements at the mid-height of the microchannel only. Three-dimensional measurements are an ongoing investigation.

We used three types of water-in-oil (w/o) emulsions for our experiments (Table I). The disperse phase consisted of a combination of deionized water and glycerol at different weight fractions. The viscosity of these mixtures was measured using a rheometer (TA Instruments ARES-G2). We note that the measured values are similar to those obtained from the literature.²⁴ The viscosity ratio, λ, is defined as the ratio of the dynamic viscosity of the disperse phase, μ_d, to the dynamic viscosity of the continuous phase, μ_c. The interfacial tension between the phases was measured by Rame-Hart 290 contact angle goniometer (Rame-Hart Co., NJ, USA). The continuous phase was a hydrofluoroether (HFE-7500, 3M) containing an ammonium salt of Krytox (DuPont) at 2% w/w as a surfactant for stabilizing the droplets against coalescence. The droplets were generated using a flow-focusing device fabricated in poly(dimethylsiloxane) (PDMS).²⁵ This channel was rendered hydrophobic by treatment with Aquapel (Pittsburgh, PA). In all our experiments, the droplet volume was 90 pl with a volume dispersity of <3%.

After the droplets were generated, they were collected and allowed to cream to the top of a syringe for 12 h. As the specific gravity of HFE-7500 is 1.6, we were able to obtain an emulsion with a high volume fraction (φ ∼ 85%) through buoyancy forces. The concentrated emulsion was then re-injected into straight microchannels for μPIV experiments. The height of the microchannels was fixed at 30 μm, less than the diameter of our droplets when spherical. The droplets had a “pancake” shape and flowed as a monolayer in the microchannel. We used straight microchannels of five different widths (w = 50, 100, 150, 200, and 250 μm, respectively). We define a transverse confinement parameter $Cp≡wD$ as the channel width w divided by the droplet diameter D (=61.8 μm). Here, D is calculated assuming the drop had a cylindrical shape with a height equal to that of the channel (=30 μm), a cross section of a circle, and was un-deformed by neighboring droplets or the side walls. The confinement parameters corresponding to the five channel widths used were Cp = 0.81, 1.62, 2.43, 3.24, and 4.05, respectively. At the volume fraction of emulsion used, these channels could fit one, two, three, four, and five rows of droplets, respectively. For Cp = 0.81, 1.62, and 2.43, we conducted experiments at three viscosity ratios (λ = 0.81, 1.08, and 13.31). For Cp = 3.24 and 4.05, we conducted experiments at a single viscosity ratio (λ = 0.81) only. Here, the purpose of conducting experiments at these two Cp values is to show that the internal flow patterns inside the drops at Cp > 2.43 can be predicted by the flow patterns at Cp = 2.43, rather than studying the effect of the viscosity ratio. As such, only a single viscosity ratio was used. The flow of the emulsion was driven by a syringe pump at fixed volumetric flow rates (5–25 μl/h) but the capillary number, Ca, was kept constant for all experiments at Ca = 10⁻³. The capillary number is defined as Ca = μ_cGr/γ where G is the strain rate in the channel, r is the droplet radius, and γ is the interfacial tension. The Reynolds number based on the width of the channel, Re, varied from 10⁻² to 10⁻¹.

We built an in-house μPIV system²⁶ to investigate the flow fields within each droplet in the concentrated emulsion. Neutrally buoyant red fluorescent particles with a diameter of 0.5 μm (R500, Microgenics Corp.) were used as tracer particles within the droplets. Particles were seeded in the disperse phase prior to the droplet generation process. We excited the fluorescent particles using a 532 nm solid state continuous wave laser with a rated total power of 500 mW. In this setup, two separate lasers of the same wavelength at 200 mW and 300 mW were combined using a polarizing beam splitter. This polarizing beam splitter transmitted p-polarization and reflected s-polarization. As the lasers were p-polarized, in order to obtain an s-polarization, a half-wave plate could be used or one laser could be rotated by 90° with respect to the other. In our experiments, we chose the latter method. This combined light source was then coupled into an optical fiber cable to provide epi-illumination through an inverted microscope. Fig. 1(a) shows a scheme of this setup.

We acquired images at the mid-height of the microchannel using a high speed camera (Phantom v341) operated in a continuous mode. A frame rate of 1000 frames per second (fps) was used to ensure that there was an average of about 3–6 pixels of particle displacements between each frame. An exposure time of 990 μs was used to minimize the particle image blurring and to obtain images with a signal-to-noise ratio of at least 2. Fig. 1(b) shows a representative instantaneous raw image acquired. A minimum of 1000 images were captured for each experiment. A 20× microscope objective (NA = 0.42) was used for image acquisition. The current imaging setup yielded a pixel scale of 2.12 pixels/μm. The depth of focus of this setup was 3.1 μm while the depth of correlation was 12.9 μm.²⁷ 

Optical distortion can arise in our imaging system due to the curvature of the droplet surface and differences in the index of refraction between the continuous phase and the disperse phase. This issue can be resolved by using fluids with matching indices of refraction^5,6,13 but at the expense of limiting the types of emulsions we could test. In our experiment, we chose to test emulsions formed between aqueous solutions and fluorinated solvent HFE-7500 because they are used in many practical droplet microfluidics applications. For this combination of liquids, errors due to optical distortion were previously shown to be negligible.⁵ 

We used a two-step cross correlation technique for image processing. The two-step process improved both the accuracy and the resolution of the technique. The coarse and fine cross correlation windows used were 32 × 32 pixels and 16 × 16 pixels, respectively. The windows had a 50% overlap which produced vectors that were spaced every 8 pixels apart. In this case, the spacing between the vectors was 3.77 μm. For an average droplet size of 60 μm, about 15 velocity vectors were available in each direction to define the flow field within a given droplet.

We performed a vector post-processing step to detect and then replace any outlier velocity vector.²⁸ A given vector was considered as outlier when its value was seven times higher than the standard deviation of its neighbor in a 3 × 3 kernel. Once an outlier was determined, it was replaced with the median value based on its neighborhood velocity vectors in the same 3 × 3 kernel.

We acquired brightfield images in addition to fluorescent images during the experiments. These brightfield images showed clear boundaries between the continuous phase and the disperse phase. From these images, we could determine the centroid of each droplet within the field-of-view. The locations of the centroid were tracked from frame-to-frame to determine their displacements. We also monitored the shape and the size of the droplets to ensure that they did not change between frames. We then used the displacement information to shift our velocity vector fields before performing ensemble-averaging.

First, we investigate the flow of an emulsion at a fixed viscosity ratio (λ = 0.81) in a narrow channel where a single row of droplets spans the entire width of the channel (i.e., Cp = 0.81). For all the figures in this paper, the flow direction was from left to right. Figs. 2(a) and 2(b) show the instantaneous and ensemble-averaged velocity vector fields within the droplets. For all velocity vectors shown, we subtracted the velocity of the centroid of each droplet from the velocity measured directly by μPIV. The colors of the velocity vectors shown in the figures represent the magnitude of the normalized total velocity U_tot^∗, which is non-dimensionalized by the droplet centroid velocity. For U_tot^∗ only, the droplet centroid velocity was not subtracted to facilitate the visualization of the magnitude of the flow. Both the instantaneous velocity vector field and the ensemble-averaged results show that the velocity inside the droplet was faster near the wall than in the middle of the droplet.

Fig. 2(c) shows the normalized vorticity computed from the ensemble-averaged results. Normalized vorticity was non-dimensionalized by the droplet diameter and droplet centroid velocity. The vorticity results appear to show that there are two co-rotating structures on each half of the droplet. The top half had clockwise (blue) rotating structures while the bottom half had counterclockwise (red) rotating structures. Here, we use “top” and “bottom” to refer to the part of the droplet or channel wall located at y^∗ > 0 and y^∗ < 0, respectively.

To further confirm the presence of a pair of co-rotating structures on each half of the droplet, we used the Q-criterion method.²⁹ The Q-criterion represents the local balance between the strain rate and vorticity magnitude and is defined in the following equation:

where S and Ω are the symmetric and antisymmetric parts of the velocity gradient tensor ∇u, respectively (i.e., $S=12∇u+∇uT$ and $Ω=12∇u−∇uT$⁠). A rotational structure is defined to be present when Q > 0. Physically, it indicates that rotation dominates over strain. The Q-criterion does not indicate the direction of rotation, however. To show the direction, we colored the iso-contours of the Q-criterion by the signs of vorticity in Fig. 2(d). The non-dimensional Q-criterion, Q^∗, was calculated based on the non-dimensional velocity components and spatial coordinates. The Q-criterion results confirm the presence of a pair of co-rotating structures on each side of the droplet. These pairs of co-rotating structures have also been shown previously in large droplets¹⁵ but would be expected to combine into a single structure as droplet size decreases.

Second, we investigate the flow of an emulsion (λ = 0.81) in channels with two or three rows of droplets spanning the entire width of the channel at Cp = 1.62 and 2.43, respectively (Figs. 3 and 4). The instantaneous and ensemble-averaged results (Figs. 3(a), 3(b), 4(a), and 4(b)) show that the internal flow pattern differed from that at Cp = 0.81. For these cases, the velocity inside a given droplet was the fastest in the region adjacent to another row of droplets.

The vorticity plots in Figs. 3(c) and 4(c) show that: (1) The row of droplets adjacent to the wall had a single rotation. The droplets by the top wall had a counterclockwise rotation, while the ones by the bottom wall had a clockwise rotation. (2) All droplets away from the wall sandwiched by other rows of drops had two counter-rotating vorticity patterns. Each droplet had a clockwise rotation in the upper half of the droplet and a counterclockwise rotation in the lower half of the droplet. Again, to confirm the presence of vortical structures inside the droplets, we analyzed the Q-criterion. The results of Q-criterion (Figs. 3(d) and 4(d)) clearly show rotational structures, and the results are consistent with the vorticity results. We observe that the values of vorticity or Q-criterion at Cp ≥ 1.62 are higher than those at Cp = 0.81, indicating the higher strength for these rotational structures at Cp ≥ 1.62.

The observed internal flow pattern extends to droplets in channels at Cp = 3.24 and 4.05 (with 4 and 5 rows of drops across the width of the channel, respectively). Fig. 8 of the  Appendix shows the Q-criterion results for these two cases. Similar to those at Cp = 2.43, the droplets adjacent to the wall had a single rotational structure and the droplets away from the wall had a pair of counter-rotating structures.

Third, we study the effect of viscosity ratio on the flow patterns inside the droplets for different channel widths. The viscosity ratios investigated ranged over three orders of magnitude from λ ∼ 10⁻¹ to 10¹ as shown in Table I. Figs. 5(a)-5(c) show the ensemble-averaged Q-criterion results for three different viscosity ratios at Cp = 0.81. For all cases, the flow patterns were the same. The top half of each droplet had two co-rotating clockwise structures while the bottom half of each droplet had two co-rotating counterclockwise structures. The results for confinement parameters Cp = 1.62 and 2.43 were similar across different viscosity ratios tested and are not shown for brevity. Overall, we found that the flow patterns inside the droplets were independent of viscosity ratios for the conditions studied. Nevertheless, the magnitude of these rotations appeared to weaken as the viscosity ratio increased. Fig. 5(d) shows the normalized, relative streamwise velocity component, $urel*$⁠, in the middle of a droplet for the three viscosity ratios tested. Consistent with the trend in the magnitude of Q-criterion, the velocity gradient decreased as the viscosity ratio increased.

In our experiments, we made three important observations: (1) The directions of rotation in a single row of concentrated emulsion appeared to oppose some previous studies of a single droplet,^{5,6,12–14,16} while they also appeared to agree with others.^15,17 (2) The rotational patterns in a concentrated emulsion flowing in a channel having a confinement parameter of Cp ≥ 1.62 (fitting two or more rows of droplets) differed from those at Cp = 0.81 (fitting a single row of droplets). (3) The viscosity ratio did not play a role in the observed flow patterns.

To understand these observations, we first consider the boundary conditions in the flow of a single droplet. We assume a no-slip boundary condition at the wall between the continuous phase and the wall, as well as at the liquid-liquid interface between the disperse phase and the continuous phase.³⁰ To the best of our knowledge, there has been no experimental or computational evidence that showed otherwise in the liquids used in this paper, although some recent work has shown slip at a liquid-liquid interface between two immiscible polymers.^31–33 The tangential velocity due to the no-slip boundary conditions can then be expressed as

where u_∥w and u_∥o are the tangential velocities in water and in oil, respectively. Another boundary condition is the matching of shear stress at the w/o interface

where μ_w are μ_o are the dynamic viscosities of water and oil, respectively, and r is the radial direction. It is apparent from Eq. (4) that a mismatch in viscosity between the oil and water phases will cause a mismatch in the velocity gradient across the w/o interface. In general, decreasing the viscosity of the disperse phase relative to the continuous phase will increase the velocity gradient in the disperse phase. A large gradient tends to create increasingly pronounced rotational patterns in the droplets.⁵ 

The boundary conditions Eqs. (2)–(4) also hold for the flow of a concentrated emulsion with multiple rows of droplets. They do not predict the number and the direction of rotational structures, however. Further consideration of the velocity in the continuous phase relative to the disperse phase is necessary to understand the details of the flow patterns inside the droplets.¹⁰ 

Our results for Cp = 0.81 show that the top and the bottom halves of a droplet had clockwise and counterclockwise rotational structures, respectively. This trend is the opposite of some reports using single droplets in the dilute limit (where the drops were separated sufficiently far apart such that there were no hydrodynamic interactions among the drops¹¹) in a similar microfluidic system.^{5,6,12–14,16} We have verified that the directions of rotation did not change in our system using droplets in the dilute limit. A possible reason for the discrepancy is the difference in the transverse droplet mobility, β_T, defined here as the ratio of the velocity of the droplet to the average velocity of the thin film of continuous phase between either the two rows of droplets or between a row of droplets and the wall as shown in the following equation:

Droplet mobility plays an important role in determining the flow patterns inside the droplets. For example, if β_T < 1, we expect the flow inside the droplet adjacent to the walls to be faster than the central region of the droplet. The flow inside the droplets, as induced by the continuous phase, will lead to a clockwise (or counterclockwise) rotational structure in the top (or bottom) half of the droplet. If β_T > 1, the opposite trend is expected. It is, therefore, possible that β_T < 1 in our case. Previously, it was reported that the droplet mobility depends on the composition and concentration of the surfactant used.³⁴ As such, even though the experimental conditions and the choice of the liquids for the disperse and continuous phases were similar in our study and the previous study,⁵ the different surfactants used could lead to different β_T, and thus the different rotation directions observed. Indeed, Wang et al., who showed flow patterns inside a droplet similar to ours, stated that β_T < 1 for their system.¹⁷ 

In our experiments with multiple rows of droplets, we found that the droplets at the wall have one rotational structure, while the droplets away from the wall have two rotational structures irrespective of the viscosity ratio. Although the flow was 3D and complex, we attempt to explain qualitatively the number of rotational structures (single or two counter-rotating), as observed at the mid-height of the channel, by considering the differences in the velocities of the continuous phase sandwiching each row of drops: (1) If the velocities of the continuous phase are the same above and below a droplet, one should expect to see a pair counter-rotating rotational structure in a droplet. (2) If the velocities of the continuous phase are different (e.g., slow above the droplet and fast below the droplet), one might expect the flow inside the droplet to be driven primarily by the side with a higher velocity and lead to a single rotational structure only.

In our experiments, for the row of droplets adjacent to the wall, because of the requirement to satisfy the no-slip boundary condition at the wall, we expect that the average velocity of the continuous phase between the droplets and the wall to be slower than that between adjacent rows of droplets. It is thus plausible that such differences in velocity lead to a single rotational region in the drops by the wall as observed. On the other hand, for the row of droplets sandwiched by other rows of droplets, there should be insignificant difference in the average velocity of the continuous phase above and below these droplets. We thus expect the flow in the droplet as induced by the continuous phase to be symmetrical about the center of the droplet. Such flow leads to two counter-rotating vortices in the row of droplets away from the wall as observed in our experiments.

The flow of small droplets in the continuous phase appeared to support that the continuous phase between rows of droplets flowed faster than the droplets. We doped the concentrated emulsion with small droplets which acted as tracers for the continuous phase. Fig. 6 shows a tracer between adjacent rows of droplets (red circle) and a tracer between two neighboring droplets within the same row (blue circle). As can be seen, the tracer circled in red overtook the droplets, indicating that the velocity of the continuous phase between adjacent rows of droplets was faster than the velocity of the droplet. On the other hand, the tracer circled in blue did not move between two droplets within the same row, indicating that the velocity of the continuous phase between neighboring droplets in the same row had the same velocity as the droplets. We note that these observations are expected because of the direction of the applied pressure gradient. The pressure gradient should be zero orthogonal to the direction of flow. The rotational structures inside the droplets are thus primarily induced by the continuous phase flowing parallel to, rather than that orthogonal to, the direction of bulk flow.

To summarize our results, Fig. 7 shows a sketch of the relative velocity profiles for experiments at confinement parameters of Cp = 0.81, 1.62, and 2.43. The velocity profiles drawn inside the droplets are based on the measured experimental results. Although we could not perform μPIV in the continuous phase due to the large seeding particle size relative to the thickness of the thin film in the continuous phase, it is possible to speculate the velocity profiles in the continuous phase based on the no-slip and matching shear stress boundary conditions at the interface. At Cp = 0.81, the velocity profile in the continuous phase could approximate the velocity profile of a flow driven by a combination of pressure (due to the applied pressure gradient) and shear (due to the moving w/o interface and the stationary oil-wall interface). At Cp ≥ 1.62, in between the row of droplets and the wall, the velocity profile in the continuous phase is the same as the one for Cp = 0.81, even though the magnitude here should be less than that for Cp = 0.81 (by comparing the velocity inside the droplet close to the wall in Figs. 2(b) and 3(b)). In between rows of droplets, the continuous phase is at a higher average velocity than the disperse phase, and we propose that the velocity profile between adjacent rows of droplets could be parabolic. We note that these velocity profiles are proposed based on the measurements in the droplets and the boundary conditions only and should be ultimately verified by methods that can measure the flow in the continuous phase or by computational methods.

Previous studies reported that the viscosity ratio determines the flow patterns inside the drops while our experiment showed otherwise across the range of values tested. The lack of dependence on the viscosity ratio here was due to the fact that the flow of the continuous phase relative to the drops was not heavily dependent on the viscosity ratio at the capillary number (Ca = 10⁻³) tested. Fig. 9 of the  Appendix shows that at Cp = 1.62, similar to the case at Cp = 2.43 and λ = 0.81 (Fig. 6), the continuous phase between rows of droplets flowed faster than the disperse phase at a viscosity ratio λ = 13.31. The lack of dependence on the viscosity ratio is further supported by the similar shapes of the velocity profiles, which govern the flow patterns inside the droplets, for all different viscosity ratios as seen in Fig. 5(d). As discussed previously, the direction and the number of the rotational structures inside the droplets are determined by the droplet mobility and the differences in velocity of the continuous phase sandwiching the droplet. If neither change significantly with the viscosity ratio, the effect of viscosity ratio on the flow patterns inside the droplets is expected to be minimal. We note that it is possible that the viscosity ratio has more pronounced effects on the flow patterns than that reported here at other capillary numbers. The investigation at an increased range of capillary numbers and the comparison of our results with prior work are under current investigation.

To conclude, the μPIV experiments here showed that rotational structures inside the droplets of a concentrated emulsion always exist and are independent of the viscosity ratio for the conditions tested at a capillary number of Ca = 10⁻³. The direction of rotation of these structures depends on the droplet mobility. The number of rotational structures depends on the difference in the velocities of the continuous phase sandwiching the droplet. Such difference, in turn, depends on the confinement parameter, i.e., the number of rows of droplets with respect to the channel width and the location of the droplets in the channel. In terms of applications, our results imply that the degree of mixing could differ in droplets within a concentrated emulsion depending on the size of the drops relative to the size of the channel, as well as the position of the drops in the channel. If uniform mixing is desired, strategies to swap the position of the drops within the emulsion might be necessary.³⁵ Work is in progress to study the 3D internal flow profiles, as well as the effects of other parameters including the capillary number, the volume fraction of the emulsion, and the choice and concentration of the surfactant.

We acknowledge support from the National Science Foundation through the NSF CAREER Award No. 1454542. S.T. acknowledges additional support from the 3M Untenured Faculty Award and the Stanford Woods Institute for the Environment. We also acknowledge Lucas Blauch for his help with viscosity measurements.

Here, we present the flow patterns for confinement parameters Cp = 3.24 and 4.05 (Fig. 8), and the comparison between the velocity of the continuous phase and that of the disperse phase for λ = 13.31 (Fig. 9). We observe that the sketch of velocity profile shown in Fig. 7 can be extended to confinement parameters Cp = 3.24 and 4.05 (Fig. 8) with the same argument discussed earlier.