# Flow and clogging of capillary droplets

Yuxuan Cheng<EMAIL_ADDRESS>Department of Physics, Yale University,
New Haven, Connecticut, 06520, USA.  Benjamin F. Lonial  Department of
Physics, Emory University, Atlanta, GA 30322, USA.  Shivnag Sista  Department
of Physics, Yale University, New Haven, Connecticut, 06520, USA.  David J.
Meer  Department of Physics, Emory University, Atlanta, GA 30322, USA.  Anisa
Hofert  Department of Physics, Emory University, Atlanta, GA 30322, USA.  Eric
R. Weeks  Department of Physics, Emory University, Atlanta, GA 30322, USA.
Mark D. Shattuck  Benjamin Levich Institute and Physics Department, The City
College of New York, New York, New York 10031, USA.  Corey S. O’Hern
<EMAIL_ADDRESS>Department of Physics, Yale University, New Haven,
Connecticut, 06520, USA.  Department of Mechanical Engineering and Materials
Science, Yale University, New Haven, Connecticut, 06520, USA.  Program in
Computational Biology and Bioinformatics, Yale University, New Haven,
Connecticut, 06520, USA.

###### Abstract

Capillary droplets form due to surface tension when two immiscible fluids are
mixed. We describe the motion of gravity-driven capillary droplets flowing
through narrow constrictions and obstacle arrays in both simulations and
experiments. Our new capillary deformable particle model recapitulates the
shape and velocity of single oil droplets in water as they pass through narrow
constrictions in microfluidic chambers. Using this experimentally validated
model, we simulate the flow and clogging of single capillary droplets in
narrow channels and obstacle arrays and find several important results. First,
the capillary droplet speed profile is nonmonotonic as the droplet exits the
narrow orifice, and we can tune the droplet properties so that the speed
overshoots the terminal speed far from the constriction. Second, in obstacle
arrays, we find that extremely deformable droplets can wrap around obstacles,
which leads to decreased average droplet speed in the continuous flow regime
and increased probability for clogging in the regime where permanent clogs
form. Third, the wrapping mechanism causes the clogging probability in
obstacle arrays to become nonmonotonic with surface tension $\Gamma$. At large
$\Gamma$, the droplets are nearly rigid and the clogging probability is large
since the droplets can not squeeze through the gaps between obstacles. With
decreasing $\Gamma$, the clogging probability decreases as the droplets become
more deformable. However, in the small-$\Gamma$ limit the clogging probability
increases, since the droplets are extremely deformable and cause clogs as they
wrap around the obstacles. The results from these studies are important for
developing a predictive understanding of capillary droplet flows through
complex and confined geometries.

††preprint: APS/123-QED

##  I Introduction

Capillary droplets, which consist of immiscible liquid droplets dispersed in
another continuous phase, are ubiquitous in food science [1, 2], the
pharmaceutical industry [3, 4], petroleum products [5, 6], and pollution
remediation [7]. Being able to control how capillary droplets flow and
interact with their surroundings is crucial for enhancing food emulsion
stability, optimizing drug production and delivery, and improving oil recovery
rates, while minimizing environmental risks associated with oil exploration.
For example, the flow rate of capillary droplets can be tuned by varying the
channel width and by designing obstacle arrays in microfluidic chambers [8, 9,
10, 11, 12, 13], cell sorting instruments [14, 15, 16], and cell encapsulation
devices [17, 18, 19]. However, we do not yet have a fundamental understanding
of the dynamics of capillary droplets flowing through complex, confined
geometries as a function of the droplet properties.

Flows of hard, frictional grains through narrow constrictions, such as
hoppers, silos, and obstacle arrays, display complex spatiotemporal dynamics
[20, 21, 22]. Over $50$ years ago, Beverloo and co-workers [23] carried out
experimental studies of hopper flows for a wide range of dry granular
materials and proposed an empirical form for the flow rate $Q$ (in the
continuous flow rate regime) versus the orifice width $w$:
$Q\sim(w/\sigma_{\rm avg}-k)^{\beta}$, where $\sigma_{\rm avg}$ is the mean
grain size, the flow arrests for $w/\sigma_{\rm avg}\lesssim k$, $k\gtrsim
1.5$ for rigid grains, and $\beta$ is the power-law scaling exponent that can
be tuned by varying the ratio of the dissipation from particle-particle
contacts and from the background fluid [24]. In addition, researchers have
found that placing an obstacle $\sim 3\sigma_{\rm avg}$ above the orifice can
reduce clogging in silo flows of hard, frictional grains [25, 26, 27, 28, 29,
30, 31]. Experimental and computational studies have also shown that the
permeability of an obstacle array to granular flow can be tuned by varying the
density, shape, and size of the obstacles [32, 33, 34, 35].

In contrast, for deformable particles, it is possible to find $Q>0$ with no
clogging even for $w/\sigma_{\rm avg}<1$. Unlike hard grains, deformable
particles can change their shapes to squeeze through narrow constrictions that
are much smaller than their undeformed diameters. However, there have been few
studies of flows of deformable particles through narrow constrictions. Many
experimental studies have instead focused on hydrogel particle flows in
hoppers and silos with $w/\sigma_{\rm avg}\sim 2$-$3$, which do not give rise
to large particle deformation. For example, recent studies have found that
hydrogel particles with smaller elastic moduli are less likely to clog [36,
37, 38] and possess flow rates that are much faster than those for hard
particles, such as glass beads [39]. These prior results for hydrogel
particles do not address how particle flow rates are affected by large
particle deformation. Thus, an important open question is: does significant
particle shape change (i.e. by more than $50\%$) still enhance flow through
narrow constrictions and other confined geometries?

The interplay between particle deformability and the geometry of obstacle
arrays can give rise to a new physical mechanism, particle wrapping, which can
strongly influence the dynamics of capillary droplets as they flow through
obstacle arrays. When colliding with an obstacle, the highly deformable
droplet can wrap around the obstacle with both ends of the droplet elongating
in the flow direction until the energy penalty from the droplet surface
tension matches the work done by the fluid driving force. At this point, both
ends of the droplet become stationary, and then one end of the droplet moves
in the opposite direction of the driving force, causing the particle to
unwrap. Particle wrapping can even lead to single-particle clogging in
obstacle arrays. In this case, the droplet is in force balance after wrapping
around an obstacle, where the fluid forces on the droplet are supported by
contact forces on the droplet from neighboring obstacles. Thus, particle
wrapping can lead to much longer transit times through obstacle arrays
compared to those for hard particles with the same undeformed diameter.
Despite its importance, the particle wrapping mechanism has not been studied
in detail in droplet flows through obstacle arrays.

A key factor that contributes to the knowledge gap concerning flows of highly
deformable particles in confined geometries is the lack of an accurate and
efficient computational model for these flows. At submicron length scales,
e.g. within microfluidic devices, inertial forces are typically much smaller
than capillary forces generated by particle surface tension, which strongly
influences the behavior of multi-phase fluid flows. Computational models for
multi-phase fluid flows include grid-based Eulerian techniques, such as the
volume-of-fluid [40] and Lattice Boltzmann methods [41, 42], and particle-
based Lagrangian techniques, such as the smoothed particle hydrodynamics
method [43, 44, 45]. Grid-based methods determine the amount of fluid and gas
in each cell and calculate the surface energy by reconstructing the fluid-gas
interface. However, this method does not strictly conserve mass, and the
surface energy can depend on the interface reconstruction algorithm [40]. For
particle-based methods, the multi-phase system is discretized into “particles”
with each particle assigned as gas or liquid with a specified mass. Therefore,
interface tracking is not needed, and the multi-phase fluid behavior is
governed by particle-particle interactions. However, it is challenging to
determine a priori the particle-particle interactions that yield the correct
contact angles at the multi-phase interfaces [46, 43]. Particle-based methods
are also computationally demanding since they simulate the droplet interiors
in addition to the interfaces. As a result, most particle-based simulations of
multi-phase fluids are limited to short time scales and small systems.

In this article, we describe combined experimental and computational studies
that enable the development of new efficient computational methodologies and
provide new insights for flows of deformable particles through narrow
constrictions and obstacle arrays. In the experimental studies, we measure the
speed and shape of capillary droplets as they flow through narrow openings as
a function of the orifice width and gravitational driving force. For the
computational studies, we develop a new deformable particle (DP) model[47, 48,
49] with surface tension (in two dimensions, 2D) that can be implemented to
efficiently model thousands of capillary droplets flowing through complex,
confined geometries. Varying surface tension in the experiments by more than a
factor of $2$ is challenging, and thus we carry out simulations of extremely
deformable particles with surface tensions that vary by more than a factor of
$10^{3}$, which accentuates the particle wrapping mechanism. We first validate
the DP simulations by comparing the droplet shape dynamics in experiments and
simulations for a single droplet flowing through narrow constrictions. We then
investigate the important physical parameters that control the flow rate and
clogging propensity of single deformable particles moving through obstacle
arrays using the experimentally validated DP simulations.

We find several important results. First, when a capillary droplet is pushed
by gravity through a narrow constriction, its speed initially decreases as it
enters the orifice and deforms. After it reaches a minimum speed, the droplet
begins to speed up as it exits the orifice. However, as the speed increases,
it can overshoot the droplet’s terminal speed. We find that the amplitude of
the overshoot decreases with increasing gravitational driving force (relative
to the capillary forces). Second, we show that highly deformable droplets flow
more slowly in obstacle arrays compared to droplets with larger surface
tension, due to the wrapping mechanism. Thus, for capillary droplets flowing
through obstacle arrays, the clogging probability, as characterized by the
average distance $\lambda$ traveled by a droplet before it clogs, is non-
monotonic as a function of surface tension. For small values of surface
tension, $\lambda$ first increases with increasing surface tension, reaches a
peak that depends on the obstacle spacing, and then decreases with continued
increases in surface tension. In contrast, models of soft particles that do
not model explicit changes in particle shape cannot capture particle
wrapping[50, 51], and thus cannot accurately model flows of capillary droplets
moving through obstacle arrays in the low surface tension limit.

The remainder of the article is organized as follows. In Section II, we
describe the computational methods including descriptions of the soft (semi-
rigid) particle model[50, 51] and deformable particle model[47, 48, 49] for
droplets and the equations of motion and boundary conditions for discrete
element method (DEM) simulations of droplet flows through narrow channels and
obstacle arrays. In Section III, we introduce the experimental setup and
method to generate capillary droplets and flows. In Section IV, we describe
the main results from the combined experimental and simulation studies. We
first characterize the droplet speed and shape dynamics as a function of the
distance $h$ from the orifice for both oil-in-water experiments and DPM
simulations of droplet flows through narrow channels. We calibrate the results
for the droplet speed and shape from the DP simulations to those from the
experiments. We then quantify the nonmonotonic behavior in the droplet speed
profile and the clogging probability of single droplets flowing through
obstacle arrays as a function of the droplet surface tension. In Section V, we
provide conclusions and propose future research directions, including studies
of multi-droplet flows through obstacle arrays with droplet breakup and
coalescence that can predict the steady-state droplet size distribution. We
also include five Appendices to supplement the discussion in the main text. In
Appendix A, we present experimental data for the near-wall droplet speed
profile. In Appendix B, we investigate how the results of the DP simulations
depend on the number of vertices $N_{v}$ in each deformable particle and the
perimeter elasticity $K_{l}$. In Appendix C, we provide additional details
concerning the definition of a clogging event and the clogging probability for
droplets flowing through obstacle arrays. In Appendix D, we identify the
surface tension regime for which the simulations can achieve continuous flows
in obstacle arrays. In Appendix E, we quantify the dependence of the steady-
state speed of the droplets on the near-wall drag coefficient.

##  II Simulation Methods

In this section, we describe the methods for simulating 2D gravity-driven
flows of a single droplet through narrow channels and obstacle arrays. We
first illustrate the geometry of the narrow channels and obstacle arrays. We
then introduce two models for describing the droplet shape and interactions
with the walls: 1) the soft particle (SP) model, which tracks the centers of
mass of the droplets and enforces excluded volume interactions using pairwise
purely repulsive forces that increase with the overlap between a droplet and
the walls and 2) the deformable particle (DP) model that treats each droplet
as a polygon with $N_{v}$ vertices, employs a shape-energy function to
penalize changes in droplet area and changes in the separations between
vertices, and enforces excluded volume interactions using pairwise purely
repulsive forces between the droplet vertices and walls. For both models, we
define the forces on the droplets that arise from the shape-energy function,
droplet-wall interactions, and dissipative forces from the background fluid,
and provide the equations of motion for the droplet trajectories. Finally, we
describe the method for initializing the droplet positions and velocities in
narrow channels and obstacle arrays.

###  II.1 Geometry for flows through narrow channels and obstacle arrays

The narrow channels in the simulations are constructed using two infinitely
long wedges, one on the top and one on the bottom that form an orifice between
them with width $w<\sigma$, where $\sigma$ is the diameter of the undeformed
droplet. (See Fig.1 (a).) The obstacle arrays are located within $L\times L$
simulation boxes with periodic boundary conditions and $L=500\sigma$. The
circular obstacles with diameter $\sigma_{\rm ob}$ are placed inside the
simulation box using a Poisson sampling procedure with separations $w>w_{\rm
ob}$. (See Fig.1 (b).)

###  II.2 Deformable particle model

When droplets squeeze through narrow constrictions, they can experience
dramatic shape changes. To efficiently describe large changes in shape that
occur during gravity-driven droplet flows, we employ the recently developed
deformable particle (DP) model [47, 48, 49]. In 2D, the droplets are modeled
as deformable polygons with $N_{v}$ vertices. The shape-energy function
$U_{s}$ for each deformable particle includes three terms:

$U_{s}=\frac{k_{a}}{2}(a-a_{0})^{2}+\frac{k_{l}N_{v}}{2}\sum_{i=1}^{N_{v}}(l_{i}-l_{0})^{2}+U_{\gamma}.$
(1)

The first term imposes a harmonic energy penalty for changes in the droplet
area $a$ from the preferred value $a_{0}$ and $k_{a}$ controls the
fluctuations in drloplet area. The second term imposes a harmonic energy
penalty for deviations in the separations $l_{i}$ between adjacent vertices
$i$ and $i+1$ from the equilibrium length $l_{0}$ and $k_{l}$ controls
fluctuations in the separations between adjacent vertices. This term ensures
that vertices are distributed evenly on the droplet surface, preventing them
from aggregating in the direction of gravity. The factor of $N_{v}$ in the
numerator of the second term of Eq. 1 makes $U_{s}$ independent of $N_{v}$. We
characterize the shape of the droplet using the dimensionless shape parameter:
${\cal A}=(\sum_{i=1}^{N_{v}}l_{i})^{2}/(4\pi a)$. The third term gives the
energy arising from line tension:

$U_{\gamma}=\gamma_{\rm so}l_{\rm so}+\gamma_{\rm sw}(L-l_{\rm
so})+\gamma_{\rm ow}l_{\rm ow},$ (2)

where $\gamma_{\rm so}$, $\gamma_{\rm sw}$, and $\gamma_{\rm ow}$ are the
solid-oil, solid-water, and oil-water line tension coefficients. (See Fig. 1
(c) for a schematic of the solid wall-oil droplet-water interface in 2D.)
$l_{\rm so}$ and $l_{\rm ow}$ are the contact lengths at the oil-solid and
oil-water interfaces, respectively. $L-l_{\rm sw}$ is the length of the solid
wall that is not in contact with the oil droplet. In the oil-in-water
experiments described here, thin layers of water exist between the oil
droplets and chamber walls, thus effectively setting $\gamma\equiv\gamma_{\rm
so}=\gamma_{\rm ow}=\gamma_{\rm sw}$.

Figure 1: Schematic of (a) a droplet (shaded red) with undeformed diameter
$\sigma$ flowing through a narrow orifice (shaded black) with opening width
$w=0.4\sigma$ and (b) a droplet (shaded red) flowing through an obstacle array
(shaded blue) with random obstacle positions, obstacle sizes $\sigma_{\rm
ob}=0.3\sigma$, and minimum separation $w_{\rm ob}=\sigma$. The small open
circles represent the droplet’s trajectory. The rightward pointing arrow
indicates the direction of gravity. (c) Illustration of the DP model in 2D
with solid-oil $\gamma_{\rm so}$, oil-water $\gamma_{\rm ow}$, and solid-water
$\gamma_{\rm sw}$ line tensions. In the oil-in-water experiments, a thin layer
of water exists between the oil droplets and chamber walls, and thus
$\gamma_{\rm so}=\gamma_{\rm ow}=\gamma_{\rm sw}$. $l_{\rm so}$ is the contact
length between the solid wall and oil droplet, $L$ is the length of the solid
wall, $l_{i}$ is the separation between adjacent droplet vertices, and $a$ is
the droplet area. Figure 2: Series of images of a single droplet flowing
through a narrow orifice with width $w=0.4\sigma$ from (a) oil-in-water
experiments with undeformed droplet diameter $\sigma\approx 3.5~{}\rm mm$ and
tilt angle $\theta\approx 28^{\circ}$ and (b) a DP simulation in 2D at
dimensionless line tension $\Gamma=\Gamma^{*}$ and near-wall drag coefficient
$b_{0}=b_{0}^{*}$. The scale bar indicates $1$ mm. Below panel (a), we provide
the distances of the droplet center of mass to the orifice $h$ at which the
images are captured. We find that $\Gamma^{*}\approx 0.16$ and
$b_{0}^{*}/b_{\infty}\approx 0.064$ minimize the deviation in the droplet’s
center of mass speed $\Delta_{v}$ between the DP simulations and experiments.
These best-fit values give $\Delta_{v}\approx 0.09$ and
$\Delta_{\cal{A}}\approx 0.01$. In panel (a), we overlay the shape of the
droplet from the DP simulations at $h=0$ (red dashed line) onto the
corresponding droplet image for the experiments (black solid line). Droplet
(c) shape parameter $\mathcal{A}$ and (d) center of mass speed in the driving
direction $v_{g}$ plotted as a function of $h/\sigma$ for both experiments
(dashed lines) and DP simulations (solid lines). We estimate the dimensionless
surface tension in the experiments to be $\Gamma_{\rm exp}\approx 0.57$.

We assume that there are purely repulsive forces between the drloplets and
hopper walls or obstacles, which are modeled using the following pairwise
interaction potential:

$U_{w}=\sum_{i=1}^{N_{v}}\frac{\epsilon_{w}}{2}\left(1-\frac{2d_{i}}{\delta}\right)^{2}\Theta\left(1-\frac{2d_{i}}{\delta}\right),$
(3)

where $\epsilon_{w}$ sets the strength of the repulsive interactions, $d_{i}$
is the minimum distance between the center of vertex $i$ of the droplet and
the wall or obstacle surface, and $\delta$ is the vertex diameter. The
Heaviside step function $\Theta(\cdot)$ ensures that the pair forces are non-
zero only between vertex $i$ overlaps a chamber wall or obstacle.

The total potential energy $U$ for the droplet is the sum of the shape-energy
function $U_{s}$, gravitational potential energy $U_{g}$, and particle-wall
interactions $U_{w}$,

$U=U_{s}+U_{g}+U_{w},$ (4)

where $U_{g}=-Mgh$, $h$ is the height of the center of mass of the droplet,
$g$ is the gravitational acceleration, $M=\rho a_{0}$ is the mass of the
droplet with areal mass density $\rho$ and area $a_{0}=\pi\sigma^{2}/4$.

We also consider dissipative forces from the background fluid acting on each
droplet vertex $i$:

$\vec{F}_{i}^{\zeta}=-\frac{\zeta_{i}}{N_{v}}{\vec{\nu}}_{i},$ (5)

where ${\vec{\nu}}_{i}$ is the velocity of vertex $i$ and $\zeta_{i}$ is the
viscous drag coefficient:

$\zeta_{i}=b_{\infty}\left(1+\frac{b_{0}}{b_{\infty}}\left(\frac{\sigma}{d_{i}}\right)^{\alpha}\right).$
(6)

In Eq. 6, $b_{\infty}$ is the drag coefficient when the droplet is far from
the walls, $d_{i}$ is the minimum distance between the center of vertex $i$
and the chamber wall (or obstacle), and $b_{0}$ is the near-wall drag
coefficient. The exponent $\alpha=4/3$ is calibrated to the experiments on
single droplet flows in narrow channels described in Appendix A.

Thus, for the DP model, the equations of motion for the position
${\vec{r}}_{i}$ of droplet vertex $i$ is

$M_{i}\frac{\partial^{2}\vec{r}_{i}}{\partial
t^{2}}=-\vec{\nabla}_{r_{i}}U+\vec{F}_{i}^{\zeta},$ (7)

where $M_{i}=M/N_{v}$ is the mass of vertex $i$.

From Eqs. 1-3, we introduce the dimensionless line tension
$\Gamma=\gamma/(g\rho\sigma^{2})$ and three additional dimensionless energy
scales for the DP model in a gravitational field in 2D:
$K_{a}=k_{a}\sigma/(g\rho)$, $K_{l}=k_{l}/(g\rho\sigma)$, and
$E_{w}=\epsilon_{w}/(g\rho\sigma^{3})$. In these studies, we vary $\Gamma$ in
the range $10^{-2}<\Gamma<10$. We set $K_{a}>10^{4}$ so that the fluctuations
in the droplet areas are small, $K_{l}/\Gamma<10^{-3}$ so that the line
tension energy dominates the perimeter elasticity (see Appendix B), and
$E_{w}=10^{4}$ to minimize the vertex-wall (and vertex-obstacle) overlaps. The
time $t$ is measured in units of $t_{0}=\sqrt{\sigma/g}$ and the equations of
motion are integrated using a modified velocity-Verlet algorithm with time
step $\Delta t=10^{-3}t_{0}$.

###  II.3 Soft particle model

We will also describe gravity-driven single-droplet flows through obstacle
arrays using the soft particle (SP) model, and compare the results from the SP
model to those of the DP model to better understand the influence of droplet
shape change on the flows. In general, there are three contributions to the
total potential energy: 1) the shape-energy function $U_{s}$, 2) the
gravitational potential energy $U_{g}$, and 3) the droplet-obstacle
interaction energy $U_{w}$. For the SP model, $U_{s}=0$ and excluded volume
interactions between the droplet and obstacles are generated by allowing
overlaps between them[50, 51, 36, 37]. The pairwise droplet-obstacle
interaction energy is

$U_{w}=\frac{\epsilon_{\rm
w}}{2}\left(1-\frac{2d}{\sigma}\right)^{2}\Theta\left(1-\frac{2d}{\sigma}\right),$
(8)

where $\sigma$ is the SP droplet diameter, $d=r-\sigma_{\rm ob}/2$ is the
distance between the center of the droplet and the obstacle surface,
$\sigma_{\rm ob}$ is the obstacle diameter, $r$ is the separation between the
center of the droplet and the center of the closest obstacle, and
$\epsilon_{\rm w}$ is the characteristic energy of the repulsive interactions.
The Heaviside step function $\Theta(\cdot)$ ensures that the pair forces are
non-zero only when the droplet overlaps with an obstacle. Thus, the total
potential energy of an SP droplet in an obstacle array is

$U=U_{g}+U_{w}.$ (9)

For the SP model, we consider following form for the viscous drag forces on
droplets moving in a background viscous fluid:

$\vec{F}_{\zeta}=-\zeta_{\rm sp}\vec{v},$ (10)

where $\zeta_{\rm sp}$ is the drag coefficient and $\vec{v}$ is the velocity
of the SP droplet. We note that, unlike Eq. 6 for the DP model, we set
$\zeta_{\rm sp}=b_{\infty}$. In contrast to the DP model, the SP model
considers droplet deformation by allowing overlaps between the droplets and
obstacles, which makes it difficult to describe distance-dependent viscous
drag. For the SP model, the equations of motion for the droplet position
${\vec{r}}_{i}$ are

$M\frac{\partial^{2}\vec{r}}{\partial
t^{2}}=-\vec{\nabla}_{r}U+\vec{F}_{\zeta}.$ (11)

We integrate Eq. 11 using a modified velocity-Verlet integration scheme with
time step $\Delta t=10^{-3}t_{0}$. For the SP model, an important
dimensionless energy scale is the ratio of the excluded volume interactions to
the gravitational potential energy, $E_{\rm w}=\epsilon_{\rm
w}/(g\rho\sigma^{3})$. Here, we vary $E_{\rm w}$ over the range
$10^{-1}<E_{\rm w}<10^{2}$.

Figure 3: Droplet (a) shape parameter $\mathcal{A}$ and (b) center of mass
speed $v_{g}$ plotted as a function of distance from the orifice $h/\sigma$
for both experiments (dashed lines) and DP simulations (solid lines) at
$w=0.9\sigma$. For (a) and (b), the dimensionless surface tension $\Gamma_{\rm
exp}\approx 5.80$, tilt angle $\theta\approx 8.5^{\circ}$, and undeformed
droplet diameter $\sigma\approx 1.98$ mm. After minimizing $\Delta_{v}$
between the experiments and DP simulations, we obtain $\Gamma^{*}\approx 1.7$,
$b_{0}^{*}/b_{\infty}\approx 0.21$, $\Delta_{v}\approx 0.08$, and
$\Delta_{\cal{A}}\approx 0.004$. The inset to (a) shows the droplet shape from
the DP simulations at maximum deformation with $\mathcal{A}_{\rm max}\approx
1.014$. (c) $\mathcal{A}$ and (d) $v_{g}$ plotted as a function of $h/\sigma$
for both experiments and DP simulations at $w=0.7\sigma$. For (c) and (d),
$\Gamma_{\rm exp}\approx 3.47$, $\theta\approx 5.0^{\circ}$, and
$\sigma\approx 3.30$ mm. After minimizing $\Delta_{v}$, we obtain
$\Gamma^{*}\approx 1.08$, $b_{0}^{*}/b_{\infty}\approx 0.07$,
$\Delta_{v}\approx 0.06$, and $\Delta_{\cal{A}}\approx 0.005$. The inset to
(c) shows the droplet shape from the DP simulations at maximum deformation
with $\mathcal{A}_{\rm max}\approx 1.093$.

###  II.4 Simulation initialization

For the DP model, we initialize the droplets in narrow channels and obstacle
arrays in the shape of regular polygons with edge lengths equal to their
equilibrium values $l_{0}=\sqrt{4a_{0}N_{v}\tan(\pi/N_{v})}/N_{v}$. For
droplet flows in narrow channels, we initially place the droplet far from the
orifice, typically at $h=-10\sigma$, so that the droplet reaches terminal
speed before touching the channel walls. For droplets in obstacle arrays, we
initially place the droplet in the center of the simulation box. However, this
placement typically causes overlaps between the droplet and nearby obstacles.
Therefore, after the initial placement of the droplet, we carry out energy
minimization (in the absence of gravity) to eliminate overlaps between the
droplet and obstacles. We then integrate the equations of motion with gravity
and the droplet begins to flow through the obstacle array until it reaches
steady state or clogs. (See Sec. IV.2.)

##  III Experimental Methods

Below, we will compare the SP and DP simulation results for gravity-driven
droplet flows in narrow channels and obstacle arrays to experimental studies
of quasi-2D gravity-driven flows of octane oil droplets in water. In this
section, we describe the details of our experimental studies. We construct the
sample chambers using laser-cut plastic film with thickness $400$ $\mu$m,
sandwiched between two glass microscope slides. The film is attached to the
slides using Norland Optical Adhesive [24]. The diameters of the droplets are
tuned so that they are larger than the chamber thickness to ensure that each
droplet is quasi-two-dimensional when confined by the microscope slides.
Obstacles in the experimental studies are produced by placing small drops of
UV adhesive between the slides, which then cure into solid cylinders spanning
the sample chamber thickness.

Oil-in-water droplets are prepared via direct injection of the oil into the
chamber. We use octane with density $\rho_{o}=0.703$ g/ml suspended in water
with density $\rho_{w}=0.997$ g/ml. Despite the absence of other oil droplets,
we find that it is still necessary to use a 0.5% Tween 20 nonionic detergent
solution to prevent sticking of the oil droplets to the microscope slide
surfaces.

We image the droplets using a $1.6\times$ lens (0.05 Numerical Aperture)
attached to a Leica DM4500 B inverted microscope. Movies are acquired by
attaching a $0.35\times$ C-mount to a ThorLabs DCC1545M camera, and recording
at a frame rate commensurate with the droplet speed. The microscope is tilted
slightly from the horizontal to cause gravitationally driven motion of the
droplet. In particular, the tilt angle of the microscope is changed by placing
a spacer underneath the front edge of the microscope, and the angle is
measured relative to the ground via a Wixley digital angle gauge.

The droplet motion is affected by viscous drag from the water phase (viscosity
$\eta_{w}\approx 1$ mPa$\cdot$s) and viscous drag between the oil and glass
plates ($\eta_{o}\approx 100$ mPa$\cdot$s). The droplets try to remain round
due to surface tension ($\gamma_{\rm ow}\approx 10$ mJ/m2, as measured by the
“Dropometer” from Droplet Labs). However, under gravity, the droplets deform
when pressed against other droplets or the sample chamber walls.

We can define a quasi-2D dimensionless surface tension for the experimental
studies:

$\Gamma_{\rm exp}=\frac{\gamma_{\rm exp}}{g_{\rm exp}\rho_{o}\sigma^{2}}.$
(12)

Here, $\gamma_{\rm exp}$ is the oil-water interfacial surface tension, $g_{\rm
exp}=\sin\theta g(\rho_{w}-\rho_{o})/\rho_{o}$ is the acceleration imposed by
oil-in-water buoyancy, $g\sim 9.8$ m/s2 is the gravitational constant,
$\sigma$ is the average undeformed diameter of the droplets, and $\theta$ is
the microscope tilt angle.

##  IV Results

In this section, we describe the results from the numerical simulations of
gravitationally-driven droplet flows through narrow channels and obstacle
arrays. We first present a method to calibrate the DP simulations to
experimental results for droplet flows through narrow channels. We then show
results for the droplet speed in the direction of the driving $\nu_{g}$ as a
function of distance from the orifice. We find that the droplet speed profile
is nonmonotonic and can even possess an overshoot, such that the droplet speed
exceeds the terminal speed $\nu_{t}$ as it exits the orifice. Using the DP and
SP models, we also investigate gravity-driven droplet flows through obstacle
arrays. For the DP simulations, we find non-monotonic behavior for the
clogging probability as a function of the dimensionless line tension $\Gamma$.
We show that droplets can clog at small $\Gamma$ when they are highly
deformable and can wrap around the obstacles. We emphasize that the SP model
does not include the wrapping mechanism since wrapping requires significant
changes in the droplet shape. As shown in many previous studies of flows of
nearly hard particles in confined geometries, the droplets can also clog when
they are nearly rigid at large $\Gamma$. At intermediate $\Gamma$, we find a
regime of continuous droplet flows through the obstacle arrays. We observe
qualitatively similar results over a wide range of obstacle densities.

###  IV.1 Calibration of DP model using experiments of a single droplet
flowing through narrow channels

We first carried out both experiments and DP simulations of a single droplet
flowing through narrow channels under the influence of gravity. To enable a
comparison between the experimental and simulation results, we define the
root-mean-square deviations in the droplet shape parameter $\Delta_{\cal{A}}$
and speed $\Delta_{\nu}$:

$\Delta_{\cal{A}}=\sqrt{\sum_{i=1}^{N_{f}}\frac{({\cal{A}}_{\rm
exp,i}-{\cal{A}}_{\rm sim,i})^{2}}{N_{f}}}$ (13)

and

$\Delta_{\nu}=\sqrt{\sum_{i=1}^{N_{f}}\frac{((\nu_{\rm exp,i}-\nu_{\rm
sim,i})/\nu_{\rm sim,i})^{2}}{N_{f}}},$ (14)

where ${\cal{A}}_{\rm exp,i}$ is the droplet shape parameter in the $i$th
frame in the experimental video, and ${\cal{A}}_{\rm sim,i}$ is the droplet
shape parameter in the DP simulations when the droplet center of mass is at
the same location as it is in the experiments at frame $i$. Similarly,
$\nu_{\rm exp,i}$ and $\nu_{\rm sim,i}$ are the droplet speed in the direction
of gravity at the $i$th frame in the experiments and simulations. The
summations in Eqs. 13 and 14 are over all $N_{f}$ frames in the experimental
videos. Despite using a small time step in the DP simulations, we might not
find a frame in the simulations where the droplet has the same position as
that in the experiments. In this case, we use linear interpolation between
frames to generate ${\cal{A}}_{\rm sim,i}$ and $\nu_{\rm sim,i}$.

We seek to minimize $\Delta_{\cal{A}}$ and $\Delta_{\nu}$ by tuning the
dimensionless line tension $\Gamma$ and near-wall drag coefficient
$b_{0}/b_{\infty}$ in the DP simulations. In practice, $\Delta_{\cal{A}}$ is
more than an order of magnitude smaller than $\Delta_{\nu}$, and thus we only
minimize $\Delta_{v}$ over $\Gamma$ and $b_{0}/b_{\infty}$. $b_{\infty}$ is
determined by matching the droplet’s terminal speed far from wall in the DP
simulations and experiments. In particular, we equate the dimensionless
terminal speed,

$V_{t}\equiv\frac{v_{t}}{\sqrt{\sigma g_{\rm
exp}}}=\frac{mg}{b_{\infty}\sqrt{\sigma g_{\rm exp}}},$ (15)

in the experiments and DP simulations.

In Fig. 2 (a) and (b), we show the droplet shapes as a function of distance
$h$ from the orifice (with width $w=0.4\sigma$) from trajectories of gravity-
driven droplet flows through narrow channels in experiments and DP simulations
using the values $\Gamma^{*}=0.16\pm 0.01$ and $b_{0}^{*}/b_{\infty}=0.064\pm
0.003$ that minimize $\Delta_{\nu}$. The error bars for $\Gamma^{*}$ and
$b_{0}^{*}$ are the standard deviations in $\Gamma^{*}$ and $b_{0}^{*}$ from
fitting the the DP simulations to eight independent experimental trials. For
the experiments in Fig. 2 (a), the terminal velocity $v_{t}\approx 7$ mm/s far
from the wall, undeformed droplet diameter $\sigma\approx 3.5~{}\rm mm$, and
tilt angle $\theta\approx 28^{\circ}$. Thus, for the data in Fig. 2 (a), we
estimate the dimensionless surface tension for the droplets to be $\Gamma_{\rm
exp}\approx 0.57$ and the dimensionless terminal velocity to be $V_{t}\approx
0.08$. In Fig. 2 (c) and (d), we plot $\mathcal{A}$ and $v_{g}$ as a function
of $h/\sigma$ for both the calibrated DP simulations and experiments for the
data presented in Fig. 2 (a) and (b). We find that the droplet shapes and
speed profiles for the DP simulations and experiments are similar over the
full range of $h$, with $\Delta_{\cal{A}}\approx 0.01$ and $\Delta_{v}\approx
0.09$. Note that despite the simplicity of the 2D DP simulations, $\Gamma^{*}$
is comparable to $\Gamma_{\rm exp}$ (with $\Gamma^{*}/\Gamma_{\rm exp}\approx
0.3$), emphasizing the predictive power of 2D DP simulations.

Figure 4: Droplet center of mass speed in the direction of the gravitational
driving (normalized by the terminal speed) $v_{g}/v_{t}$ plotted as a function
of $h/\sigma$ for droplet flows in narrow channels with orifice width
$w=0.7\sigma$ in the experiments (solid lines) for several values of the tilt
angle, $2.7^{\circ}$ (black), $4.5^{\circ}$ (blue), and $6.0^{\circ}$ (red),
and calibrated DP simulations (dashed lines) with $\Gamma^{**}=1.85$ (black),
$1.11$ (blue), and $0.83$ (red). The downward arrow indicates increasing tilt
angle $\theta$.

In Fig. 3, we compare the results for the calibrated DP simulations and
experiments of gravity-driven droplet flows in narrow channels for two
additional orifice widths and tilt angles. In Fig. 3 (a) and (b), we show
${\cal A}$ and $v_{g}$ versus $h/\sigma$ for $w/\sigma=0.90$, $\theta\approx
8.5^{\circ}$, $\sigma\approx 1.98$ mm, $v_{t}\approx 0.97$ mm/s, $V_{t}\approx
0.028$, and $\Gamma_{\rm exp}\approx 5.80$. In Fig. 3 (c) and (d), we show
${\cal A}$ and $v_{g}$ versus $h/\sigma$ for $w/\sigma=0.70$, $\theta\approx
5.0^{\circ}$, $\sigma\approx 3.30$ mm, $v_{t}\approx 1.1$ mm/s, $V_{t}\approx
0.032$, and $\Gamma_{\rm exp}\approx 3.47$. We can obtain the dimensionless
line tension $\Gamma^{**}$ that minimizes $\Delta_{\nu}$ for the experimental
data in Fig. 3 (with tilt angle $\theta^{**}$ and undeformed droplet diameter
$\sigma^{**}$) using the values of $\Gamma^{*}$, $\theta^{*}$, and
$\sigma^{*}$ from the calibrated experiments and DP simulations in Fig. 2
using
$\Gamma^{**}/\Gamma^{*}=(\sigma^{*})^{2}\sin\theta^{*}/((\sigma^{**})^{2}\sin\theta^{**})$.
We find $\Gamma^{**}=1.7\pm 0.1$ for $\theta^{**}=8.5^{\circ}$ and
$\sigma^{**}\approx 1.98$ mm (with $\Gamma_{\rm exp}\approx 5.8$) and
$\Gamma^{**}=1.08\pm 0.09$ for $\theta^{**}=5.0^{\circ}$ and
$\sigma^{**}=3.30$ mm (with $\Gamma_{\rm exp}\approx 3.47$). For all tilt
angles and droplet diameters studied, we find that the ratio of the
dimensionless line tension that minimizes $\Delta_{v}$ to the dimensionless
surface tension in experiments satisfies $\Gamma^{**}/\Gamma_{\rm exp}\approx
0.3$, which likely stems from the fact that we compare values of the
dimensionless line tension in 2D to values of the dimensionless surface
tension in 3D. We find that the near-wall drag coefficients that minimize
$\Delta_{v}$ do not depend strongly on droplet diameter and tilt angle,
$b_{0}^{*}/b_{\infty}\approx 0.06$-$0.2$.

A key feature of the droplet speed profile $v_{g}$ is that it can possess
nonmonotonic behavior and $v_{g}$ can even exceed the terminal speed $v_{t}$
as the droplet exits the orifice. For example, in Fig. 3 (d), $\nu_{g}$
overshoots $v_{t}$ near $h/\sigma\approx 0.3$. Similar phenomena have been
observed in pressure-driven droplet flows in obstacle arrays, where the
pressure drop across the orifice changes from positive to negative, which
indicates that the droplet speed exceeds the background fluid velocity[42]. To
further investigate the effects of the gravitational driving on the
nonmonotonic speed profile, we show $\nu_{g}$ versus $h/\sigma$ for fixed
$w/\sigma=0.7$ for several values of the tilt angle: $\theta=2.7^{\circ}$,
$4.5^{\circ}$, and $6^{\circ}$ in Fig. 4. The nonmonotonic droplet speed
profile can be understood by analyzing energy transfer while the droplet moves
through the narrow channel. The droplet first impacts the chamber wall at
$h/\sigma<0$ and deforms while squeezing through narrow orifice. During this
process, the kinetic and gravitational energy of the droplet transfers into
surface energy of the deformed droplet, which corresponds to the first speed
minimum near $h/\sigma\approx-0.25$. The large surface energy is transferred
into kinetic energy as the droplet continues to move through the orifice,
which can give rise to speeds that exceed the terminal velocity near
$h/\sigma\approx 0.3$. However, the excess kinetic energy is quickly
dissipated by viscous drag, and $\nu_{g}$ decreases again, reaching a second
local minimum near $h/\sigma\approx 0.5$. Finally, the speed increases again,
approaching $\nu_{t}$ as the droplet moves far from the orifice since the drag
coefficient $\zeta$ decreases toward $b_{\infty}$ for $h/\sigma\gg 1$. (See
Eq. 6.)

As shown in Fig. 4, the amplitude of the overshoot in the droplet speed
increases with decreasing $g$ and $\theta$. The excess shape energy $\Delta
U_{s}$ gained during droplet deformation is largely determined by the width of
the orifice, $w/\sigma$. In contrast, the droplet kinetic energy at the
terminal speed $E_{t}$ decreases with increasing $g$. Thus, for $h>0$, the
excess surface energy $\Delta U_{s}$ transferred into kinetic energy
$E_{k}\sim\Delta U_{S}$, generates enhanced droplet speed overshoots
$\nu_{g}\sim\sqrt{E_{k}/E_{t}}$ as $g$ or $\theta$ decreases.

Figure 5: Oil-in-water droplets flowing through an obstacle array can wrap
around obstacles in (a) experiments with $\Gamma_{\rm exp}\approx 0.88$, tilt
angle $\theta\approx 90^{\circ}$, undeformed droplet diameter $\sigma\approx
2$ mm, obstacle diameter $\sigma_{\rm ob}/\sigma\approx 0.3$-$0.4$, and
minimum obstacle separation $w_{\rm ob}/\sigma\approx 0.2$ and (b) calibrated
DP simulations with $\Gamma\approx 0.3$ and a similar obstacle array to that
used in (a). In experiments, the driving force is directed from left to right
and in the DP simulations gravity is directed from left to right. The droplets
can also squeeze through narrow constrictions formed by adjacent obstacles in
a similar set of (c) experiments and (d) DP simulations as those in (a) and
(b).

###  IV.2 Clogging of single droplets in obstacle arrays

Droplet shape dynamics plays an important role in determining the motion and
clogging of droplets flowing through obstacle arrays. For example, in Fig. 5
(a) and (c), we show that droplets can wrap around obstacles and squeeze
through narrow constrictions formed by adjacent obstacles. These shape change
mechanisms do not occur for nearly rigid particles and can be accurately
modeled using the calibrated DP simulations with surface tension described in
Sec. IV.1, as shown in Fig. 5 (b) and (d).

For minimum obstacle separation $w_{\rm ob}/\sigma\gg 1$, single droplets flow
continuously through obstacle arrays without clogging. However, when $w_{\rm
ob}/\sigma<1$, single droplets can become permanently clogged within the
obstacle array. We assume that as the droplet moves through the obstacle
array, it successively interacts with each obstacle, with each interaction
contributing to the probability of clogging. This scenario suggests that
particle clogging is a Poisson process, confirmed by numerous studies prior
studies[52, 53, 54]. This means that the droplet can move a distance $r>0$
without clogging with cumulative probability distribution function,

$P(r)=e^{-r/\lambda},$ (16)

where $\lambda$ is clogging decay length such that $P(\lambda)=e^{-1}$.
Additional details concerning the measurements of $P(r)$ in the DP simulations
are included in Appendix C.

Several previous studies of clogging of soft particles in hoppers have shown
that softer particles (i.e. with smaller elastic moduli) are less likely to
form clogs [36, 37, 38], and thus would possess larger values of the clogging
decay length $\lambda$. Does this result also hold for the highly deformable
droplets considered in this study? To address this question, we carry out
numerical simulations of single droplets (using the SP and DP models) moving
through obstacle arrays in the clogging regime with $w_{\rm ob}/\sigma\leq
0.5$. We measure $\lambda$ as a function of $\Gamma$ (for the DP model),
$K_{\rm w}$ (for the SP model), and the deformability $\delta_{s}/\sigma$,
given by the fraction of the droplet diameter that deforms under its own
weight. We find that the wrapping mechanism plays an important role in
determining the clogging statistics for single droplets moving through
obstacle arrays. In particular, clogs can form when droplets wrap around an
obstacle as shown in Fig. 6 (a). The droplet is mainly supported by the
obstacle around which the droplet has wrapped, with additional stabilizing
forces from adjacent obstacles. For less deformable particles, clogs form when
droplets cannot squeeze through narrow channels between adjacent obstacles as
shown in Fig. 6 (b).

In Fig. 6 (c), we plot $\lambda/\sigma$ as a function of the line tension
$\Gamma$ from DP model simulations of single droplets moving through obstacle
arrays for $\sigma_{\rm ob}/\sigma=0.3$ and $w_{\rm ob}/\sigma=0.2$, $0.3$,
$0.4$, and $0.5$. We find that $\lambda/\sigma$ versus $\Gamma$ is non-
monotonic. At large $\Gamma$, where droplet shape deformation is small,
decreasing $\Gamma$ increases $\lambda$, which implies that softer particles
are less likely to clog. However, at small $\Gamma$, where droplet shape
deformation is large, we find the opposite behavior. In this regime,
decreasing $\Gamma$ decreases $\lambda$, which implies that softer particles
are more likely to clog. The nonmonotonic behavior of $\lambda(\Gamma)$ is
caused by the two different mechanisms for clogging (wrapping versus
squeezing) as shown in Fig. 6 (a) and (b). At large $\Gamma$, droplets become
clogged when squeezing between obstacles. It is easier for nearly rigid
droplets to squeeze through the gaps between obstacles as $\Gamma$ decreases,
which gives rise to decreased clogging probability and increased $\lambda$. In
contrast, at small $\Gamma$, highly deformable droplets can easily squeeze
through narrow gaps between obstacles, but they can wrap around an obstacle
and clog. Thus, increased deformability enhances the wrapping mechanism,
increasing the clogging probability, and decreasing $\lambda$. In Fig. 6 (c),
we also show that $\lambda$ increases with $w_{\rm ob}$ as expected, since
wider gaps between obstacles reduces clogging.

Figure 6: Images of a single droplet (shaded red) from DP simulations moving
through obstacle arrays (shaded blue) with random positions. Clogs form via
the (a) wrapping and (b) squeezing mechanisms. (c) Clogging decay length
$\lambda/\sigma$ plotted as a function of line tension $\Gamma$ for the
obstacle arrays in (a) and (b) with $\sigma_{\rm ob}/\sigma=0.3$ and $w_{\rm
ob}/\sigma=0.2$ (red), $0.3$ (green), $0.4$ (blue), and $0.5$ (black).

We also carry out SP simulations of single droplets moving through the same
obstacle arrays in Fig. 6 to understand the effects of droplet stiffness
versus deformability on clogging. In Fig. 7 (a), we plot $\lambda/\sigma$
versus the dimensionless stiffness $E_{\rm w}$ for the SP model;
$\lambda/\sigma$ decreases monotonically with increasing $E_{\rm w}$ over the
full range of $E_{\rm w}$, and as expected, $\lambda/\sigma$ increases with
$w_{\rm obs}/\sigma$. These results emphasize that the SP model can only
generate clogs via the squeezing mechanism, not the wrapping mechanism.

To directly compare the results for single droplet motion through obstacle
arrays using the SP and DP models, we consider the droplet deformation
$\delta_{s}$ caused by its own weight [36]. (See Fig. 7 (b).) To calculate
$\delta_{s}$, we initially turn off gravity and place an SP droplet or a DP
droplet with diameter $\sigma$ next to a wall. We then turn on gravity and
measure the distance $\delta_{s}$ that the droplet sinks. In Fig. 7 (c), we
plot $\lambda/\sigma$ versus $\delta_{s}/\sigma$ for both the SP and DP models
of single droplet motion through obstacle arrays. $\lambda/\sigma$ versus
$\delta_{s}$ is similar for the SP and DP models for $\delta_{s}/\sigma<0.25$,
but the results differ significantly for $\delta_{s}/\sigma>0.25$. At large
$\delta_{s}/\sigma$, overlaps between the obstacles and SP droplets are not
strongly penalized, clogs are highly unlikely, and
$\lambda/\sigma\rightarrow\infty$. However, at large $\delta_{s}/\sigma$, DP
droplets are extremely deformable, can wrap around obstacles, and thus
$\lambda/\sigma\rightarrow 0$. Comparing the results for the SP and DP models
allows us to assess the importance of shape deformability, particle softness,
and volume conservation in determining the clogging probability.

Figure 7: (a) Clogging decay length $\lambda/\sigma$ plotted as a function of
normalized stiffness $E_{\rm w}$ from SP model simulations of single droplets
moving through obstacle arrays with $\sigma_{\rm ob}/\sigma=0.3$ and $w_{\rm
ob}/\sigma=0.2$ (red), $0.3$ (green), $0.4$ (blue), and $0.5$ (black). The
arrow indicates increasing $w_{\rm ob}/\sigma$. (b) Schematic of the
deformation $\delta_{s}$ that a droplet undergoes under its own weight for the
(left) SP and (right) DP models. The solid (dashed) outlines of the droplet
represent the droplet before (after) gravity is added. (c) $\lambda/\sigma$
plotted as a function of $\delta_{s}/\sigma$ for both SP (black) and DP (red)
models for obstacle arrays with $w_{\rm ob}/\sigma=0.2$. Figure 8: (a) Average
speed of the center mass of the droplet in the direction of gravity
$\langle\nu_{g}\rangle/\nu_{t}$ plotted as a function of the dimensionless
line tension $\Gamma$ for DP simulations of continuous single-droplet flows
through obstacle arrays with minimum gap sizes $w_{\rm ob}/\sigma=1.0$ (red),
$1.1$ (green), $1.2$ (blue), and $1.3$ (black). The arrow indicates increasing
$w_{\rm ob}/\sigma$. (b) $\langle\nu_{g}\rangle/\nu_{t}$ plotted versus
$\Gamma$ for DP simulations of continuous single-droplet flows through
obstacle arrays with $w_{\rm ob}/\sigma=1.0$ and $\sigma_{\rm ob}/\sigma=0.3$
(red), $0.6$ (green), and $0.9$ (blue). The arrow indicates increasing
$\sigma_{\rm ob}$.

###  IV.3 Continuous single-droplet flows in obstacle arrays

We now investigate continuous flows of single droplets through obstacle arrays
(using the SP and DP models). We calculate the average droplet speed
$\langle\nu_{g}\rangle$ (over each droplet trajectory and multiple random
obstacle arrays, see Appendix D) as a function of the size and placement of
the obstacles. For simplicity, we set $V_{t}=0.03$ and $b_{0}=0$, but the
results also hold for non-zero $b_{0}$. (See Appendix E.) Intuitively, one
might expect that increasing the line tension $\Gamma$ would lead to increased
resistance to droplet shape changes, and thus decreases in
$\langle\nu_{g}\rangle$. However, in Fig. 8 (a), we find that
$\langle\nu_{g}\rangle$ increases with $\Gamma$ for $w_{\rm ob}/\sigma\geq 1$.
The decrease in $\langle\nu_{g}\rangle$ with decreasing $\Gamma$ occurs
because droplets with small surface tension can easily wrap around an
obstacle. When droplets wrap around an obstacle, the two ends of the droplet
extend in the flow direction until the droplet surface energy matches the work
done by the driving force. After the ends of droplet stop moving in the flow
direction, one end of the droplet moves in the opposite direction of the
driving force, which allows the droplet to unwrap. This wrapping and
unwrapping process significantly slows droplet movement and decreases
$\langle\nu_{g}\rangle$.

In Fig. 8 (b), we quantify continuous flows after fixing $w_{\rm
ob}/\sigma=1.0$ and varying $\sigma_{\rm ob}/\sigma$. We plot
$\langle\nu_{g}\rangle/\nu_{t}$ versus $\Gamma$ for $\sigma_{\rm
ob}/\sigma=0.3$, $0.6$, and $0.9$. We find that increasing $\sigma_{\rm
ob}/\sigma$ increases $\langle\nu_{g}\rangle$ at small $\Gamma$, indicating
that it is more difficult for droplets to wrap around larger obstacles. In the
large-$\Gamma$ limit, the single-droplet flows for all $\sigma_{\rm
obs}/\sigma$ approach the same $w_{\rm obs}$-dependent plateau value in
$\langle\nu_{g}\rangle$. Droplet wrapping does not occur in this regime since
droplet deformation would require large surface energy.

In Fig. 9 (a), we show $\langle\nu_{g}\rangle/\nu_{t}$ versus the droplet-
obstacle stiffness $E_{\rm w}$ using the same obstacle arrays in Fig. 8 for
the SP model. In contrast to the DP simulations, we find that
$\langle\nu_{g}\rangle/\nu_{t}$ decreases with increasing $E_{\rm w}$, which
indicates that softer particles flow more rapidly for the SP model. This
result stems from the fact that droplet softness in the SP model is
represented by larger allowed overlaps between the droplet and obstacles.
Therefore, in the $E_{\rm w}\rightarrow 0$ limit, droplets would flow through
obstacles without slowing down, resulting in
$\langle\nu_{g}\rangle/\nu_{t}\rightarrow 1$. In contrast, in the
$\Gamma\rightarrow 0$ limit in DP simulations, the droplets wrap around each
obstacle that they encounter. This difference in the flow behavior of single
SP and DP droplets emphasizes the importance of modeling explicit
deformability in droplet flows through obstacle arrays. In Fig. 9 (b), we plot
$\langle\nu_{g}\rangle/\nu_{t}$ versus droplet deformation $\delta_{s}/\sigma$
at $w_{\rm ob}/\sigma=1.1$ for both the DP and SP models. We find that
$\langle\nu_{g}\rangle/\nu_{t}$ for the SP and DP models converge to same
value in the rigid-droplet limit ($\delta_{s}/\sigma\rightarrow 0$), but they
possess diverging behavior for large droplet deformation,
$\delta_{s}/\sigma\rightarrow 1$.

Figure 9: (a) Average speed of the center mass of the droplet in the direction
of gravity $\langle\nu_{g}\rangle/\nu_{t}$ plotted as a function of $E_{\rm
w}$ for SP simulations of continuous single-droplet flows through obstacle
arrays with minimum gap sizes $w_{\rm ob}/\sigma=1.0$ (red), $1.1$ (green),
$1.2$ (blue), and $1.3$ (black). The arrow indicates increasing $w_{\rm
ob}/\sigma$. (b) $\langle\nu_{g}\rangle/\nu_{t}$ plotted as a function of
$\delta_{s}/\sigma$ for single-droplet flows through obstacle arrays with
$w_{\rm ob}/\sigma=1.1$ for both SP (black) and DP (red) models.

##  V Discussion and Conclusions

In this article, we carried out coordinated experiments and numerical
simulations of gravity-driven flows of single droplets in narrow channels and
obstacle arrays. We first developed deformable particle (DP) simulations with
line tension in 2D and calibrated them to the experiments. The DP simulations
recapitulate the experimental results for the droplet shape and speed for
flows within narrow channels with an accuracy of less than $10\%$. Using the
calibrated DP simulations, we find several important results concerning
single-droplet flows in narrow channels and obstacle arrays. First, the
droplet speed possesses nonmonotonic behavior as the droplet traverses the
narrow channels. The droplet speed can even overshoot the terminal speed and
the overshoot is enhanced as the gravitational driving is decreased. Second,
we showed that the clogging probability is nonmonotonic with the line tension
for droplets flowing through obstacle arrays. The nonmonotonic behavior is
caused by two distinct clogging mechanisms (droplet wrapping and squeezing) in
obstacle arrays. Stiffer droplets become trapped when they cannot squeeze in
between obstacles, while highly deformable droplets clog when they wrap around
obstacles, resulting in high clogging probabilities for droplets with both
large and small deformabilities. We also compared the DP simulation results
for flows in obstacle arrays to those for the frequently used soft particle
(SP) model, which does not include explicit shape deformability. The DP and SP
simulations of single-droplet flows in obstacle arrays are similar in the
rigid-droplet limit, but diverge for deformable droplets since the SP model
can only model clogging via the squeezing mechanism (and not the wrapping
mechanism).

These results suggest several promising future research directions. First, the
current numerical simulations employ the DP model with line tension in 2D.
When calibrating the 2D DP model to the experimental results, we obtain values
for the dimensionless line tension $\Gamma^{*}$ that are a factor of $\sim 3$
smaller than the dimensionless surface tension $\Gamma_{\rm exp}$ in
experiments. Thus, in future work we will develop the 3D DP model [49] with
surface tension to describe quasi-2D flows of droplets in narrow channels and
obstacle arrays. Second, in the current work, we do not quantitatively model
the background fluid. However, the hydrodynamics of the background fluid can
influence the droplet dynamics, such as causing the droplets to follow stream
lines, which can prevent droplet-obstacle contact. Therefore, in the future,
we can include the effects of Stokesian dynamics of the background fluid on
the droplet motion.

Third, we observe in experiments that droplets undergoing sufficiently large
deformations break up into smaller droplets. However, our current numerical
simulations do not account for this phenomenon. In future studies, we will
develop a quantitative framework to predict the conditions under which droplet
breakup occurs and use it to determine the size distribution of daughter
droplets emerging from obstacle arrays. This improved understanding of droplet
breakup will be valuable in microfluidic applications [55, 56], where
generating droplets of a specific size is required [57]. We also plan to
investigate multi-droplet flows through obstacle arrays. The efficiency of the
DP model will enable us to simulate systems consisting of thousands of
droplets, and explore under what conditions droplet interactions give rise to
coalescence, where two droplets in contact merge to form a larger droplet
[58]. Droplet breakup and coalescence play competing roles in determining the
equilibrium size distribution of droplets moving through obstacle arrays [59].
In future studies, we will investigate their combined effects to
quantitatively predict the droplet size distribution exiting from the obstacle
arrays.

Figure 10: The ratio $\nu_{g}/(\nu_{t}-\nu_{g})$ plotted versus distance
$d/\sigma$ from the chamber wall for the experiments. The droplet has diameter
$\sigma=0.98~{}\rm mm$, tilt angle $\theta=5^{\circ}$, and terminal velocity
far away from the wall, $v_{t}=0.65~{}\rm mm/s$. An image from the experiment
is shown in the inset. The dashed line has slope $4/3$.

## Appendix A

In this Appendix, we show experimental measurements of the power-law exponent
$\alpha$ that appears in Eq. 6 for the distance-dependent drag coefficient. In
the low Reynolds number (overdamped) regime, Eq. 6 can be rewritten as:

$\nu_{g}\zeta=\nu_{g}b_{\infty}\left(1+\frac{b_{0}}{b_{\infty}}\left(\frac{\sigma}{d}\right)^{\alpha}\right)\approx-
Mg=\nu_{t}b_{\infty}.$ (17)

After rearranging Eq. 17, we find

$b_{0}\left(\frac{\sigma}{d}\right)^{\alpha}\approx\left(\frac{\nu_{t}}{\nu_{g}}-1\right)b_{\infty},$
(18)
$\alpha\log_{10}\left(\frac{d}{\sigma}\right)\approx\log_{10}\left(\frac{\nu_{g}}{\nu_{t}-\nu_{g}}\right)+\log_{10}\left(\frac{b_{\infty}}{b_{0}}\right).$
(19)

In Fig. 10, we plot $v_{g}/(v_{t}-v_{g})$ versus $d/\sigma$ on a
$\log_{10}$-$\log_{10}$ scale, where $d$ is the minimum distance between
chamber wall and the droplet. The droplet shown in the inset has diameter
$\sigma=0.98~{}\rm mm$, tilt angle $\theta=5^{\circ}$, and terminal velocity
far from wall $v_{t}=0.65~{}\rm mm/s$. We find that the experimental data in
Fig. 10 is well-fit by a power-law exponent $\alpha=4/3$.

## Appendix B

In this Appendix, we investigate the effects of the number of vertices $N_{v}$
and perimeter elasticity $K_{l}$ on the results for droplet flows through
narrow channels using the DP simulations. The DP model with
$N_{v}\rightarrow\infty$ and $K_{l}\rightarrow 0$ is the most accurate for
describing droplet flows. However, to limit computational time, we choose
$N_{v}$ such that it is sufficiently large that further increases do not
significantly change the results. In Fig. 11, we show the droplet speed
profile, $\nu_{g}/\nu_{t}$ versus $h/\sigma$, for $N_{v}=32$, $64$, $128$,
$256$, $512$, and $1024$. We find that for $N_{v}\leq 64$, the velocity
profile is quite noisy. For $N_{v}\geq 128$, the velocity profile is smooth on
a linear scale. For the DP simulations presented here, we use $N_{v}=256$.

Figure 11: Droplet center of mass speed $\nu_{g}/\nu_{t}$ plotted as a
function of distance $h/\sigma$ from the orifice for DP simulations of a
droplet flowing through a narrow channel with $w=0.7\sigma$, $\Gamma=1.08$,
and $N_{v}=32$, $64$, $128$, $256$, $512$, and $1024$. The inset is a close-up
of $\nu_{g}/\nu_{t}$ near the first minimum in the speed at
$h/\sigma\approx-0.25$ to highlight the results for varying $N_{v}$, from $32$
(blue) to $1024$ (red).

For the DP simulations to accurately model line tension, it is important that
$K_{l}/\Gamma\rightarrow 0$. However, the perimeter elasticity in Eq. 1 is
important for ensuring the stability of the DP simulations. If the DP model
only includes line tension without perimeter elasticity, the droplet vertices
can aggregate in the direction of gravity, which gives rise to an uneven
distribution of vertices around the droplet perimeter and can cause the
droplet to overlap the side walls. Thus, we need to select a non-zero value of
$K_{l}$ to uniformly distribute vertices on the droplet surface, but small
enough such that the results for the DP model are dominated by line tension
and do not depend on $K_{l}$. In Fig. 12, we show the root-mean-square
deviations in the droplet speed $\Delta_{v}$ between the DP simulations and
experiments of flow in narrow channels plotted as a function of
$K_{l}/\Gamma$. $\Delta_{v}\sim 0.08$ reaches a plateau for
$K_{l}/\Gamma\lesssim 10^{-2}$. In the present work, we set
$K_{l}/\Gamma=10^{-3}$.

Figure 12: Root-mean-square deviation in the droplet speed $\Delta_{\nu}$
between the DP simulations and experiments for a droplet flowing through a
narrow orifice with width $w/\sigma=0.9$ plotted as a function of
$K_{l}/\Gamma$. $\Delta_{v}\sim 0.08$ reaches a plateau for
$K_{l}/\Gamma<10^{-2}$.

## Appendix C

In this Appendix, we include additional details concerning the definition of a
permanent clog and the clogging decay length $\lambda$ for droplets moving in
obstacle arrays. In the simulations, a clog occurs when the droplet kinetic
energy $E_{k}$ (normalized by the kinetic energy at the terminal speed) decays
below a threshold $E_{k}/E_{t}<10^{-17}$. In Fig. 13 (a), we show examples of
$E_{k}/E_{t}$ versus time $t/t_{0}$ for cases when the droplet clogs ($w_{\rm
ob}/\sigma=0.3$) and does not clog ($w_{\rm ob}/\sigma=1.1$).

To determine the clogging decay length $\lambda$, we first measure the
cumulative probability $P$ that the droplet does not clog as a function of the
droplet displacement $r$ in the obstacle array. We run $N_{s}=50$ simulations
for a total time $T=2\times 10^{5}t_{0}$ with different randomized initial
droplet positions and randomized obstacle arrays. $P(r)=1-N_{\rm
clog}(r)/N_{s}$, where $N_{\rm clog}$ is the number of times that the droplet
clogs before it moves a distance $r$. In Fig. 13 (b), we plot $P$ versus
$r/\sigma$ for DP simulations with $\Gamma\approx 1.09$, $\sigma_{\rm
ob}/\sigma=0.3$, and $w_{\rm ob}/\sigma=0.2$. At the beginning of each
simulation, $r=0$ and $P=1$. As each droplet traverses the obstacle array, the
cumulative probability $P$ of not clogging decreases exponentially with
increasing $r$.

Figure 13: (a) Droplet kinetic energy $E_{k}$, normalized by the kinetic
energy at the terminal speed $E_{t}$, plotted as a function of time $t/t_{0}$
for droplets moving in obstacle arrays with $\Gamma=0.5$, $\sigma_{\rm
ob}/\sigma=0.3$, and both $w_{\rm ob}/\sigma=1.1$ (red solid line) and $0.3$
(black solid line). For $w_{\rm ob}/\sigma=0.3$, $E_{k}/E_{t}$ decays
exponentially toward the clogging threshold $10^{-17}$. For $w_{\rm
ob}/\sigma=1.1$, the droplet flows continuously without clogging. (b)
Cumulative probability $P(r)$ of the droplet not forming a clog plotted as a
function of the displacement in the obstacle array $r/\sigma$ using DP
simulations with $\Gamma\approx 1.09$, $\sigma_{\rm ob}/\sigma=0.3$, and
$w_{\rm ob}/\sigma=0.2$. The red solid curve represents $P=\exp(-r/\lambda)$,
where $\lambda/\sigma=24\pm 1$ is the clogging decay length.

## Appendix D

In this Appendix, we show that the DP simulations of single droplet flows in
obstacle arrays can reach a continuous flow regime (with no permanent clogs)
for minimum obstacle spaceing $w_{\rm ob}/\sigma>1$. In Fig. 14, we plot the
droplet displacement $r$ versus time $t$ for DP simulations with
$\Gamma=0.04$-$1$. We determine the steady state speed $\nu_{g}$ for a single-
droplet trajectory in a given obstacle array by measuring the slope of
$r/\sigma$ versus $t/t_{0}$. $\langle\nu_{g}\rangle$ is averaged over $10$
random obstacle arrays.

Figure 14: Droplet displacement $r$ in the direction of gravity $r/\sigma$
plotted as a function of time $t/t_{0}$ as the droplet traverses an obstacle
array with $w_{\rm ob}/\sigma=1.1$ and $\sigma_{\rm ob}/\sigma=0.3$. The arrow
indicates increasing $\Gamma$ from $0.04$ (blue) to $1$ (red).

## Appendix E

In this Appendix, we show that the choice of the near-wall drag coefficient
$b_{0}$ does not qualitatively affect the results in Sec.IV.3. As discussed
previously, the SP model does not include distance-dependent drag forces
(Eq.6). In Fig. 9 (c) and Fig.7 (b), we set $b_{0}=0$ in the DP simulations so
that when we compare the DP and SP simulations, both have the same form for
the drag force. In Fig. 15, we plot $\langle\nu_{g}\rangle/\nu_{t}$ versus
$\Gamma$ for single droplets flowing through obstacle arrays with $w_{\rm
ob}/\sigma=1.0$. All other parameters are the same as those in Fig.8 (a),
except the near-wall drag coefficient varies: $b_{0}/b_{\infty}=0$, $0.2$, and
$0.4$. We find that increasing $b_{0}/b_{\infty}$ decreases
$\langle\nu_{g}\rangle/\nu_{t}$, but it monotonically increases with $\Gamma$
for all $b_{0}/b_{\infty}$.

Figure 15: Average speed of the center mass of the droplet in the direction of
gravity $\langle\nu_{g}\rangle$ normalized by the terminal velocity $\nu_{t}$
plotted as a function of the dimensionless line tension $\Gamma$ for obstacle
arrays with $w_{\rm ob}/\sigma=1.0$. All other parameters are the same as
those used in Fig.8 (a), except we vary the near-wall drag coefficient:
$b_{0}/b_{\infty}=0$ (red), $0.2$ (blue), and $0.4$ (black). The arrow
indicates increasing $b_{0}/b_{\infty}$.

## Acknowledgments

We acknowledge support from NSF Grants No. CBET-2002782 (Y. C., S.S., and C.
S. O.), No. CBET-2002815 (B. L., P. H., and E. R. W.), No. CBET-2306371 (D.
M.), and No. CBET-2002797 (M. D. S.). This work was also supported by the High
Performance Computing facilities operated by Yale’s Center for Research
Computing.

## References

  * Chung and McClements [2014] C. Chung and D. J. McClements, _Food Structure_ , 2014, 1, 106–126.
  * Bai _et al._ [2021] L. Bai, S. Huan, O. J. Rojas and D. J. McClements, _Journal of Agricultural and Food Chemistry_ , 2021, 69, 8944–8963.
  * Freitas _et al._ [2006] S. Freitas, G. Hielscher, H. P. Merkle and B. Gander, _Ultrasonics Sonochemistry_ , 2006, 13, 76–85.
  * Albert _et al._ [2019] C. Albert, M. Beladjine, N. Tsapis, E. Fattal, F. Agnely and N. Huang, _Journal of Controlled Release_ , 2019, 309, 302–332.
  * Krawczyk _et al._ [1991] M. A. Krawczyk, D. T. Wasan and C. Shetty, _Industrial & Engineering Chemistry Research_, 1991, 30, 367–375.
  * Umar _et al._ [2018] A. A. Umar, I. B. M. Saaid, A. A. Sulaimon and R. B. M. Pilus, _Journal of Petroleum Science and Engineering_ , 2018, 165, 673–690.
  * Cerff _et al._ [2021] B. Cerff, D. Key and B. Bladergroen, _Sustainability_ , 2021, 13, 12339.
  * Bick _et al._ [2021] A. D. Bick, J. W. Khor, Y. Gai and S. K. Y. Tang, _Proceedings of the National Academy of Sciences_ , 2021, 118, e2017822118.
  * Sohrabi _et al._ [2020] S. Sohrabi, N. kassir and M. Keshavarz Moraveji, _RSC Advances_ , 2020, 10, 27560–27574.
  * Guo _et al._ [2012] M. T. Guo, A. Rotem, J. A. Heyman and D. A. Weitz, _Lab Chip_ , 2012, 12, 2146–2155.
  * Casadevall i Solvas and deMello [2011] X. Casadevall i Solvas and A. deMello, _Chemical Communications_ , 2011, 47, 1936–1942.
  * Ding _et al._ [2020] Y. Ding, P. D. Howes and A. J. deMello, _Analytical Chemistry_ , 2020, 92, 132–149.
  * Schneider _et al._ [2013] T. Schneider, J. Kreutz and D. T. Chiu, _Analytical Chemistry_ , 2013, 85, 3476–3482.
  * Xue _et al._ [2015] C. Xue, J. Wang, Y. Zhao, D. Chen, W. Yue and J. Chen, _Micromachines_ , 2015, 6, 1794–1804.
  * Baret _et al._ [2009] J.-C. Baret, O. J. Miller, V. Taly, M. Ryckelynck, A. El-Harrak, L. Frenz, C. Rick, M. L. Samuels, J. B. Hutchison, J. J. Agresti, D. R. Link, D. A. Weitz and A. D. Griffiths, _Lab Chip_ , 2009, 9, 1850–1858.
  * Beattie _et al._ [2014] W. Beattie, X. Qin, L. Wang and H. Ma, _Lab Chip_ , 2014, 14, 2657–2665.
  * Joensson and Andersson Svahn [2012] H. N. Joensson and H. Andersson Svahn, _Angewandte Chemie International Edition_ , 2012, 51, 12176–12192.
  * Kaminski _et al._ [2016] T. S. Kaminski, O. Scheler and P. Garstecki, _Lab Chip_ , 2016, 16, 2168–2187.
  * Matuła _et al._ [2020] K. Matuła, F. Rivello and W. T. S. Huck, _Advanced Biosystems_ , 2020, 4, 1900188.
  * Thomas and Durian [2016] C. C. Thomas and D. J. Durian, _Physical Review E_ , 2016, 94, 022901.
  * Pascot _et al._ [2020] A. Pascot, N. Gaudel, S. Antonyuk, J. Bianchin and S. Kiesgen De Richter, _Chemical Engineering Science_ , 2020, 224, 115749.
  * Tang and Behringer [2011] J. Tang and R. P. Behringer, _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , 2011, 21, 041107.
  * Beverloo _et al._ [1961] W. Beverloo, H. Leniger and J. van de Velde, _Chemical Engineering Science_ , 1961, 15, 260–269.
  * Cheng _et al._ [2022] Y. Cheng, J. D. Treado, B. F. Lonial, P. Habdas, E. R. Weeks, M. D. Shattuck and C. S. O’Hern, _Soft Matter_ , 2022, 18, 8071–8086.
  * Zuriguel _et al._ [2011] I. Zuriguel, A. Janda, A. Garcimartín, C. Lozano, R. Arévalo and D. Maza, _Physical Review Letters_ , 2011, 107, 278001.
  * Lozano _et al._ [2012] C. Lozano, A. Janda, A. Garcimartín, D. Maza and I. Zuriguel, _Physical Review E_ , 2012, 86, 031306.
  * Gao _et al._ [2021] G.-J. J. Gao, F.-L. Yang, M. C. Holcomb and J. Blawzdziewicz, _Physical Review E_ , 2021, 103, 062904.
  * Alessio [2021] S. D. Alessio, _European Journal of Physics_ , 2021, 42, 065808.
  * Gai _et al._ [2022] Y. Gai, A. Montessori, S. Succi and S. K. Y. Tang, _Physical Review Fluids_ , 2022, 7, 080501.
  * Gao _et al._ [2019] G.-J. J. Gao, J. Blawzdziewicz, M. C. Holcomb and S. Ogata, _Granular Matter_ , 2019, 21, 25.
  * Bick _et al._ [2021] A. D. Bick, J. W. Khor, Y. Gai and S. K. Y. Tang, _Proceedings of the National Academy of Sciences_ , 2021, 118, e2017822118.
  * Reichhardt and Reichhardt [2018] C. Reichhardt and C. J. O. Reichhardt, _Physical Review Letters_ , 2018, 121, 068001.
  * Nguyen _et al._ [2017] H. T. Nguyen, C. Reichhardt and C. J. O. Reichhardt, _Physical Review E_ , 2017, 95, 030902.
  * Reichhardt and Reichhardt [2021] C. Reichhardt and C. J. O. Reichhardt, _Soft Matter_ , 2021, 17, 1548–1557.
  * Péter _et al._ [2018] H. Péter, A. Libál, C. Reichhardt and C. J. O. Reichhardt, _Scientific Reports_ , 2018, 8, 10252.
  * Hong _et al._ [2017] X. Hong, M. Kohne, M. Morrell, H. Wang and E. R. Weeks, _Physical Review E_ , 2017, 96, 062605.
  * Tao _et al._ [2021] R. Tao, M. Wilson and E. R. Weeks, _Physical Review E_ , 2021, 104, 044909.
  * Alborzi _et al._ [2022] S. Alborzi, B. G. Clark and S. M. Hashmi, _Soft Matter_ , 2022, 18, 4127–4135.
  * Wang _et al._ [2021] J. Wang, B. Fan, T. Pongó, K. Harth, T. Trittel, R. Stannarius, M. Illig, T. Börzsönyi and R. Cruz Hidalgo, _Soft Matter_ , 2021, 17, 4282–4295.
  * Hirt and Nichols [1981] C. Hirt and B. Nichols, _Journal of Computational Physics_ , 1981, 39, 201–225.
  * Silva _et al._ [2024] D. P. F. Silva, R. C. V. Coelho, I. Pagonabarraga, S. Succi, M. M. Telo da Gama and N. A. M. Araújo, _Soft Matter_ , 2024, 20, 2419–2441.
  * Coelho _et al._ [2023] R. C. V. Coelho, D. P. F. Silva, A. M. R. Maschio, M. M. Telo da Gama and N. A. M. Araújo, _Physics of Fluids_ , 2023, 35, 013304.
  * Breinlinger _et al._ [2013] T. Breinlinger, P. Polfer, A. Hashibon and T. Kraft, _Journal of Computational Physics_ , 2013, 243, 14–27.
  * Gingold and Monaghan [1977] R. A. Gingold and J. J. Monaghan, _Monthly Notices of the Royal Astronomical Society_ , 1977, 181, 375–389.
  * Monaghan [2005] J. J. Monaghan, _Reports on Progress in Physics_ , 2005, 68, 1703–1759.
  * Morris [2000] J. P. Morris, _International Journal for Numerical Methods in Fluids_ , 2000, 33, 333–353.
  * Boromand _et al._ [2018] A. Boromand, A. Signoriello, F. Ye, C. S. O’Hern and M. D. Shattuck, _Physical Review Letters_ , 2018, 121, 248003.
  * Boromand _et al._ [2019] A. Boromand, A. Signoriello, J. Lowensohn, C. S. Orellana, E. R. Weeks, F. Ye, M. D. Shattuck and C. S. O’Hern, _Soft Matter_ , 2019, 15, 5854–5865.
  * Wang _et al._ [2021] D. Wang, J. D. Treado, A. Boromand, B. Norwick, M. P. Murrell, M. D. Shattuck and C. S. O’Hern, _Soft Matter_ , 2021, 17, 9901–9915.
  * Durian [1995] D. J. Durian, _Physical Review Letters_ , 1995, 75, 4780–4783.
  * Durian [1997] D. J. Durian, _Physical Review E_ , 1997, 55, 1739–1751.
  * Hong _et al._ [2022] X. Hong, K. W. Desmond, D. Chen and E. R. Weeks, _Physical Review E_ , 2022, 105, 014603.
  * Thomas and Durian [2015] C. C. Thomas and D. J. Durian, _Physical Review Letters_ , 2015, 114, 178001.
  * Marin _et al._ [2018] A. Marin, H. Lhuissier, M. Rossi and C. J. Kähler, _Physical Review E_ , 2018, 97, 021102.
  * Link _et al._ [2004] D. R. Link, S. L. Anna, D. A. Weitz and H. A. Stone, _Physical Review Letters_ , 2004, 92, 054503.
  * Chung _et al._ [2009] C. Chung, K. H. Ahn and S. J. Lee, _Journal of Non-Newtonian Fluid Mechanics_ , 2009, 162, 38–44.
  * Tazikeh Lemeski _et al._ [2024] A. Tazikeh Lemeski, S. M. Seyyedi, M. Hashemi-Tilehnoee and A. S. Naeimi, _Nature Scientific Reports_ , 2024, 14, 13324.
  * Dolganov _et al._ [2021] P. V. Dolganov, A. S. Zverev, K. D. Baklanova and V. K. Dolganov, _Physical Review E_ , 2021, 104, 014702.
  * Barai and Mandal [2019] N. Barai and N. Mandal, _International Journal of Multiphase Flow_ , 2019, 111, 1–15.