# Dynamics of a droplet migration in oscillatory and pulsating microchannel
flows and prediction and uncertainty quantification of its lateral equilibrium
position using Multi Fidelity Gaussian processes

Ali Lafzi  Department of Agricultural and Biological Engineering, Purdue
University, West Lafayette, Indiana 47907, USA  Sadegh Dabiri  Department of
Agricultural and Biological Engineering, Purdue University, West Lafayette,
Indiana 47907, USA

###### Abstract

Dynamics of a droplet in oscillatory and pulsating flows of a Newtonian fluid
in a microchannel has been studied numerically. The effects of oscillation
frequency, surface tension, and channel flow rate have been explored by
simulating the drop within a microchannel. These types of flows introduce new
equilibrium positions for the drop compared to steady flows with similar
conditions. The simulation results are very sensitive to the grid resolution
due to the unsteady behavior of the base flow. Therefore, a set of fine grids
have been used in this study to capture the physics of this problem more
accurately. However, these fine grids make the computations significantly
expensive. Therefore, a Multi Fidelity Gaussian processes method with two
levels of fidelity has been used to predict the results of the remaining fine-
grid simulations along with their uncertainties based on their correlations
with those of the coarse-grid cases over a wide range of input parameters.

##  I Introduction

Transient dynamics of droplets and their lateral equilibrium positions in
microchannels depend on several parameters including the rheology of the
carrier fluid, the geometry of channel, the channel flow rate, and the drop
characteristics such as shape, size, interfacial tension, etc, which can be
deployed for numerous biological and biomedical applications Gossett  _et al._
(2010); Toner and Irimia (2005); Gossett  _et al._ (2012); Karimi  _et al._
(2013); Di Carlo  _et al._ (2007).

Most of the effective parameters underlying the physics of particle motion in
microchannels and label-free approaches are inherent to the system, particle
properties, and are out of human’s complete control. Therefore, adding an
extra parameter as a tool for more direct control over the dynamics of
particles and better manipulating them can be very beneficial; especially,
coming up with equilibrium positions between the channel center and the wall
can be very useful for separation purposes Marson  _et al._ (2018). The other
restriction in inertial microfluidics is the limitation of working with
particle sizes of a few microns or larger only Mutlu  _et al._ (2018). This is
because of the inverse relationship between the particle size and its
traveling distance until focusing, which means that the separation of smaller
particles requires longer channels Lan and Khismatullin (2012).

Exploiting an oscillatory flow in the channel has been recommended as a
solution to address the aforementioned issues Mutlu  _et al._ (2018); Lafzi
_et al._ (2020). The oscillation frequency as the new extra parameter does not
have the disadvantage of adversely affecting the biological properties of
cells like most of the active methods Erlandsson and Robinson (2011). Besides,
this type of flow as an alternative to the traditional microfluidics has
attracted researchers in recent years. Chaudhury  _et al._ Chaudhury  _et al._
(2016) have observed a complicated trajectory of the droplet in the lateral
direction suspended in an oscillatory flow as opposed to the smooth cross-
stream migration under steady flow conditions. Pawłowska  _et al._ Pawłowska
_et al._ (2017) have reported the use of hydrogel nanofilaments as an
appropriate substitute for the long and deformable macromolecules and studied
their dynamics in an oscillatory microchannel flow. They have shown that the
final position of these particles fluctuates around the flow axis Pawłowska
(2018). Sarkar and Schowalter Sarkar and Schowalter (2000) have investigated
the dynamics of a viscoelastic drop in steady and oscillating extensional
flows and have shown that although their deformation behaviors are naturally
dissimilar, their maximum values are close in the long-time limit.

The presented study investigates different aspects of the droplet migration in
oscillatory and pulsating microchannel flows numerically and compares them
with those in the steady ones. The main challenge of performing this numerical
work is using a sufficient grid resolution to resolve the underlying physics
more accurately. Since the background flow is unsteady, the mesh has to be
fine enough to capture the velocity gradients and lift forces acting on the
drop correctly. Consequently, the transient dynamics of the drop and its final
equilibrium position, which is the main goal of the presented study, are all
affected by this choice. Therefore, after carrying out the mesh dependence
study, a very fine grid has been selected to enable getting the realistic
results we are interested in. However, running simulations with these grids
can take a few months in some cases. Hence, it is impossible to get all the
data in the wide range of input parameters we are looking for by just running
all these simulations. An alternative approach is to produce some new data
using a coarser grid resolution and implement the Multi Fidelity Gaussian
processes (MFGP) algorithm to predict the outputs of the finer mesh.
Perdikaris  _et al._ Perdikaris  _et al._ (2017) have proposed a probabilistic
framework and a recursive version of MFGP. This method has been tested in
several benchmark problems involving both synthetic and real multi-fidelity
datasets such as the one employed by Babaee  _et al._ Babaee  _et al._ (2016).
The implementation of the recursive MFGP and the obtained results in this
study are elaborated in the following sections.

##  II Methodology

###  II.1 The fluid dynamics simulations

A single Newtonian droplet has been placed in a laminar flow of an
incompressible Newtonian fluid in a rectangular microchannel with a square
cross section. A schematic of the configuration is illustrated in Fig. 1. The
density and viscosity ratios ($\eta$ and $\lambda$, respectively) are set to 1
for most of the simulations. The front-tracking method Unverdi and Tryggvason
(1992) is used to update the interface position. In this method, the main
governing equations are solved in the fixed Eulerian grid, and the obtained
information is used to update the properties across the droplet surface
containing thousands of moving Lagrangian elements. The governing equations to
be solved in the entire computational domain are the following:

$\nabla\cdot\textbf{u}=0,$ (1) $\frac{\partial(\rho\textbf{u})}{\partial
t}+\nabla\cdot(\rho\textbf{uu})=-\nabla
p+\nabla\cdot\bm{\tau}+\iint{}\gamma\kappa\delta(\bm{x}-\bm{x_{i}})\textbf{n}dA,$
(2)

where ${\rho}$ is the density of the fluid, $p$ represents the pressure, u is
the velocity vector, $t$ is the time,
$\bm{\tau}=\mu(\nabla\textbf{u}+\nabla\textbf{u}^{T})$ is the stress tensor in
which $\mu$ is the fluid viscosity, $\kappa$ is the curvature at the
interface, $\gamma$ is the interfacial tension, $\delta$ is the Dirac delta
function, $x$ is an arbitrary location in the whole computational domain,
$x_{i}$ is such position on the drop interface, and n is the unit normal
vector to a point on the interface. The given delta function is defined as:

$\delta(x)=\tilde{D}(x)\tilde{D}(y)\tilde{D}(z),$ (3)
$\tilde{D}(x)=\frac{1}{4\Delta}(1+\cos(\frac{\pi}{2\Delta}(x))),|x|\leq
2\Delta,$ (4)

where $\Delta$ is the constant Eulerian grid size.

To generate the oscillatory flow, a cosine wave of pressure gradient with a
constant amplitude (in the form of ${P_{0}}cos({\omega}t)$) is applied along
the channel ($x$ direction) to change the direction of the flow symmetrically.
For the pulsating case, this pressure gradient has the shape of
${P_{0}}(a+bcos({\omega}t))$, where $a$ and $b$ are the weights of the steady
and oscillatory components, respectively, and $a+b=1$. The periodic boundary
condition is applied in the $x$ direction, and the no-slip condition is
applied on the walls in the $y$ and $z$ directions. $W$ and $U_{0}$ (channel
centerline velocity corresponding to the steady case) are used as the
characteristic length and velocity, respectively. In other words,
$x^{*}={\frac{x}{W}}$, $u^{*}={\frac{u}{U_{0}}}$,
$t^{*}={\frac{t}{{\frac{W}{U_{0}}}}}$,
$P^{*}={\frac{P}{\mu{\frac{U_{0}}{W}}}}$,
$T^{*}={\frac{T}{{\frac{W}{U_{0}}}}}$ (where $T$ is the period), and
$\omega^{*}={\frac{2\pi}{T^{*}}}$. Three dimensionless parameters describe and
affect the motion of the drop: (i) Reynolds number
$Re={\frac{{\rho}U_{0}2W}{\mu}}$, expressing the ratio between inertial forces
to viscous ones (ii) Capillary number $Ca={\frac{\mu U_{0}}{\gamma}}$, which
denotes the ratio of viscous stress to the interfacial tension, where high
$Ca$ corresponds to a highly deformable drop (iii) The dimensionless
oscillation frequency ($\omega^{*}$). The effect of the frequency can also be
embedded in the Womersley number $Wo=(\frac{\omega
W^{2}}{\nu})^{\frac{1}{2}}$, where $\nu$ is the kinematic viscosity of the
fluid. This number compares the transient inertial effects to viscous forces
Dincau  _et al._ (2020). The blockage ratio of the drop (${\frac{a}{W}}$) is
constant and equals to $0.3$ for most of the cases studied here. The drop is
assumed to have a spherical initial shape and is released at
${\frac{y}{W}}=0.55$ and ${\frac{z}{W}}=1.07$. The initial location of the
drop is arbitrary since it does not alter its equilibrium position Pan  _et
al._ (2016); Lan and Khismatullin (2012); Chaudhury  _et al._ (2016); Razi and
Pourghasemi (2017). The axes of symmetry have been avoided. A fine grid of
$196\times 114\times 114$ (high fidelity level) and a coarse grid of
$128\times 76\times 76$ (low fidelity level) in the $x$, $y$, and $z$
directions, respectively, is used for the simulations having $Re=10$. The
simulations with higher $Re$ numbers require even a finer grid resolution
($256\times 152\times 152$). Since this adds another level of fidelity to the
problem, we have omitted the effect of $Re$ in the second part of the paper
for more simplicity. Therefore, the input parameters that affect the output of
interest, which is the distance of the equilibrium position of the drop from
the channel center, are reduced to the Capillary number ($Ca$) and frequency
($\omega^{*}$).

Figure 1: Schematic of the problem setup

###  II.2 Multi fidelity GP

Gaussian processes (GP) is an example of a continuous stochastic process in
which, a multi-dimensional Gaussian distribution is assigned to the function
values at the input points as the prior knowledge. In other words, the value
of the function that we are aiming to predict follows a Gaussian probability
density at each point:

$\mathbf{f}_{1:n}|\mathbf{x}_{1:n}\sim\mathcal{N}\left(\mathbf{m}(\mathbf{x}_{1:n}),\mathbf{K}(\mathbf{x}_{1:n},\mathbf{x}_{1:n})\right),$
$\mathbf{m}(\mathbf{x}_{1:n})=\left(m(\mathbf{x}_{1}),\dots,m(\mathbf{x}_{n})\right),$
$\mathbf{K}(\mathbf{x}_{1:n},\mathbf{x}_{1:n})=\left(\begin{array}[]{ccc}k(\mathbf{x}_{1},\mathbf{x}_{1})&\dots&k(\mathbf{x}_{1},\mathbf{x}_{n})\\\
\vdots&\ddots&\vdots\\\
k(\mathbf{x}_{n},\mathbf{x}_{1})&\dots&k(\mathbf{x}_{n},\mathbf{x}_{n})\end{array}\right).$

Where $\mathbf{m}$ is the mean vector, and $\mathbf{K}$ is the covariance
matrix whose elements denote the correlations between the input points. These
correlations are quantified by a covariance (kernel) function that encodes our
prior belief and knowledge about the targeted function like its continuity,
differentiability, etc. The choice of this function is very important for a
good prediction; otherwise, we need to have sufficient amount of data to
compensate this lack of knowledge. One of the most typical examples for this
function is the squared exponential (SE) kernel given below:

$k(\mathbf{x},\mathbf{x}^{\prime})=v\exp\left\\{-\frac{1}{2}\sum_{i=1}^{d}\frac{(x_{i}-x_{i}^{\prime})^{2}}{\ell_{i}^{2}}\right\\},$

In which, $v$ is known as the signal strength, and $\ell_{i}$ is the length
scale of the $i$-th input dimension of the GP. In this study, we also use the
SE function as we believe the function of interest is continuous and
infinitely differentiable.

If we have $n$ observations (training set), consisting of
$\mathbf{x}_{1:n}=(\mathbf{x}_{1},\dots,\mathbf{x}_{n})$ and
$\mathbf{y}_{1:n}=(y_{1},\dots,y_{n})$, and $n^{*}$ test points
($\mathbf{x}^{*}_{1:n^{*}}$) that we would like to predict the function values
at ($\mathbf{f}^{*}_{1:n^{*}}$), we use the following model according to the
definition of GP as our likelihood:

$p(\mathbf{f}_{1:n},\mathbf{f}^{*}_{1:n^{*}}|\mathbf{x}_{1:n},\mathbf{x}^{*}_{1:n^{*}})=\mathcal{N}\left(\left(\begin{array}[]{c}\mathbf{f}_{1:n}\\\
\mathbf{f}^{*}_{1:n^{*}}\end{array}\right)\middle|\left(\begin{array}[]{c}\mathbf{m}(\mathbf{x}_{1:n})\\\
\mathbf{m}(\mathbf{x}^{*}_{1:n^{*}})\end{array}\right),\left(\begin{array}[]{cc}\mathbf{K}(\mathbf{x}_{1:n},\mathbf{x}_{1:n})&\mathbf{K}(\mathbf{x}_{1:n},\mathbf{x}^{*}_{1:n^{*}})\\\
\mathbf{K}(\mathbf{x}^{*}_{1:n^{*}},\mathbf{x}_{1:n})&\mathbf{K}(\mathbf{x}^{*}_{1:n^{*}},\mathbf{x}_{1:n^{*}})\end{array}\right)\right),$

After performing the Bayes’ rule and assuming the observations are not noisy
($\mathbf{y}_{1:n}=\mathbf{f}_{1:n}$) and
$\mathcal{D}=(\mathbf{x}_{1:n},\mathbf{y}_{1:n})$ being the observed data, we
get the posterior distribution for $\mathbf{f}^{*}$ as below:

$p(\mathbf{f}^{*}_{1:n^{*}}|\mathbf{x}^{*}_{1:n^{*}},\mathcal{D})=\mathcal{N}\left(\mathbf{f}^{*}_{1:n^{*}}\middle|\mathbf{m}_{n}(\mathbf{x}^{*}_{1:n^{*}}),\mathbf{K}_{n}(\mathbf{x}^{*}_{1:n^{*}},\mathbf{x}^{*}_{1:n^{*}})\right),$

Where the posterior mean function is:

$m_{n}(x)=m(x)+\mathbf{k}(x,\mathbf{x}_{1:n})\left(\mathbf{K}(\mathbf{x}_{1:n},\mathbf{x}_{1:n})\right)^{-1}\left(\mathbf{y}_{1:n}-\mathbf{m}(\mathbf{x}_{1:n})\right),$

And the posterior covariance function is:

$k_{n}(x,x^{\prime})=k(x,x^{\prime})-\mathbf{k}(x,\mathbf{x}_{1:n})\left(\mathbf{K}(\mathbf{x}_{1:n},\mathbf{x}_{1:n})\right)^{-1}\mathbf{k}^{T}(x,\mathbf{x}_{1:n}),$

With
$\mathbf{k}(x,\mathbf{x}_{1:n})=\left(k(x,\mathbf{x}_{1}),\dots,k(x,\mathbf{x}_{n})\right)$
being the cross-covariance vector.

Multi-fidelity modelling enables accurate inference of quantities of interest
by synergistically combining realizations of low-cost/low-fidelity models with
a small set of high-fidelity observations. This is particularly effective when
the low and high-fidelity models exhibit strong correlations and can lead to
significant computational gains over approaches that solely rely on high-
fidelity models Perdikaris  _et al._ (2017). The MFGP method introduced by
Perdikaris  _et al._ Perdikaris  _et al._ (2017) is fundamentally very similar
to that of the Kennedy and O’Hagan Kennedy and O’Hagan (2000), but it is also
capable of capturing the more complex nonlinear cross-correlations between the
data. Assuming that we have $s$ fidelity levels (which is 2 in our case),
based on this method, the GP has to be done $s$ times recursively. In other
words, the GP is done in its standard, usual way on the first level, and for
the next levels, the outputs of the previous level $t-1$ are the inputs for
the level $t$. This is done in the way described below Perdikaris  _et al._
(2017):

$f_{t}(\mathbf{x})=g_{t}(\mathbf{x},f_{*t-1}(\mathbf{x})),$

Where $f_{*t-1}(\mathbf{x})$ is the posterior distribution of the previous
level $t-1$ evaluated at the $\mathbf{x}$, being the inputs of the current
level $t$. The unknown function $g_{t}$ follows the description of GP:

$g_{t}\sim\mathcal{G}\mathcal{P}(f_{t}|0,k_{t}((\mathbf{x},f_{*t-1}(\mathbf{x})),(\mathbf{x}^{\prime},f_{*t-1}(\mathbf{x}^{\prime}));\theta_{t}))$

According to Le Gratiet and Garnier Le Gratiet and Garnier (2014), this scheme
has the same posterior distribution predicted by the fully coupled scheme of
Kennedy and O’Hagan Kennedy and O’Hagan (2000). This procedure implies two
conditions: (i) The training sets need to have a nested structure (i.e.
$\mathbf{x}_{t}\subseteq\mathbf{x}_{t-1}$) (ii) The Markov property:

$cov\left\\{f_{t}(\mathbf{x}),f_{t-1}(\mathbf{x}^{\prime})|f_{t-1}(\mathbf{x})\right\\}=0,\forall\mathbf{x}\neq\mathbf{x}^{\prime}$

which translates into assuming that given the nearest point
$f_{t-1}(\mathbf{x})$, we can learn nothing more about $f_{t}(\mathbf{x})$
from any other model output $f_{t-1}(\mathbf{x}^{\prime})$, for any
$\mathbf{x}\neq\mathbf{x}^{\prime}$.

Since $f_{*t-1}$ and $\mathbf{x}$ belong to inherently 2 different spaces, a
more structured kernel function, coupling the elements from the same space
together, is the following Perdikaris  _et al._ (2017):

$k_{t_{g}}=k_{t_{\rho}}(\mathbf{x},\mathbf{x}^{\prime};\theta_{t_{\rho}})\cdot
k_{t_{f}}(f_{*t-1}(\mathbf{x}),f_{*t-1}(\mathbf{x}^{\prime});\theta_{t_{f}})+k_{t_{\delta}}(\mathbf{x},\mathbf{x}^{\prime};\theta_{t_{\delta}}),$

where $k_{t_{\rho}}$, $k_{t_{f}}$, and $k_{t_{\delta}}$ are valid covariance
functions, and $\theta_{t_{\rho}}$, $\theta_{t_{f}}$, and
$\theta_{t_{\delta}}$ are their hyperparameters, respectively. For our
application, we have chosen the SE kernel function:

$k(x,x^{\prime};\theta_{t})=\sigma_{t}^{2}\exp\left\\{-\frac{1}{2}\sum_{i=1}^{d}\mathcal{W}_{i,t}(x_{i}-x_{i}^{\prime})^{2}\right\\},$

where $\sigma_{t}^{2}$ is a variance parameter, and
$\bigl{\\{}\mathcal{W}_{i,t}\bigr{\\}}_{i=1}^{d}$ are the Automatic Relevance
Determination (ARD) weights corresponding to the fidelity level $t$.

As previously discussed, the posterior distribution of the first level is
obtained by the normal GP. Therefore, it is Gaussian. However, this is not
generally the case for the subsequent levels, except for the case of $t=2$.
Therefore, for all the cases with $t\geq 2$, we have to perform predictions
given uncertain inputs, where the uncertainty is propagated along each
recursive step. Thus, the posterior distribution is given by Perdikaris  _et
al._ (2017):

$p(f_{*t}(\mathbf{x}_{*})):=p(f_{t}(\mathbf{x}_{*},f_{*t-1}(\mathbf{x}_{*}))|f_{*t-1},\mathbf{x}_{*},\mathbf{y}_{t},\mathbf{x}_{t})=\int
p(f_{t}(\mathbf{x}_{*},f_{*t-1}(\mathbf{x}_{*}))|\mathbf{y}_{t},\mathbf{x}_{t},\mathbf{x}_{*})p(f_{*t-1}(\mathbf{x}_{*}))d\mathbf{x}_{*}$

The predictive mean and variance of all posteriors
$p(f_{*t}(\mathbf{x}_{*})),t\geq 2$ is calculated using Monte Carlo
integration of this equation. We can then sample from this distribution.

##  III Results and Discussion

###  III.1 Oscillatory flow

Dynamics of the droplet in an oscillatory flow in the microchannel has been
studied, and the effects of $Re$, $Ca$, and $\omega^{*}$ have been
investigated. $Re$ ranges between $10$ and $100$, $Ca$ ranges between $0.09$
and $10$, and $\omega^{*}$ values are chosen such that for a channel with a
square cross-section of $100\mu m$ and water as the working fluid at room
temperature, the frequency values range between 2Hz and 1600Hz (or a
corresponding $Wo$ number between $0.2$ and $6.3$). Although a maximum
frequency of 200Hz is mostly reported in the literature Dincau  _et al._
(2020), recent works have claimed of generating frequencies of around 1KHz
Vishwanathan and Juarez (2020). The equilibrium position of the drop is a
result of the competition between the lateral lift forces acting on it,
including the wall effect and the deformation-induced lift forces, both acting
towards the channel center, and the shear gradient force acting towards the
wall Fay  _et al._ (2016); Zhang  _et al._ (2016); Rivero-Rodriguez and Scheid
(2018). Magnus and Saffman lift forces are often very small compared to the
other mentioned components and can be neglected Zhang  _et al._ (2016); Martel
and Toner (2014); Stoecklein and Di Carlo (2018). The boundary wall causes the
particle to have rotational and translational velocities different from those
of the adjacent fluid, which is caused by an uneven distribution of
vorticities around the particle Ho and Leal (1974); Bagchi and Balachandar
(2002). This induces a higher pressure in the gap between the particle and the
wall, which repels the particle away from the wall Zhang  _et al._ (2016). The
existing curvature of the fluid velocity profile makes the magnitude of the
velocity of the fluid on the wall side much higher than the channel center
side from the particle frame of reference. This inequality causes a low
pressure on the wall side leading to a shear gradient lift force that pushes
the particle towards the wall Zhang  _et al._ (2016). Following the analytical
results of Chan and Leal Chan and Leal (1979), the deformability-induced lift
force for droplets or bubbles that have a distance higher than their diameter
from the wall, which is the case in our simulations, is given by Stan  _et
al._ (2013):

$F_{L,deformation}=Ca_{p}\mu V_{avg}a(\frac{a}{W})^{2}\frac{d}{W}f(\lambda),$
(5)
$f(\lambda)=\frac{128\pi}{(\lambda+1)^{3}}\left\\{\frac{11\lambda+10}{140}(3\lambda^{2}-\lambda+8)-\frac{3(19\lambda+16)}{14(3\lambda+2)}(2\lambda^{2}+\lambda-1)\right\\},$
(6)

Where $Ca_{p}={\frac{\mu U_{0}}{\gamma}}{\frac{a}{W}}$ is the drop capillary
number, $V_{avg}$ is the average velocity of the carrier fluid across the
channel, $d$ is the distance of the drop from the channel center, and
$\lambda$ is the viscosity ratio between the inner and outer fluids.

The code has been validated by comparing the drop deformations at $Ca=0.2$ and
different Deborah numbers with those of the Aggarwal and Sarkar Aggarwal and
Sarkar (2007). The results are in good agreement with a maximum error of
$0.72\%$. To validate the inertial effects, the focal points of the drop at
$Re=8.25$, $Ca=0.18$, $\frac{a}{W}=0.2$ and $Re=21$, $Ca=0.14$,
$\frac{a}{W}=0.3$ have been compared with those presented by Marson  _et al._
Marson  _et al._ (2018). Our obtained focal points lie within their
corresponding uncertainty bands. Furthermore, we have shown that the numerical
results are independent of the distance between $2$ consecutive drops in an
infinite domain in the flow direction. This has been done by comparing the
drop trajectory at $Re=10$, $Ca=1$, and $\omega^{*}=0.1$ for three different
channel lengths of $4W$, $6W$, and $8W$ in our simulation setup. The maximum
difference between the drop trajectories for $L=4W$ and $L=6W$ is $0.0003W$,
and the one between those of $L=4W$ and $L=8W$ is $0.0005W$. The results have
also been shown to be grid independent, by comparing the equilibrium positions
for the case of $Re=10$, $Ca=1$, and $\omega^{*}=0.1$ with two different grids
of $196\times 114\times 114$ and $256\times 152\times 152$. The difference
between their focal distances from the center is $0.0009W$.

Figure 2 shows the distance of the droplet focal point from the channel center
($d^{*}$) at different values of $Wo$, $Ca$, and $Re$. The drop focal point in
the steady flow moves towards the center by increasing the $Ca$ due to the
increase in the deformation force Pan  _et al._ (2016); Lan and Khismatullin
(2012); Hadikhani  _et al._ (2018) and shifts towards the wall as $Re$
increases because of the improvement in the strength of the shear gradient
force Mortazavi and Tryggvason (2000); Hadikhani  _et al._ (2018); Di Carlo
_et al._ (2009). Since $W$ and $\nu$ are constant in this work, $\omega^{*}$
and $Wo$ are directly related to each other. Hence, we can use them
interchangeably. The droplet travels at locations far from the wall; so the
wall lift can be neglected in our study Zhang  _et al._ (2016); Ho and Leal
(1974). Both deformation and shear gradient lift forces, as the remaining
active forces, depend on $V_{avg}$ and $d^{*}$ that vary as the simulations
proceed. The dependence of deformation force on these 2 parameters is apparent
from equation 5. The parameter $d^{*}$ and the flow velocity at the drop
location affect the magnitude of the difference between the velocities on the
wall and center sides from the drop frame of reference. Hence, both $d^{*}$
and $V_{avg}$ determine the magnitude of the shear gradient force. For non-
steady flows, including oscillatory and pulsating ones, $V_{avg}$ is time-
dependent. The average of this $V_{avg}$ in each corresponding periodic cycle
decreases as the $\omega^{*}$ (or $Wo$) increases. Consequently, the averages
of both forces in one periodic cycle change by changing the $Wo$ value keeping
other parameters fixed, leading to different equilibrium positions as we can
see in Fig. 2. A complete explanation of the relationship between the focal
point and parameters like $\omega^{*}$, $Ca$, and $Re$ can be found in our
previous work Lafzi  _et al._ (2020). According to this figure, the focal
point is closest to the channel center at the highest $Wo$ except for $Re=10$
and $Ca=1$ and $Re=10$ and $Ca=1.67$. This can be explained based on the shape
of the flow velocity profile elaborated below.

Figure 2: Distance of the equilibrium position of the droplet from the channel
center as a function of $Wo$ for different cases

Stokes number in the oscillatory and pulsating flows is defined as
$St=W\sqrt{\frac{\omega}{\nu}}$. O’Brien O’Brien (1975) has solved these types
of flows in rectangular channels analytically and quantified the shape of the
velocity profile by calculating the ratio between the velocity at the center
and its average across the cross-section. Based on the values in table 1 of
this paper, the profile maintains its parabolic shape up to $St=3$ for a
square channel. Above this approximate value, the profile starts to become
more like a flat, plug-like profile O’Brien (1975). Figure 3 illustrates the
shapes of the dimensionless averaged velocity along the flow direction near
the drop focal points at $St=3$ ($Wo=3.15$) and $St=6$ ($Wo=6.3$). These
shapes are consistent with the findings of O’Brien (1975) and Karbaschi  _et
al._ (2014). At $Wo=3.15$ ($\omega^{*}=2$) the velocity shape is still
parabolic and is similar to those of the lower $Wo$ values. However, this
shape changes to plug like at $Wo=6.3$ ($\omega^{*}=8$). The average of
$V_{avg}$ at this $Wo$ is the lowest because it has the highest frequency
among others. Therefore, the value of deformation force on average is very
small according to equation 5. However, due to the shape of the velocity
profile, the relative flow velocity from the drop frame of reference is very
small near the center and very large near the wall. Thus, there is a strong
shear gradient force although the average of $V_{avg}$ is small. Consequently,
the focal point at $Wo=6.3$ does not follow the trend observed for the lower
$Wo$ numbers and is pushed towards the wall.

Furthermore, by taking a closer look at Fig. 3(c), we can see that the
velocity has an opposite sign near the wall. Due to the present hysteresis in
this type of flow, there is a lag in the response of fluid to the change of
flow direction Noguchi (2010); Chaudhury  _et al._ (2016).

(a)

(b)

(c)

Figure 3: Normalized flow direction averaged velocity profiles on the $xy$ and
$xz$ planes at $Re=10$, $Ca=1.67$, and (a) $\omega^{*}=2$ $(Wo=3.15,St=3)$,
(b) $\omega^{*}=8$ $(Wo=6.3,St=6)$, and (c) their shape comparison on the $xy$
plane with the same constant $z$

The Taylor deformation parameter of the drop is defined as:

$D=\frac{L-B}{L+B}$ (7)

In which, $L$ is the principal major axis, and $B$ is the principal minor axis
of an equivalent ellipsoidal particle. The parameter $D$ is oscillatory for
all the flows in this study except for the steady one. This parameter is
proportional to the shear rate, which depends on the flow velocity. As a
result, the average of $D$ is lower at higher frequencies Zhao and Bagchi
(2011). This is reflected in Fig. 4 by visualizing the average of $D$ in the
corresponding periodic cycle as a function of time. This trend also implies
that the amount of oscillations in the deformation is lower at higher
frequencies since the minimum deformation in each cycle is zero. Moreover, it
is apparent that the case with $\omega^{*}=8$ has a deformation of close to 0,
which confirms that it has a very low $V_{avg}$. In fact, the droplet in this
case travels about only $0.02W$ along the flow direction and remains almost
spherical though the $Ca=1$ corresponds to a very deformable drop.
Nevertheless, it migrates around $0.08W$ from the initial location to its
focal point.

Previous works have reported that drops in the flow regimes of high $Ca$ are
elongated significantly leading to their break up Lan and Khismatullin (2012);
Pan  _et al._ (2016); Marson  _et al._ (2018). This is the reason for the
absence of any data for $Ca=10$ and $Re=10$ for steady and lower frequency
flows in Fig. 2. The drop undergoes a very large deformation in these cases,
which is not what it experiences at higher frequencies. A similar argument
holds for $Re=37.8$ and $Ca=1.67$ since the drop is able to deform more easily
at higher $Re$ Inamuro  _et al._ (2003).

Figure 4: Averaged deformation parameter vs averaged dimensionless time at
$Re=10$, $Ca=1$, and different frequencies

As it was discussed previously, the velocity in the presented type of flows is
oscillatory. Consequently, all the active forces also experience fluctuations
leading to oscillations in the trajectory of the drop except for the steady
flows Chaudhury  _et al._ (2016). This can be seen in Fig. 5 where the whole
lateral migration patterns at $Re=10$, $Ca=1$, and all frequencies,
$\omega^{*}=1$, and $\omega^{*}=0.01$, showing the oscillations in the
trajectories are illustrated in figures 5(a), 5(b), and 5(c), respectively,
and the amplitudes of oscillations after focusing ($A^{*}$) for different
cases are depicted in Fig. 5(d). In these figures, $d^{*}$ is the
dimensionless, time-dependent distance of the drop from the channel center,
and the insets of figures 5(b) and 5(c) show the trajectory at the last $2$
periodic cycles. The plots depicted in the latter $2$ figures qualitatively
agree with those of a previous study and become more like a helical, spiral
pathway as the frequency increases Chaudhury  _et al._ (2016). The minimum
velocity in each periodic cycle is zero, occurring when the flow direction
changes. The higher average velocity at lower frequencies implies a higher
maximum velocity in the corresponding cycle. The higher the difference between
the maximum and minimum velocities, the higher is the oscillations amplitude
in the forces and in the trajectory. Therefore, similar to the aforementioned
discussion of deformation oscillations, the parameter $A^{*}$ decreases as
$\omega^{*}$ or $Wo$ increases (Fig. 5(d)). This is also the case for the
oscillations amplitude along the flow direction (the comparison between
figures 5(b) and 5(c)). The existence of these oscillations at all $Wo$ values
is noteworthy Chaudhury  _et al._ (2016). It is vital to mention that all the
values reported in Fig. 2 are the averages of $d^{*}$ in the last periodic
cycle.

(a)

(b)

(c)

(d)

Figure 5: (a) Migration patterns at $Re=10$, $Ca=1$, and various frequencies,
(b) Drop trajectory versus the flow direction at $Re=10$, $Ca=1$, and
$\omega^{*}=1$, and (c) $\omega^{*}=0.01$, and (d) Amplitude of oscillations
around the equilibrium point after focusing for different cases

Focusing time can be considered as an important factor in the design and
performance of the microfluidic system. Nevertheless, it strongly depends on
where the particle is initially released. Therefore, the average migration
velocity is a better parameter for a more meaningful comparison under
different circumstances. The average migration velocity ($v^{*}$) for
different cases is shown in Fig. 6. This velocity is computed by calculating
the distance between the initial and equilibrium positions and dividing it by
the focusing time. The corresponding focusing time is determined when the
trajectory reaches within $0.015W$ of the focal point. This figure expresses
that the average migration velocity decreases as $Wo$ increase. This pattern
is also observed in the average velocity along the flow direction. Besides,
the average migration velocity in the steady flows increases by increasing the
$Ca$, which is in agreement with the findings of Alghalibi  _et al._ Alghalibi
_et al._ (2019).

Figure 6: The average migration velocity vs $Wo$

Viscosity ratio in the range of $1\leq\lambda<13.2$ has a weak effect on the
drop migration Marson  _et al._ (2018); Pan  _et al._ (2016). The effect of
density ratio ($\eta$) on the drop focal point is even less Mortazavi and
Tryggvason (2000). Therefore, we have limited our study on the effect of these
two parameters on the drop migration only to one case as shown in Fig. 7. The
drop with a $\lambda$ higher than $1$ focuses closer to the wall for the $Ca$
we are studying here Marson  _et al._ (2018); Pan  _et al._ (2016); Mortazavi
and Tryggvason (2000). We can also see that the change in the $\eta$ and
$\lambda$ does not change the distance between the focal points in the steady
and oscillatory flows significantly.

Figure 7: The effect of density and viscosity ratios on the distance of the
drop focal point from the channel center at $Re=10$ and $Ca=1.67$

###  III.2 Pulsating flow

Changing the steady flow to oscillatory type helps achieve different focal
points, which can have potential applications in cell sorting and separation
as a great advantage. However, the oscillatory flow has a zero net throughput,
which can be counted as a disadvantage considering many microfluidic
applications depending on high-throughput systems. One solution to fix this
issue is to combine both steady and oscillatory parts in the pressure gradient
and make the flow regime to be pulsating Dincau  _et al._ (2020). A pulsating
flow has the advantages of having non-zero net throughput as well as
introducing a new equilibrium point for the drop. The latter occurs because
the drop dynamics in the pulsating flow is very similar to the one in a pure
oscillatory flow with an equivalent frequency between $0$ and that of the
oscillatory portion of the pulsating pressure gradient. This claim is valid
since the weights of steady and oscillatory portions add up to 1. We can see
this equivalent frequency in Fig. 8(a) and the inset of Fig. 8(c) where the
trajectories in the pulsating flows have frequencies of almost half of those
of the oscillatory ones with the same frequencies. This feature also enables
the existence of cases with high $Ca$ or $Re$ and low $\omega^{*}$ that are
absent in the figures of the previous section. This is because the equivalent
frequency of the pulsating flow makes the average deformation lower compared
to the pure steady flow or oscillatory flow with a lower frequency. It is
crucial to mention that the highest weights for the steady portions of the
pulsating cases depicted in Fig. 8 denote flows in which the drop experiences
the highest feasible deformation without breaking up or being significantly
elongated. This has been done to show the highest possible changes in the
focal points.

In addition, the directionality of the focal points and averages of
deformation values obey the expected trend at each combination of $Re$ and
$Ca$ in Fig. 8. For instance, the pulsating flow at $Re=37.8$ and $Ca=1.67$
has an equivalent frequency less than $0.1$. Therefore, we can see a focal
point closer to the channel center and a higher average deformation according
to Fig. 8(a) and Fig. 8(b), respectively. Similarly, we can observe a focal
point closer to the center at $Re=10$ and $Ca=10$ for the pulsating case with
the larger steady portion. This is because this case has an equivalent
frequency between $0.1$ and $0.5$ and less than the one with the lower steady
portion. Fig. 8(c) reflects this pattern. The deformation behavior in Fig.
8(d) is also qualitatively similar to that shown in Fig. 8(b). It is also
momentous to pay attention to the $0.012W$ difference between the focal points
of the pulsating cases at $Re=10$ and $Ca=10$ although the difference between
the weights of their steady portions is only $0.01$.

(a)

(b)

(c)

(d)

Figure 8: Evolution of the oscillatory and pulsating dynamics of the drop at
$Re=37.8$ and $Ca=1.67$ by illustrating its (a) trajectory, (b) deformation,
and at $Re=10$ and $Ca=10$ by visualizing the corresponding (c) trajectory,
and (d) deformation

Tables 2 and 2 quantify the average migration velocities for the flow regimes
discussed above. The observed trend in these tables is consistent with the
information provided in the discussion of Fig. 6. The pulsating case at
$Re=37.8$ and $Ca=1.67$ has the lowest equivalent frequency. Hence, it has the
highest average migration velocity among others. Similarly, the $v^{*}$ at
$Re=10$ and $Ca=10$ decreases as the equivalent frequency increases.

Pressure gradient form |  Average migration velocity
---|---
$P_{0}(0.93+0.07cos(0.1t))$ |  0.00191
$P_{0}cos(0.1t)$ |  0.00141
$P_{0}cos(0.5t)$ |  0.00012
Table 1: Average migration velocities for different flow regimes at $Re=37.8$
and $Ca=1.67$

Pressure gradient form |  Average migration velocity
---|---
$P_{0}cos(0.1t)$ |  0.000389
$P_{0}(0.17+0.83cos(0.5t))$ |  0.000297
$P_{0}(0.07+0.93cos(0.5t))$ |  0.000213
$P_{0}cos(0.5t)$ |  0.000205
Table 2: Average migration velocities for different flow regimes at $Re=10$
and $Ca=10$

The effect of droplet size on its focal point in the oscillatory and pulsating
flows at $Re=10$ and $Ca=10$ is summarized in Table 3. We observe that
reducing the drop size pushes its equilibrium location towards the wall Pan
_et al._ (2016); Mortazavi and Tryggvason (2000); Di Carlo  _et al._ (2009);
Wang  _et al._ (2017); Bazaz  _et al._ (2020). Furthermore, this size
reduction appears to enhance the change in the focal point at different values
of equivalent frequency.

Pressure gradient form | Equilibrium distance from center
---|---
|  $\mathbf{\frac{a}{W}=0.3}$ |  $\mathbf{\frac{a}{W}=0.2}$
$P_{0}cos(0.1t)$ |  0.138 |  0.193
$P_{0}(0.17+0.83cos(0.5t))$ |  0.166 |  0.241
$P_{0}cos(0.5t)$ |  0.183 |  0.266

Table 3: The effect of drop size on its focal point at $Re=10$ and $Ca=10$

###  III.3 MFGP

In this section, we evaluate the MFGP performance by having a dataset
consisting of 29 low-fidelity and 22 high-fidelity observations. Fig. 9 shows
all these observations together. It can be seen that the required nested
structure, as mentioned earlier, is satisfied. In other words, for any data
point in the high-fidelity level, there is a corresponding point in the low-
fidelity level.

Figure 9: Observations of the high and low fidelity levels

Fig. 10 illustrates the predictions on the distance of the drop equilibrium
position from the channel center over the range of 0 to 1 for the frequency
and 0.09 to 1.67 for the Capillary number. The algorithm is trained over the
entire available data. A few sanity checks have been done to ensure that these
results make sense. First, Fig. 10(a) denotes that at any fixed value of the
Capillary number, there is a global optimum frequency for which the
equilibrium distance from the center (the $z$ axis) is an extremum. This is
also previously observed in the simulations for the input values in the
dataset (Fig. 2). Secondly, the predictions made by the high-fidelity response
are slightly higher than those of the low-fidelity response, which is also
compatible with the simulations outcomes. Lastly, we know from the underlying
physics that in the steady flow (at a frequency of 0), the distance decreases
by increasing the Capillary number Pan  _et al._ (2016); Lan and Khismatullin
(2012). This can be seen in both contour plots of the mean low and high
responses (figures 10(b) and 10(c)) as well as being confirmed quantitatively
in the code. Fig. 10(d) illustrates a relatively low variance for the high-
fidelity response, being our main goal, in the entire domain. The red-colored
region in this plot corresponds to places where the density of data points is
less, and hence, the predicted outputs have more uncertainty.

(a)

(b)

(c)

(d)

Figure 10: Prediction of the MFGP algorithm by illustrating the (a) surface
plot of the high-fidelity response colored by red and low-fidelity response
colored by green, (b) mean response of the low-fidelity level, (c) mean
response of the high-fidelity level, and (d) variance of the high-fidelity
response

Since there is no analytical solution to compare the predictions with, we
decided to split the high-fidelity data into training and test sets in this
section. This helps us evaluate the performance of the implemented MFGP. We
assign 15 training and 7 test points at the high-fidelity level randomly.
Then, we train the MFGP on the whole low-fidelity data (29 points) and only
the training high-fidelity data, and evaluate the predictive responses on the
high-fidelity test points, since those are our main targets. We do this entire
procedure 500 times to eliminate the dependence of the results on the test
points choice. This especially helps us examine the algorithm performance at
the regions with less amount of data. After this bootstrapping, the average of
the mean squared error (MSE) was 0.00015 and the average of the $R^{2}$ score
was 0.9858. The successful reproduction of these results is also checked. Fig.
11 shows this evaluation at the last (500th) test set. Fig. 11(a) expresses
the predictions of $d^{*}$ along with their 95% credible intervals and denotes
that the observed data lies within the shaded uncertainty bands. In this
figure, the x axis label (i) denotes the index of each test point. Figure
11(b) visualizes the true or known value of $d^{*}$ versus its predicted value
at each test point. This figure illustrates the strength of the model as the
plotted points are very close to the line of $y=x$, which is an indication of
the agreement between predictions and observations. The exact and predicted
correlations between the high and low responses of $d^{*}$ are also very close
to each other according to Fig. 11(c). All of these are evidence for the
strong and successful performance of the implemented MFGP algorithm.

(a)

(b)

(c)

Figure 11: Evaluation of the MFGP performance on a random test set by
illustrating the (a) observations and predictions with their uncertainties,
(b) predictions versus observations, and (c) exact versus predicted
correlations

##  IV Conclusions

Determination and control of the particle’s equilibrium position in the
microchannels are extremely crucial as it can help in a variety of
microfluidics applications. This importance, as well as the need to overcome
the issue of designing impractically long channels to work with sub-micron
particles, led us to do some simulations to capture the dynamics of a single
droplet suspended in an oscillatory flow within the channel. The drawback of
the zero net throughputs of the oscillatory flow has been addressed by
modifying it further to become a pulsating flow. Both types of flows bring new
equilibrium locations to the system. They also enable the presence of droplets
at high $Ca$ or $Re$ that could break up in the steady or a very low-frequency
regime. Moreover, fluctuations in the trajectory of the drop have been
observed. It has been shown that the amplitude of these oscillations, the
average of the oscillatory deformation, and the average migration velocity all
decrease by increasing the frequency. The dependence of the drop focal point
on the shape of the velocity profile has been investigated as well. It has
been explored that this equilibrium position moves towards the wall in a plug-
like profile, which is the case in very high $Wo$ numbers. Due to the
significant cost of these simulations, a recursive version of the Multi
Fidelity Gaussian processes has been used to replace the numerous high-
fidelity simulations that cannot be afforded numerically. The MFGP algorithm
is used to predict the equilibrium distance of the drop from the channel
center for a given range of the interplaying input parameters, namely the
Capillary number and frequency, assuming a constant Reynolds number. In
addition, its performance was evaluated by randomly shuffling the high-
fidelity data 500 times and assigning 31.8% of it as the test set for an
accurate quantitative comparison each time. The algorithm outputs high
statistical scores, which is an indication of its reasonably accurate
performance.

## References

  * Gossett  _et al._ (2010) D. R. Gossett, W. M. Weaver, A. J. Mach, S. C. Hur, H. T. K. Tse, W. Lee, H. Amini, and D. Di Carlo, Analytical and bioanalytical chemistry 397, 3249 (2010).
  * Toner and Irimia (2005) M. Toner and D. Irimia, Annu. Rev. Biomed. Eng. 7, 77 (2005).
  * Gossett  _et al._ (2012) D. R. Gossett, H. T. K. Tse, J. S. Dudani, K. Goda, T. A. Woods, S. W. Graves, and D. Di Carlo, Small 8, 2757 (2012).
  * Karimi  _et al._ (2013) A. Karimi, S. Yazdi, and A. Ardekani, Biomicrofluidics 7, 021501 (2013).
  * Di Carlo  _et al._ (2007) D. Di Carlo, D. Irimia, R. G. Tompkins, and M. Toner, Proceedings of the National Academy of Sciences 104, 18892 (2007).
  * Marson  _et al._ (2018) R. L. Marson, Y. Huang, M. Huang, T. Fu, and R. G. Larson, Soft matter 14, 2267 (2018).
  * Mutlu  _et al._ (2018) B. R. Mutlu, J. F. Edd, and M. Toner, Proceedings of the National Academy of Sciences 115, 7682 (2018).
  * Lan and Khismatullin (2012) H. Lan and D. B. Khismatullin, International journal of multiphase flow 47, 73 (2012).
  * Lafzi  _et al._ (2020) A. Lafzi, A. H. Raffiee, and S. Dabiri, Physical Review E 102, 063110 (2020).
  * Erlandsson and Robinson (2011) P. G. Erlandsson and N. D. Robinson, Electrophoresis 32, 784 (2011).
  * Chaudhury  _et al._ (2016) K. Chaudhury, S. Mandal, and S. Chakraborty, Physical Review E 93, 023106 (2016).
  * Pawłowska  _et al._ (2017) S. Pawłowska, P. Nakielski, F. Pierini, I. K. Piechocka, K. Zembrzycki, and T. A. Kowalewski, PloS one 12 (2017).
  * Pawłowska (2018) S. Pawłowska, IPPT Reports on Fundamental Technological Research , 1 (2018).
  * Sarkar and Schowalter (2000) K. Sarkar and W. R. Schowalter, Journal of non-newtonian fluid mechanics 95, 315 (2000).
  * Perdikaris  _et al._ (2017) P. Perdikaris, M. Raissi, A. Damianou, N. Lawrence, and G. E. Karniadakis, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 473, 20160751 (2017).
  * Babaee  _et al._ (2016) H. Babaee, P. Perdikaris, C. Chryssostomidis, and G. Karniadakis, Journal of Fluid Mechanics 809, 895 (2016).
  * Unverdi and Tryggvason (1992) S. O. Unverdi and G. Tryggvason, Journal of computational physics 100, 25 (1992).
  * Dincau  _et al._ (2020) B. Dincau, E. Dressaire, and A. Sauret, Small 16, 1904032 (2020).
  * Pan  _et al._ (2016) D.-y. Pan, Y.-q. Lin, L.-x. Zhang, and X.-m. Shao, Journal of Hydrodynamics 28, 702 (2016).
  * Razi and Pourghasemi (2017) M. Razi and M. Pourghasemi, Computers & Fluids 154, 200 (2017).
  * Kennedy and O’Hagan (2000) M. C. Kennedy and A. O’Hagan, Biometrika 87, 1 (2000).
  * Le Gratiet and Garnier (2014) L. Le Gratiet and J. Garnier, International Journal for Uncertainty Quantification 4 (2014).
  * Vishwanathan and Juarez (2020) G. Vishwanathan and G. Juarez, Microfluidics and Nanofluidics 24, 1 (2020).
  * Fay  _et al._ (2016) M. E. Fay, D. R. Myers, A. Kumar, C. T. Turbyfield, R. Byler, K. Crawford, R. G. Mannino, A. Laohapant, E. A. Tyburski, Y. Sakurai, _et al._ , Proceedings of the National Academy of Sciences 113, 1987 (2016).
  * Zhang  _et al._ (2016) J. Zhang, S. Yan, D. Yuan, G. Alici, N.-T. Nguyen, M. E. Warkiani, and W. Li, Lab on a Chip 16, 10 (2016).
  * Rivero-Rodriguez and Scheid (2018) J. Rivero-Rodriguez and B. Scheid, Journal of Fluid Mechanics 855, 1242 (2018).
  * Martel and Toner (2014) J. M. Martel and M. Toner, Annual review of biomedical engineering 16, 371 (2014).
  * Stoecklein and Di Carlo (2018) D. Stoecklein and D. Di Carlo, Analytical chemistry 91, 296 (2018).
  * Ho and Leal (1974) B. Ho and L. Leal, Journal of fluid mechanics 65, 365 (1974).
  * Bagchi and Balachandar (2002) P. Bagchi and S. Balachandar, Journal of Fluid Mechanics 473, 379 (2002).
  * Chan and Leal (1979) P.-H. Chan and L. Leal, Journal of Fluid Mechanics 92, 131 (1979).
  * Stan  _et al._ (2013) C. A. Stan, A. K. Ellerbee, L. Guglielmini, H. A. Stone, and G. M. Whitesides, Lab on a Chip 13, 365 (2013).
  * Aggarwal and Sarkar (2007) N. Aggarwal and K. Sarkar, Journal of Fluid Mechanics 584, 1 (2007).
  * Hadikhani  _et al._ (2018) P. Hadikhani, S. M. H. Hashemi, G. Balestra, L. Zhu, M. A. Modestino, F. Gallaire, and D. Psaltis, Lab on a Chip 18, 1035 (2018).
  * Mortazavi and Tryggvason (2000) S. Mortazavi and G. Tryggvason, Journal of Fluid Mechanics 411, 325 (2000).
  * Di Carlo  _et al._ (2009) D. Di Carlo, J. F. Edd, K. J. Humphry, H. A. Stone, and M. Toner, Physical review letters 102, 094503 (2009).
  * O’Brien (1975) V. O’Brien, Journal of the Franklin Institute 300, 225 (1975).
  * Karbaschi  _et al._ (2014) M. Karbaschi, A. Javadi, D. Bastani, and R. Miller, Colloids and Surfaces A: Physicochemical and Engineering Aspects 460, 355 (2014).
  * Noguchi (2010) H. Noguchi, Physical Review E 81, 061920 (2010).
  * Zhao and Bagchi (2011) M. Zhao and P. Bagchi, Physics of Fluids 23, 111901 (2011).
  * Inamuro  _et al._ (2003) T. Inamuro, R. Tomita, and F. Ogino, International Journal of Modern Physics B 17, 21 (2003).
  * Alghalibi  _et al._ (2019) D. Alghalibi, M. E. Rosti, and L. Brandt, Physical Review Fluids 4, 104201 (2019).
  * Wang  _et al._ (2017) Q. Wang, D. Yuan, and W. Li, Micromachines 8, 197 (2017).
  * Bazaz  _et al._ (2020) S. R. Bazaz, A. Mashhadian, A. Ehsani, S. C. Saha, T. Krüger, and M. E. Warkiani, Lab on a Chip 20, 1023 (2020).