Collective behavior of crowded drops in microfluidic systems

Droplet microfluidics, in which microdroplets serve as individual reactors, has enabled a wide range of high-throughput biochemical processes. Unlike solid wells typically used in current biochemical assays, droplets are subject to instability and can undergo breakup, especially under fast flow conditions. Although the physics of single drops has been studied extensively, the flow of crowded drops or concentrated emulsions—where droplet volume fraction exceeds $\sim 80\%$—is relatively unexplored in microfluidics. In this paper and the related invited lecture from the 74th Annual Meeting of the American Physical Society's Division of Fluid Dynamics, we describe the collective behavior of drops in a concentrated emulsion by tracking the dynamics and the fate of individual drops within the emulsion. At the slow flow limit of the concentrated emulsion, we observe an unexpected order, where the velocity of individual drops in the emulsion exhibits spatiotemporal periodicity. This periodicity is surprising from both fluid and solid mechanics points of view. We show the phenomenon can be explained by treating the emulsion as a soft crystal undergoing plasticity, in a nanoscale system comprising thousands of atoms as modeled by droplets. Our results represent a type of collective order which can have practical use in on-chip droplet manipulation. As the flow rate increases, the emulsion transitions from a solidlike to a liquidlike material, and the spatiotemporal order in the flow is lost. At the fast flow limit, droplet breakup starts to occur. We show that droplet breakup within the emulsion follows a probability distribution, in stark contrast to the deterministic behavior in classical single-drop studies. In addition to capillary number and viscosity ratio, breakup probability is governed by a confinement factor that measures drop size relative to a characteristic channel length. The breakup probability arises from the time-varying packing configuration of the drops. Replacing surfactant with nanoparticles as droplet stabilizers suppresses breakup and increases the throughput of droplet processing by $>300\%$. Strategic placement of an obstacle suppresses breakup by $>{10}^3$-fold. Finally, we discuss recent progress in computation methods for recapitulating the flow of concentrated emulsions.

Figure 1

Summary of our work on concentrated emulsions.

Figure 2

(a) Microparticle image velocimetry (μPIV) experimental setup for flow measurements of droplets within a concentrated emulsion in a microfluidic channel. (b) A representative instantaneous fluorescent raw image was acquired using the μPIV technique. Representative three rows of droplets: (c) Instantaneous and (d) ensemble-averaged velocity vector fields. The droplet centroid velocity is subtracted from all the velocity vectors shown. The velocity vectors are colored by the magnitude of normalized total velocity $U_{\mathrm{tot}}^{*}$, which is normalized by the droplet centroid velocity. (e) Color contours of normalized vorticity field ω* computed using the ensemble-averaged results. Vorticity is normalized by the droplet diameter and droplet centroid velocity. Red and blue indicate counterclockwise and clockwise vorticity, respectively. (f) Color contours of normalized $Q$ criterion $Q^{*}$ computed using the ensemble-averaged results. Colors are shown based on the results of vorticity. The spatial coordinates of $x^{*}$ and $y^{*}$ are normalized by the width of a single-drop-wide channel ($=50\mu \mathrm{m}$). Reproduced from Ref. [74], with the permission of AIP Publishing.

Figure 3

Sketch illustrating flow inside a single row, two rows, and three rows of droplets within a concentrated emulsion spanning the entire width of the channel. The velocity profiles inside the droplets (black lines) are based on the experimental results, while those in the continuous phase (shown in the insets) are proposed based on the no-slip and matching shear stress boundary conditions at the water-oil interface (gray lines). All velocity profiles shown are relative to the droplet velocity. Red and blue arrows show the counterclockwise and clockwise vorticity, respectively. Reproduced from Ref. [74], with the permission of AIP Publishing.

Figure 4

(a) T1 droplet rearrangement process. Droplets labeled #1 and #3, which were not in contact initially, converge at the end of the T1 process. Droplets labeled #2 and #4, which were in contact initially, diverge at the end of the T1 process. (b) A representative instantaneous fluorescent raw image was acquired using the microparticle image velocimetry (μPIV) technique. (c)–(h) Phase-averaged flow field in the rearrangement zone $N=3–2$. (c)–(e) Velocity vector field and (f)–(h) $Q$ criterion for three time points within one T1 event at $t=0$, 0.5, and 1.0T, where T is the period of a T1 event. The T1 event occurs at an angle of $+{60}^{\circ }$ with respect to the flow direction. Droplets involved in the T1 event are outlined by solid lines while the ones not involved are outlined by dotted lines. Flow is from left to right. Velocity vector and its color are shown relative to the droplet velocity and the color also represents the phase-averaged normalized total velocity, ${\tilde{U}}_{\mathrm{tot}}^{*}$. Phase-averaged normalized $Q$ criterion ${\tilde{Q}}^{*}$ is colored by the phase-averaged normalized vorticity ${\tilde{w}}^{*}$. (i) Nondimensional circulation Γ* vs nondimensional time $t_{\mathrm{N}32}^{*}$ for the four droplets involved in the rearrangement process. Here, circulation is nondimensionalized by circulation inside the droplets at the constriction, and time is nondimensionalized by the duration of the T1 event for this rearrangement zone (i.e., from $N=3$ to 2). Every other data point is shown for clarity purposes. Solid lines are curve fits to the data points and are guides for the eyes only. The error in the data points estimated using the propagation of uncertainty analysis is about $\mathrm{\Delta }{\mathrm{\Gamma }}^{*}=\pm 0.37$. (j) Nondimensional average velocity inside a droplet $U_{\mathrm{drop}}^{*}$ vs nondimensional time $t_{\mathrm{N}32}^{*}$ for the four droplets involved in the rearrangement process. Here, the average velocity inside a droplet is nondimensionalized by the droplet centroid velocity at the constriction. Reproduced from Ref. [78], with the permission of AIP Publishing.

Figure 5

(a) Microscope image of the concentrated emulsion in the tapered microfluidic channel. (b) The time-averaged total velocity of drops at various locations in the channel. The velocity was normalized to the maximum total velocity measured at the constriction. The position of each marker represents the position where the droplet velocities were averaged. The color of the marker represents the magnitude of droplet velocity at that location. Note that the markers overlap and appear as continuous lines which reflect the trajectories of the drops. The rearrangement zones are indicated by blue arrows. (c) and (d) Time-averaged $x$ component of droplet velocity $\overline{u_x}$ across the width of the channel measured at the locations marked by the blue and green boxes in (b). The velocities shown were normalized to ${\overline{u_x}}_{,\max }$, the maximum $x$ component of velocity measured in the constriction. (e) and (f) Kymographs of instantaneous droplet velocities in the rearrangement zone $(N=7–6)$ and nonrearrangement zone $\left(N=7\right)$, respectively. $N$ is the number of rows of drops across channel width. The position ($t$, $y$) of each marker represents the time $t$ and the $y$ position of a drop in the channel when the velocity of the drop was measured. The color of the marker represents the magnitude of the total velocity of the drop normalized to the maximum velocity measured at the constriction. (g) Instantaneous velocity vectors of the drops in the channel. Reproduced from Ref. [101].

Figure 6

(a) A T1 event involves one pair of diverging drops (#2 and #4) and one pair of converging drops (#1 and #3). Snapshots of the T1 cascade which is equivalent to a dislocation glide. Three T1 events were highlighted: event $1(t=0–0.13\mathrm{s})$, where the diverging (converging) drops were colored in red (green); event $2(t=0.13–0.25\mathrm{s})$, where the diverging (converging) drops were colored in purple (yellow); event $3(t=0.25–0.38\mathrm{s})$, where the diverging (converging) drops were colored in red (green). (b) and (c) Burgers circuits (pink lines) of the dislocations in the crystal at $t=0.19\mathrm{s}$ [also indicated by the red box in (a)], and at $t=0.58\mathrm{s}$, respectively. The pink arrow indicates the Burgers vector b. The green line indicates the slip plane, and the green arrow indicates the direction of motion of the dislocation. The black arrows indicate the stress acting on the crystal. Spatial and temporal periodicity of dislocation dynamics: (d) $x/a$ positions of T1s in the channel as a function of time. The mean spacing $s$ between the seven sets of slips shown is $s=5.5$. $a$ is the diameter of one drop. Each trace represents the $x/a$ positions of the T1s within one rearrangement zone where $N$ rows of drops reduces to $N-1$. (e) Fluctuation in the $x/a$ position of the dislocations as a function of channel width expressed in terms of the number of rows of drops $N$. The height of the gray error bar represents the maximum fluctuation $x$/$a$ measured (also plotted in the inset). (f) Scaled $y$ position ($y^{\prime }$) of the dislocations as a function of time. For each trace, the maximum (or minimum) $y^{\prime }$ value represents the top (or bottom) of the channel wall within one rearrangement zone where $N$ rows of drops reduces to $N-1$. The $y$ positions of the T1s oscillate with period $T$. For both (a) and (c), the positions plotted are that of the drops with five neighbors for convenience. (g) The period $T$ scales linearly with ${\left(N\text{–}\frac{1}{2}\right)}^2$. The periods are calculated using the fast Fourier transform of the data in (c). Gray markers represent data from four independent experiments at the same flow conditions. Black markers represent the average from these experiments. The black line is the best linear fit to the average. Reproduced from Ref. [101].

Figure 7

Voronoi tessellation during a T1 event at (a) $Ca=4.9\times {10}^{-7}$ and (b) $Ca=2.1\times {10}^{-3}$. The green lines represent the Voronoi cell edges. The red line represents the growing edge shared by converging drops during the T1 events. The blue dots represent the centers of 5-edged Voronoi cells. (c) Log-log plot showing $T^{\prime }$ as a function of $Ca$. Each data point consists of a batch of $T^{\prime }$ measurements of $>200$ T1 events within one rearrangement zone at a given flow rate $Q$. The vertical and the horizontal error bars represent the standard deviation in $T^{\prime }$ and $Ca$ for all $>200$ measurements within one rearrangement zone at one $Q$ value, respectively (see Fig. 5). The three dashed lines that follow the data points have logarithmic slopes of −0.14, −0.61, and −0.37. The two vertical dashed lines indicate the transition between between Regimes 1 and 2 at $C{a}_{1-2}$, and between Regimes 2 and 3 at $C{a}_{2-3}$. (d) $d{T}^{\prime }$, $T^{\prime }$, and $De$ as a function of $Ca$. $C{a}_{\mathrm{S}-\mathrm{L}}$ represent the solidlike-to-liquidlike transition when $De=1$. Reproduced from Ref. [92].

Figure 8

Trajectories of 2000 droplets advected through the channel for (a) $1.1\times {10}^{-7}<Ca<1.6\times {10}^{-6}$ (within Regime 1), (b) $8.5\times {10}^{-6}<Ca<6.3\times {10}^{-5}$ (within Regime 2), and (c) $7.9\times {10}^{-4}<Ca<8.5\times {10}^{-3}$ (within Regime 3). Each marker represents the instantaneous position of the centroid of one droplet. The insets are time averaged, $x$ component of droplet velocity ($\overline{U_x}$) across the width of the channel calculated at the locations marked by the red and green boxes, respectively. $\overline{U_x}$ is normalized to ${\overline{U_x}}_{\max }$, the maximum $x$ component of velocity observed in the constriction. The rearrangement zones are indicated by blue arrows in (a) and (b). Snapshots of successive T1 events at (d) $\langle Ca\rangle =2.5\times {10}^{-7}$ (Regime 1), (e) $\langle Ca\rangle =1.2\times {10}^{-5}$ (Regime 2), and (f) $\langle Ca\rangle =1.4\times {10}^{-3}$ (Regime 3). The blue arrows indicate the direction along which successive T1 events occur. For T1 events along the 60° slip plane, the converging (diverging) drops are colored in green (red). For T1 events along the direction of imposed flow, the converging (diverging) drops are colored in yellow (purple). For both (d) and (e), each sequence shows 5 successive T1 events along the slip plane. For (f), the sequence shows 5 successive T1 events along the slip plane and the flow direction, respectively. Reproduced from Ref. [92].

Figure 9

(a) Images of the flow of a concentrated emulsion consisting of 40 pL drops in a channel with a constriction. The average velocities of the flow at the constriction were (i) $v=6.2\mathrm{cm}/\mathrm{s}$ and (ii) $v=64.8\mathrm{cm}/\mathrm{s}$. The blue and red boxes correspond to regions where we tracked drops upstream and downstream. (b) Histograms showing the frequency of occurrence of droplet sizes upstream (blue) and downstream (red) of the constriction, with plots (i) and (ii) corresponding to images (i) and (ii) in (a). (c) Breakup fraction as a function of capillary number at different drop sizes and constriction geometries. (d) Breakup fraction as a function of capillary number for emulsions with different viscosity ratios. A1–A8 in (c) represent experiments with different combinations of droplet size and channel constriction confinement (defined by constriction width and height), and B1–B5 in (d) represent experiments with various viscosity ratios. (e) Breakup fraction as a function of the product of capillary number, viscosity ratio, and confinement factor. The dashed lines are for visual guides only. $Ca=\frac{{\mu }_{c}Gr}{\sigma }$, where ${\mu }_c$ is the dynamic viscosity of the continuous phase, and $G$ is the strain rate in the constriction. See our previous work for detail [20, 79]. Reproduced from Refs. [20, 79] with permission from the Royal Society of Chemistry.

Figure 10

(a) Scheme of the microchannel used in our experiments with two possible outcomes of a pair of droplets entering the constriction. The initial offset $\mathrm{\Delta }x$ determines if breakup occurs. Only 3 droplets of the concentrated emulsion are shown. (b) Breakup regime map of normalized droplet pair offset as a function of $Ca$. (c) Snapshots of drop pairs in (i) Region I, (ii) Region II (break), (iii) Region II (no break), and (iv) Region III at $Ca=0.014$. The leading-edge offset is (i) 3.54 μm, (ii) 28.32 μm, (iii) 28.32 μm, and (iv) 34.22 μm, respectively. Reproduced from Ref. [102], with the permission of AIP Publishing.

Figure 11

(a) Nanoparticle (NP) droplet breakup fraction as a function of the product of capillary number ($Ca$) and confinement factor ($cf$). The black dashed line is a visual guide. The gray dashed line is the visual guide adapted from surfactant droplets [see Fig. 9]. (b) Effect of viscosity ratio on breakup of NP emulsion droplets. Breakup fraction as a function of capillary number at different viscosity ratios λ: A5 ($\lambda =0.78$), B1 ($\lambda =1.38$), B2 ($\lambda =2.90$), and B3 ($\lambda =17.04$). The dashed line is for visual guide only. Reproduced from Ref. [10], with the permission of AIP Publishing.

Figure 12

(a) Schematic diagram of the channel geometry. (b) Microscopic image of the emulsion flowing in the channel. The image was converted to a binary image to increase the contrast of droplet borders. (c) Breakup fraction as a function of volume fraction at different $Q_{\mathrm{total}}$ or $Ca$, and (d) as a function of $Q_{\mathrm{total}}$ or $Ca$ at various volume fractions. The curves are fits from the data. In (c) and (d), each data point and the corresponding error bar are obtained from $N=4000–21000$ drops. (e) Droplet throughput as a function of volume fraction at different breakup fractions. See our work for a detailed description of data processing [107]. Reproduced from Ref. [107], with the permission of AIP Publishing.

Figure 13

(a) Schematic diagram of the channel geometry. (b) Drop breakup fraction in the presence of an obstacle normalized to that in the absence of the obstacle (β) as a function of normalized obstacle-to-constriction distance, $x$/$D_d$ for various obstacle diameters ($D_o$) and capillary numbers $Ca$. Datasets F#.5 (where $\#=1–6$), G#.5 (where $\#=1–10$), I#.5 (where $\#=1–10$) have $D_o=90$, 112.5, and 150 μm, respectively, but the same $Ca=10.8\times {10}^{-3}\pm 0.5\times {10}^{-3}$. Datasets G#.3 (where $\#=4–7$), G#.4 (where $\#=1–9$), and G#.5 (where $\#=1–10$) have the same $D_o=112.5\mu \mathrm{m}$ but different $Ca=7.6\times {10}^{-3}\pm 0.6\times {10}^{-3},9.3\times {10}^{-3}\pm 0.6\times {10}^{-3}$, and $10.8\times {10}^{-3}\pm 0.5\times {10}^{-3}$, respectively. See Ref. [118] for details of the experiments. (c) Regime map of β at various $l_{\mathrm{pow}}/D_d$ and $l_{\mathrm{ow}}/D_d$ values. All data shown here were measured from channels with a half angle $\theta ={30}^{\circ }$ and $Ca$ ranging from $5.46\times {10}^{-3}$ to $13.1\times {10}^{-3}$. Regions 0 to V are defined in the map. The red, black, and green markers indicate cases with $\beta >1.0$, $0.1<\beta \le 1.0$, and $\beta \le 0.1$, respectively. The green + marker indicates the case with minimum $\beta =1.2\times {10}^{-3}$. The light gray area indicates geometries that were inaccessible experimentally as $l_{\mathrm{ow}}$ cannot exceed $l_{\mathrm{pow}}$. (d) Snapshots of the emulsion flowing in representative channel geometries in regions 0 to V. The circle filled with blue hashed lines is the obstacle. The scale bar is 100 μm. (e) Time series showing that drops alternate in an orderly manner and enter the constriction one at a time in a channel with an optimally placed obstacle. Reproduced from Ref. [118].

Figure 14

(a) Double emulsion formation through hierarchical flow-focusing microchannel obtained via the volume-of-fluid (VOF)-finite volume approach [163]. (b) Numerical simulation of droplets assembly in a divergent channel. Simulation made with aphros software [143]. (c) Boundary integral method (BIM) simulation of droplets passing through a constriction [137]. (d) Droplets generation in a microfluidic T-junction via VOF simulations [141].

Figure 15

(a) A snapshot of the droplets in a concentrated emulsion, flowing from left to right. The upper panel corresponds to experiment, and lower to lattice Boltzmann method (LBM) simulation [164]. (b) Soft flowing crystals formation within a microfluidic focuser: LBM color-gradient simulations vs experiments [121]. (c) Deformation and breakup dynamics of droplets within a tapered channel: LBM color-gradient simulations vs experiments [161]. (d) Jet-breakup in soft granular materials flowing in a microfluidic focuser: LBM color-gradient simulations vs experiments [166]. (e) Free-energy LBM simulations of dense emulsions flowing in a straight channel [167].