Droplet breakup against an isolated obstacle

Abstract

We describe combined experiments and simulations of single droplet breakup during flow-driven interactions with a circular obstacle in a quasi-two-dimensional microfluidic chamber. Due to a lack of in-plane confinement, the droplets can also slip past the obstacle without breaking. Droplets are more likely to break when they have a higher flow velocity, larger size (relative to the obstacle radius R), smaller surface tension, and for head-on collisions with the obstacle. We also observe that droplet–obstacle collisions are more likely to result in breakup when the height of the sample chamber is increased. We define a nondimensional breakup number Bk ∼ Ca that accounts for changes in the likelihood of droplet break up with variations in these parameters, where Ca is the Capillary number. As Bk increases, we find in both experiments and discrete element method (DEM) simulations of the deformable particle model that the behavior changes from droplets never breaking (Bk ≪ 1) to always breaking for Bk ≫ 1, with a rapid change in the probability of droplet breakup near Bk = 1. We also find that Bk ∼ S^4/3, where S characterizes the symmetry of the collision, which implies that the minimum symmetry required for breakup is controlled by a characteristic distance h ∼ R.

---

1 Introduction

Droplet formation is key to mixing two immiscible liquids to form an emulsion,¹ spread of some transmissible diseases via airborne droplets,² and inkjet printing.³ Furthermore, microfluidic devices are used to form droplets for lab-on-a-chip applications.⁴ Making droplets often involves starting with two liquids, adding energy by shaking, stirring, or otherwise flowing the two liquids, and thus mixing the fluids into large droplets of one fluid mixed into the other. Previous work has studied how fluid flow patterns at the individual droplet scale cause large droplets to break into smaller droplets. The simplest case is an isolated droplet in a shear flow.⁵ The surface tension of the droplet tries to minimize its surface area, and thus acts to maintain a spherical shape. Competing with surface tension, viscous stresses caused by the fluid shear flow try to stretch the droplet. For sufficiently fast flows, these viscous forces make the droplets deform or even tear themselves apart into smaller droplets. The capillary number is a dimensionless parameter that measures the ratio of viscous forces to surface tension forceswhere γ is the surface tension and the viscous forces are quantified by the continuous phase dynamic viscosity μ and a characteristic flow velocity v. For Ca > Ca_c, droplet breakup occurs, where the critical value Ca_c depends on the specific geometry of the fluid flow.^5–7 Prior work has studied droplet breakup in relatively simple microfluidic geometries, for example, in T-junctions,^6–12 constrictions where droplets drip from a nozzle,¹³ and narrow channels with an obstacle in the middle where the droplet wraps around both sides of the obstacle and then breaks.¹⁴ In these prior experiments, at large Ca viscous effects dominate causing increased droplet breakup. At small Ca, the droplets can deform, but they do not breakup.

Droplet breakup is less understood in more complex geometries, such as porous media. The flow of two immiscible fluids through porous media is important for numerous applications, such as petroleum extraction,^15,16 pharmaceutical manufacturing,^17,18 and agricultural and food production.¹⁹ Droplets moving through porous media are also crucial for understanding the flow of groundwater pollutants, such as PFAS.^20,21 One key feature of porous media is that droplets can be found in channels larger than their diameter, allowing them to assume complex shapes not observed in more confined geometries.

Many previous studies have considered flows through fully wetted porous media (e.g. oil fills the pore space) that is invaded by an immiscible fluid (water).^22–25 The effect of surface tension is enhanced if one of the fluids forms droplets, thus greatly increasing the interfacial area between the two fluids. A stream of droplets moves differently through a porous medium compared to invasion of a continuous phase fluid into a porous medium,²⁶ mainly due to the increase in interfacial area. Previous studies of droplets flowing through porous media have in many cases not considered deformable droplets and droplets that can break up.^27–29

In this article, we seek to understand droplet flow through a simplified porous medium, i.e., a single droplet interacting with a single obstacle, which is droplet flow through a porous medium in the limit of small droplet and obstacle densities (see Video S1). In Fig. 1, we show that the droplets can either move around the obstacle or wrap around it and break into two smaller droplets due to the flow. Several parameters determine whether the droplet breaks up or not. The first parameter influencing breakup is Ca. Faster flows have larger viscous forces that push the droplet against the obstacle, while surface tension prevents the droplet from deforming. Second, larger droplets are easier to deform and break, so the droplet size is an important parameter. A third parameter is the symmetry of the collision between the droplet and obstacle. Head-on collisions of the droplets with the obstacle lead to a higher probability of breakup; for example, compare Fig. 1(a), (b) with (e), (f). The symmetry of the collision was not relevant in prior experimental studies.^7,8,14,30 Droplets can also slide around the obstacle, requiring only small droplet deformations as shown in Fig. 1(e) and (f).

Fig. 1 (left) Images of droplets moving from left to right through the array of obstacles with collisions highlighted by red circles in each image pair. (a) and (b) A head-on collision (with symmetry S = 1.0) causes a large droplet with velocity v = 1.8 mm s⁻¹ to break up; (c) and (d) an asymmetric collision with S = 0.5 between a large droplet with v = 1.4 mm s⁻¹ causes droplet break up; and (e) and (f) an asymmetric collision with S = 0.5 between a small droplet with v = 1.1 mm s⁻¹ does not lead to break up. The time interval Δt between the left and right images in each row is determined by the velocity v so that the distance traveled is fixed, vΔt = 400 µm. (right) Microfluidics design of the sample chamber, which is ≈22 mm long. Droplets form at the middle top region and exit through the central channel into the wider region below. In the wider region, droplets collide with small obstacles before exiting the chamber at the bottom outlet.

We also develop mesoscale simulations to model droplet shape evolution and breakup due to stresses arising from interactions between the droplets and the obstacles and continuous phase fluid. For the simulations, we will employ the deformable particle model with surface tension³¹ supplemented by a geometric criterion for the onset of droplet breakup. An advantage of the deformable particle model is that it includes only a small number of physics-informed parameters that can be calibrated to the experimental results. The deformable particle simulations will allow us to map the regions of parameter space, such as the continuous phase fluid viscosity, droplet surface tension, droplet-to-obstacle size ratio, and collision symmetry, where droplet break up does and does not occur. Another advantage of the deformable particle simulations is that they enable exploration of regions of parameter space that are difficult to access experimentally. In particular, the simulations allow independent variation of the surface tension and continuous phase fluid viscosity, so that we can isolate and quantify their individual effects on droplet deformation and breakup. In contrast, in the experiments, the propensity for droplet breakup is tuned indirectly through changes in the flow conditions, effectively varying the capillary number.

The remainder of this article is organized as follows. In Section 2, we describe the experimental setup and method to generate, control, track, and quantify droplet–obstacle collisions. In Section 3, we describe mesoscale simulations of fluid-flow driven droplet–obstacle interactions using the deformable particle model with a geometric criterion for droplet breakup. In Section 4, we quantify when droplet break up occurs in the experiments and simulations by defining a nondimensional “breakup number” Bk, which is Ca multiplied by the ratio of the area of the droplet to that of the obstacle, and other geometric factors. Thus, droplet breakup is more likely at higher Ca, for droplets that are larger relative to the obstacle, and droplets that incur head-on collisions with the obstacle. Section 5 concludes with a summary of the results and promising future research directions. We include three appendices that provide details of the microfluidics device design (Appendix A), validation of the continuous phase fluid model used in the simulations (Appendix B), and independent variation of the surface tension and viscosity in the simulations to determine their effects on droplet breakup (Appendix C).

2 Experimental methods

The experimental flow cell, as shown in Fig. 1, consists of droplets, a continuous phase fluid flow, and obstacle pillars. The droplets are composed of water, a rhodamine dye added to saturation, and 1% tween-20 by mass, well above the critical micelle concentration. We did not observe depletion forces, since the surfactant is in the droplet phase and there was a constant flow of fresh silicon oil into the cell. The results were replicated using a grocery-store food coloring instead of rhodamine, suggesting that the properties of the specific dye are unimportant. The surfactant is included to provide a barrier to coalescence, although occasionally coalescence events are observed and these events are excluded from the data analysis. The continuous phase fluid phase is one of two silicon oils with kinematic viscosity ν_oil = 50 cSt or 100 cSt. Both oils have density ρ_oil = 960 kg m⁻³. These values yield a dynamic viscosity of μ = νρ = 480 or 960 Pa s depending on which silicon oil is used. The surface tension between the continuous phase fluid and droplets is γ ≈ 20 mN m⁻¹. The silicon oil is injected at a flow rate of 20–60 µL min⁻¹, driving the droplets with a velocity in the range of v = 0.5 to 2 mm s⁻¹. The water is injected at a flow rate in the range of 15 to 120 µL h⁻¹, which creates droplets using the pinch-off effect.¹³ The droplet diameters are in the range D₀ ≈ 100 to 600 µm, and the obstacle radii are in the range R = 60 to 120 µm. Based on these values, the Reynolds number, defined as Re = ρ_oilvR/μ_oil is ≤O(10⁻²) for all experiments, indicating that inertial effects are negligible.

The microfluidic chamber is made from polydimethylsiloxane (PDMS) created by pouring a degassed mixture of unsolidified polymer and curing agent (7.5 : 1 ratio by mass) onto a silicon wafer with the desired pattern, shown in Fig. 1, etched into it. This etching is performed using photolithography on an SU-8 surface, with a photomask ordered from ARTNET Pro, Inc. A profilometer measured the depth of the etching on the silicon wafer to be z = 85 µm ± 5 µm, which sets the sample chamber thickness z. Another sample chamber with z = 45 µm ± 5 µm is used for a subset of experiments to test the influence of z on the results. The PDMS is then allowed to solidify on this etching, either over a weekend or overnight with 70 °C heating. After the PDMS chambers are prepared and cut from the wafer, they are bonded to a microscope slide covered in PDMS using oxygen plasma cleaning, which allows it to serve as the “floor” of the chamber.³²

To avoid the need for 3D imaging, we created quasi-2D water-in-oil droplets with volume V that satisfies . Thus, the droplets are “pancake” shaped with small out-of-plane curvature, which is observable in Fig. 1 as a dark “border ring.” The outline of the droplets aids in detecting when one region contains a single concave droplet versus two convex droplets.

In untreated sample chambers, we observe that water droplets can stick to the PDMS chambers. To prevent adhesive forces, we coat the chambers with Aquapel which increases the hydrophobicity of the surfaces. The process works best when Aquapel flows through and is heated to 70 °C³³ for at least 20 minutes. We also find that fresh Aquapel creates a less hydrophobic surface than Aquapel stored in a degassed syringe for at least 24 hours. In addition, Aquapel spoils, visibly changes color, and loses its hydrophobicity after ≈2 weeks. Thus, we always use Aquapel within a few days after it arrives at the laboratory, but after storing overnight in the syringe.

The microfluidic chamber experiments are recorded using a LEICA DMIRB microscope at 60 frames per second with a ThorLabs DCC1645C – USB 2.0 CMOS camera at 640 × 512 resolution. We used a 1.6× objective (0.05 numerical aperture air lens) and a 1.5× zoom, resulting in videos with a scale of 5.28 µm per pixel. These videos are processed by separating the images into three regions: the continuous phase which is ignored, the pink centers of the droplets, and the black border of each droplet, which is assigned to the pink region that the border encircles. We can detect when two or more droplets are in contact, and we reject droplets that are in contact with other droplets during collisions with an obstacle. After segmentation and particle tracking,³⁴ we have data on the droplet area, velocity, position relative to obstacles, whether the droplet broke, and if so, the sizes of the daughter droplets after the collision for 5056 droplet–obstacle collisions.

The droplet–obstacle collisions are obtained under several experimental conditions. The standard parameters are the following: obstacle radius R = 85 µm, continuous phase dynamic viscosity μ_oil = 480 Pa s, and sample chamber thickness z = 85 µm. In addition to the standard set of parameters, we also vary each parameter one at a time, investigating μ_oil = 960 Pa s by changing the silicon oil, R = 60 µm and 120 µm, and z = 45 µm. Each experimental movie (Supplementary Video) contains multiple instances of droplet–obstacle collisions both with and without break up, depending on the experimental conditions. The velocity of droplets is modified by changing the continuous phase fluid pump rate, and the size of droplets is modified by controlling the relative flow rates between the continuous phase and droplet fluid.¹³ The angle of collision (later defined as the symmetry parameter in Section 4) is spontaneously varied by droplets within each experiment as they move through the arrays since neighboring droplets slightly modify each others' flow paths.

3 Simulation methods

3.1 Deformable particle model

We performed simulations of a single droplet colliding with a single obstacle using the deformable particle model, which can accurately model large deformations of capillary droplets flowing through confined geometries.³¹ In two-dimensions, the droplet is defined as a deformable polygon with N_v vertices, whose positions and velocities are the degrees of freedom of the system (see Fig. 2). The mass of the droplet is uniformly distributed among the vertices, and the motion of the vertices is determined by the droplet shape-energy function:

Fig. 2 Schematic of a droplet (shaded pink) with area A interacting with an obstacle (shaded gray) with radius R. The inset highlights the vertices that define the droplet surface in the deformable particle model. We also define vertex center-to-center distance r_ij and the arc length s_ij between vertices i and j. The droplet neck thickness, d_neck, is identified using the method described in Section 3.3.

The first term in eqn (2) imposes a harmonic energy penalty for changes in the droplet area A from the equilibrium value A_eq and k_a controls the fluctuations in the droplet area. This term represents the analog of the bulk modulus of the droplet in 2D. The second term in the shape-energy function imposes a harmonic energy penalty for deviations in the separations l_i between adjacent vertices i and i + 1 from the equilibrium length l_eq (which is also the diameter of each of the vertices) and k_l controls fluctuations in l_i. This term ensures that the vertices are evenly distributed on the droplet surface, preventing them from clumping when the droplet interacts with the obstacle. The factor of N_v in the numerator of the second term of eqn (2) makes U_s independent of N_v. The third term is the energy arising from line tension. We observe in the experiments that the droplet and the obstacle are coated by a thin layer of the continuous phase fluid (i.e. oil) and hence

where γ_2D ∼ γz is the line tension corresponding to the oil–water interface, and P is the droplet perimeter.

To prevent overlap between the droplet and the obstacle, we assume that the droplet interacts with the obstacle via pairwise, purely repulsive spring interactions between the obstacle and each of the droplet vertices:

where d_i is the distance between the center of the vertex i and the surface of the obstacle and ε_w sets the scale of the repulsive interactions. The Heaviside step function Θ(·) ensures that the force is non-zero only when vertex i overlaps with the obstacle. The total potential energy U of the droplet is given by the sum of the shape-energy function and droplet–wall interaction energy:

U = U_w + U_s. (5)

3.2 Modeling the effect of the continuous phase fluid

To mimic the experiments, the fluid flow in the simulations is pressure-driven. We neglect the effect of the droplet on the continuous phase fluid profile, but we include the drag force on each droplet vertex i from the fluid flow,where v⃑_i is the velocity of vertex i, v⃑_f(r⃑_i) is the velocity of the fluid at vertex i, μ is the fluid viscosity, and is the diameter of the undeformed droplet. In eqn (6), the factor of 1/N_v ensures that the drag force on the droplet is independent of the number of vertices. To model the flow field, we usewhich enforces no-slip boundary conditions on the surface of the obstacle for the radial v_f,r and angular v_f,θ components of v⃑_f. v_f∞ is the velocity of the fluid far from the obstacle. The coordinate system is defined so that the origin is at the center of the obstacle, the horizontal axis is aligned with the fluid flow, r is the distance from the origin, and θ is the angle relative to the horizontal axis. Eqn (7) and (8) represent Stokes' solution for creeping flow around a circular obstacle, which provide an accurate description of the fluid flow close to the obstacle. We observe that when the droplets are much smaller than the obstacle, this choice for the flow field yields an accurate trajectory for the droplet around the obstacle, since the droplets effectively act as passive tracers, advected by the continuous-phase flow with no appreciable disturbance to the flow field. For droplets that are much larger than the obstacle, the droplet can distort the flow field. Nevertheless, we find that the deformation of large droplets is insensitive to the form of the flow field in these studies that are in the low Reynolds number regime, provided that the no-slip boundary condition at the obstacle is satisfied (see Appendix B).

3.3 Mesoscale modeling of droplet breakup

When a droplet in a shallow microfluidic channel encounters an obstacle, the confinement forces it to deform around the obstacle and produce a neck that thins as the in-plane deformation increases. The neck thickness decreases through a combination of viscous drainage and capillary pressure gradients set by the channel height. Once the neck reaches a critical thickness at which capillary stresses can no longer sustain a connected interface across the confined gap, small perturbations at the interface grow and the neck ruptures, resulting in break up of the droplet. In our 2D deformable particle simulations, we use a simple, geometric criterion to determine when a droplet breaks up.

We define the neck thickness d_neck as the smallest of the center-to-center distances r_ij between every pair of vertices i and j subject to the constraint that the length s_ij of the shortest arc joining them (as measured along the perimeter of the droplet) is larger than a threshold d_min-sep = 0.28D₀ (see Fig. 2). Minimizing r_ij subject to this constraint on s_ij prevents unphysical breakup events. Specifically, in the limit d_min-sep → 0, daughter droplets form with an unrealistically small number of vertices. Conversely, as d_min-sep → D₀, it becomes increasingly difficult to identify a vertex pair whose shortest connecting arc length exceeds the threshold. Consequently, d_min-sep must fall within a small range. If the neck thickness falls below a smaller threshold d_neck < d_min-neck at any point during the simulation, we break the droplet into two daughter droplets along the line defining the neck. The value of d_min-neck can be calibrated to experimental results for the likelihood of droplet breakup. We find that setting d_min-neck = 0.17D₀ results in the best match to the experimental results.

3.4 Definitions of the model parameters

In the simulations, we use D₀ as the characteristic length scale and t₀ = D₀/v_f∞ as the characteristic time scale. Using these along with the droplet mass M, we define the dimensionless viscosity [small mu, Greek, tilde] = μt₀D₀/M, line tension [small gamma, Greek, tilde]_2D = γ_2Dt₀/(Mv_f∞), edge-length spring constant k̃_l = k_lD₀/γ_2D, and area spring constant k̃_a = k_aD₀⁴/(Mv_f∞²). We impose fluid incompressibility of the droplet by setting k̃_a > 10⁴ and fix k̃_l/[small gamma, Greek, tilde]_2D < 0.05 so that the line tension energy dominates the perimeter spring energy. In a Hele–Shaw geometry³⁵ such as that used in the experiments, the deformation of a pancake-shaped droplet is governed by a balance between viscous pressure variations induced by confinement and the restoring capillary pressure associated with interfacial curvature. Because the flow is pressure driven, the dominant viscous forcing acting on the droplet originates from lubrication pressure within the thin wetting film that separates the droplet interface from the confining plates. This lubrication pressure p satisfies³⁶where x denotes the flow-direction and y is the direction normal to the flow. If the characteristic thickness of the lubrication layer is denoted by δ, a simple scaling argument gives uswhere v is the droplet–velocity. The restoring pressure scale associated with surface tension is set by the inverse radius of curvature of the interface,

The degree of droplet deformation is therefore controlled by the ratio of these two pressure scales,

where Ca is the experimentally measured capillary number. Classical lubrication theory predicts that, for a Hele–Shaw gap of height z, the lubrication film thickness scales as δ ∼ zCa^2/3.^37–39 In the present experiments, Ca = O(10⁻³), which implies

Since D₀ and z are of the same order-of-magnitude in our setup, it follows that

This result is consistent with the strong amplification of the viscous stresses induced by confinement. In the simulations, we define a capillary number Ca_sim ≡ [small mu, Greek, tilde]v_f∞/[small gamma, Greek, tilde]_2D ∼ (μv_f∞/γ)(D₀/z). Here, we used the fact that γ_2D ∼ γz, with the additional length scale D₀ arising because the line tension has different physical dimensions than the surface tension. The viscous force in our simulation acts on each vertex and therefore, the viscosity coefficient being used here corresponds to the local viscous forcing that gives rise to the amplification of viscous stresses induced by confinement between two plates. Thus, we can write:

since z ∼ D₀ are the same order of magnitude. Thus, the capillary number in the 2D simulations corresponds to the effective capillary number Ca_eff (not Ca), and a direct comparison between experiments and simulations requires the conversion

Determining the precise value of the conversion factor (δ/D₀)² is nontrivial, as it depends on details of the lubrication film that are not directly accessible in the experiments. We therefore treat this value as a conversion factor that relates the simulation capillary number Ca_sim and the experimentally measured Ca. Its value is determined by requiring that the onset of droplet breakup occurs at the same value of the capillary number in the simulations and experiments, with all other parameters held fixed. This procedure gives

Ca = 9.8 × 10⁻⁵Ca_sim, (17)

which is consistent with eqn (14). Henceforth, when discussing the simulation results, we report Ca from eqn (17).

3.5 Equations of motion

The equations of motion for vertex i of the deformable particle arewhere m = M/N_v is the mass of each vertex and M is the total mass of the droplet. We integrate eqn (18) using a modified velocity-Verlet numerical integration scheme with time step Δt = 10⁻⁴t₀. The droplet is initialized as a regular polygon of N_v sides with area A_eq. We then set the edge lengths equal to their equilibrium values . At the start of the simulation, we place the droplet at rest at a distance of 5D₀ from the center of the obstacle to allow it to reach v_f∞ before it collides with the obstacle.

4 Results

We seek to understand the key physical properties that determine droplet breakup as a droplet collides with a circular obstacle. All of the data concerning the droplet (i.e. the area, symmetry, and velocity) is collected at the frame of the first contact between the droplet and the obstacle. An important parameter governing droplet breakup is the offset between the droplet trajectory and the obstacle. Collisions are commonly characterized using the impact parameterwhere v⃑ is the velocity of the center of mass of the droplet and r⃑_cc is the vector joining the center of the obstacle to the center of mass of the droplet. This definition does not capture the droplet shape when it collides with the obstacle. Thus, we define the symmetry parameter S that better captures the shape of the droplet while retaining the information contained in the impact parameter. S is defined using the observed droplet area relative to a centerline at the time of the collision as shown in Fig. 3. The centerline passes through the center of the obstacle and is parallel to each droplet's center of mass velocity when that droplet first makes contact with the obstacle. The symmetry parameter is defined using the two subareas A_large and A_small above and below the centerline, where the total area is A = A_large + A_small:

Fig. 3 Illustration of the definition of the symmetry parameter S. When a droplet first contacts an obstacle (shaded gray), a dividing line (dashed line) is drawn through the droplet to form two regions with areas A_large (shaded blue) and A_small (shaded yellow). The dividing line is parallel to the center of mass velocity vector of the droplet (large white arrow) and passes through the center of the obstacle.

S = 1 indicates a perfectly symmetric collision, and S = 0 indicates a collision where all of the droplet is on one side of the obstacle. Droplet breakup requires S > 0, otherwise the droplet slides around the obstacle. S has an uncertainty of ±0.088 due to noise in measuring the areas. We emphasize that S is measured at the moment the droplet first contacts the obstacle and we do not focus on the subsequent evolution of A_large and A_small.

In Section 4.1, we will show that we can identify distinct regions in the parameter space where droplets break up versus where droplets do not break up using the symmetry parameter S and the combined control parameter CaÃz̃, where à = A/R² is the nondimensional droplet area and z̃ = z/R is the nondimensional chamber thickness. One can also nondimensionalize the geometrical parameters using the diameter of the undeformed droplet D₀, resulting in the nondimensional quantities: z/D₀ and A/D₀². However, since the droplet area is approximately constant during the experiments (until it breaks), A ≈ πD₀²/4, which makes A/D₀² ≈ π/4. In this case, the important nondimensional parameters become z/D₀ and R/D₀ (i.e. the obstacle size relative to the droplet size). We prefer to non-dimensionalize using R instead of D₀ because the droplets are often deformed slightly from circular at the moment they first touch an obstacle, making direct measurements of A easier in experiments, rather than D₀.

We find that the boundary separating the breakup and no breakup regimes follows a power-law relation in S and CaÃz̃. In Section 4.2, we use the power-law scaling relation to construct a dimensionless breakup number Bk, where Bk ≫ 1 indicates that the droplet will break up and Bk ≪ 1 indicates that the droplet will not break up. In Section 4.3, we focus on droplets that undergo breakup and determine how the ratio of the daughter droplet areas depends on S. Finally, in Section 4.4, we measure the minimum neck thickness that a droplet can sustain without breaking up and use these measurements to validate the break up model employed in the deformable particle model simulations.

4.1 Parameter regimes for droplet breakup

In Fig. 4, we separate the experimental data for droplet–obstacle collisions into two sets in the parameter space of S versus CaÃz̃: (left) droplets that break up and (right) droplets that do not break up. These results show that droplet break up involves a tradeoff between S and CaÃz̃; for example, a droplet that collides with the obstacle off-center (S < 1) can be made to break up by increasing the velocity, since Ca ∼ v. The two clouds of experimental data are best separated by a power-law scaling of the form:

S_c = α(CaÃz̃)^β, (21)

where the prefactor α ≈ 0.083 and the power-law scaling exponent β ≈ −0.74. Since S_c scales with Ca, droplets are more likely to break up when they have a large velocity, are immersed in a fluid with large viscosity, and have small surface tension.

Fig. 4 The symmetry parameter S plotted versus CaÃz̃ for all droplet collisions separated into those for which the droplets (left) break up, (right) do not break up, and both overlaid (middle). The separating line (black dashed line) is given in eqn (21) with power-law scaling exponent β = −0.74.

Eqn (21) shows that larger droplets are also easier to break, i.e. a larger object feels more pressure when pushed against a smaller object; the smaller object is effectively sharper for larger Ã. Our data have à ∼ O(1). When à ≪ 1, we expect different behavior than that shown in Fig. 4 because the obstacle will appear more like a wall, which is a situation that has been previously studied.⁷

Additionally, droplets in thicker chambers are more likely to break up. Between the parallel plates of the sample chamber, the droplet interface curves between the top and bottom surface with a radius of curvature ∼z/2.⁴⁰ This small curvature increases the Laplace pressure inside the droplet, ΔP ∼ γ/z. Thus, thicker droplets are “softer” and more easily deform and break.

In the experiments, the surface tension is fixed and the viscosity is varied between two values, 50 cSt and 100 cSt. Consequently, in experiments we varied the capillary number primarily by changing the droplet velocity. In contrast, the simulations allow independent changes of the surface tension and viscosity. By varying these parameters independently, we show in Appendix C that S_c ∼ ([small mu, Greek, tilde]/[small gamma, Greek, tilde])^β with β ≈ −0.72, which is in agreement with the experimental results as shown in Fig. 5. Note that because the simulations are two-dimensional, the chamber thickness z̃ cannot be varied. Hence, we set z̃ = 1 for the simulation data in Fig. 5. We performed a sweep over a range of values of S and CaÃz̃, noting whether the droplet breaks up or not for each (S,CaÃz̃) pair. We then hold CaÃz̃ fixed and increase S until we find the largest value S_c for which the droplet does not break up. Repeating this for all of the CaÃz̃ values in our simulation gives us a set of points that form the separation boundary between regions of no-breakup and break up. We find that this separation curve is best described by the power law

S_c = 0.089(CaÃz̃)^−0.72, (22)

which is in excellent agreement with the separation curve obtained from experiment, as can be seen in Fig. 5 and by direct comparison to eqn (21).

Fig. 5 Results from the deformable particle model (DPM) simulations (stars) showing droplets that break up (upper right; green) and do not break up (lower left; violet) as a function of S and CaÃz̃ overlaid on the experimental data from Fig. 4. We set z̃ = 1 for the simulation data. The best fit lines that separate the droplets that break up and do not break up have slopes −0.72 (black dashed; simulations) and −0.74 (orange solid; experiments).

4.2 Breakup number

To determine the power-law exponent in eqn (21), we use the maximum likelihood method.⁴¹ We start by assuming an equation for the dividing line of the form:

S_c = α(CaA^ε₁R^ε₂z^ε₃)^β, (23)

where β is the power-law exponent. Given that A, R, and z all have dimensions related to length, we have the constraint 2ε₁ + ε₂ + ε₃ = 0 such that eqn (23) is nondimensional. We next define the “breakup number”

Bk = α^1/βCaA^ε₁R^ε₂z^ε₃/S^1/β, (24)

to quantify the distance from the dividing line S_c. With these definitions, droplets that break up will have large values of Bk, and droplets that do not break up will have Bk → 0.

We apply the maximum likelihood method for separating the experimental data for which droplets break up and droplets do not break up because the experimental data are noisy near S_c. As shown in Fig. 4, some droplet–obstacle collisions are observed on the “wrong” side of the dividing line. There are examples in the data set close to S_c where for similar experimental parameters, some droplets break and others do not. One reason for this behavior could be small uncertainties in the experimental measurements of velocity or S stemming from the image analysis. However, these quantities are defined at the moment a droplet first contacts the obstacle, and a more likely concern is that droplets can and do change their speed and direction of motion as they interact with the obstacle. In addition, other droplets near the given droplet–obstacle collision may influence the flow of oil around the droplet, again changing the behavior of the droplet as it interacts with the obstacle. Therefore, for the experimental data, we consider breakup as a probabilistic process, and Bk = 1 corresponds to the case where droplets are equally likely to break or not.

To mathematically implement the breakup probability, we define the breakup characteristic k as:

We then define the probability of observing outcome k, using x = (Bk_α,β)^1/w (where the subscripts indicate that Bk is a function of α and β) and the function

which means that for Bk_α,β ≫ 1 break up is likely [P(k = 1) → 1] and for Bk_α,β ≪ 1 no break up is likely [P(k = 0) → 1]. P(x,k) is a sigmoid function of ln(Bk), where w is the width of the sigmoid. We then define the likelihood L of observing the data by a product over all droplet–obstacle collision events i as:

Terms in the product are close to 0 for data points with a low probability of the actual outcome, and close to 1 when the prediction matches the actual outcome. Therefore, incorrectly chosen α and β dramatically reduce the likelihood due to the contributions from many incorrectly assigned points. In contrast, optimal α and β will create many more matches and a much larger total likelihood. w accounts for the width of the region where droplets have an intermediate chance of break up. Maximizing the logarithm of the likelihood⁴¹ yields power-law exponent β = −0.74 and α ≈ 0.083.

We also use the maximum log-likelihood method to calculate the exponents ε₁, ε₂, and ε₃ in eqn (24) that are subject to the constraint 2ε₁ + ε₂ + ε₃ = 0. We find that S_c ∼ A¹z¹R⁻³ as shown in Fig. 4. As an additional check, we allowed the power-law scaling exponent ε_μ for μ in eqn (24) to vary. Maximizing the log-likelihood returned ε_μ = 1, confirming the Ca dependence for S_c. We determine the uncertainty of α, β, and w using the bootstrapping method (with 100 samples of 2528 randomly selected data points). We find α = 0.083 ± 0.006, β = −0.74 ± 0.03, and w = 0.098 ± 0.006.

We now present an argument to interpret the observed value of the power-law exponent β ≈ −0.74 in eqn (21). We will show that an exponent of β = −3/4 corresponds to a conserved geometric quantity during the collision, the minimum height of A_small. Consider droplet–obstacle collisions with fixed velocity, geometry, and fluid properties, but varying A. In this case, droplet–obstacle collisions on the dividing line in Fig. 4 are described by S_c = (A/Λ)^β, where all of the other parameters are subsumed into Λ with units of area. We can then use eqn (20) to rewrite S_c in terms of the two droplet subareas:

Solving for A_small gives:

A comparison of droplet–obstacle collisions for two droplet areas is shown in Fig. 6, where the yellow-shaded circular segments have area A_small. The figure defines the length scale h and droplet diameter D₀, and in the limit h ≪ D₀ one can show A_small ∼ D₀^1/2h^3/2. We can relate the diameter and area of the droplet using D₀ ∼ A^1/2 and express eqn (29) for A_small in terms of h:

h ∼ A_small^2/3A^−1/6. (30)

Fig. 6 Image of two droplets colliding with an obstacle (shaded gray). The droplets have equal velocities (v = 1.0 mm s⁻¹), identical chamber geometries, the same fluid parameters, but one droplet is twice the radius of the other. The smaller droplet has A = 0.02 mm². The two droplets have Bk = 1, and the smaller area below the symmetry line (in yellow) of the larger droplet is 1.4 times that of the smaller droplet. The droplet radius D₀/2 and height of the smaller circular section h are indicated.

Next, consider the case h ∼ A⁰, where h has a fixed value that is independent of A for a droplet–obstacle collision at S_c. With these assumptions, we substitute eqn (29) into eqn (30) to obtain

A⁰ ∼ (A^1+β)^2/3A^−1/6. (31)

To ensure that h is independent of A, we must set β = −3/4. The power-law exponents β = −0.74 obtained in the experiments and β = −0.72 obtained in the simulations are in excellent agreement with this scaling analysis. Using eqn (29) and (30), we can also show how the small circular segment height h scales with the obstacle radius R, sample thickness z, and Ca:

Thus, our results in Fig. 4 and 5 suggest that the separating curve S_c depends on the droplet area A such that h is constant, which yields the scaling behavior in eqn (32).

We now consider experimental measurements of Bk. Using eqn (24), we can write the breakup number as:

Bk = 28CaÃz̃S^4/3, (33)

where Bk < 1 indicates that the droplets are not likely to break up, and Bk > 1 indicates that the droplets are likely to break up. In Fig. 7, we plot the fraction N_break/N_total of droplet–obstacle collisions that yield droplet break up as a function of Bk. The probability is a sigmoid function with a narrow width: the range of Bk over which the probability grows from 5% to 95% is only a factor of ≈2 in Bk. The results in Fig. 7, which extend over six orders of magnitude in Bk, illustrate that we have identified the key parameters that determine droplet break up during single droplet–obstacle collisions.

Fig. 7 The fraction N_break/N_total of droplet–obstacle collisions in the experiments (black solid line) that yield droplet break up plotted versus Bk. Eqn (26) (violet solid line) and the step function Θ(Bk) (pink solid line) are also shown. (inset) A close-up of N_break/N_total in the main panel near Bk = 1.

4.3 Area ratio of daughter droplets

As noted above, all of the experiments are in the low Reynolds number regime, Re ≤ O(10⁻²). In this regime, we find that nearly all droplet break up events result in the formation of only two daughter droplets. If we denote the area of the smaller daughter droplet as A₁ and that of the larger daughter droplet as A₂, we can define the daughter droplet area ratio 0 < A₁/A₂ < 1, which depends on Ca and à of the original droplet. We contrast these areas with A_small and A_large, which are defined for each droplet at the first contact between the droplet and the obstacle. In off-centered collisions (S < 1), we find that A₁ < A_small and A₂ > A_large (with A₁ + A₂ = A_small + A_large from mass conservation), because the larger lobe partially slides along the obstacle prior to breakup, dragging the smaller lobe with it, and thereby redistributing area before break up. However, when Ca ≫ 1, we expect that the timescale for mass redistribution is much larger than the timescale for droplet deformation (and subsequent break up). In this limit, A₁ → A_small and A₂ → A_large, and thus

In Fig. 8, we plot the daughter droplet area ratio A₁/A₂ versus S for droplet–obstacle collisions that yield droplet break up from experiments and simulations. We observe that all of the simulation data and nearly all of the experimental data occur below this limiting curve. However, given the presence of continuous phase fluid flow driven by other nearby droplets, it is possible for some of the droplet fluid to move into the smaller lobe such that A₁/A₂ exceeds the limiting prediction. Other confounding influences on A₁/A₂ in the experiment are droplets that change velocity while they are in contact with the obstacle, and droplets that are pre-deformed upon the initial contact with the obstacle (due to the motion of other nearby droplets).

Fig. 8 The daughter droplet area ratio A₁/A₂ for droplet–obstacle collisions that yield droplet break up plotted verses the symmetry S of the collision in experiments (blue filled circles). We also show eqn (34) (dotted black line), where the droplet breaks along the collision axis with no fluid exchange, A₁ = A_small, and A₂ = A_large. The simulation data (solid lines) are shaded according to the value of Ca from Ca = 0.0006 to 0.0045. Two linear clusters of data are observed; these correspond to experiments where droplets collide reproducibly over a narrow range of S, leading to A₁/A₂ ∼ S.

We also show A₁/A₂ versus S in Fig. 8 from the simulations for a range of Ca. For the simulations, we find that A₁/A₂ approaches S/(2 − S) in the large-Ca limit in agreement with eqn (34). In particular, A₁/A₂ = 1 when S = 1. In this case, the two lobes are identical and there is no fluid redistribution across the lobes during break up for all Ca. For A₁/A₂ to reach S/(2 − S) for S < 1, Ca must progressively increase as S decreases.

4.4 Droplet neck thickness

We also examine the minimum neck thickness d_neck (as defined in Section 3.3) attained during droplet–obstacle interactions, considering both droplets that undergo breakup and those that do not. We display the evolution of the neck of a breaking droplet from experiments in Fig. 9.

Fig. 9 Experimental snapshots of a droplet breaking up. The images are taken at intervals of 1/60 s to show the neck thinning and breaking. The image analysis places the exact moment of breakup between the 3rd and 4th frames from the left.

The neck is experimentally measured in the single video frame immediately before the frame where the droplet has broken into two daughter droplets (for example, the middle image of Fig. 9). Image analysis identifies the perimeter pixels of the droplet. We assign each pixel an index i ordered sequentially along the perimeter. We then compute the distance between each pixel pair (i,j). In the space (i,j) the line i = j forms the global minimum as the distance from each pixel to itself is zero. There is also a local minimum away from the diagonal that corresponds to the minimum distance between two pixels with indices i and j far apart (along the contour of the perimeter), which identifies the neck thickness and the two specific pixels i and j defining the neck.

In Fig. 10, we show the frequency distribution of the thinnest neck d_neck observed in the experiments. We find that droplets which undergo breakup never appear in their final recorded frame with neck thicknesses larger than d_neck ≈ 70 µm (green histogram), while droplets that do not break always have d_neck ≳ 140 µm (violet histogram). Based on studies of instability-mediated fluid filament breakup, we expect that there is a critical neck thickness d_min-neck below which the droplet will always break up.⁴² The histogram of neck thicknesses for non-breaking droplets sets a ceiling on d_min-neck: droplets with d_neck < 140 µm always break up. It is possible if we observed a much larger number of droplets, we would find rarer cases of non-breaking droplets with slightly thinner necks that would reduce the ceiling. The value of d_neck-thin is expected to depend on factors such as the droplet size, chamber thickness, and possibly the surface tension; we have insufficient data to determine this dependence so we simply estimate d_min-neck ≈ 140 µm.

Fig. 10 Frequency distribution for the minimum neck thickness d_neck of droplets that (green) do and (violet) do not break up in experiments. (insets) Example droplets (i)–(iii) for regions (i)–(iii) of d_neck marked on the histograms, with fluid flow from right to left, obstacle radius R = 85 µm (shaded pink), and neck outline (orange solid line). The data shown are for the standard experimental parameters given in Section 2.

The green histogram (corresponding to droplets that break) depends on the frame rate of the camera (60 frames per second). Droplet neck thinning likely has a rapid stage at the end⁴² but nonetheless is expected to be a continuous process, so a higher speed camera will reduce the observed d_neck as seen one video frame before breakup. The green histogram can thus be interpreted as an observation that some droplets with a neck thickness up to ∼70 µm can thin to zero on a time scale of 1/60 s or faster. The location and breadth of this histogram for breaking droplets reflects the finite frame rate of the camera, which captures droplets at different stages of neck-thinning prior to breakup. As a result, we expect that increasing the frame rate would decrease the width and location of the green histogram.

In the simulations, we set d_min-neck = 0.167D₀. Using the experimentally observed median diameter of droplets that undergo breakup, D₀ = 481 µm, this corresponds to d_min-neck ≈ 80.3 µm, which is consistent with the experimental bound.

5 Conclusions

We carried out coordinated experiments and simulations of quasi-two-dimensional deformable droplets flowing in microfluidic chambers that collide with cylindrical pillars. For some conditions, droplet–obstacle collisions give rise to droplet break up. For others, droplets collide with the obstacle, but slide around it, and do not break up. Break up is influenced by the interplay between the viscous surface stresses that deform the droplet and the surface tension that resists deformation, which is quantified by the capillary number Ca in eqn (1). We find that droplet break up also depends on several geometrical parameters. For example, larger droplets relative to the obstacle are more likely to break up. Droplets in thinner sample chambers are less likely to break up, as a result of the larger internal Laplace pressure that resists droplet deformation required for break up. In addition, we show that the symmetry of the droplet trajectory relative to the center of the obstacle influences droplet break up, i.e. head-on collisions with S = 1 maximize the likelihood of break up. The results for droplet break up can be collapsed using the nondimensional breakup number Bk, which is proportional to CaÃz̃S^4/3, and à and z̃ are the droplet area and sample chamber thickness normalized by the obstacle radius R. Bk accurately predicts the likelihood of droplet breakup in experiments over six orders of magnitude of Bk, with a narrow region of Bk near Bk = 1 where the droplet break up probability transitions from zero to one with increasing Bk.

The experimental results are also verified through discrete element method simulations using the deformable particle model with line tension in 2D. We demonstrate that incorporating a geometric criterion for droplet break up related to the neck thickness is sufficient to capture the droplet break up statistics. In particular, we find that the characteristic symmetry for which a droplet undergoes breakup is S_c = α(CaÃz̃)^β, which is the same functional form observed in experiments, and the power-law exponent β ≈ −0.72 is in close agreement with experiments. Further, the geometric break up criterion implemented in the simulations is consistent with the experimental observations that droplets whose neck thickness is above a threshold value do not break up.

Our work gives additional insight into the daughter droplets formed after break up. We observe that single obstacle-induced droplet breakup events predominantly produce two daughter droplets, and the resulting daughter droplet area ratio is constrained by the initial symmetry of the collision. For off-centered collisions, interfacial sliding and lobe interactions prior to break up redistribute area, causing the daughter droplet areas to deviate from those defined at first contact and the daughter droplet area ratio to be below the upper bound A₁/A₂ = S/(2 − S). In the limits of large capillary number and large droplet area, area redistribution is negligible and the daughter droplet area ratio approaches the upper bound, which is verified by the simulations.

These results suggest several promising future research directions. First, we can extend our studies of droplet break up to three dimensions, where thin droplet necks are completely unstable due to surface tension. Second, droplets can coalesce, as well as break up.^43,44 In future studies, we will investigate droplet coalescence in a microfluidic porous medium. With a fundamental understanding of both droplet break up and coalescence, we will be able to predict the resulting droplet size distributions as they move through the medium.⁴⁵ As droplets flow though a microfluidic porous medium composed of obstacles, we expect that a steady-state size distribution will be reached where break up events are balanced by coalescence events.⁴⁶ We can also consider the problem of fluid wetting; in real porous media flows, droplets are not perfectly dewetted to the media.⁴⁷ This effect may have significant effects on fluid-flow induced droplet break up and coalescence.

Author contributions

Conceptualization: all authors; methodology: David J. Meer (experimental). Shivnag Sista (computational); investigation: David J. Meer (experimental). Shivnag Sista (computational); formal analysis: David J. Meer, Eric R. Weeks (experimental). Shivnag Sista, Mark D. Shattuck, Corey S. O'Hern (computational); data curation: David J. Meer (experimental). Shivnag Sista (computational); writing – original draft: David J. Meer (experimental). Shivnag Sista (computational); writing – review & editing: all authors; visualization: all authors; supervision: Mark D. Shattuck, Corey S. O'Hern, Eric R. Weeks; funding acquisition: Mark D. Shattuck, Corey S. O'Hern, Eric R. Weeks.

Conflicts of interest

There are no conflicts to declare.

Data availability

Data for this article, including experimental microscopy videos and corresponding image analyses, computational data and code to recreate the figures are available in the UNC dataverse https://doi.org/10.15139/S3/CZM4AN.

Supplementary information (SI) is available, containing a movie from the experiment. See DOI: https://doi.org/10.1039/d5sm01266j.

Appendices

Appendix A: Microfluidic chamber design

In this appendix, we describe individual sections of the microfluidics chamber pattern labeled 1–7 in Fig. 11:

Fig. 11 We label regions of the microfluidics chamber pattern. The design was converted into a photomask after it was designed in openSCAD.

(1) (yellow) Inputs for the continuous oil-based phase (left) and the droplet water-based phase (right). The design within the input regions catches dust and debris brought in with the fluid.

(2) (light green) The central region where droplet formation occurs via pinch-off of the water phase.¹³ The droplets then continue down the central channel, and the oil travels through all three channels.

(3) (green) Extra oil input regions to modify the symmetry of the droplet–obstacle collisions. These inputs were not used for any of the experiments in this article, and once they were filled with oil, they did not affect the flow of the droplets.

(4) (light blue) The oil phase expands in cross-sectional area, and slows down. This main expansion occurs before the droplet phase enters the main chamber, which ensures a minimal velocity gradient between the main observing region and central droplet channel.

(5) (dark blue) Structural obstacles to ensure minimal wall effects, a consistent cross-sectional area of the flow, and structural stability to prevent the chamber from collapsing.

(6) (violet) The main imaging area, where obstacles have radius R.

(7) (orange) The exit region, where the wall separation decreases and flow velocity increases. Therefore, no droplet behavior is recorded in this region, though the obstacles still have radius R.

Appendix B: Comparing droplet shape and motion from simulations and experiments

In this appendix, we show that the deformable particle model simulations (described in Section 3.2) accurately recapitulate the flow trajectories of the droplets and the droplet shapes as they interact with the obstacles, especially in the limit of small droplets relative to the obstacles D₀/R ≪ 1. In Fig. 12(a)–(d), we provide experimental images of a droplet with D₀/R = 0.78 that does not break up with the simulation results overlaid on the experimental results. The bottom panel shows the trajectory of the droplet around the obstacle for both the experiments and simulations. In particular, we plot the y-coordinate of the droplet center of mass y_CM versus the x-coordinate x_CM, where the droplet flow is in the x-direction. Note that for the simulations, the center-of-mass coordinates are obtained from the x- and y-coordinates of the vertices x_i and y_i:

Fig. 12 (a)–(d) Snapshots of a small droplet (D₀/R = 0.78) flowing past an obstacle (R = 85 µm) in experiments with snapshots from the deformable particle model simulations overlaid in red. The droplet moves from panels (a)–(d) in ∼0.43 seconds. The bottom panel shows the y-component of the droplet center of mass y_CM plotted versus the x-component x_CM for the simulations (solid line) and experiments (filled circles).

In the bottom panel of Fig. 12, we find strong qualitative agreement between the droplet trajectories in experiments and simulations, which indicates that the simplified fluid model is appropriate for the experimental studies.

As the droplet size increases relative to the obstacle size, its motion induces perturbations to the surrounding flow field, leading to distortions of the flow streamlines and higher-order hydrodynamic effects on the droplet dynamics. Resolving these effects would require two-way coupled computational fluid dynamics simulations with moving deformable interfaces. Given that we are focusing on the low Reynolds number regime, viscous dissipation dominates and inertial effects are suppressed, which limits the extent to which droplet-induced flow perturbations affect the droplet motion.

We can quantify the deformation of the droplet using the dimensionless shape parameter: , where is the droplet perimeter, for a circle, and indicates shape-deformation due to interactions of the droplet with the obstacle. In Fig. 13(a)–(f), (Video S2) we show experimental images of a droplet that is larger than the obstacle with D₀/R = 2.5 as it interacts with the obstacle. The bottom panel shows a plot of the shape parameter versus x_CM of the droplet for both experiments and simulations. Once again, we find qualitative agreement between the simulations and experiments. This level of agreement indicates that, despite its simplifying assumptions, the fluid model in the simulations can capture droplet shape evolution with sufficient accuracy to reliably predict the transition between break up and no-breakup regimes.

Fig. 13 (a)–(f) Snapshots of a large droplet (D₀/R = 2.5) flowing past an obstacle (R = 85 µm) in experiments with snapshots from the deformable particle model simulations overlaid in red. The droplet moves from panels (a)–(f) in ∼1.05 seconds. The bottom panel shows the droplet shape parameter A plotted versus the x-component of the droplet center of mass x_CM for the simulations (solid line) and experiments (filled circles).

Appendix C: Verifying the role of capillary number in droplet breakup

In this appendix, we describe results from simulations of a single droplet colliding with a single obstacle, while tuning the normalized droplet line tension [small gamma, Greek, tilde] and fluid viscosity [small mu, Greek, tilde] separately. For each [small gamma, Greek, tilde], we vary S and [small mu, Greek, tilde] to determine S_c([small mu, Greek, tilde]) that separates droplet break up from no break up behavior. In the top panel of Fig. 14, we show that the separating curve scales as S_c ∼ [small mu, Greek, tilde]^β with β ≈ −0.72 for each value of [small gamma, Greek, tilde]. In the bottom panel, we show that all of the S_c curves collapse when the horizontal axis is scaled as [small mu, Greek, tilde]/[small gamma, Greek, tilde]. The collapse demonstrates that the combined influence of fluid viscosity and droplet surface tension on droplet break up is governed by the capillary number Ca = vμ/γ.

Fig. 14 (top) The boundary in the symmetry parameter S_c that separates break up and non-breakup outcomes for droplet–obstacle collisions as a function of the normalized fluid viscosity [small mu, Greek, tilde] for several values of [small gamma, Greek, tilde]. The dashed lines have slope β = −0.72. (bottom) The same data in the top panel except the horizontal axis is rescaled as [small mu, Greek, tilde]/[small gamma, Greek, tilde], showing that S_c collapses onto a single curve: S_c ∼ Ca^β (dashed line) with β = −0.72, which matches the experimental observations of Section 4.1.

Acknowledgements

We thank the laboratories of J. Burton, C. Roth, and S. Urazdhin for use of their equipment in microfluidics fabrication. We thank T. Brzinski for useful discussions. We acknowledge support from NSF Grant No. CBET-2333224 (D. J. M. and E. R. W.), No. CBET-2333222 (S. S. and C. S. O.), and No. CBET-2333223 (M. D. S.).

References

1. J. Bibette, L. Calderon and P. Poulin, Rep. Prog. Phys., 1999, 62, 969–1033.
2. W. C. K. Poon, A. T. Brown, S. O. L. Direito, D. J. M. Hodgson, L. L. Nagard, A. Lips, C. E. MacPhee, D. Marenduzzo, J. R. Royer, A. F. Silva, J. H. J. Thijssen and S. Titmuss, Soft Matter, 2020, 16, 8310–8324.
3. P. Tabeling, Introduction to Microfluidics, Oxford University Press, Oxford, Reprint edn, 2010.
4. S. Selimovic, F. Gobeaux and S. Fraden, Lab Chip, 2010, 10, 1696–1699.
5. J. J. M. Janssen, A. Boon and W. G. M. Agterof, Fluid Mech. Transp. Phenom., 1994, 40, 1929–1939.
6. D. R. Link, S. L. Anna, D. A. Weitz and H. A. Stone, Phys. Rev. Lett., 2004, 92, 054503.
7. T. Fu, Y. Ma, D. Funfschilling and H. Z. Li, Chem. Eng. Sci., 2011, 66, 4184–4195.
8. W. L. Cheng, R. Sadr, J. Dai and A. Han, Biomed. Microdevices, 2018, 20, 72.
9. D. H. Wankawala and P. K. Mondal, Sādhanā Acad. Proc. Eng. Sci., 2023, 48, 274.
10. D. H. Wankawala, A. Q. Shen and P. K. Mondal, Proc. R. Soc. A, 2025, 481, 20240594.
11. A. Tiribocchi, M. Durve, M. Lauricella, A. Montessori, D. Marenduzzo and S. Succi, Phys. Fluids, 2023, 35, 063106.
12. S. S. Dhruvkumar, H. Wankawala and P. K. Mondal, Europhys. Lett., 2025, 150, 23001.
13. A. S. Utada, A. Fernandez-Nieves, H. A. Stone and D. A. Weitz, Phys. Rev. Lett., 2007, 99, 094502.
14. S. Protière, M. Z. Bazant, D. A. Weitz and H. A. Stone, Europhys. Lett., 2010, 92, 54002.
15. S. Betancur, L. Quevedo and C. M. Olmos, Can. J. Chem. Eng., 2024, 102, 2583–2607.
16. M. G. Gerritsen and L. J. Durlofsky, Annu. Rev. Fluid Mech., 2005, 37, 211–238.
17. K. Saitoh, K. Koichi, F. Yabiku, Y. Noda, M. D. Porter and M. Shibukawa, J. Chromatogr. A, 2008, 1180, 66–72.
18. C. A. Browne and S. S. Datta, Proc. Natl. Acad. Sci. U. S. A., 2024, 121, e2320962121.
19. K. Schroen, C. Berton-Carabin, D. Renard, M. Marquis, A. Boire, R. Cochereau, C. Amine and S. Marze, Micromachines, 2021, 12, 863.
20. S. Chen and B. Guo, Water Resour. Res., 2023, 59, e2023WR034664.
21. J.-W. Liu, K.-H. Wei, S.-W. Xu, J. Cui, J. Ma, X.-L. Xiao, B.-D. Xi and X.-S. He, Sci. Total Environ., 2008, 756, 144142.
22. S. S. Datta, H. Chiang, T. S. Ramakrishnan and D. A. Weitz, Phys. Rev. Lett., 2013, 111, 064501.
23. S. S. Datta and D. A. Weitz, Europhys. Lett., 2013, 101, 14002.
24. S. S. Datta, J.-B. Dupin and D. A. Weitz, Phys. Fluids, 2014, 26, 062004.
25. S. S. Datta, T. S. Ramakrishnan and D. A. Weitz, Phys. Fluids, 2014, 26, 022002.
26. A. Parrenin, S. Kooij, C. J. M. Van Rijn and D. Bonn, Phys. Fluids, 2024, 36, 103323.
27. M. Izaguirre and S. Parsa, Soft Matter, 2024, 20, 3585–3592.
28. E. Benet, A. Badran, J. Pellegrino and F. Vernerey, J. Membr. Sci., 2017, 535, 10–19.
29. E. Benet, G. Lostec, J. Pellegrino and F. Vernerey, Phys. Rev. E, 2018, 97, 042607.
30. J. K. Nunes, S. S. H. Tsai, J. Wan and H. A. Stone, J. Phys. D: Appl. Phys., 2013, 46, 114002.
31. Y. Cheng, B. F. Lonial, S. Sista, D. J. Meer, A. Hofert, E. R. Weeks, M. D. Shattuck and C. S. O'Hern, Soft Matter, 2024, 20, 8036–8051.
32. B. K. Barnes, H. Ouro-Koura, J. Derickson, S. Lebarty, J. Omidokun, N. Bane, O. Suleiman, E. Omagamre, M. J. Fotouhi, A. Ogunmolasuyi, A. Dominguez, L. Gonick and K. S. Das, Am. J. Phys., 2021, 89, 372–382.
33. S. Nawar, J. K. Stolaroff, C. Ye, H. Wu, D. T. Nguyen, F. Xin and D. A. Weitz, Lab Chip, 2020, 20, 147–154.
34. J. C. Crocker and D. G. Grier, J. Colloid Interface Sci., 1996, 179, 298–310.
35. Y. Ling, J.-M. Fullana, S. Popinet and C. Josserand, Phys. Fluids, 2016, 28, 062001.
36. H. Schlichting, Boundary-Layer Theory, McGraw-Hill, New York, 6th edn, 1968.
37. I. Shukla, N. Kofman, G. Balestra, L. Zhu and F. Gallaire, J. Fluid Mech., 2019, 874, 1021–1040.
38. F. P. Bretherton, J. Fluid Mech., 1961, 10, 166–188.
39. G. Taylor and P. G. Saffman, Q. J. Mech. Appl. Math., 1959, 12, 265–279.
40. K. W. Desmond, P. J. Young, D. Chen and E. R. Weeks, Soft Matter, 2013, 9, 3424–3436.
41. P. Nelson, Physical Models of Living Systems, Chiliagon Science, 2nd edn, 2022.
42. X. Shi, M. P. Brenner and S. R. Nagel, Science, 1994, 265, 219–222.
43. N. Bremond, A. R. Thiam and J. Bibette, Phys. Rev. Lett., 2008, 100, 024501.
44. A. Lai, N. Bremond and H. A. Stone, J. Fluid Mech., 2009, 632, 97–107.
45. A. G. Yiotis, L. Talong and D. Salin, Phys. Rev. E:Stat., Nonlinear, Soft Matter Phys., 2013, 87, 033001.
46. E. Amstad, S. S. Datta and D. A. Weitz, Lab Chip, 2014, 14, 705–709.
47. S. K. Woche, M.-O. Goebel, R. Mikutta, C. Schurig, M. Kaestner, G. Guggenberger and J. Bachmann, Sci. Rep., 2017, 7, 42877.

---

Footnote

† Contributed equally and share first authorship.

---