# A reduced model for droplet dynamics in shear flows at finite capillary
numbers

Diego Taglienti
Department of Physics & INFN, University of Rome “Tor Vergata”,
Via della Ricerca Scientifica 1, 00133, Rome, Italy.
<EMAIL_ADDRESS>
&Fabio Guglietta
Helmholtz Institute Erlangen-Nürnberg for Renewable Energy (IEK-11),
Forschungszentrum Jülich GmbH, Cauerstraße 1, 91058 Erlangen, Germany
&Mauro Sbragaglia
Department of Physics & INFN, University of Rome “Tor Vergata”,
Via della Ricerca Scientifica 1, 00133, Rome, Italy.

###### Abstract

We propose an extension of the Maffettone-Minale (MM) model to predict the
dynamics of an ellipsoidal droplet in a shear flow. The parameters of the MM
model are traditionally retrieved in the framework of the perturbation theory
for small deformations, i.e., small capillary numbers ($\text{Ca}\ll 1$)
applied to Stokes equations. In this work, we take a novel route, in that we
determine the model parameters at finite capillary numbers
($\text{Ca}\sim{\cal O}(1)$) without relying on perturbation theory results,
while retaining a realistic representation in creeping time and steady
deformation attained by the droplet for different realizations of the
viscosity ratio $\lambda$ between the inner and the outer fluids. This
extended MM (EMM) model hinges on an independent characterization of the
process of droplet deformation via numerical simulations of Stokes equations
employing the Immersed Boundary - Lattice Boltzmann (IB-LB) numerical
techniques. Issues on droplet breakup are also addressed and discussed within
the EMM model.

##  1 Introduction

Understanding the deformation of droplets induced by hydrodynamic flows has
been a central topic in rheology for almost a century. In 1934, G. I. Taylor
faced this problem by investigating the deformation of a pure droplet using
the approximation of small deformations [1, 2]: in his study, Einstein’s work
for suspensions of solid spheres [3] was extended to the case of droplets with
radius $R$ and viscosity $\lambda\mu$ deformed by the action of an external
flow with shear intensity $\dot{\gamma}$ and viscosity $\mu$ (see Fig. 1).
Droplet deformation is quantified via the dimensionless capillary number
$\text{Ca}=\mu R\dot{\gamma}/\sigma$, weighting the importance of the external
viscous forces against the surface tension $\sigma$. Taylor’s work, whose
results are valid in the perturbative limit of $\text{Ca}\ll 1$, proved to be
a fundamental approach to the problem, paving the way for several other
studies which refine and extend the investigation in many directions, e.g. by
considering higher orders in the perturbation theory [4, 5], analyzing the
time-dependent deformation properties [6, 7], introducing the effects of
different types of flows [8, 9], plugging additional complexities in the
suspension system [10, 11]. Various reviews have been written on the topic,
covering experimental, theoretical and numerical aspects [12, 13, 14, 15].
Droplet dynamics is a very complex multiscale problem and, in many instances,
it is desirable to possess a simpler phenomenological description via reduced
models. Maffettone & Minale (MM) [16] proposed an equation for the evolution
of a droplet which is assumed to be ellipsoidal at all times. The droplet
shape is described via a second order tensor $\bm{\bm{S}}$ evolving with the
equation

$\displaystyle\frac{d\bm{\bm{S}}}{dt^{\prime}}-\text{Ca}\,\left(\bm{\Omega}\cdot\bm{S}-\bm{S}\cdot\bm{\Omega}\right)=-f_{1}^{(\mbox{MM})}(\lambda)\left[\bm{S}-g\left(\bm{S}\right)\bm{I}\right]+\text{Ca}\,f_{2}^{(\mbox{MM})}(\lambda)\left(\bm{E}\cdot\bm{\bm{S}}+\bm{\bm{S}}\cdot\bm{E}\right)\
,$ (1)

where $\bm{E}$ and $\bm{\Omega}$ represent the symmetric and asymmetric
velocity gradient tensors respectively, $\bm{I}$ the unity tensor and
$g\left(\bm{\bm{S}}\right)$ a nonlinear function of the $\bm{\bm{S}}$ tensor
[16]. Time $t^{\prime}={t}/{\tau}$ is made dimensionless via the droplet
characteristic time $\tau=\mu R/\sigma$. The model parameter
$f_{1}^{(\mbox{MM})}(\lambda)$ can be linked to the creeping time of the
droplet $\tau_{\mbox{creep}}$, i.e., the time the droplet takes to reach a
stationary state, while $f_{2}^{(\mbox{MM})}(\lambda)$ is linked to the
deformation of the droplet induced by the flow [16]. The dependence of these
two model parameters on the viscosity ratio $\lambda$ is fixed by requiring
that in the limit of small Ca the model recovers the known results of
perturbation theory [17], in agreement with Taylor’s theory of small
deformation [1, 2]. The MM model gained great success thanks to its ability to
reproduce experimental data [16, 18, 19, 20] and it has been extended in many
directions, e.g. to account for non-Newtonian effects based on experimental
and theoretical inputs [21, 22, 10, 23], to investigate the deformation of a
drop in confined flows [24, 25] or to account for the tank-treading motion of
the membrane in a soft suspension [26]. Based on experimental observations of
droplet deformation up to subcritical regimes [27, 20, 28], some models were
also proposed to predict the non-linear response at high values of capillary
number [28, 29, 30]. An extensive review on the variety of reduced models for
droplet dynamics can be found in [31].
When designing a reduced model for droplet dynamics, it is a standard practice
to determine the model coefficients by requiring a theoretical matching with
the asymptotic perturbative expansion in the capillary number Ca coming from
Stokes equations. In this paper we take a different route, developing an
alternative strategy to construct the model parameters via direct matching
with numerical simulations of Stokes equations at $\text{Ca}\sim{\cal O}(1)$.
Specifically, we will study the process of droplet deformation subject to a
simple shear flow and construct an extended MM (EMM) model with parameters
$f_{1}^{(\mbox{EMM})}(\lambda,\text{Ca})$ and
$f_{2}^{(\mbox{EMM})}(\lambda,\text{Ca})$ capable of reproducing creeping
times and non-linear rheology at $\text{Ca}\sim{\cal O}(1)$. Thus, we retain
the attractiveness and ease of a single second-order tensorial model, plugging
a functional dependency on Ca in the model coefficients. The EMM model is then
studied at changing $\lambda$ up to the subcritical regime [27, 20, 32]. The
numerical solutions of Stokes equations are obtained via numerical simulations
employing the Immersed Boundary - Lattice Boltzmann (IB-LB) method [33, 34,
35, 36], a hybrid numerical method widely employed for simulating the dynamics
of deformable soft interfaces [37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48,
49, 50, 51]. Our study is also instrumental to further strengthen the
computational applicability of this numerical technique in situations aimed at
capturing the features of droplet deformations in subcritical regimes.
The paper is organized as follows: in Sec. 2 we review the basic ingredients
of the MM model and summarize relevant model predictions that will be relevant
for our analysis; in Sec. 3 we describe the IB-LB method used in this work to
simulate the dynamics of a single droplet in simple shear flow; results and
discussions are presented in Sec. 4; conclusions and a summary of our results
will be presented in Sec. 5.

Figure 1: Shear plane section of a droplet with viscosity $\lambda\mu$
deforming under the effect of a simple shear flow with intensity
$\dot{\gamma}$ applied in an outer fluid with visocisty $\mu$. The droplet has
an initial spherical shape with radius $R$ and deforms into a non-spherical
shape due to the action of the shear flow. The shape of the droplet is
described via a second-order tensor $\bm{\bm{S}}(t)$. The major ($L(t)$) and
minor ($B(t)$) axes in the shear plane are used to compute the time dependent
Taylor index $D(t)$, while $\theta(t)$ measures the inclination angle with
respect to the streamflow direction ($x-$axis).

##  2 Maffettone-Minale (MM) Model and proposed extension (EMM)

The Maffettone-Minale (MM) model [16] is structured by assuming that the
droplet retains an ellipsoidal shape at all times and the second-order tensor
$\bm{\bm{S}}$ used to describe its shape is assumed to evolve due to two main
competing effects: the drag exerted by the inner and outer fluid and the
surface tension tendency to restore the spherical geometry; while being
deformed, the droplet is assumed to further conserve the volume. The
dimensionless equation describing the time evolution of $\bm{\bm{S}}$ proposed
by MM follows as:

$\displaystyle\frac{d\bm{\bm{S}}}{dt^{\prime}}-\text{Ca}\,\left(\bm{\Omega}\cdot\bm{S}-\bm{S}\cdot\bm{\Omega}\right)=-f_{1}\left[\bm{S}-g\left(\bm{S}\right)\bm{I}\right]+\text{Ca}\,f_{2}\left(\bm{E}\cdot\bm{\bm{S}}+\bm{\bm{S}}\cdot\bm{E}\right)\
,$ (2)

where
$\bm{E}=\frac{1}{2}\left(\bm{\nabla}\bm{\bm{u}}+\bm{\nabla}\bm{\bm{u}}^{T}\right)$
and
$\bm{\Omega}=\frac{1}{2}\left(\bm{\nabla}\bm{\bm{u}}-\bm{\nabla}\bm{\bm{u}}^{T}\right)$
are the symmetric and asymmetric parts of the velocity gradient
$\bm{\nabla}\bm{\bm{u}}$, the LHS contains the Jaumann derivative,
$g\left(\bm{\bm{S}}\right)=\frac{3III_{S}}{II_{S}}$ is an additional nonlinear
function plugged in to satisfy the constraint of volume conservation, with:

$II_{S}=\frac{1}{2}\left[\mbox{tr}\left(\bm{S}\right)^{2}-\mbox{tr}\left(\bm{S}^{2}\right)\right]\
,\qquad III_{S}=\det(\bm{S})\ ,$ (3)

and Ca is the capillary number:

$\text{Ca}=\frac{\mu R\dot{\gamma}}{\sigma}\ ,$ (4)

where $\mu$ is the viscosity of the fluid, $R$ is the radius of the droplet at
rest and $\sigma$ is the surface tension.
In order to evaluate the coefficients $f_{1}$ and $f_{2}$, MM expand the shape
of the drop $\bm{\bm{S}}$ in the capillary number Ca as follows:

$\bm{\bm{S}}=\bm{I}+\text{Ca}\bm{\bm{F}}\ ,$ (5)

with $\bm{\bm{F}}$ capturing the deformations “beyond sphericity" whose
magnitude is expressed by the capillary number Ca, and they substitute it in
Eq.(2), obtaining:

$\frac{d\bm{\bm{F}}}{dt}=-f_{1}\bm{\bm{F}}+2f_{2}\bm{\bm{E}}\ .$ (6)

They finally compare the latter with the expansion to the first order in Ca
(see refs.[16, 17] and references therein for further details on the
perturbative approaches at the problem), finding:

$\displaystyle
f_{1}=f_{1}^{(\mbox{MM})}(\lambda)=\frac{40\left(\lambda+1\right)}{\left(2\lambda+3\right)\left(19\lambda+16\right)}\
,$ $\displaystyle f_{2}=f_{2}^{(\mbox{MM})}(\lambda)=\frac{5}{2\lambda+3}\ .$
(7)

The information on the deformation can be retrieved from the eigenvectors of
$\bm{\bm{S}}$, whose eigenvalues correspond to the square of the ellipsoid
semiaxis $L^{2},B^{2},W^{2}$. The major ($L$) and minor ($B$) axes in the
shear plane are used to compute the Taylor index $D=\frac{L-B}{L+B}$ as well
as the inclination angle $\theta$ (see Fig. 1). We are interested at the
deformation $D$ and the inclination angle $\theta$ in simple shear flow, whose
matrices $\bm{\bm{E}}$ and $\bm{\bm{\Omega}}$ are given by:

$\bm{E}=\frac{1}{2}\begin{pmatrix}0&1&0\\\ 1&0&0\\\ 0&0&0\\\
\end{pmatrix}\qquad\bm{\Omega}=\frac{1}{2}\begin{pmatrix}0&1&0\\\ -1&0&0\\\
0&0&0\end{pmatrix}\ .$ (8)

By substituting them in Eq. (2) and imposing
$\frac{d\bm{\bm{S}}}{dt}=\bm{\bm{0}}$, the model provides the following
analytical formulae for both $D$ and $\theta$ at the stead-state [16]:

$D=\frac{\sqrt{f_{1}^{2}+\text{Ca}^{2}}-\sqrt{f_{1}^{2}+\text{Ca}^{2}-f_{2}^{2}\text{Ca}^{2}}}{f_{2}\text{Ca}}\
,$ (9) $\theta=\frac{1}{2}\arctan\left(\frac{f_{1}}{\text{Ca}}\right).$ (10)

The model parameter $f_{1}$ physically represents the inverse creeping time
$\tau_{\mbox{{creep}}}$, i.e., the time the droplet takes to reach the
stationary state. As remarked above, however, this coefficient is fixed by
requiring a theoretical matching with the perturbative results coming from
Stokes equation, hence the validity of Eq. (10) remains restricted to the
small Ca regime. Nonetheless, $f_{1}^{-1}$ is a measurable time, e.g. it can
be measured from a direct analysis of the time evolution of the Taylor index
$D$. One could therefore think of modifying $f_{1}$ by directly measuring the
creeping time via dedicated numerical simulations; one can then use the steady
state results from droplet deformation to modify $f_{2}$ in such a way that
the steady deformation predicted by the MM model in Eq. (9) matches the one
that is measured in numerical simulations. In this way, one can define an
“Extended Maffettone-Minale" model (hereafter denoted with EMM) with model
coefficients $f_{1}^{(\mbox{EMM})}$ and $f_{2}^{(\mbox{EMM})}$ that naturally
acquire a dependency on Ca

$f_{1}=f_{1}^{(\mbox{EMM})}(\lambda,\text{Ca})\ ,\hskip
14.45377ptf_{2}=f_{2}^{(\mbox{EMM})}(\lambda,\text{Ca})\ ,$ (11)

and will reduce to the original MM coefficients when $\text{Ca}\rightarrow 0$

$\lim_{\text{Ca}\rightarrow
0}f_{1,2}^{(\mbox{EMM})}(\lambda,\text{Ca})=f^{(\mbox{MM})}_{1,2}(\lambda).$
(12)

The resulting EMM model will then possess (by construction) more realistic
creeping times and a more realistic steady deformation at finite Ca. This is
the strategy that we will pursue to extend the MM model.

##  3 Numerical Simulations: Immersed Boundary - Lattice Boltzmann (IB-LB)
Method

Numerical simulations are needed for the characterization of the droplet
deformation process under shear flow. In general, various numerical methods
have been developed to simulate bulk viscous flows in presence of deformable
soft suspensions, like boundary element methods [52, 53, 32, 54], volume-of-
fluid methods [55] and Immersed Boundary-Lattice Boltzmann (IB-LB) methods
[36]. In this work we use the IB-LB method. This method has already been used
in previous works for investigating the dynamics of droplets and viscoelastic
capsules [37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 56, 47, 48, 49, 50, 51],
hence we only recall the relevant features of the method (the interested
reader can found more details in the aforementioned works and references
therein).
The IB-LB method hinges on the Immersed Boundary (IB) model that couples the
interface of the droplet to the fluid, and on the lattice Boltzmann method
(LB) that simulates the bulk viscous flows. The LB method is a kinetic
approach to simulate hydrodynamics [36, 57]: we will briefly illustrate its
working principles taking into consideration the bulk fluid in the outer
droplet region, with macroscopic density $\rho$, velocity $\bm{\bm{u}}$ and
viscosity $\mu$ (the same reasoning applies to the bulk fluid inside the
droplet). The target equations for the LB method are the continuity equation
and the Navier-Stokes equations

$\frac{\partial\rho}{\partial t}+\bm{\nabla}\cdot(\rho\bm{\bm{u}})=0\ ,$ (13)
$\rho\left(\frac{\partial\bm{\bm{u}}}{\partial
t}+(\bm{\bm{u}}\cdot\boldsymbol{\nabla})\bm{\bm{u}}\right)=-\boldsymbol{\nabla}p+\mu\boldsymbol{\nabla}^{2}\bm{\bm{u}}+\bm{\bm{F}}_{{\mbox{ext}}}\
,$ (14)

where $\bm{\bm{F}}_{{\mbox{ext}}}$ accounts for external forces. Instead of
solving the hydrodynamic equations directly, the LB method evolves in time the
probability distribution functions that can stream along a finite set of
directions, according to the finite set of kinetic velocities
$\bm{\bm{c}}_{i}$ ($i=0...Q-1$). We implement the so-called D3Q19 velocity
scheme, featuring $Q=19$ directions [36, 57], with kinetic velocities

$\bm{\bm{c}}_{i}=\begin{cases}(0,0,0)&i=0\\\ (\pm\Delta x/\Delta
t,0,0),(0,\pm\Delta x/\Delta t,0),(0,0,\pm\Delta x/\Delta t)&i=1-6\\\
(\pm\Delta x/\Delta t,\pm\Delta x/\Delta t,0),(0,\pm\Delta x/\Delta
t,\pm\Delta x/\Delta t),(\pm\Delta x/\Delta t,0,\pm\Delta x/\Delta t)&i=7-18\
,\end{cases}$ (15)

where $\Delta x$ is the lattice spacing and $\Delta t$ is the discretised time
step. The probability distribution function on the fluid (Eulerian) node with
coordinates $\bm{\bm{x}}$ moving with discrete velocity $\bm{\bm{c}}_{i}$ at
time $t$ is represented by $\mbox{n}_{i}(\bm{\bm{x}},t)$, and the LB equation
reads:

$\mbox{n}_{i}(\bm{\bm{x}}+\bm{\bm{c}}_{i}\Delta t,t+\Delta
t)-\mbox{n}_{i}(\bm{\bm{x}},t)=\Delta
t\left[\Omega_{i}(\bm{\bm{x}},t)+S_{i}(\bm{\bm{x}},t)\right]\ ,$ (16)

where the left-hand side represents the streaming along the direction $i$,
while in the right-hand side we find the collision term:

$\Omega_{i}(\bm{\bm{x}},t)=-\frac{1}{\tau_{\mbox{LB}}}[\mbox{n}_{i}(\bm{\bm{x}},t)-\mbox{n}_{i}^{(\mbox{eq})}(\bm{\bm{x}},t)]\
,$ (17)

expressing the relaxation towards a local equilibrium
$\mbox{n}_{i}^{(\mbox{eq})}(\bm{\bm{x}},t)$ with characteristic relaxation
time $\tau_{\mbox{LB}}$. The relaxation process is supplemented by the
presence of external forces implemented in the LB via the source term
$S_{i}(\bm{\bm{x}},t)$ according to the so-called Guo scheme [58]:

$S_{i}(\bm{\bm{x}},t)=\left(1-\frac{\Delta
t}{2\tau_{\mbox{LB}}}\right)\frac{w_{i}}{c_{s}^{2}}\left[\left(\frac{\bm{\bm{c}}_{i}\cdot\bm{\bm{u}}}{c_{s}^{2}}+1\right)\bm{\bm{c}}_{i}-\bm{\bm{u}}\right]\cdot\bm{\bm{F}}_{{\mbox{ext}}}\
,$ (18)

where $w_{i}=w(|\bm{\bm{c}}_{i}|^{2})$ are suitable weights, such that
$w_{0}=1/3$, $w_{1-6}=1/18$, $w_{7-18}=1/36$. The populations $\mbox{n}_{i}$
are used to compute the hydrodynamic density and momentum fields:

$\displaystyle\rho(\bm{\bm{x}},t)=\sum_{i}\mbox{n}_{i}(\bm{\bm{x}},t)\;,$
(19a)
$\displaystyle\rho\bm{\bm{u}}(\bm{\bm{x}},t)=\sum_{i}\bm{\bm{c}}_{i}\mbox{n}_{i}(\bm{\bm{x}},t)\
.$ (19b)

Finally, the local equilibrium depends on $\bm{\bm{x}}$ and $t$ via the
density $\rho(\bm{\bm{x}},t)$ and velocity field $\bm{\bm{u}}(\bm{\bm{x}},t)$:

$\mbox{n}_{i}^{(\mbox{eq})}(\bm{\bm{x}},t)=w_{i}\rho\left(1+\frac{\bm{\bm{u}}\cdot\bm{\bm{c}}_{i}}{c_{s}^{2}}+\frac{(\bm{\bm{u}}\cdot\bm{\bm{c}}_{i})^{2}}{2c_{s}^{4}}-\frac{\bm{\bm{u}}\cdot\bm{\bm{u}}}{c_{s}^{2}}\right)\
,$ (20)

where $c_{\mbox{s}}=\Delta x/\Delta t\sqrt{3}$ is the speed of sound. Details
on the computation of the force $\bm{\bm{F}}_{{\mbox{ext}}}$ are given below.
The LB equation recovers the continuity and Navier-Stokes equations (see Eq.
(13) and Eq. (14), respectively) when fluctuations around the local
equilibrium are small. In such a case, hydrodynamic description is granted
with viscosity $\mu=\rho c_{\mbox{s}}^{2}\left(\tau_{\mbox{LB}}-\frac{\Delta
t}{2}\right)$ and pressure $p=c_{s}^{2}\rho$.

The interface of the droplet is represented by a finite set of Lagrangian
points linked together to build a 3D triangular mesh: on each triangular
element, the stress is given by the surface tension
$\bm{\bm{\tau}}_{\sigma}^{(\bm{\bm{S}})}=\sigma\bm{I}$. Details on the
computation of the force $\bm{\bm{\varphi}}_{i}$ on the $i-$th Lagrangian
point from the stress tensor are given in [39]. The force
$\bm{\bm{\varphi}}_{i}$ is thus communicated from the Lagrangian point with
coordinates $\bm{\bm{r}}_{i}$ to the surrounding Eulerian (fluid) nodes, while
the velocity of the Lagrangian point $\dot{\bm{\bm{r}}}_{i}(t)$ is retrieved
from the fluid velocity computed on the neighbouring Eulerian nodes (see Eq.
(19b)). Such a two-way coupling is handled via an interpolation according to
the IB method [59, 36]:

$\displaystyle\bm{\bm{F}}_{{\mbox{ext}}}(\bm{\bm{x}},t)=\sum_{i}\bm{\bm{\varphi}}_{i}(t)\Delta(\bm{\bm{r}}_{i}-\bm{\bm{x}})\
,$ (21a)
$\displaystyle\dot{\bm{\bm{r}}}_{i}(t)=\sum_{\bm{\bm{x}}}\bm{\bm{u}}(\bm{\bm{x}},t)\Delta(\bm{\bm{r}}_{i}-\bm{\bm{x}})\Delta
x^{3}\ ,$ (21b)

where the $\Delta$ function is the discretised Dirac delta (see [36]).
IB-LB numerical simulations are performed in a 3D box with size
$L_{\mbox{x}}\times L_{\mbox{y}}\times L_{\mbox{z}}=120\times 200\times 120$
lattice units. The linear shear flow is generated via the bounce-back method
for moving walls [36] that are located in $y=\pm L_{\mbox{y}}/2$ with velocity
$\bm{\bm{U}}_{\mbox{w}}=(\pm\dot{\gamma}L_{\mbox{y}}/2,0,0)$, respectively,
while periodic boundary conditions are set along the $x$ and $z$ directions.
The channel width $L_{\mbox{y}}$ is chosen to be large enough to avoid effects
of confinement on the droplet deformation [60]. Once the linear shear flow is
developed, the droplet is placed in the centre of the 3D box. The interface of
the droplet is discretised with a 3D triangular mesh made of 81920 faces and
with radius $R=19$ lattice units. We further recover Stokes equation and the
incompressibility equations by choosing a very small Reynolds number
$\text{Re}=\rho Ru/\mu$ and allowing for small density fluctuations.
Simulations are performed in the limit of small Reynolds number (i.e.,
$\text{Re}<10^{-3}$). The information on the major and minor axes is retrieved
from the eigenvalues $\mathcal{I}_{i}$ of the inertia tensor
$\bm{\bm{\mathcal{I}}}$, that is,
$L=2\sqrt{5(\mathcal{I}_{2}+\mathcal{I}_{3}-\mathcal{I}_{1})/(2M)}$ and
$B=2\sqrt{5(\mathcal{I}_{2}+\mathcal{I}_{3}-\mathcal{I}_{1})/(2M)}$, where M
is the mass of the drop [37]. To compute the inclination angle $\theta$, we
use again the coordinates of the major eigenvector.

Figure 2: Time evolution of the Taylor index $D$ at increasing Ca (lighter to
darker colors). Time is made dimensionless via droplet time $\tau=\mu
R/\sigma$. Numerical simulations data (points) are fitted via Eq. (22) (solid
lines). Different viscosity ratios $\lambda$ are considered. Panel (a):
$\lambda=0.2$ and $\text{Ca}=0.1$ ( ), $\text{Ca}=0.2$ ( ),$\text{Ca}=0.3$ (
), $\text{Ca}=0.4$ ( ); Panel (b): $\lambda=1$ and $\text{Ca}=0.1$ ( ),
$\text{Ca}=0.2$ ( ), $\text{Ca}=0.3$ ( ), $\text{Ca}=0.4$ ( ); Panel (c):
$\lambda=2$ and $\text{Ca}=0.15$ ( ), $\text{Ca}=0.3$ ( ), $\text{Ca}=0.45$ (
), $\text{Ca}=0.55$ ( ). Figure 3: Comparison between numerical simulations
data (points), a stretched exponential fit (solid lines, Eq. (22) with
$\delta\neq 1$) and pure exponential fit (dashed lines, Eq. (22) with
$\delta=1$) for the time evolution of the Taylor index $D$ for different
$\lambda$ and selected values of Ca. Time is made dimensionless via droplet
time $\tau=\mu R/\sigma$. Panel (a): $\lambda=0.2$, $\text{Ca}=0.4$ ( ),
$\delta=0.742$; Panel (b): $\lambda=1$, $\text{Ca}=0.4$ ( ), $\delta=0.723$;
Panel (c): $\lambda=2$, $\text{Ca}=0.55$ ( ), $\delta=0.796$. Figure 4: Panel
(a): characteristic dimensionless creeping frequency $f_{1}^{(\mbox{EMM})}$ as
a function of the capillary number Ca (see Eq. (22)) for $\lambda=0.2$ ( ),
$\lambda=1$ ( ), $\lambda=2$ ( ). The horizontal dashed lines indicate the
values of $f_{1}^{(\mbox{MM})}(\lambda)$ (see Eq. (2)), while solid lines
refer to a polynomial fit for $f_{1}^{(\mbox{EMM})}$ (see Eq. (23) and Table
1). Panel (b): stretching factor $\delta$ (see Eq. (22)) as a function of Ca
for different values of $\lambda$. The horizontal dashed line is drawn in
correspondence of $\delta=1$.

##  4 Results and discussions

We start our investigation by analyzing the process of droplet deformation for
different viscosity ratios $\lambda$. We choose $\lambda$ to span one decade,
$\lambda\in[0.2,2]$, which is a realistic range [61, 62] where we can
characterize the droplet deformation process from the small deformations up to
the breakup occurring at moderate values of capillary number Ca [5, 27, 32].
In Fig. 2, we report the evolution of the Taylor index (see Fig. 1) for
different values of the capillary number Ca. Time is made dimensionless with
respect to the droplet time $\tau=\mu R/\sigma$. The linearization of the MM
model would imply the deformation to exponentially increase with a single
characteristic time equal to $f_{1}^{-1}$ up to the steady deformation value
$D_{\infty}$, i.e., $D(t)=D_{\infty}(1-e^{-f_{1}t/\tau})$: this behaviour is
in fact in good agreement with the numerical data only for small deformations,
while for larger deformations we find that the following stretched exponential
function

$D(t)=D_{\infty}\left[1-e^{-(f_{1}t/\tau)^{\delta}}\right]\ ,$ (22)

better adapts to the data, where $\delta$ is the stretching factor used as a
fitting parameter. The improved agreement brought by the stretched exponential
function is highlighted in Fig. 3. The stretched exponential function is
mostly used to describe non-linear relaxation processes [63, 64] and its basis
can be understood in the context of the relaxation of homogeneous glasses and
glass-forming liquids [65, 66], i.e., physical systems that are known to
exhibit multiple relaxation times [67, 68, 69]. For $\delta=1$, we reduce to a
linear relaxation process with a single characteristic time, while $\delta<1$
signals the emergence of a spectrum of characteristic times wherein the value
$\tau_{\mbox{{creep}}}=\left[f_{1}^{(\mbox{EMM})}\right]^{-1}$ is a
representative one. The values of $\delta$ and $f_{1}^{(\mbox{EMM})}$ that we
extracted from this fitting procedure are displayed in Fig. 4, showing that an
increase in Ca results in a decrease of both $\delta$ and
$f_{1}^{(\mbox{EMM})}$. As the droplet approaches large deformations, it gets
closer to the critical point of breakup and the characteristic time for the
deformation process $\tau_{\mbox{{creep}}}$ is expected to get larger, i.e.,
the characteristic creeping frequency $f_{1}^{(\mbox{EMM})}$ gets smaller [70,
32]. To characterize the approach towards criticality in a more quantitative
way, we analyse the subcritical scaling of the creeping time
$\tau_{\mbox{{creep}}}$, expected to follow a power law
$\tau_{\mbox{{creep}}}\sim(\text{Ca}_{\mbox{cr}}-\text{Ca})^{-1/2}$ as showed
in [70]. In Fig. 5, we show that the creeping time $\tau_{\mbox{{creep}}}$
follows the above mentioned scaling in the subcritical regime, which has been
used to evaluate the critical capillary number for breakup
$\text{Ca}_{\mbox{cr}}$. We evaluate $\text{Ca}_{\mbox{cr}}\sim 0.49$ for
$\lambda=0.2$, $\text{Ca}_{\mbox{cr}}\sim 0.435$ for $\lambda=1$ and
$\text{Ca}_{\mbox{cr}}\sim 0.62$ for $\lambda=2$: these values are in good
agreement with other evaluations of $\text{Ca}_{\mbox{cr}}$ reported in the
literature for a droplet in simple shear flow [27, 32].

Figure 5: Subcritical scaling for the dimensionless characteristic creeping
time $\tau_{\mbox{{creep}}}=[f_{1}^{(\mbox{EMM})}]^{-1}$ for $\lambda=0.2$ (
), $\lambda=1$ ( ) and $\lambda=2$ ( ). Dashed lines indicate the scaling
prediction with exponent $-1/2$ [70]. Figure 6: Stationary inclination angles
for different $\lambda$. We report results of numerical simulations (points)
for $\lambda=0.2$ ( , Panel (a)), $\lambda=1$ ( , Panel (b)), $\lambda=2$ ( ,
Panel (c)), as well as the prediction of Eq. (10) with
$f_{1}=f_{1}^{(\mbox{MM})}$ (dashed lines) and with
$f_{1}=f_{1}^{(\mbox{EMM})}$ (solid lines). The colored regions refer to the
evaluated breakup regime.  Figure 7: Stationary droplet deformations for
different $\lambda$. We report results of numerical simulations (points) for
$\lambda=0.2$ ( , Panel (a)), $\lambda=1$ ( , Panel (b)), $\lambda=2$ ( ,
Panel (c)), as well as the prediction of Eq. (9) with
$f_{1,2}=f_{1,2}^{(\mbox{MM})}$ (dashed lines) and the one with the values
$f_{1}=f_{1}^{(\mbox{EMM})}$ and the optimal $f_{2}=f_{2}^{(\mbox{EMM})}$ that
allow to match the numerical results (solid lines). The colored regions refer
to the evaluated breakup regime.

We now consider the prediction of the MM model for the steady inclination
angle $\theta_{\infty}$ in Eq. (10). In Fig. 6, we report a comparison between
the steady values of the inclination angle $\theta_{\infty}$ extracted from
numerical simulations (filled points) and those predicted via Eq. (10) with
the coefficient $f_{1}^{(\mbox{EMM})}$ being estimated from the creeping times
(solid lines); the prediction based on Eq. (10) with $f_{1}^{(\mbox{MM})}$ is
also reported (dashed lines). The latter well adapts to the results from
numerical simulations for small values of Ca, while a mismatch appears for
larger values of Ca. The prediction based on the EMM model coefficient
$f_{1}^{(\mbox{EMM})}$, instead, shows a good match with the numerical data,
with only a slight mismatch that is observed for the largest value of
viscosity ratio $\lambda$ considered. This is due to the development of an
overshooting behaviour in the transient dynamics at high Ca that is not caught
via a monotonic function such as Eq. (22). As shown in Fig. 9, the overshoot
grows in $\lambda$, driving the slight mismatch observed in Fig. 6. We have
verified that by promoting the fitting function given in Eq. (22) to a non
monotonous function would improve the accuracy of the creeping time
evaluation. Overall, the agreement that we observe between the results of
numerical simulations and steady state predictions from the EMM model is
remarkable. We believe this result is highly non trivial, in that we are able
to predict the steady state inclination angle via an independent measurement
of the creeping time.

We now want to consider the prediction for the steady deformation in Eq. (9)
and show the strategy adopted to compute $f_{2}$. Once we know the value of
$f_{1}^{(\mbox{EMM})}$, we can look for an optimal value of
$f_{2}=f_{2}^{(\mbox{EMM})}$ such that the prediction of the deformation in
Eq. (9) matches the data coming from numerical simulations for finite Ca. This
matching procedure is illustrated in Fig. 7, where we report the results for
the steady deformation extracted from the numerical simulations and the
prediction based on Eq. (9) using the MM model parameters
$f_{1,2}^{(\mbox{MM})}$ (dashed lines). As already observed for the steady
inclination angle (see Fig. 6), the MM prediction in Eq. (9) with
$f_{1,2}=f_{1,2}^{(\mbox{MM})}$ well adapts to the data for small Ca, while a
mismatch emerges at $\text{Ca}\sim{\cal O}(1)$. Our working strategy is then
to identify, for each value of Ca, the value of $f_{2}^{(\mbox{EMM})}$ that
allows Eq. (9) to perfectly match with the simulation data. We find that this
matching is possible for the values of $\lambda$ that we considered (solid
lines). The values of $f_{2}^{(\mbox{EMM})}$ that we extracted with this
procedure are reported in Fig. 8 and also compared with an heuristic
$\text{Ca}-$dependent formula proposed by MM (Eq.(34) in [16]). The latter
compares very well with our proposed extensions only for the case $\lambda>1$,
where we report a very similar increase in $f_{2}$ at increasing Ca.
Qualitatively, the model parameter $f_{2}$ is expected to acquire a dependence
on Ca, in that higher values of Ca would imply a droplet deformation process
mainly driven by the viscous forces, while surface tension effects become
negligible in such situation. More quantitatively, the model parameter $f_{2}$
can be viewed as an “affinity" parameter [71] accounting for an additional
surface slip which effectively removes part of the droplet elongation, whereas
the full rigid body rotation is retained coherently with the rotating flow
[72]. When $\lambda\neq 1$, a common behaviour is shown in the approaching of
the affine limit $f_{2}=1$ at higher Ca: the affine motion of the interface
changes with respect to the magnitude of the capillary number and the value
$f_{2}=1$ is approached near or at the corresponding critical configuration
for the droplet, signaling that an affine motion is ultimately required for
the breakup [73, 74]. We also notice that for the case $\lambda=1$ we have
$f_{2}^{(\mbox{MM})}=1$, so that in the case of equiviscous fluids the
effective slip at the interface is lacking and the droplet is already in the
affine configuration at small Ca. We have fitted the coefficients
$f_{1,2}^{(\mbox{EMM})}$ with a polynomial expansion in the capillary number
Ca:

$f_{1,2}^{(\mbox{EMM})}(\lambda,\text{Ca})=\sum_{i}a_{i}^{(1,2)}(\lambda)\text{Ca}^{i}\
,$ (23)

in Tab. 1 we report details for a forth-order fit for $\lambda=0.2,1,2$.

| $\lambda=0.2$ | $\lambda=1$ | $\lambda=2$
---|---|---|---
$i$ | $a_{i}^{(1)}$ | $a_{i}^{(2)}$ | $a_{i}^{(1)}$ | $a_{i}^{(2)}$ | $a_{i}^{(1)}$ | $a_{i}^{(2)}$
0 | 0.713 | 1.470 | 0.457 | 1 | 0.317 | 0.714
1 | 0.416 | 0.347 | 0.235 | 0 | 0.270 | -0.028
2 | -5.333 | -6.040 | -4.546 | 0 | -2.283 | 3.996
3 | 16.776 | 20.021 | 20.536 | 0 | 8.338 | -9.188
4 | -24.339 | -27.703 | -34.798 | 0 | -9.91 | 6.611
Table 1: Values of the interpolated coefficients $a_{i}^{(1,2)}(\lambda)$ for
$f_{1,2}^{(\mbox{EMM})}$ (Eq. (23), with forth order interpolation). Figure 8:
We report the coefficient $f_{2}^{(\mbox{EMM})}$ as a function of the
capillary number Ca that allows to match the prediction of the steady
deformation in Eq. (9) with numerical simulations data (points) (see Fig. 7)
for $\lambda=0.2$ ( ), $\lambda=1$ ( ), $\lambda=2$ ( ). The colored regions
refer to a range of Ca where breakup occurs. The heuristic prediction reported
in Eq. (34) of [16] is also displayed (colored dashed lines). The limit where
the affine deformation takes place ($f_{2}=1$) is indicated with a black
dashed line; solid lines refer to a polynomial fit for $f_{2}^{(\mbox{EMM})}$
as in Table 1.

Lastly, in Fig. 9, we report the transient values of the Taylor index for
different viscosity ratios $\lambda$, comparing the EMM (Eq. (2) with
$f_{1,2}=f_{1,2}^{(\mbox{EMM})}$) and the original MM model (Eq. (2) with
$f_{1,2}=f_{1,2}^{(\mbox{MM})}$). When comparing model predictions and data
from numerical simulations in the transient region, the EMM model matches
better than the MM one. When either the viscosity ratio $\lambda$ or the
capillary number Ca increase, the Taylor index $D$ predicted by both MM and
EMM models shows an overshoot111Numerical simulations also show overshoots,
but they are tiny and cannot be well appreciated in Fig. 9.. In particular,
this overshoot causes the EMM model to overestimate the data of the numerical
simulations during the transient dynamics; however, the steady value of the
deformation predicted by the EMM model, on the contrary of the MM model,
matches (by construction) the one observed in numerical simulations. In order
to delve deeper into the origin of the mismatch observed in Fig. 9 between
model predictions and simulation data, we analyze in Fig. 10 the shape of the
droplet in correspondence of some selected times for the highest Ca analyzed
in Fig. 9 (see arrows therein). We compare the shapes predicted by the EMM
model (colored ellipses) with those coming from the numerical simulations
(filled points). Dashed lines represent ellipses with both the major and minor
axes ($L$ and $B$) and inclination angle $\theta$ coming from the inertia
tensor computed in the numerical simulations (see Sec. 3). Filled points and
dashed lines in Fig. 10 show a good agreement, suggesting that the shape of
the droplet remains ellipsoidal during the whole numerical simulation. As
$\lambda$ grows, Panels (c,f,i) highlight that the rotation of the droplet
induced by $\bm{\bm{\Omega}}$ leads to a slight mismatch between model
predictions and numerical simulations. More quantitatively, the discrepancy
for $L$ and $B$ that we observe for the highest $\lambda$ in Fig. 10 is of the
order of $10\%$ and $15\%$ in Panels (c and f), respectively; the error
reduces to about $5\%$ in Panel (i). Concerning the inclination angle, the
mismatch is of the order of $15\%$ in Panel (c), $5\%$ in Panel (f) and up to
$30\%$ in Panel (i). In the transient equiviscous case (Panel (e)), a mismatch
for $L$, $B$ and $\theta$ of the same order of Panels (c and f), is observed.
We hasten to remark that discrepancies between numerical simulations and EMM
model predictions show up only when we move close to the critical point, while
being very mitigated for smaller - but still finite Ca.

Figure 9: Time evolution of the Taylor index for different $\lambda$ and Ca.
Numerical simulations data (points) are compared with the original MM model as
in Eq. (2) with $f_{1,2}=f_{1,2}^{(\mbox{MM})}$. We also report the
predictions of Eq. (2) with $f_{1,2}=f_{1,2}^{(\mbox{EMM})}$ that we evaluated
from the creeping times and the steady deformations (see Figs. 4 and 8). Panel
(a): $\lambda=0.2$ and $\text{Ca}=0.1$ ( ), $\text{Ca}=0.2$ (
),$\text{Ca}=0.3$ ( ), $\text{Ca}=0.4$ ( ); Panel (b): $\lambda=1$ and
$\text{Ca}=0.1$ ( ), $\text{Ca}=0.2$ ( ), $\text{Ca}=0.3$ ( ), $\text{Ca}=0.4$
( ); Panel (c): $\lambda=2$ and $\text{Ca}=0.15$ ( ), $\text{Ca}=0.3$ ( ),
$\text{Ca}=0.45$ ( ), $\text{Ca}=0.55$ ( ). Figure 10: Droplet shapes in the
shear plane for selected times indicated by arrows in Fig. 9. Numerical
results for $\lambda=0.2$ ( ), $\lambda=1$ ( ) $\lambda=2$ ( ) are compared
with the ellipsoidal shape (dashed lines) obtained via the major/minor axes
$L,B$ and inclination angle extracted from the inertia tensor in numerical
simulations (see Sec. 2 for details). Colored areas indicate the ellipsoidal
shapes obtained via the use of Eq. (2) with $f_{1,2}=f_{1,2}^{(\mbox{EMM})}$.

##  5 Summary and Conclusions

We have studied the deformation process of a droplet under simple shear flow
at finite capillary numbers ($\text{Ca}\sim{\cal O}(1)$) and extracted
relevant physical information to extend the popular Maffettone-Minale (MM)
model [16]. In the latter, the shape of the droplet is described by the
second-order tensor $\bm{\bm{S}}$, while the flow properties are expressed by
the symmetric ($\bm{\bm{E}}$) and asymmetric ($\bm{\bm{\Omega}}$) parts of the
velocity gradient $\bm{\nabla}\bm{\bm{u}}$ (see Sec. 2). The MM model hinges
on two parameters ($f_{1}(\lambda)$ and $f_{2}(\lambda)$) that are linked to
the creeping time of the droplet and the deformation induced by the flow,
respectively (see Eq. (2)). In the MM model, such parameters are found by
requiring a theoretical matching with results from perturbation theory at
small Ca [1, 7, 17]: for this reason, the MM model is not appropriate for a
quantitative prediction of the deformation $D$ and the inclination angle
$\theta$ for finite values of the capillary number $\text{Ca}\sim{\cal O}(1)$
(see Figs. 6-7). In this paper, we provide an extension of the MM model
(referred as EMM model) via a direct evaluation of both the creeping time and
the steady deformation at $\text{Ca}\sim{\cal O}(1)$ with independent
numerical simulations. The latter have been performed with the Immersed
Boundary - Lattice Boltzmann (IB-LB) method [36, 38, 39], that allowed precise
characterization of the droplet dynamics from small deformations up to
subcritical regimes. The proposed EMM model is given by:

$\frac{d\bm{\bm{S}}}{dt^{\prime}}-\text{Ca}\,\left(\bm{\Omega}\cdot\bm{S}-\bm{S}\cdot\bm{\Omega}\right)=-f_{1}^{(\mbox{EMM})}(\lambda,\text{Ca})\left[\bm{S}-g\left(\bm{S}\right)\bm{I}\right]+\text{Ca}\,f_{2}^{(\mbox{EMM})}(\lambda,\text{Ca})\left(\bm{E}\cdot\bm{\bm{S}}+\bm{\bm{S}}\cdot\bm{E}\right)\
,$ (24)

where the model coefficients $f_{1,2}^{(\mbox{EMM})}(\lambda,\text{Ca})$ have
been characterized as a function of both $\lambda$ and Ca (see Figs. 4, 5, 8
and Eq. (23)). The EMM model embeds a realistic creeping time as well as a
realistic deformation at $\text{Ca}\sim{\cal O}(1)$, and its predictions on
the steady deformation angle are found to be in good agreement with the
results of numerical simulations (see Fig. 6). This result is not trivial,
since we are able to predict the steady deformation angle via an independent
measure of the creeping time; moreover, it can suggest purposeful experimental
setups aimed at the evaluation of stationary state properties. Improved
agreement is also observed during the whole transient process of droplet
deformation (see Fig. 9), although some enhanced overshoots are observed with
respect to the numerical simulations, suggesting that a further model
extension is needed to better capture the full dynamics of droplet
deformation. Investigating whether these extensions would be possible by
retaining a second order tensorial description or by including an higher order
expansion is surely an interesting topic. Also, the results in Fig. 5 pave the
way for the determination of $\text{Ca}_{\mbox{cr}}$ via the EMM model,
including a detailed analysis for the characterization of the droplet shapes
very close to criticality. Another prospective for future research would be to
investigate the degree of universality of our results with respect to a change
in the topology of the driving flows. This could help, for example, to extend
earlier studies on the statistics of droplet deformation and breakup in a
generic time-dependent flow [75, 76, 77], in that one can study how such
statistics is affected by the interplay between the characteristic times of
the driving flow and the characteristic droplet time that depends on the local
flow intensity.

## Acknowledgments

We wish to acknowledge Fabio Bonaccorso for support. We also acknowledge
useful discussion with Luca Biferale. This work received funding from the
European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation programme (grant agreement No 882340).

## References

  * [1] Geoffrey Ingram Taylor  “The viscosity of a fluid containing small drops of another fluid”  In _Proceedings of the Royal Society of London. Series A_ 138.834 The Royal Society London, 1932, pp. 41–48
  * [2] Geoffrey Ingram Taylor  “The formation of emulsions in definable fields of flow”  In _Proceedings of the Royal Society of London. Series A_ 146.858 The Royal Society London, 1934, pp. 501–523
  * [3] Albert Einstein  “A new determination of molecular dimensions”  In _Ann. Phys._ 19, 1906, pp. 289–306
  * [4] Charles E Chaffey and Howard Brenner  “A second-order theory for shear deformation of drops”  In _Journal of Colloid and Interface Science_ 24.2 Elsevier, 1967, pp. 258–269
  * [5] D Barthes-Biesel and Andreas Acrivos  “Deformation and burst of a liquid droplet freely suspended in a linear shear field”  In _Journal of Fluid Mechanics_ 61.1 Cambridge University Press, 1973, pp. 1–22
  * [6] RG Cox  “The deformation of a drop in a general time-dependent fluid flow”  In _Journal of Fluid Mechanics_ 37.3 Cambridge University Press, 1969, pp. 601–623
  * [7] NA Frankel and Andreas Acrivos  “The constitutive equation for a dilute emulsion”  In _Journal of Fluid Mechanics_ 44.1 Cambridge University Press, 1970, pp. 65–78
  * [8] FS Hakimi and WR Schowalter  “The effects of shear and vorticity on deformation of a drop”  In _Journal of Fluid Mechanics_ 98.3 Cambridge University Press, 1980, pp. 635–645
  * [9] GK Youngren and A Acrivos  “On the shape of a gas bubble in a viscous extensional flow”  In _Journal of Fluid Mechanics_ 76.3 Cambridge University Press, 1976, pp. 433–442
  * [10] Francesco Greco  “Drop deformation for non-Newtonian fluids in slow flows”  In _Journal of Non-Newtonian Fluid Mechanics_ 107.1-3 Elsevier, 2002, pp. 111–131
  * [11] Anja Vananroye, Peter Van Puyvelde and Paula Moldenaers  “Effect of confinement on droplet breakup in sheared emulsions”  In _Langmuir_ 22.9 ACS Publications, 2006, pp. 3972–3974
  * [12] JM Rallison  “The deformation of small viscous drops and bubbles in shear flows”  In _Annual Review of Fluid Mechanics_ 16.1 Annual Reviews 4139, 1984, pp. 45–66
  * [13] Howard A Stone  “Dynamics of drop deformation and breakup in viscous fluids”  In _Annual Review of Fluid Mechanics_ 26.1 Annual Reviews 4139, 1994, pp. 65–102
  * [14] Peter Fischer and Philipp Erni  “Emulsion drops in external flow fields—the role of liquid interfaces”  In _Current Opinion in Colloid & Interface Science_ 12.4-5 Elsevier, 2007, pp. 196–205
  * [15] Vittorio Cristini and Yung-Chieh Tan  “Theory and numerical simulation of droplet dynamics in complex flows—a review”  In _Lab on a Chip_ 4.4 Royal Society of Chemistry, 2004, pp. 257–264
  * [16] PL Maffettone and Mario Minale  “Equation of change for ellipsoidal drops in viscous flow”  In _Journal of Non-Newtonian Fluid Mechanics_ 78.2-3 Elsevier, 1998, pp. 227–241
  * [17] JM Rallison  “Note on the time-dependent deformation of a viscous drop which is almost spherical”  In _Journal of Fluid Mechanics_ 98.3 Cambridge University Press, 1980, pp. 625–633
  * [18] S Torza, RG Cox and SG Mason  “Particle motions in sheared suspensions XXVII. Transient and steady deformation and burst of liquid drops”  In _Journal of Colloid and Interface Science_ 38.2 Elsevier, 1972, pp. 395–411
  * [19] Stefano Guido and Marco Villone  “Three-dimensional shape of a drop under simple shear flow”  In _Journal of Rheology_ 42.2 The Society of Rheology, 1998, pp. 395–415
  * [20] BJ Bentley and L Gary Leal  “An experimental investigation of drop deformation and breakup in steady, two-dimensional linear flows”  In _Journal of Fluid Mechanics_ 167 Cambridge University Press, 1986, pp. 241–283
  * [21] Pier Luca Maffettone and Francesco Greco  “Ellipsoidal drop model for single drop dynamics with non-Newtonian fluids”  In _Journal of Rheology_ 48.1 The Society of Rheology, 2004, pp. 83–100
  * [22] Mario Minale  “Deformation of a non-Newtonian ellipsoidal drop in a non-Newtonian matrix: extension of Maffettone–Minale model”  In _Journal of Non-Newtonian Fluid Mechanics_ 123.2-3 Elsevier, 2004, pp. 151–160
  * [23] S Guido, M Simeone and F Greco  “Deformation of a Newtonian drop in a viscoelastic matrix under steady shear flow: experimental validation of slow flow theory”  In _Journal of Non-Newtonian Fluid Mechanics_ 114.1 Elsevier, 2003, pp. 65–82
  * [24] Mario Minale  “A phenomenological model for wall effects on the deformation of an ellipsoidal drop in viscous flow”  In _Rheologica acta_ 47.5-6 Springer, 2008, pp. 667–675
  * [25] Mario Minale, Sergio Caserta and Stefano Guido  “Microconfined shear deformation of a droplet in an equiviscous non-Newtonian immiscible fluid: Experiments and modeling”  In _Langmuir_ 26.1 ACS Publications, 2010, pp. 126–132
  * [26] Dhruv Arora, Marek Behr and Matteo Pasquali  “A tensor-based measure for estimating blood damage.”  In _Artificial organs_ 28 11, 2004, pp. 1002–15
  * [27] Harold P Grace  “Dispersion phenomena in high viscosity immiscible fluid systems and application of static mixers as dispersion devices in such systems”  In _Chemical Engineering Communications_ 14.3-6 Taylor & Francis, 1982, pp. 225–277
  * [28] Abdulwahab S Almusallam, Ronald G Larson and Michael J Solomon  “A constitutive model for the prediction of ellipsoidal droplet shapes and stresses in immiscible blends”  In _Journal of Rheology_ 44.5 The Society of Rheology, 2000, pp. 1055–1083
  * [29] Nancy E Jackson and Charles L Tucker III  “A model for large deformation of an ellipsoidal droplet with interfacial tension”  In _Journal of Rheology_ 47.3 The Society of Rheology, 2003, pp. 659–682
  * [30] Wei Yu and Mosto Bousmina  “Ellipsoidal model for droplet deformation in emulsions”  In _Journal of Rheology_ 47.4 The Society of Rheology, 2003, pp. 1011–1039
  * [31] Mario Minale  “Models for the deformation of a single ellipsoidal drop: a review”  In _Rheologica acta_ 49.8 Springer, 2010, pp. 789–806
  * [32] V Cristini et al.  “Drop breakup and fragment size distribution in shear flow”  In _Journal of Rheology_ 47.5 The Society of Rheology, 2003, pp. 1283–1298
  * [33] Zhi-Gang Feng and Efstathios E Michaelides  “The immersed boundary-lattice Boltzmann method for solving fluid–particles interaction problems”  In _Journal of Computational Physics_ 195.2 Elsevier, 2004, pp. 602–628
  * [34] Junfeng Zhang, Paul C Johnson and Aleksander S Popel  “An Immersed Boundary Lattice Boltzmann Approach to Simulate Deformable Liquid Capsules and Its Application to Microscopic Blood Flows”  In _Physical biology_ 4.4, 2007, pp. 285–295
  * [35] Michael M Dupin et al.  “Modeling the flow of dense suspensions of deformable particles in three dimensions”  In _Physical Review E_ 75.6 APS, 2007, pp. 066707
  * [36] Timm Krüger et al.  “The Lattice Boltzmann Method - Principles and Practice”, 2016
  * [37] T. Krüger  “Computer Simulation Study of Collective Phenomena in Dense Suspensions of Red Blood Cells under Shear”  Springer, 2012
  * [38] Ping Li and Junfeng Zhang  “A finite difference method with subsampling for immersed boundary simulations of the capsule dynamics with viscoelastic membranes”  In _International Journal for Numerical Methods in Biomedical Engineering_ 35.6 Wiley Online Library, 2019, pp. e3200
  * [39] Fabio Guglietta et al.  “On the effects of membrane viscosity on transient red blood cell dynamics”  In _Soft Matter_ 16.26 Royal Society of Chemistry, 2020, pp. 6191–6205
  * [40] Ali Rezghi and Junfeng Zhang  “Tank-treading dynamics of red blood cells in shear flow: On the membrane viscosity rheology”  In _Biophysical Journal_ 121.18 Elsevier, 2022, pp. 3393–3410
  * [41] Ali Rezghi, Ping Li and Junfeng Zhang  “Lateral migration of viscoelastic capsules in tube flow”  In _Physics of Fluids_ 34.1 AIP Publishing LLC, 2022, pp. 011906
  * [42] Ping Li and Junfeng Zhang  “Similar but Distinct Roles of Membrane and Interior Fluid Viscosities in Capsule Dynamics in Shear Flows”  In _Cardiovascular Engineering and Technology_ 12.2 Springer, 2021, pp. 232–249
  * [43] Timm Krüger, Markus Gross, Dierk Raabe and Fathollah Varnik  “Crossover from tumbling to tank-treading-like motion in dense simulated suspensions of red blood cells”  In _Soft Matter_ 9.37 Royal Society of Chemistry, 2013, pp. 9008–9015
  * [44] Timm Krüger, David Holmes and Peter V Coveney  “Deformability-based red blood cell separation in deterministic lateral displacement devices—A simulation study”  In _Biomicrofluidics_ 8.5 AIP Publishing, 2014, pp. 054114
  * [45] Alessandro Coclite, Sergio Ranaldo, Giuseppe Pascazio and Marco D Tullio  “A Lattice Boltzmann dynamic-immersed boundary scheme for the transport of deformable inertial capsules in low-Re flows”  In _Computers & Mathematics with Applications_ 80.12 Elsevier, 2020, pp. 2860–2876
  * [46] A Coclite et al.  “Kinematic and dynamic forcing strategies for predicting the transport of inertial capsules via a combined lattice Boltzmann–immersed boundary method”  In _Computers & Fluids_ 180 Elsevier, 2019, pp. 41–53
  * [47] Christian Bächer and Stephan Gekle  “Computational modeling of active deformable membranes embedded in three-dimensional flows”  In _Physical Review E_ 99.6 APS, 2019, pp. 062418
  * [48] Christian Bächer et al.  “Antimargination of microparticles and platelets in the vicinity of branching vessels”  In _Biophysical journal_ 115.2 Elsevier, 2018, pp. 411–425
  * [49] Christian Bächer, Lukas Schrack and Stephan Gekle  “Clustering of microscopic particles in constricted blood flow”  In _Physical Review Fluids_ 2.1 APS, 2017, pp. 013102
  * [50] Stephan Gekle  “Strongly accelerated margination of active particles in blood flow”  In _Biophysical journal_ 110.2 Elsevier, 2016, pp. 514–520
  * [51] Moritz Lehmann, Sebastian Johannes Müller and Stephan Gekle  “Efficient viscosity contrast calculation for blood flow simulations using the lattice Boltzmann method”  In _International Journal for Numerical Methods in Fluids_ 92.11 Wiley Online Library, 2020, pp. 1463–1477
  * [52] JM Rallison and Andrias Acrivos  “A numerical study of the deformation and burst of a viscous drop in an extensional flow”  In _Journal of Fluid Mechanics_ 89.1 Cambridge University Press, 1978, pp. 191–200
  * [53] Constantine Pozrikidis  “Boundary integral and singularity methods for linearized viscous flow”  Cambridge university press, 1992
  * [54] J. Gounley, G. Boedec, M. Jaeger and M. Leonetti  “Influence of surface viscosity on droplets in shear flow”  In _Journal of Fluid Mechanics_ 791 Cambridge University Press, 2016, pp. 464–494
  * [55] Jie Li, Yuriko Y Renardy and Michael Renardy  “Numerical simulation of breakup of a viscous drop in simple shear flow through a volume-of-fluid method”  In _Physics of Fluids_ 12.2 American Institute of Physics, 2000, pp. 269–282
  * [56] Alessandro Coclite, Marco Donato Tullio, Giuseppe Pascazio and Paolo Decuzzi  “A combined Lattice Boltzmann and Immersed boundary approach for predicting the vascular transport of differently shaped particles”  In _Computers & Fluids_ 136 Elsevier, 2016, pp. 260–271
  * [57] Sauro Succi  “The lattice Boltzmann equation: for complex states of flowing matter”  Oxford University Press, 2018
  * [58] Zhaoli Guo, Chuguang Zheng and Baochang Shi  “Discrete lattice effects on the forcing term in the lattice Boltzmann method”  In _Physical review E_ 65.4 APS, 2002, pp. 046308
  * [59] Charles S Peskin  “The immersed boundary method”  In _Acta numerica_ 11 Cambridge University Press, 2002, pp. 479–517
  * [60] Stefano Guido  “Shear-induced droplet deformation: Effects of confined geometry and viscoelasticity”  In _Current opinion in colloid & interface science_ 16.1 Elsevier, 2011, pp. 61–70
  * [61] Richard R Lessard and Greg DeMarco  “The significance of oil spill dispersants”  In _Spill Science & Technology Bulletin_ 6.1 Elsevier, 2000, pp. 59–68
  * [62] Fernando Kassis Carvalho et al.  “Viscosity, surface tension and droplet size of sprays of different formulations of insecticides and fungicides”  In _Crop Protection_ 101 Elsevier, 2017, pp. 19–23
  * [63] JC Phillips  “Stretched exponential relaxation in molecular and electronic glasses”  In _Reports on Progress in Physics_ 59.9 IOP Publishing, 1996, pp. 1133
  * [64] Dmitry A Fedosov  “Multiscale modeling of blood flow and soft matter”  Citeseer, 2010
  * [65] Richard G Palmer, Daniel L Stein, Elihu Abrahams and Philip Warren Anderson  “Models of hierarchically constrained dynamics for glassy relaxation”  In _Physical Review Letters_ 53.10 APS, 1984, pp. 958
  * [66] Marcel Potuzak, Roger C Welch and John C Mauro  “Topological origin of stretched exponential relaxation in glass”  In _The Journal of Chemical Physics_ 135.21 American Institute of Physics, 2011, pp. 214502
  * [67] A Jurlewicz and K Weron  “A relationship between asymmetric Lévy-stable distributions and the dielectric susceptibility”  In _Journal of Statistical Physics_ 73.1 Springer, 1993, pp. 69–81
  * [68] DC Johnston  “Stretched exponential relaxation arising from a continuous sum of exponential decays”  In _Physical Review B_ 74.18 APS, 2006, pp. 184430
  * [69] Daniel C Elton  “Stretched exponential relaxation”  In _arXiv preprint arXiv:1808.00881_ , 2018
  * [70] Jerzy Bławzdziewicz, Vittorio Cristini and Michael Loewenberg  “Critical behavior of drops in linear flows. I. Phenomenological theory for drop dynamics near critical stationary states”  In _Physics of Fluids_ 14.8 American Institute of Physics, 2002, pp. 2709–2718
  * [71] Erich J Windhab et al.  “Emulsion processing—from single-drop deformation to design of complex processes and products”  In _Chemical Engineering Science_ 60.8-9 Elsevier, 2005, pp. 2101–2113
  * [72] RJ Gordon and WR Schowalter  “Anisotropic fluid theory: a different approach to the dumbbell theory of dilute polymer solutions”  In _Transactions of the Society of Rheology_ 16.1 The Society of Rheology, 1972, pp. 79–97
  * [73] PHM Elemans, HL Bos, JMH Janssen and HEH Meijer  “Transient phenomena in dispersive mixing”  In _Chemical engineering science_ 48.2 Elsevier, 1993, pp. 267–276
  * [74] Han EH Meijer, Jos MH Janssen and Patrick D Anderson  “Mixing of immiscible liquids”  In _Mixing and compounding of polymers: theory and practice_ Hanser, New York, 1994, pp. 85–148
  * [75] Luca Biferale, Charles Meneveau and Roberto Verzicco  “Deformation statistics of sub-Kolmogorov-scale ellipsoidal neutrally buoyant drops in isotropic turbulence”  In _Journal of Fluid Mechanics_ 754 Cambridge University Press, 2014, pp. 184–207
  * [76] Said Elghobashi  “Direct numerical simulation of turbulent flows laden with droplets or bubbles”  In _Annual Review of Fluid Mechanics_ 51.1, 2019, pp. 217–244
  * [77] Vamsi Spandan, Detlef Lohse and Roberto Verzicco  “Deformation and orientation statistics of neutrally buoyant sub-Kolmogorov ellipsoidal droplets in turbulent Taylor–Couette flow”  In _Journal of Fluid Mechanics_ 809 Cambridge University Press, 2016, pp. 480–501