# Droplet dynamics in homogeneous isotropic turbulence with the immersed
boundary-lattice Boltzmann method

Diego Taglienti1&Fabio Guglietta1,†&Mauro Sbragaglia1
1Department of Physics & INFN, Tor Vergata University of Rome, Via della
Ricerca Scientifica 1, 00133, Rome, Italy
<EMAIL_ADDRESS>
###### Abstract

We develop a numerical method for simulating the dynamics of a droplet
immersed in a generic time-dependent velocity gradient field. This approach is
grounded on the hybrid coupling between the lattice Boltzmann (LB) method,
employed for the flow simulation, and the immersed boundary (IB) method,
utilized to couple the droplet with the surrounding fluid. We show how to
enrich the numerical scheme with a mesh regularization technique, allowing
droplets to sustain large deformations. The resulting methodology is adapted
to simulate the dynamics of droplets in homogeneous and isotropic turbulence,
with the characteristic size of the droplet being smaller than the
characteristic Kolmogorov scale of the outer turbulent flow. We report on
statistical results for droplet deformation and orientation, collected from an
ensemble of turbulent trajectories, as well as comparisons with theoretical
models in the limit of small deformation.

##  1 Introduction

Understanding the dynamics of droplets in turbulence is fundamentally
important [1] and practically relevant for a broad spectrum of applications
including mixing and blending processes [2, 3, 4], petroleum industry [5, 6],
environmental remediation [7, 8], food industry [9, 10], and hemodynamic
engineering [11]. The characterization of droplet features such as elongation
and orientation in such scenarios is challenging: the flow inside and outside
the droplet couple in a non-linear way at the droplet interface, and this
coupling is further enriched by the complexity brought by time-dependent
strain rates exerted by the flow surrounding the droplet. Analytical
approaches are only possible in the limiting case of small deformations and/or
with limiting assumptions on the properties of the flow [12, 13, 14, 15]. The
use of numerical simulations is therefore of paramount importance to
understand the dynamics of such droplets [16, 17, 18, 19, 20, 21, 22, 23].

The dynamics of droplets in turbulent flows primarily depends on the ratio
between the droplet characteristic scale (e.g. the radius of the droplet at
rest, cf. Fig. 1) and the Kolmogorov dissipation length $\eta$, expressing the
scale at which the energy spectrum of the turbulent flow is cut off [24, 25,
26]. Droplets larger than the Kolmogorov length scale are subjected to
inertia-dominated hydrodynamic stresses, whereas viscous stresses dominate for
droplets smaller than the Kolmogorov scale (sub-Kolmogorov droplets). Here we
focus on the latter regime, where the flow deforming the droplet can be
described by the Stokes equations matched with the turbulent far field [27,
28] (cf. Sec. 2). This is a crucial assumption for having analytical reference
models to compare with, in the limit of small deformations.

The dynamics of sub-Kolmogorov droplets has been investigated with various
numerical approaches. The boundary integral approach was used to capture
critical sub-Kolmogorov droplet size and breakup rates in [27], while the
lattice Boltzmann (LB) method based on the diffuse-interface multiphase
approach was validated for sub-Kolmogorov investigations and employed to model
high deformations in Ref. [28]; statistical analysis on droplets dynamics was
performed in Refs. [29, 30, 28] based on phenomenological models for
ellipsoidal droplet dynamics [31]. In this work, we employ a novel route based
on the hybrid coupling between the immersed boundary (IB) method [32] and the
LB method [33, 34]. In contrast to multiphase/multicomponent LB approaches
[35, 36], the hybrid IBLB method could have several advantages: it inherently
preserves the sharp-interface limit of hydrodynamics without creating spurious
currents [37] and it allows to easily introduce various interfacial
properties, like the surface viscosity [38, 39] or generic constitutive laws
[40, 41, 42, 43, 44], while having the possibility of using finer meshes to
model them [45]. Furthermore, accurate modeling of the outer flow via
linearization of turbulent velocity gradients can be easily achieved with the
IBLB approach employing elementary LB schemes [46]. Overall, if from one side
the use of the IBLB method for the simulation of sub-Kolmogorov droplet
dynamics in turbulent flows seems appropriate, the applicability of the method
in this context has never been quantitatively assessed so far. This paper aims
to fill this gap.

The paper is organized as follows: in Sec. 2, we review the continuum
equations and control parameters of our system; in Sec. 3, we describe the
IBLB method used to simulate the dynamics of a single sub-Kolmogorov droplet
in homogeneous isotropic turbulence; simulation results and comparisons with
theoretical models in the limit of small deformation are presented in Sec. 4;
conclusions and a summary of our results will be presented in Sec. 5.

Figure 1: System and numerical setup description: trajectory and velocity
gradients $\nabla\bm{u}$ of a droplet are probed from direct numerical
simulations (DNS) [47] (panel (a)). The turbulent velocity gradients are then
imposed at the boundaries of the numerical domain used to perform immersed
boundary-lattice Boltzmann (IBLB) simulations (panel (b)). Relevant droplet
features are indicated in panel (c) (more details are given in the text).
Evolution in the IBLB simulation is sketched in panels (d-f) for three
different times: flow streamlines are reported with colormap following
velocity gradient magnitude (cf. panel (a)). A complete visualization of the
evolution of the droplet is shown in the Ancillary Files.

##  2  Problem Statement

We examine a system where a droplet with rest radius $\displaystyle R$ and
surface tension $\displaystyle\sigma$ is immersed in homogeneous isotropic
turbulence (cf. Fig. 1). The dynamic viscosity of the carrier fluid is
$\displaystyle\mu$, while $\displaystyle\chi\mu$ is the dynamic viscosity of
the droplet phase, so that $\displaystyle\chi$ denotes the system viscosity
ratio. In our analysis, we investigate the deformation of droplets in the
centered frame reference moving along a given trajectory. In this description,
we assume to be in the sub-Kolmogorov regime, i.e., $R/\eta\ll 1$, with
$\displaystyle\eta=(\nu^{3}/\epsilon)^{1/4}$ being the Kolmogorov length
scale, $\nu=\mu/\rho$ the kinematic viscosity of the turbulent flow with
$\rho$ the fluid density, $\tau_{\eta}=(\nu/\epsilon)^{1/2}$ the Kolmogorov
time scale and
$\displaystyle\epsilon=\nu\left\langle\nabla\bm{u}\cdot\nabla\bm{u}\right\rangle$
the turbulent dissipation rate [24, 25, 26], with
$\displaystyle\langle...\rangle$ indicating the ensemble mean. Within the sub-
Kolmogorov assumption, the droplet Reynolds number:

$\text{Re}_{\mbox{drop}}=\frac{R^{2}}{\nu\tau_{\eta}}=\left(\frac{R}{\eta}\right)^{2}\
,$ (1)

is inherently small, indicating the disparity in the inertial stress
$\Pi_{\rho}=\rho R^{2}/\tau_{\eta}^{2}$ and viscous stress
$\Pi_{\mu}=\mu/\tau_{\eta}$, with the latter dominating the overall
hydrodynamic stress. This is crucial to assume the dynamics is given by the
Stokes equations for an incompressible flow:

$\displaystyle\bm{\nabla}\cdot\bm{\Pi^{\prime}}$
$\displaystyle=\bm{\nabla}\cdot\bm{u}^{\prime}=0$ inside the droplet, (2a)
$\displaystyle\bm{\nabla}\cdot\bm{\Pi}$
$\displaystyle=\bm{\nabla}\cdot\bm{u}=0$ outside the droplet, (2b)

with:

$\displaystyle\bm{\Pi}^{\prime}=-p^{\prime}\bm{I}+\chi\mu\left(\bm{\nabla}\bm{u}^{\prime}+\bm{\nabla}\bm{u}^{\prime\dagger}\right)\
,$ (3a)
$\displaystyle\bm{\Pi}=-p\bm{I}+\mu\left(\bm{\nabla}\bm{u}+\bm{\nabla}\bm{u}^{\dagger}\right)\
,$ (3b)

being the inner and outer hydrodynamic stress, with $\bm{I}$ the 3D identity
matrix and $\displaystyle p^{\prime},p$ the inner and outer pressure,
respectively. Denoting with $\bm{x_{c}}$ the droplet center, the far-field
flow surrounding the droplet matches a linear expansion constructed from the
turbulent velocity gradients $\displaystyle\nabla\bm{u}(t)$:

$\bm{u}(\bm{x},t)=\left(\bm{x}-\bm{x_{c}}\right)\cdot\nabla\bm{u}(t).$ (4)

Denoting with $\tau_{\sigma}=\mu R/\sigma$ the droplet characteristic time set
by the balance between viscous effects and surface tension, the overall
droplet deformation state is controlled by the capillary number:

$\text{Ca}=\frac{\tau_{\sigma}}{\tau_{\eta}}\ .$ (5)

For small values of Ca, a small deformation regime is expected, wherein the
time scale for the droplet is much smaller than the characteristic time of the
turbulent flow, and the droplet tends to retract very soon to any deformation
imposed by the outer flow; on the other hand, for large values of Ca, the
droplet can be exposed for longer times to the straining effects of the outer
flow and non-linear effects emerge with larger deformations as well as
misalignment with the underlying flow [27]. In our analysis, the state of the
droplet is assessed by looking at both its elongation and orientation with
respect to the straining direction at all times (cf. Fig. 1). Information on
the deformation is retrieved from the largest eigenvalue
$\mathcal{\lambda}_{\mbox{max}}$ (corresponding to the eigenvector
$\bm{\hat{v}}_{\mbox{max}}$) of the morphological tensor $\bm{S}$, i.e., the
second order tensor describing the droplet shape [48, 14]:

$\bm{S}\bm{\hat{v}}_{\mbox{max}}=\lambda_{\mbox{max}}\bm{\hat{v}}_{\mbox{max}}\
.$ (6)

Assuming an ellipsoidal shape, the square root of
$\mathcal{\lambda}_{\mbox{max}}$ coincides with the major semi-axis of the
droplet (cf. Fig. 1). The instantaneous droplet orientation is studied via the
orientation parameter $\beta$ [27]:

$\displaystyle\beta=\frac{\bm{\hat{v}}_{\mbox{max}}\cdot\bm{E}\cdot\bm{\hat{v}}_{\mbox{max}}}{||\bm{E}||}\
,$ (7)

where $\bm{E}=\left(\nabla\bm{u}+\nabla\bm{u}^{\dagger}\right)/2$ is the
symmetric velocity gradient matrix representing the straining component of the
flow, with $||\bm{E}||$ being its Frobenius norm. The orientation parameter
$\displaystyle\beta$ weights the effects of the instantaneous strain
$\bm{E}(t)$ on the major elongation direction attained by the droplet at time
$t$: $\beta$ depends on the inclination angle $\theta$ between the droplet and
the straining direction of the flow (cf. Fig. 1), it is maximum when
$\theta=0$ (complete alignment), and minimum when $\theta=\pi/2$ (complete
misalignment). For example, in stationary shear flow $\beta=1/\sqrt{2}$ for a
droplet aligned with the strain direction and $\beta=-1/\sqrt{2}$ when the
droplet elongation is orthogonal to it.
In Sec. 4, we compare our numerical results with the linear theory in [49, 14]
adapted for turbulent gradients:

$\frac{\partial\bm{S}(t)}{\partial
t}=-\frac{1}{\tau_{\sigma}}\frac{40\left(\chi+1\right)}{\left(2\chi+3\right)\left(19\chi+16\right)}\left(\bm{S}(t)-\bm{I}R^{2}\right)+\\\
+2R^{2}\text{Ca}\frac{5}{2\chi+3}\bm{E}(t)\ .$ (8)

This equation is found by assuming an expansion of the morphological tensor
$\bm{S}$ in power series of Ca plugged into the Stokes equations (2) and
retaining only the linear contributions in Ca [50, 51]. We will use this
equation to get predictions for the droplet deformation and orientation when
$\text{Ca}\longrightarrow 0$ and validate the results of IBLB numerical
simulations in such a limit.

##  3 Numerical Method

###  3.1 Immersed boundary-lattice Boltzmann (IBLB) model

Various numerical methods have been developed for simulating bulk viscous
flows with deformable suspensions, such as boundary element methods [52, 53,
27, 54], volume-of-fluid methods [55], and lattice Boltzmann methods [46]. In
our work, we employ the immersed boundary-lattice Boltzmann (IBLB) method,
which has been extensively used in previous works for investigating
viscoelastic capsules [56, 42, 38, 39, 57, 58, 59, 60, 61, 62]. The method
combines the IB method, used to couple the sharp droplet interface with the
fluid, with LB method, which simulates bulk viscous flows. Few works
specifically focused on using such a method for simulating droplet dynamics in
stationary flows [38, 39, 63, 44]: here, we take a step further, adapting the
method for the simulation of droplet dynamics in a complex turbulent flow.

The LB method is a kinetic approach to simulate hydrodynamics, evolving
probability distribution functions in discrete directions. In this work, we
use the D3Q19 velocity scheme with a set of 19 discrete velocities
$\displaystyle\bm{c}_{i}$ [46, 34], which are indicated in Fig. 2, with
$\bm{c}_{0}$ related to the resting population. The LB method recovers both
the continuity and Navier-Stokes equations when fluctuations around the local
equilibrium are small [46, 34]. In the following, all quantities are given in
lattice Boltzmann units (LBU), in which we assume both the lattice and time
spacing to be unitary. In the chosen lattice scheme, the evolution of
probability distribution functions $\displaystyle f_{i}$ with discrete
velocities $\bm{c}_{i}$ at coordinate $\displaystyle\bm{x}$ and time
$\displaystyle t$ is regulated by the LB equation:

$f_{i}(\bm{x}+\bm{c}_{i},t+1)-f_{i}(\bm{x},t)=\left[\Omega_{i}(\bm{x},t)+S_{i}(\bm{x},t)\right]\
,$ (9)

including the streaming along direction $i$ on the l.h.s., the collision term
$\Omega_{i}$ and the force source $S_{i}$ on the r.h.s. Here, we implemented
the standard Bhatnagar-Gross-Krook (BGK) collision operator [64, 46]:

$\Omega_{i}(\bm{x},t)=-\frac{1}{\tau_{\mbox{f}}}[f_{i}(\bm{x},t)-f_{i}^{(\mbox{eq})}(\bm{x},t)]\
,$ (10)

stating that distribution functions relax towards equilibrium
$f_{i}^{(\mbox{eq})}(\bm{x},t)$ with characteristic time $\tau_{\mbox{f}}$.
The relaxation time is directly linked to the kinematic viscosity of the
system $\nu$ [33, 34]:

$\nu=c_{\mbox{s}}^{2}\left(\tau_{\mbox{f}}-\frac{1}{2}\right)\ ,$ (11)

where $c_{\mbox{s}}=1/\sqrt{3}$ is the speed of sound. The local equilibrium
dependency on $\bm{x}$ and $t$ is due to the density $\rho(\bm{x},t)$ and
velocity field $\bm{u}(\bm{x},t)$ [64]:

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

where $w_{i}=w(|\bm{c}_{i}|^{2})$ are suitable weights, such that $w_{0}=1/3$,
$w_{1-6}=1/18$, $w_{7-18}=1/36$. During the relaxation process, external
forces $\bm{F}_{{\mbox{ext}}}$ can be added in the LB equation via the source
term $S_{i}(\bm{x},t)$ according to the Guo scheme [65]:

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

The populations $\displaystyle f_{i}$ can then be used to retrieve the density
and momentum fields used to compute the equilibrium and source terms:

$\displaystyle\rho(\bm{x},t)=\sum_{i}f_{i}(\bm{x},t)\;,$ (14a)
$\displaystyle\rho\bm{u}(\bm{x},t)=\sum_{i}\bm{c}_{i}f_{i}(\bm{x},t)+\frac{\bm{F}_{{\mbox{ext}}}}{2}\
,$ (14b)

with the half-force correction appearing in the momentum [46, 65].
The interface of the droplet is represented by a finite set of Lagrangian
points forming a 3D triangular mesh with zero thickness. On each triangle, the
stress $\bm{\Pi}_{\text{drop}}$ given by the surface tension is:

$\bm{\Pi}_{\text{drop}}=\sigma\bm{I}_{2}\ ,$ (15)

where $\bm{I}_{2}$ is the 2D identity matrix. Details on the computation of
the force $\bm{\varphi}_{i}$ on the $i-$th Lagrangian point from the stress
tensor are given in [39].

The two-way coupling between Lagrangian points and fluid nodes is handled via
interpolation according to the IB method [32, 46]:

$\displaystyle\bm{F}_{{\mbox{ext}}}(\bm{x},t)=\sum_{i}\bm{\varphi}_{i}(t)\Delta(\bm{r}_{i}-\bm{x})\
,$ (16a) $\displaystyle\frac{\partial\bm{r}_{i}(t)}{\partial
t}=\sum_{\bm{x}}\bm{u}(\bm{x},t)\Delta(\bm{r}_{i}-\bm{x})\ ,$ (16b)

where the $\Delta$ function is the discretised Dirac delta and $\bm{r}_{i}(t)$
represents the coordinates of the $i$-th Lagrangian node. The $\Delta$
function can be factorised in three interpolation stencils
$\Delta(\bm{x})=\phi(x)\phi(y)\phi(z)$ which can be chosen and tweaked for
higher accuracy in the two-way coupling; here we use the 4-points
interpolation stencil [46]:

$\phi(x)=\begin{cases}\frac{1}{8}\left(3-2|x|+\sqrt{1+4|x|-4x^{2}}\right)&0\leq|x|<1\
,\\\ \frac{1}{8}\left(5-2|x|-\sqrt{-7+12|x|-4x^{2}}\right)&1\leq|x|<2\ ,\\\
0&2\leq|x|\;.\end{cases}$ (17)

Once the velocity of the Lagrangian points has been computed (cf. Eq. (16)),
the position $\bm{r}_{i}(t)$ is updated via an explicit forward Euler method:

$\bm{r}_{i}(t+1)=\bm{r}_{i}(t)+\frac{\partial\bm{r}_{i}(t)}{\partial t}\ .$
(18)

Figure 2: Representation of LB scheme applied for a specific edge (on the
left) and corner (on the right). Numbers refer to the LB discrete velocities
$\bm{c}_{i}$. Buried links are highlighted in red. Two examples of vectors
$\bm{n}^{(i)}$ normal to the planes are given on the top, for both the edge
and the corner.

Concerning the boundary conditions, the flux on the boundary domain is fixed
at each time $\displaystyle t$ by the turbulent gradients $\nabla\bm{u}(t)$
(cf. Eq. (4)): the linearization of the velocity fields is centered within the
droplet center of mass (cf. Fig. 1). In particular, we set:

$\displaystyle u_{x}=\frac{\partial u_{x}}{\partial
x}\left(l_{x}-\frac{L_{x}}{2}\right)+\frac{\partial u_{x}}{\partial
y}\left(l_{y}-\frac{L_{y}}{2}\right)+\frac{\partial u_{x}}{\partial
z}\left(l_{z}-\frac{L_{z}}{2}\right)\ ,$ (19a) $\displaystyle
u_{y}=\frac{\partial u_{y}}{\partial
x}\left(l_{x}-\frac{L_{x}}{2}\right)+\frac{\partial u_{y}}{\partial
y}\left(l_{y}-\frac{L_{y}}{2}\right)+\frac{\partial u_{y}}{\partial
z}\left(l_{z}-\frac{L_{z}}{2}\right)\ ,$ (19b) $\displaystyle
u_{z}=\frac{\partial u_{z}}{\partial
x}\left(l_{x}-\frac{L_{x}}{2}\right)+\frac{\partial u_{z}}{\partial
y}\left(l_{y}-\frac{L_{y}}{2}\right)+\frac{\partial u_{z}}{\partial
z}\left(l_{z}-\frac{L_{z}}{2}\right)\ ,$ (19c)

where $L_{x},L_{y}$ and $L_{z}$ represent the side lengths of the
computational domain, and $l_{i}\in[0,L_{i}]$ with $i=x,y,z$. Each boundary
plane is uniquely determined by fixing one of the three indices $l_{i}$ at a
value $0$ or $L_{i}$, while varying the other two $l_{j}$ in the range
$\left[0,L_{j}\right]$ (with $j\neq i$). This is required to impose a matching
of the system with the turbulent far-field, and the velocity in the bulk
domain is then naturally developed via the LB dynamics. When performing the
streaming step in the LB method (cf. Eq. (9)) on the nodes at the boundary,
special treatment is required for those populations pointing out of the
domain, whereas those in the opposite direction are somewhat undefined as
there is no lattice node where they would come from. To implement boundary
conditions for the outer flow, we adopt the Zou and He velocity boundary
conditions [66]. In particular, since the original proposed scheme is for the
D2Q9 scheme, we employ the 3D generalization for generic flow directions
presented in Ref. [67]. The implementation of such boundary conditions allows
to impose generic time dependent flows at the system boundaries regardless of
their instantaneous straining direction: this is achieved by supporting a
regular bounce back rule with a counter slip in the reflected populations at
the boundary nodes, acting as a transversal momentum correction. Among all the
populations at the boundary, we distinguish between three different cases: (a)
populations that are neither in the corners nor in the edges, (b) populations
that are in the edges but not in the corner, and (c) populations that are in
the corners. Let us introduce the tangential vectors
$\displaystyle\bm{t}_{i}=\bm{c}_{i}-\left(\bm{c}_{i}\cdot\bm{n}\right)\bm{n}$,
where $\bm{n}$ represents the normal vector to a boundary plane. We indicate
with $f_{i}$ populations corresponding to the $i$-th discrete velocity
$\bm{c}_{i}$ streaming outside of the domain, and with $f_{-i}$ those that get
reflected inside the domain, whose discrete velocity is $-\bm{c}_{i}$. For
case (a), the reflected populations $\displaystyle f_{-i}$ read [67]:

$f_{-i}=f_{i}-\left(2-|\bm{t}_{i}|\right)\frac{\rho}{6}\left(\bm{c}_{i}\cdot\bm{u}\right)-\frac{\rho}{3}\left(\bm{t}_{i}\cdot\bm{u}\right)+\\\
+\frac{1}{2}\sum_{j=0}^{18}f_{j}\left(\bm{t}_{i}\cdot\bm{c}_{j}\right)\left(1-|\bm{c}_{j}\cdot\bm{n}|\right)\
.$ (20)

As can be seen for those populations pointing perpendicular w.r.t. the
boundary plane (i.e., $\bm{t}_{i}=\bm{0}$), we remark that the factor
$2-|\bm{t}_{i}|$ must be included in the second term of the r.h.s. to retrieve
the standard bounce-back rule [46]:

$f_{-i}=f_{i}-\frac{\rho}{3}\left(\bm{c}_{i}\cdot\bm{u}\right)\ .$ (21)

For the case (b), two different planes with normal vectors $\bm{n}^{(1)}$ and
$\bm{n}^{(2)}$ must be taken in consideration (cf. Fig. 2), with the reflected
populations that read:

$f_{-i}=f_{i}-\frac{1}{4}\sum_{j=0}^{18}f_{j}\left(\bm{t}_{i}\cdot\bm{c}_{j}\right)\\\
\cdot\left(1-|\bm{c}_{j}\cdot\bm{n}^{(1)}|\right)\left(1-|\bm{c}_{j}\cdot\bm{n}^{(2)}|\right)\
.$ (22)

Among these populations, we define $f^{(1,2)}$ those associated with the two
lattice vectors $\bm{c}^{(1,2)}=\pm\left(\bm{n}^{(1)}-\bm{n}^{(2)}\right)$.
The latter constitute the so-called “buried links”, i.e., those lattice
vectors along which both $f_{i}$ and $f_{-i}$ points outside the domain (see
red arrows in Fig. 2). For such populations, Eq. (22) does not hold, and an
additional prescription is required [67]:

$f^{(1,2)}=\frac{1}{22}\sum_{i=1}^{18}f_{i}\Bigg{\\{}1-\left(1-\left|\bm{c}_{i}\cdot\left[\bm{n}^{(1)}\times\bm{n}^{(2)}\right]\right|\right)\\\
\cdot\left(1-\left|\bm{c}_{i}\cdot\left[\frac{\bm{n}^{(1)}+\bm{n}^{(2)}}{|\bm{c}_{i}|^{2}}\right]\right|\right)\Bigg{\\}}\
.$ (23)

On each edge, there is exactly one population $f^{(1)}$ streaming along
$\bm{c}^{(1)}$ with the corresponding opposite one $f^{(2)}$ streaming along
$\bm{c}^{(2)}$. As an example, in Fig. 2, we consider the edge in $y=0$ and
$z=L_{z}$: the two normals are $\bm{n}^{(1)}=(0,0,-1)$ and
$\bm{n}^{(2)}=(0,1,0)$. The buried links are therefore identified by
$\bm{c}^{(1,2)}=\pm(0,-1,-1)$, that correspond to the directions $\bm{c}_{11}$
and $\bm{c}_{13}$ (red arrows in Fig. 2). The resting population $f_{0}$ must
be corrected as $\displaystyle f_{0}=12f^{(1,2)}$ [67]. Lastly, we consider
case (c). For each corner, we have three different planes with normals
$\displaystyle\bm{n}^{(1)},\bm{n}^{(2)}$ and $\bm{n}^{(3)}$ which identify 6
buried links (cf. Fig. 2) whose corresponding lattice vectors are
$\bm{c}^{(1,\dots,6)}$. The buried vectors are such that
$\bm{c}^{(i)}\cdot(\bm{n}^{(1)}+\bm{n}^{(2)}+\bm{n}^{(3)})=0$. The
corresponding distributions assigned to these six buried vectors are:

$f^{(1,\dots,6)}=\frac{1}{18}\sum_{i=1}^{18}f_{i}\left|\bm{c}_{i}\cdot\frac{\bm{n}^{(1)}+\bm{n}^{(2)}+\bm{n}^{(3)}}{|\bm{c}_{i}|^{2}}\right|\
.$ (24)

Again, the resting population needs a special consideration, and it is set to
$f_{0}=12f^{(1,\dots,6)}$ [67].

###  3.2 Outer Flow and matching between DNSU and LBU

We consider the dynamics of sub-Kolmogorov droplets with $\displaystyle\chi=1$
in an ensemble of turbulent trajectories taken from the TURB-Lagr open-source
database [47]. This database contains over $3\times 10^{5}$ particle
trajectories obtained for the case of homogeneous isotropic turbulence via
direct numerical simulations (DNS) of the incompressible Navier-Stokes
equations using a pseudo-spectral approach [68]. The Reynolds number
calculated based on the Taylor scale is $\displaystyle\text{Re}\sim 130$. In
what we refer to as direct numerical simulation units (DNSU), the Kolmogorov
length scale is $\eta=4.2\times 10^{-3}\ \mbox{DNSU}$, and the Kolmogorov time
scale is $\displaystyle\tau_{\eta}=2.3\times 10^{-2}\ \mbox{DNSU}$. The
kinematic viscosity of the system is $\displaystyle\nu=8\times 10^{-4}\
\mbox{DNSU}$, and the dissipation rate is $\epsilon=1.4\ \mbox{DNSU}$. Each
trajectory lasts for roughly an eddy turnover time
$\displaystyle\tau_{\mbox{eddy}}=4.5\ \mbox{DNSU}\sim 200\tau_{\eta}$ so that,
by subjecting droplets to a linearization of the turbulent field as in Eq.
(4), all the relevant turbulent features are probed along each trajectory. The
local velocity gradients are probed every $\displaystyle\Delta t=2.25\times
10^{-3}\ \mbox{DNSU}$ for a total of $\displaystyle 2\times 10^{3}$ turbulent
realizations for each trajectory.

In order to match the turbulent outer flow with the IBLB simulation, we need
to apply a unit conversion from the quantities of the DNS to those used in the
IBLB simulations. Hereafter, a generic quantity $A$ is assumed to be in DNSU,
whereas the corresponding quantity in LB units (LBU) is represented as
${A}_{\text{{LBU}}}$. Therefore, the conversion factor $C_{A}$ in DNSU is such
that $A=A_{\text{{LBU}}}\ C_{A}$. We start by fixing the conversion factor for
the length $C_{x}$. By fixing a desired value of $R/\eta$, and knowing the
Kolmogorov lengthscale $\eta$, the radius of the droplet is known in DNSU and
the conversion factor $C_{x}$ can be computed as:

$C_{x}=\frac{R}{R_{\text{{LBU}}}}\ .$ (25)

Then, we consider the unit conversion for time, which can be retrieved from
the viscosity. Indeed, since $\nu=\nu_{\text{{LBU}}}C_{x}^{2}/C_{t}$, we have:

$C_{t}=C_{x}^{2}\frac{\nu_{\text{{LBU}}}}{\nu}\ .$ (26)

Within this prescription, at every time frame $\Delta t$, we read the velocity
gradient in DNSU and convert it in LBU as:

$\nabla\bm{u}_{\text{{LBU}}}=C_{t}\nabla\bm{u}\ .$ (27)

For each realization of the velocity gradient, we then perform $N$ IBLB time
steps given by:

$N=\frac{\Delta t}{C_{t}}\ .$ (28)

Note that, by using Eq. (25) and Eq. (26), we obtain that the number of IBLB
steps needed to match the dynamics of the outer flow scales with the inverse
square of $R/\eta$:

$N=\Delta
t\frac{\nu}{\nu_{\text{{LBU}}}}\left(\frac{R_{\text{{LBU}}}}{R}\right)^{2}\sim\left(\frac{\eta}{R}\right)^{2}\
.$ (29)

For a given realization of $\displaystyle\nabla\bm{u}(t)$, convergence to the
Stokes dynamics given by Eq. (2) is expected when $R/\eta\rightarrow 0$. As
discussed in Sec. 4, a good convergence is already observed when
$R/\eta\approx 1/2$, corresponding to $N\approx 500$. We remark that by
increasing the LB relaxation time $\displaystyle\tau_{\mbox{f}}$ (and
consequently $\displaystyle\nu_{\mbox{LBU}}$), the number of steps $N$
required for LB to reproduce sub-Kolmogorov droplets dynamics can be reduced
without altering the simulation results, albeit attention should be used to
not induce over relaxation [46].
Once the unit conversion between the two systems has been set, one only needs
to compute the surface tension of the drop $\sigma$ in LBU. Since we know the
value of $\tau_{\eta}$ in DNSU, we can compute $\tau_{\sigma}$ in DNSU for a
given capillary number by using its definition given in Eq. (5), and convert
it in LBU:

$\tau_{\sigma,\text{{LBU}}}=\frac{\tau_{\sigma}}{C_{t}}\ .$ (30)

The value of the surface tension in LBU is then computed as:

$\sigma_{\text{{LBU}}}=\frac{\mu_{\text{{LBU}}}R\text{{LBU}}}{\tau_{\sigma,\text{{LBU}}}}\
.\\\ $ (31)

Once the parameters $R/\eta,\text{Ca},\tau_{\mbox{f}}$, $R_{\text{{LBU}}}$ and
the computational domain size $L_{x}\times L_{y}\times L_{z}$ are set, one
needs to apply the following steps:

  1. 1.

Consider the turbulent gradient matrix ${\nabla\bm{u}}(t)$ at time $t$;

  2. 2.

Apply the conversion in Eq. (27) to initialize the boundary planes via Eq.
(19);

  3. 3.

Perform $N$ IBLB steps with the prescriptions in Sec. 3 employing the surface
tension given in Eq. (31);

  4. 4.

Consider the turbulent gradient matrix ${\nabla\bm{u}}(t+\Delta t)$ at time
$t+\Delta t$ and re-iterate steps 2-4.

In our IBLB setup, we employed a 3D box with size $\displaystyle
L_{x}=L_{y}=L_{z}=200\ \mbox{LBU}$ for simulating droplets with rest radius
$\displaystyle R_{\text{{LBU}}}=19\ \mbox{LBU}$ possessing a structured mesh
composed of 20480 triangles: these numerical parameters where chosen after
performing LB simulations to test the convergence of the LB to the
hydrodynamic equations [34, 46]. Furthermore, the chosen parameters need to be
tailored to minimize confinement effects [69, 42] and discrete curvature
effects [38, 42], and a very refined mesh is required for the latter. The
value for the LB relaxation time is $\tau_{\mbox{f}}=1.5$ LBU. Simulations ran
on Nvidia Ampere “A100” Graphic Processing Units (GPUs), with a simulation
time for each trajectory lasting roughly 3 hours.

Figure 3: Sketch of the mesh regularization employed. Panel (a) shows the same
trajectory for a not-regularized (solid line) and regularized droplet (dashed
line); panels (b) and (c) represent the Lagrangian mesh right before the time
where the not-regularized simulation destabilizes. Panel (d) sketches the
regularization process for a specific node and its neighbours (detailed
equations are given in the text).

###  3.3 Mesh regularization

The IBLB method, while inherently preserving the sharp-interface limit of
hydrodynamics, presents challenges in tracking the evolution of mesh elements,
especially in the presence of irregular distribution of surface nodes [70,
71]. As shown in [72], droplets subjected to a prolonged shearing strain rate
can incur an accumulation of mesh nodes and very loose mesh elements stretched
in the major strain direction. Both issues must be avoided when using the IBLB
method to couple fluid and interface: indeed, an extreme clustering of the
mesh nodes can cause them to overlap, whereas triangles that get too loose can
compromise the two-way coupling with the inner/outer fluid penetrating the
mesh. In Fig. 3, panel (a), we report the largest normalized eigenvalue
$\displaystyle\lambda_{\mbox{max}}/R^{2}$ as a function of the normalized time
$t/\tau_{\mbox{eddy}}$. The solid red line corresponds to an IBLB simulation
without mesh regularization, while the dashed blue line represents one that
has been regularized using the algorithm described in this section (details
provided below). In panel (b), we report the triangular mesh of the drop
without regularization right before the simulation became unstable, which
shows both loosening and clustering of triangles. In panel (c), the
regularized mesh is reported. Concerning the algorithms employed to obtain
regular meshes, one can ideally split them into two big families: the first
one involves changes in the mesh topology (e.g., creation/destruction of
nodes/triangles); the second one employees smoothing techniques, i.e., without
any change in the mesh topology, its vertices are moved in order to obtain a
regular mesh. In the specific context of IBLB simulations of drops, a hybrid
version has been implemented in Ref. [72]: by using an iterative algorithm,
the mesh nodes are moved until reaching a given tolerance; then, a “flipping
algorithm” causing a change in the connections of the mesh nodes is applied.
We note that topology changes may generate quite delicate issues since the
topology of the drop at rest is essential for the computation of the nodal
forces starting from the stress on each triangle (cf. Eq. (15)) [39, 38, 73].
Hence, we preferred not to change the topology but instead use the extension
of the Laplacian smoothing technique proposed by Ref. [74]. This technique
acts globally on each node with an iterative process: after having updated all
the mesh nodes via the IBLB step (cf. Eq. (18)), we end up with a set of mesh
coordinates $\\{\bm{x}_{o}\\}$. Before starting the Laplacian smoothing, we
initialize a set of coordinates $\\{\bm{x}_{p}\\}=\\{\bm{x}_{o}\\}$. We then
apply the following iterative algorithm:

  1. 1.

we initialise the auxiliary set of coordinates
$\\{\bm{x}_{q}\\}=\\{\bm{x}_{p}\\}$. This set is necessary to perform a
simultaneous update and represents the mesh coordinates at the previous step
of this iterative algorithm.

  2. 2.

We update each node in the set $\\{\bm{x}_{p}\\}$:

$\bm{x}_{p}^{i}=\alpha\bm{x}_{o}^{i}+\frac{1-\alpha}{N_{i}}\sum_{j\in{\mbox{neigh}}(i)}\bm{x}_{q}^{j}\
,$ (32)

with $\displaystyle\alpha\in[0,1]$ being a user-defined parameter regulating
the influence of the starting position of the node [74]. The sum runs over the
$N_{i}$ neighbors of the $i$-th node. This step represents the classic
Laplacian smoothing [74].

  3. 3.

Once all the coordinates $\bm{x}_{p}^{i}$ have been updated, one ends up with
some nodes that might have been moved inward or outward w.r.t. the coordinates
at the previous iteration $\bm{x}_{q}^{i}$ (cf. Fig. 3, panel (d)). This may
cause the creation of artificial valleys and ridges, thus undermining the
local curvature of the droplet and the overall volume conservation [74]. In
order to preserve curvature regularity, we employ the correction to the
classic Laplacian smoothing introduced in Ref. [74] (cf. Fig. 3, panel (d)).
We therefore correct the coordinates $\bm{x}_{p}^{i}$ in the following way:

$\qquad\qquad\bm{x}_{p}^{i}\mapsto\bm{x}_{p}^{i}-\left(\gamma\bm{x}_{b}^{i}-\frac{1-\gamma}{N_{i}}\sum_{j\in{\mbox{neigh}}(i)}\bm{x}_{b}^{j}\right)\
,$ (33)

where

${\bm{x}_{b}^{i}}=\bm{x}_{p}^{i}-\left[\alpha\bm{x}_{o}^{i}+(1-\alpha)\bm{x}_{q}^{i}\right]\
.$ (34)

The parameter $\displaystyle\gamma\in[0,1]$ is an additional control
parameter.

  4. 4.

We finally need a condition to exit the iterative algorithm. For each node
$\bm{x}_{p}^{i}$ in the set $\\{\bm{x}_{p}\\}$, we compute the difference
w.r.t the set of coordinates at the previous iteration $\\{\bm{x}_{q}\\}$. We
exit the loop if the $\max_{i}|\bm{x}_{p}^{i}-\bm{x}_{q}^{i}|$ is less than a
certain threshold. Otherwise, we repeat the present algorithm from step 1.

In our simulations, we find $\displaystyle\alpha=0.1,\gamma=0.51$ to be a
suitable set of parameters, although other options can be taken into
consideration [75]. Concerning the condition to break the iterative algorithm,
the threshold has been set to 0.01 LBU.

Figure 4: Largest normalized eigenvalue
$\displaystyle\lambda_{\mbox{max}}/R^{2}$ (panel (a)) and orientation
parameter $\displaystyle\beta$ (panel (b)) for a representative trajectory as
a function of time for $\text{Ca}=0.05$ and different values of $R/\eta$. Time
is made dimensionless with respect to the eddy turnover time
$\displaystyle\tau_{\mbox{eddy}}$ of the outer turbulent flow. Linear theory
results obtained from the Stokes equations (cf. Eq. (8)) are also reported.
Arrows indicate an increase in $R/\eta$.

##  4 Results

We start our numerical investigations by verifying the convergence towards the
Stokes equations (2) at changing the relative importance of the droplet radius
$R$ with respect to the Kolmogorov scale $\eta$. To this aim, we perform
numerical simulations for a fixed $\text{Ca}=0.05$ at changing the ratio
$R/\eta$. Results for the largest eigenvalue
$\displaystyle\lambda_{\mbox{max}}$ (normalized with $\displaystyle R^{2}$)
and the orientation parameter $\displaystyle\beta$ for a representative
trajectory are reported in Fig. 4. The value of Ca is chosen small enough to
ensure that the linear theory (cf. Eq. (8)) holds when $R/\eta\rightarrow 0$,
making it possible to have analytical predictions to compare with the IBLB
numerical simulations. Results for linear theory are obtained by integrating
Eq. (8) with a standard explicit Runge-Kutta method [76]. By decreasing
$R/\eta$, a very good convergence of the IBLB numerical simulations to the
linear theory predictions is observed. In the simulated time range,
$\displaystyle\lambda_{\mbox{max}}$ is shown to be more sensible to variations
in $\displaystyle R/\eta$, with an oscillating behavior emerging when
increasing $\displaystyle R/\eta$ around unity. These oscillations emerge
alongside both enhanced and delayed elongations, as well as faster
retractions, indicating an inability to adapt to the linearized Stokes
solutions when $R/\eta\geq 1$. Milder oscillations are shown for the
orientation parameter $\beta$, where departure from the linear theory is less
evident at increasing $\displaystyle R/\eta$. By further decreasing
$\displaystyle R/\eta$ the computational cost increases (cf. Eq. (29)) but the
convergence to the linear theory results does not improve sensibly (see the
comparison between the reported case $\displaystyle R/\eta=0.472$ and
$\displaystyle R/\eta=0.15$ in Fig. 4).

Figure 5: Largest normalized eigenvalue
$\displaystyle\lambda_{\mbox{max}}/R^{2}$ and orientation parameter
$\displaystyle\beta$ for a representative trajectory as a function of time for
$R/\eta=0.472$ at changing Ca: $\text{Ca}=0.05$ (panels (a-b)),
$\text{Ca}=0.15$ (panels (c-d)) and $\text{Ca}=0.25$ (panels (e-f)). Time is
made dimensionless with respect to the eddy turnover time
$\displaystyle\tau_{\mbox{eddy}}$ of the outer turbulent flow. Linear theory
results obtained from the Stokes equations (cf. Eq. (8)) are also reported.
For a complete visualization of the IBLB simulation for the higher Ca
reported, see Ancillary Files. Figure 6: PDF (probability density function)
for the largest normalized eigenvalue $\lambda_{\mbox{max}}/R^{2}$ and the
orientation parameter $\displaystyle\beta$, with $R/\eta=0.472$ and
$\text{Ca}=0.05,0.15,0.25$. Both lin-lin plots (panels (a, c)) and log-log
plots (panels (b, d)) are reported.

In Fig. 5, we complement the findings reported in Fig. 4 by comparing results
of IBLB numerical simulations and the predictions of linear theory for
$R/\eta=0.472$ at changing Ca. For growing values of Ca, agreement with the
linear theory valid for $\text{Ca}\ll 1$ is not achieved, as expected. The
value of $\displaystyle\lambda_{\mbox{max}}$ obtained from IBLB exceeds the
linear theory with enhanced peaks. Notice that a similar behavior is found
when studying the stationary deformation at finite Ca [77, 78, 79, 55, 54,
63]. Large variations in $\displaystyle\lambda_{\mbox{max}}$ observed in
panels (c) and (e) are clearly related to the intermittency in the strain
rates, and they become more evident and more persistent in time at increasing
Ca. Good agreement in $\beta$ is shown between IBLB numerical simulations and
predictions of linear theory up to $\text{Ca}=0.15$, with an overall positive
value of $\beta$, pointing to the alignment of the droplet with the straining
direction of the outer flow. When Ca is increasing, persistence of large
variations in $\displaystyle\lambda_{\mbox{max}}$ (panel (e)) goes together
with droplet misalignment with respect to the straining direction of the flow
[80]: for such events, e.g., the one shown at around
$t/\tau_{\mbox{eddy}}\approx 0.15$ in panel (f), the value of
$\displaystyle\beta$ is shown to locally reach negative values, suggesting the
inability of the droplet to adapt with the fast varying strain rates. A 3D
visualization of the IBLB simulation reported in panels (e-f) is provided as
Ancillary Files.

To delve deeper into the analysis, we performed a statistical characterization
of the PDF (probability density function) of both
$\displaystyle\lambda_{\mbox{max}}$ and $\displaystyle\beta$. To this aim, we
considered 1000 independent turbulent trajectories [47] and performed IBLB
numerical simulations for the values of $R/\eta=0.472$ and
$\text{Ca}=0.05,0.15,0.25$. Along the same trajectories, we also integrated
the linear theory derived from the Stokes equations (cf. Eq. (8)) and
extracted the corresponding statistics. The two statistics are compared in
Fig. 6. Data are displayed on both a linear scale (panels (a,c)) and a
logarithmic scale (panels (b,d)). In a linear scale (panel (a)), the PDF of
$\displaystyle\lambda_{\mbox{max}}/R^{2}$ shows a good agreement between IBLB
results and the results of the linear theory for $\text{Ca}\rightarrow 0$; at
increasing Ca, the IBLB PDF departs from linear theory predictions. In panel
(b), the tails of the PDF of $\displaystyle\lambda_{\mbox{max}}/R^{2}$ are
analyzed, showing that the IBLB PDF tail overestimates the linear theory
prediction, and this overestimation grows with Ca. More statistics are
probably needed for a more quantitative assessment of the distribution tails.
In panel (c), we report the PDF of $\displaystyle\beta$, and we observe that
the IBLB simulation results well converge to the linear theory predictions
when Ca is small. The same holds for the PDF tails in panel (d). We notice
that within the framework of linear theory, results for deformation are more
sensitive to a variation in Ca (panels (a,b)) in comparison to the results for
orientation (panels (c,d)). This is not unexpected if one thinks of the linear
theory stationary results for deformation and orientation in simple flows
[13], where the deformation changes more than the orientation angle between
the droplet and the flow direction. Overall, Fig. 6 well supports the
statement that results from IBLB simulations match the linear theory
predictions in the limit $\text{Ca}\rightarrow 0$.

##  5 Conclusions

We have investigated the application of the hybrid immersed boundary-lattice
Boltzmann (IBLB) method for the simulation of the dynamics of a droplet
subjected to a time-dependent velocity gradient. We have specialized the
method to simulate the dynamics of a droplet in a turbulent flow, with the
size of the droplet being smaller than the Kolmogorov scale of the turbulent
flow. The proposed methodology hinges on two essential ingredients: first, the
implementation of a generalized Zou-He scheme (cf. Sec. 3.1) to accommodate a
generic strain matrix at the boundaries of the IBLB simulation domain; second,
an extended Laplacian smoothing technique (cf. Sec. 3.3) to avoid
deterioration of the triangular mesh resulting in numerical instabilities. To
verify the correctness of the method, we compared the results of numerical
simulations with the results of the linear theory that can be analytically
derived from the Stokes equations in the limit of small droplet deformations
(i.e., small capillary numbers Ca). Specifically, we have shown how the
statistics of droplet deformation and orientation estimated over thousands of
eddy turn-over times of the turbulent flow well match the results of the
theory in the limit of small deformations.

Various future perspectives will open up at this stage. In this paper, we
mainly focused on sub-Kolmogorov droplet dynamics and relatively small values
of Ca. The assumption of sub-Kolmogorov droplet dynamics is crucial to have
some reference theory to compare with, but in principle the LB method can be
pushed towards situations where inertial contributions appear. Moreover,
simulations with larger Ca are possible, but larger deformations and droplet
breakup need to be suitably accounted for within the numerical simulations
[27, 29, 30] to allow a detailed characterization of breakup statistics. We
also notice that we focused on the dynamics of simple droplets with surface
tension at the interface. The IBLB method is flexible enough to allow the
adaptation of additional interfacial complexities to model soft suspensions
like vesicles or capsules, realistically employed in industrial and medical
devices [60, 81, 82, 83]. This potentiality could shed light on the dynamics
of such soft suspensions in time-dependent flows. Finally, we also remark that
the deformation of soft suspensions in complex flows has also been
investigated with phenomenological models [31, 84, 85, 86] based on the idea
that the suspended particle retains an ellipsoidal shape at all times. Earlier
studies used these models also to investigate the dynamics of droplets in
turbulent flows [29, 87, 30]: while the predictions in the small deformation
limit coincide with the linear theory that we have used in this paper,
deviations are expected to emerge at large Ca. It could then be of interest to
use the IBLB method to design improved phenomenological models better working
at large Ca, in the same spirit of Ref. [63] for droplets in simple shear
flow.

### Acknowledgments

We wish to acknowledge Fabio Bonaccorso for support. We also acknowledge
useful discussion with Luca Biferale, Michele Buzzicotti and Chiara
Calascibetta. 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). This work was supported by the Italian
Ministry of University and Research (MUR) under the FARE programme (No.
R2045J8XAW), project “Smart-HEART”. MS acknowledges the support of the
National Center for HPC, Big Data and Quantum Computing, Project CN_00000013 -
CUP E83C22003230001, Mission 4 Component 2 Investment 1.4, funded by the
European Union - NextGenerationEU.

## References

  * [1] F. Mashayek and R.V.R. Pandya  “Analytical description of particle/droplet-laden turbulent flows”  In _Progress in Energy and Combustion Science_ 29.4 Elsevier BV, 2003, pp. 329–378  DOI: 10.1016/s0360-1285(03)00029-7
  * [2] Albert Sanfeld and Annie Steinchen  “Emulsions stability, from dilute to dense emulsions — Role of drops deformation”  In _Advances in Colloid and Interface Science_ 140.1 Elsevier BV, 2008, pp. 1–65  DOI: 10.1016/j.cis.2007.12.005
  * [3] K. Tang and A. Gomez  “Generation by electrospray of monodisperse water droplets for targeted drug delivery by inhalation”  In _Journal of Aerosol Science_ 25.6 Elsevier BV, 1994, pp. 1237–1249  DOI: 10.1016/0021-8502(94)90212-7
  * [4] M.. Yukuyama, D… Ghisleni, T… Pinto and N.. Bou-Chacra  “Nanoemulsion: process selection and application in cosmetics – a review”  In _International Journal of Cosmetic Science_ 38.1 Wiley, 2015, pp. 13–24  DOI: 10.1111/ics.12260
  * [5] E.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 BV, 2005, pp. 2101–2113  DOI: 10.1016/j.ces.2004.12.003
  * [6] Narendra Kumar, Amit Verma and Ajay Mandal  “Formation, characteristics and oil industry applications of nanoemulsions: A review”  In _Journal of Petroleum Science and Engineering_ 206 Elsevier BV, 2021, pp. 109042  DOI: 10.1016/j.petrol.2021.109042
  * [7] Ming Li and Chris Garrett  “The relationship between oil droplet size and upper ocean turbulence”  In _Marine Pollution Bulletin_ 36.12 Elsevier BV, 1998, pp. 961–970  DOI: 10.1016/s0025-326x(98)00096-4
  * [8] Luky Hendraningrat and Ole Torsæter  “Metal oxide-based nanoparticles: revealing their potential to enhance oil recovery in different wettability systems”  In _Applied Nanoscience_ 5.2 Springer ScienceBusiness Media LLC, 2014, pp. 181–199  DOI: 10.1007/s13204-014-0305-6
  * [9] Gipsy Tabilo-Munizaga and Gustavo V. Barbosa-Cánovas  “Rheology for the food industry”  In _Journal of Food Engineering_ 67.1–2 Elsevier BV, 2005, pp. 147–156  DOI: 10.1016/j.jfoodeng.2004.05.062
  * [10] Laura Lee and Ian T. Norton  “Comparing droplet breakup for a high-pressure valve homogeniser and a Microfluidizer for the potential production of food-grade nanoemulsions”  In _Journal of Food Engineering_ 114.2 Elsevier BV, 2013, pp. 158–163  DOI: 10.1016/j.jfoodeng.2012.08.009
  * [11] Oguz K. Baskurt and Herbert J. Meiselman  “Blood rheology and hemodynamics”  In _Seminars in Thrombosis and Hemostasis_ 29.5 Georg Thieme Verlag KG, 2003, pp. 435–450  DOI: 10.1055/s-2003-44551
  * [12] Geoffrey Ingram Taylor  “The viscosity of a fluid containing small drops of another fluid”  In _Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character_ 138.834 The Royal Society London, 1932, pp. 41–48  DOI: 10.1098/rspa.1932.0169
  * [13] Geoffrey Ingram Taylor  “The formation of emulsions in definable fields of flow”  In _Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character_ 146.858 The Royal Society London, 1934, pp. 501–523  DOI: 10.1098/rspa.1934.0169
  * [14] 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  DOI: 10.1017/s0022112080000316
  * [15] Kausik Sarkar and William R. Schowalter  “Deformation of a two-dimensional viscous drop in time-periodic extensional flows: analytical treatment”  In _Journal of Fluid Mechanics_ 436 Cambridge University Press (CUP), 2001, pp. 207–230  DOI: 10.1017/s0022112001004013
  * [16] F. Milan, M. Sbragaglia, L. Biferale and F. Toschi  “Lattice Boltzmann simulations of droplet dynamics in time-dependent flows”  In _The European Physical Journal E_ 41.6 Springer ScienceBusiness Media LLC, 2018, pp. 6  DOI: 10.1140/epje/i2018-11613-0
  * [17] Daniel Fuster et al.  “Numerical simulation of droplets, bubbles and waves: state of the art”  In _Fluid Dynamics Research_ 41.6 IOP Publishing, 2009, pp. 065001  DOI: 10.1088/0169-5983/41/6/065001
  * [18] 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  DOI: 10.1146/annurev-fluid-010518-040401
  * [19] Alessio Roccon, Marco De Paoli, Francesco Zonta and Alfredo Soldati  “Viscosity-modulated breakup and coalescence of large drops in bounded turbulence”  In _Physical Review Fluids_ 2.8 American Physical Society (APS), 2017, pp. 083603  DOI: 10.1103/physrevfluids.2.083603
  * [20] Siddhartha Mukherjee et al.  “Droplet–turbulence interactions and quasi-equilibrium dynamics in turbulent emulsions”  In _Journal of Fluid Mechanics_ 878 Cambridge University Press (CUP), 2019, pp. 221–276  DOI: 10.1017/jfm.2019.654
  * [21] Marco Crialesi-Esposito, Marco Edoardo Rosti, Sergio Chibbaro and Luca Brandt  “Modulation of homogeneous and isotropic turbulence in emulsions”  In _Journal of Fluid Mechanics_ 940 Cambridge University Press, 2022, pp. A19  DOI: 10.1017/jfm.2022.179
  * [22] Marco Crialesi-Esposito, Sergio Chibbaro and Luca Brandt  “The interaction of droplet dynamics and turbulence cascade”  In _Communications Physics_ 6.1 Nature Publishing Group UK London, 2023, pp. 5  DOI: 10.1038/s42005-022-01122-8
  * [23] Camille Deberne, Victor Chéron, Alexandre Poux and Jorge César Brändle de Motta  “Breakup prediction of oscillating droplets under turbulent flow”  In _International Journal of Multiphase Flow_ 173 Elsevier BV, 2024, pp. 104731  DOI: 10.1016/j.ijmultiphaseflow.2024.104731
  * [24] Uriel Frisch  “Turbulence: The Legacy of A.N. Kolmogorov”  Cambridge University Press, 1995  DOI: 10.1017/cbo9781139170666
  * [25] Pierre Sagaut and Claude Cambon  “Homogeneous Turbulence Dynamics”  Cambridge University Press, 2008  DOI: 10.1017/cbo9780511546099
  * [26] Roberto Benzi and Federico Toschi  “Lectures on turbulence”  In _Physics Reports_ 1021 Elsevier BV, 2023, pp. 1–106  DOI: 10.1016/j.physrep.2023.05.001
  * [27] Vittorio Cristini, J. Bławzdziewicz, Michael Loewenberg and Lance R. Collins  “Breakup in stochastic Stokes flows: sub-Kolmogorov drops in isotropic turbulence”  In _Journal of Fluid Mechanics_ 492 Cambridge University Press (CUP), 2003, pp. 231–250  DOI: 10.1017/s0022112003005561
  * [28] Felix Milan, Luca Biferale, Mauro Sbragaglia and Federico Toschi  “Sub-Kolmogorov droplet dynamics in isotropic turbulence using a multiscale lattice Boltzmann scheme”  In _Journal of Computational Science_ 45 Elsevier, 2020, pp. 101178  DOI: 10.1016/j.jocs.2020.101178
  * [29] 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  DOI: 10.1017/jfm.2014.366
  * [30] Samriddhi Sankar Ray and Dario Vincenzi  “Droplets in isotropic turbulence: deformation and breakup statistics”  In _Journal of Fluid Mechanics_ 852 Cambridge University Press (CUP), 2018, pp. 313–328  DOI: 10.1017/jfm.2018.453
  * [31] 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  DOI: 10.1016/s0377-0257(98)00065-2
  * [32] Charles S Peskin  “The immersed boundary method”  In _Acta numerica_ 11 Cambridge University Press, 2002, pp. 479–517  DOI: 10.1017/s0962492902000077
  * [33] Roberto Benzi, Sauro Succi and Massimo Vergassola  “The lattice Boltzmann equation: theory and applications”  In _Physics Reports_ 222.3 Elsevier, 1992, pp. 145–197  DOI: 10.1016/0370-1573(92)90090-m
  * [34] Sauro Succi  “The lattice Boltzmann equation: for fluid dynamics and beyond”  Oxford university press, 2001  URL: http://global.oup.com/academic/product/the-lattice-boltzmann-equation-9780198503989?cc=it&lang=en&
  * [35] Xiaowen Shan and Hudong Chen  “Lattice Boltzmann model for simulating flows with multiple phases and components”  In _Physical review E_ 47.3 APS, 1993, pp. 1815  DOI: 10.1103/physreve.47.1815
  * [36] Li Chen et al.  “A critical review of the pseudopotential multiphase lattice Boltzmann model: Methods and applications”  In _International Journal of Heat and Mass Transfer_ 76 Elsevier BV, 2014, pp. 210–236  DOI: 10.1016/j.ijheatmasstransfer.2014.04.032
  * [37] Kevin Connington and Taehun Lee  “A review of spurious currents in the lattice Boltzmann method for multiphase flows”  In _Journal of Mechanical Science and Technology_ 26.12 Springer ScienceBusiness Media LLC, 2012, pp. 3857–3863  DOI: 10.1007/s12206-012-1011-5
  * [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, 2019, pp. e3200  DOI: 10.1002/cnm.3200
  * [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  DOI: 10.1039/d0sm00587h
  * [40] 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  DOI: 10.1016/j.jcp.2003.10.013
  * [41] 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 IOP Publishing, 2007, pp. 285  DOI: 10.1088/1478-3975/4/4/005
  * [42] T. Krüger, F. Varnik and D. Raabe  “Efficient and accurate simulations of deformable particles immersed in a fluid using a combined immersed boundary lattice Boltzmann finite element method”  In _Computers & Mathematics with Applications_ 61.12 Elsevier BV, 2011, pp. 3485–3505  DOI: 10.1016/j.camwa.2010.03.057
  * [43] Sebastian J Müller et al.  “A hyperelastic model for simulating cells in flow”  In _Biomechanics and modeling in mechanobiology_ 20.2 Springer, 2021, pp. 509–520  DOI: 10.1007/s10237-020-01397-2
  * [44] F Pelusi et al.  “A sharp interface approach for wetting dynamics of coated droplets and soft particles”  In _Physics of Fluids_ 35.8, 2023, pp. 082126  DOI: 10.1063/5.0160096
  * [45] Y. Peng et al.  “Application of multi-block approach in the immersed boundary–lattice Boltzmann method for viscous fluid flows”  In _Journal of Computational Physics_ 218.2 Elsevier BV, 2006, pp. 460–478  DOI: 10.1016/j.jcp.2006.02.017
  * [46] Timm Krüger et al.  “The Lattice Boltzmann Method - Principles and Practice”  Springer Cham, 2016  DOI: 10.1007/978-3-319-44649-3
  * [47] Luca Biferale, Fabio Bonaccorso, Michele Buzzicotti and Chiara Calascibetta  “TURB-Lagr. A database of 3D Lagrangian trajectories in homogeneous and isotropic turbulence”  In _arXiv:2303.08662_ , 2023  URL: https://api.semanticscholar.org/CorpusID:257532470
  * [48] 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  DOI: 10.1017/s0022112070001696
  * [49] 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  DOI: 10.1017/s0022112073000534
  * [50] H. Lamb  “Hydrodynamics”  Cambridge University Press, 1932  URL: https://www.cambridge.org/it/universitypress/subjects/mathematics/fluid-dynamics-and-solid-mechanics/hydrodynamics-6th-edition?format=PB&isbn=9780521458689
  * [51] 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  DOI: 10.1017/s0022112069000759
  * [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  DOI: 10.1017/s0022112078002530
  * [53] C. Pozrikidis  “Boundary Integral and Singularity Methods for Linearized Viscous Flow”  Cambridge University Press, 1992  DOI: 10.1017/cbo9780511624124
  * [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  DOI: 10.1017/jfm.2016.39
  * [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  DOI: 10.1063/1.870305
  * [56] T. Krüger  “Computer Simulation Study of Collective Phenomena in Dense Suspensions of Red Blood Cells under Shear”  Vieweg+Teubner Verlag Wiesbaden, 2012  DOI: https://doi.org/10.1007/978-3-8348-2376-2
  * [57] Fabio Guglietta et al.  “Suspensions of viscoelastic capsules: effect of membrane viscosity on transient dynamics”  In _Journal of Fluid Mechanics_ 971, 2023, pp. A13  DOI: 10.1017/jfm.2023.694
  * [58] 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  DOI: 10.1007/s13239-020-00517-4
  * [59] 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  DOI: 10.1039/c3sm51645h
  * [60] 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, 2014, pp. 054114  DOI: 10.1063/1.4897913
  * [61] Stephan Gekle  “Strongly accelerated margination of active particles in blood flow”  In _Biophysical journal_ 110.2 Elsevier, 2016, pp. 514–520  DOI: 10.1016/j.bpj.2015.12.005
  * [62] Danilo P.. Silva et al.  “Lattice Boltzmann simulation of deformable fluid-filled bodies: progress and perspectives”  In _Soft Matter, DOI:10.1039/D3SM01648J_ The Royal Society of Chemistry, 2024  DOI: 10.1039/D3SM01648J
  * [63] Diego Taglienti, Fabio Guglietta and Mauro Sbragaglia  “Reduced model for droplet dynamics in shear flows at finite capillary numbers”  In _Physical Review Fluids_ 8.1 APS, 2023, pp. 013603  DOI: 10.1103/physrevfluids.8.013603
  * [64] Yue-Hong Qian, Dominique d’Humières and Pierre Lallemand  “Lattice BGK models for Navier-Stokes equation”  In _EPL (Europhysics Letters)_ 17.6 IOP Publishing, 1992, pp. 479  DOI: 10.1209/0295-5075/17/6/001
  * [65] Zhaoli Guo, Chuguang Zheng and Baochang Shi  “Discrete lattice effects on the forcing term in the lattice Boltzmann method”  In _Phys. Rev. E Stat. Nonlin. Soft Matter Phys._ 65.4 Pt 2B American Physical Society (APS), 2002, pp. 046308  DOI: 10.1103/physreve.65.046308
  * [66] Qisu Zou and Xiaoyi He  “On pressure and velocity boundary conditions for the lattice Boltzmann BGK model”  In _Physics of Fluids_ 9.6 AIP Publishing, 1997, pp. 1591–1598  DOI: 10.1063/1.869307
  * [67] Martin Hecht and Jens Harting  “Implementation of on-site velocity boundary conditions for D3Q19 lattice Boltzmann simulations”  In _Journal of Statistical Mechanics: Theory and Experiment_ 2010.01 IOP Publishing, 2010, pp. P01018  DOI: 10.1088/1742-5468/2010/01/p01018
  * [68] Michele Buzzicotti et al.  “Lagrangian statistics for Navier–Stokes turbulence under Fourier-mode reduction: fractal and homogeneous decimations”  In _New Journal of Physics_ 18.11 IOP Publishing, 2016, pp. 113047  DOI: 10.1088/1367-2630/18/11/113047
  * [69] M. Shapira and S. Haber  “Low Reynolds Number Motion of a Droplet in Shear Flow Including Wall Effects”  In _International Journal of Multiphase Flow_ 16.2, 1990, pp. 305–321  DOI: 10.1016/0301-9322(90)90061-m
  * [70] Alexander Z. Zinchenko, Michael A. Rother and Robert H. Davis  “A novel boundary-integral algorithm for viscous interaction of deformable drops”  In _Physics of Fluids_ 9.6 AIP Publishing, 1997, pp. 1493–1511  DOI: 10.1063/1.869275
  * [71] David P. Schmidt, Meizhong Dai, Haoshu Wang and J. Perot  “DIRECT INTERFACE TRACKING OF DROPLET DEFORMATION”  In _Atomization and Sprays_ 12.5–6 Begell House, 2002, pp. 721–736  DOI: 10.1615/atomizspr.v12.i56.110
  * [72] Haobo Hua et al.  “Effect of Confinement on Droplet Deformation in Shear Flow”  In _International Journal of Computational Fluid Dynamics_ 27.8-10, 2013, pp. 317–331  DOI: 10.1080/10618562.2013.857406
  * [73] Alireza Yazdani and Prosenjit Bagchi  “Influence of membrane viscosity on capsule dynamics in shear flow”  In _Journal of Fluid Mechanics_ 718 Cambridge University Press, 2013, pp. 569–595  DOI: 10.1017/jfm.2012.637
  * [74] J. Vollmer, R. Mencl and H. Müller  “Improved Laplacian Smoothing of Noisy Surface Meshes”  In _Computer Graphics Forum_ 18.3, 1999, pp. 131–138  DOI: https://doi.org/10.1111/1467-8659.00334
  * [75] Antonio Bacciaglia, Alessandro Ceruti and Alfredo Liverani  “Surface smoothing for topological optimized 3D models”  In _Structural and Multidisciplinary Optimization_ 64.6 Springer ScienceBusiness Media LLC, 2021, pp. 3453–3472  DOI: 10.1007/s00158-021-03027-6
  * [76] J.C. Butcher  “A history of Runge-Kutta methods”  In _Applied Numerical Mathematics_ 20.3 Elsevier BV, 1996, pp. 247–260  DOI: 10.1016/0168-9274(95)00108-5
  * [77] 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  DOI: 10.1016/0021-9797(72)90255-x
  * [78] 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  DOI: 10.1122/1.550942
  * [79] 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  DOI: 10.1122/1.1562152
  * [80] M. Tjahjadi and J.. Ottino  “Stretching and breakup of droplets in chaotic flows”  In _Journal of Fluid Mechanics_ 232.1 Cambridge University Press (CUP), 1991, pp. 191  DOI: 10.1017/s0022112091003671
  * [81] Sebastian J Müller, Ben Fabry and Stephan Gekle  “Predicting cell stress and strain during extrusion bioprinting”  In _Physical Review Applied_ 19.6 APS, 2023, pp. 064061  DOI: 10.1103/physrevapplied.19.064061
  * [82] Benjamin Owen et al.  “Lattice-Boltzmann modelling for inertial particle microfluidics applications - a tutorial review”  In _Advances in Physics: X_ 8.1 Informa UK Limited, 2023, pp. 1–76  DOI: 10.1080/23746149.2023.2246704
  * [83] Dario De Marinis, Alberto Mantegazza, Alessandro Coclite and Marco Donato Tullio  “A fluid-structure interaction method for soft particle transport in curved microchannels”  In _Computer Methods in Applied Mechanics and Engineering_ 418 Elsevier, 2024, pp. 116592  DOI: 10.1016/j.cma.2023.116592
  * [84] Mario Minale  “Models for the deformation of a single ellipsoidal drop: a review”  In _Rheologica acta_ 49.8 Springer, 2010, pp. 789–806  DOI: 10.1007/s00397-010-0442-0
  * [85] Dhruv Arora, Marek Behr and Matteo Pasquali  “A tensor-based measure for estimating blood damage.”  In _Artificial organs_ 28 11, 2004, pp. 1002–15  DOI: 10.1111/j.1525-1594.2004.00072.x
  * [86] Nico Dirkes, Fabian Key and Marek Behr  “Eulerian Formulation of the Tensor-Based Morphology Equations for Strain-Based Blood Damage Modeling”  In _arXiv preprint arXiv:2402.09319_ , 2024  URL: https://arxiv.org/abs/2402.09319
  * [87] 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 (CUP), 2016, pp. 480–501  DOI: 10.1017/jfm.2016.694