# A Lattice Boltzmann dynamic-Immersed Boundary scheme for the transport of
deformable inertial capsules in low-Re flows

Alessandro Coclite<EMAIL_ADDRESS>Sergio Ranaldo
<EMAIL_ADDRESS>Giuseppe Pascazio<EMAIL_ADDRESS>Marco
D. de Tullio<EMAIL_ADDRESS>Scuola di Ingegneria, Università
degli Studi della Basilicata, Viale dell’Ateneo Lucano 10 – 85100 Potenza,
Italy  Dipartimento di Meccanica, Matematica e Management, Politecnico di
Bari, Via Re David 200 – 70125 Bari, Italy

###### Abstract

In this work, a dynamic-Immersed–Boundary method combined with a BGK-
Lattice–Boltzmann technique is developed and critically discussed. The fluid
evolution is obtained on a three-dimensional lattice with 19 reticular
velocities (D3Q19 computational molecule) while the immersed body surface is
modeled as a collection of Lagrangian points responding to an elastic
potential and a bending resistance. A moving least squares reconstruction is
used to accurately interpolate flow quantities and the forcing field needed to
enforce the boundary condition on immersed bodies. The proposed model is
widely validated against well known benchmark data for rigid and deformable
objects. Rigid transport is validated by computing the settling of a sphere
under gravity for five different conditions. Then, the tumbling of inertial
particles with different shape is considered, recovering the Jefferey orbit
for a prolate spheroid. Moreover, the revolution period for an oblate spheroid
and for a disk-like particle is obtained as a function of the Reynolds number.
The existence of a critical Reynolds number is demonstrated for both cases
above which revolution is inhibited. The transport of deformable objects is
also considered. The steady deformation of a membrane under shear for three
different mechanical stiffness is assessed. Then, the tumbling of a weakly-
deformable spheroid under shear is systematically analyzed as a function
strain stiffness, bending resistance and membrane mass.

###### keywords:

Dynamic forcing , Moving Least Squares , Fluid-Structure Interaction ,
Capsules under shear , Deforming particle

††journal: Computer & Mathematics with applications

## Introduction

By definition a capsule is composed of a membrane, with prescribed
permeability and mechanical properties, enclosing some internal medium and
controlling the exchange (e.g., fluid and momentum) between the environment
and the enclosed medium [1, 2]. Examples of capsules can be found in nature
(cells, bacteria, eggs, eye-bulbs) as well as in the man-made universe
(medicines, cosmetics, food- and oil-industries). Specifically, the transport
of micro- and nano-membranes within established flows represents one of the
physical phenomena of great relevance for a wide range of engineering
applications. The scientific community is investing its efforts in general and
reliable numerical models able to describe micro- and nano-structures
navigating complex geometries, such as: transport of cells and biological
particles in the circulatory system and in in-vitro devices [3, 4, 5, 6],
systemic injection of drugs [7, 8, 9, 10, 11, 12]; rheology of active
colloidal suspensions [13, 14, 15]; bacteria, micro-swimmers and micro-robots
[16, 17, 18, 19, 20]; microfluidics devices for oil industrial applications
with fixed or variable geometries [21, 22, 23, 24]. Indeed, all of these
applications require the the detailed description of the complex mechanism
driving such particles’ transport. Reliable and accurate fluid-mechanical
models are thus fundamental for the rational design and optimization of
capsules abilities.

In this perspective, the Immersed Boundary (IB) method emerged a very
attractive technique for the treatment of rigid and deformable structures
immersed into the fluid stream [25]. Within this method, the fluid and the
structure solutions are obtained and coupled together at each time step via a
transfer function. Specifically, the fluid evolution is obtained on an
Eulerian grid not conformed to the immersed object, the latter being
represented by a collection of Lagrangian markers responding to some
constitutive models. The choice of the adopted constitutive models is
essentially arbitrary as well as the forcing procedure (either kinematic or
dynamic) transferring boundary conditions imposed by the presence of the body
onto the neighboring Eulerian fluid cells [26]. Recognizing its absolute
flexibility, several scientists employed the IB technique adopting different
constitutive models, fluid/structure coupling strategies, and forcing
procedures for the fluid boundary conditions. [27, 28, 29, 30, 31, 32].

In this work, a dynamic-Immersed–Boundary method is combined with a BGK-
Lattice– Boltzmann (BGK-LB) technique for describing the evolution of rigid
and soft inertial capsules transportded in incompressible flows. This method
was extensively studied and validated by the authors for two-dimensional
problems [9, 33, 34, 35]. Here, the authors propose an extension of such
scheme to three-dimensional manifolds. Specifically, the BGK-LB equation is
discretized and solved onto a three-dimensional lattice with 19 reticular
velocities (D3Q19 computational molecule) describing the fluid evolution with
second order accuracy. External boundaries are treated by imposing Dirichlet
boundary conditions for the disctribution functions as formulated by Zou-He
[36]. On the other hand, the argument by Guo et al. [37] for the projection of
the forcing term onto the D3Q19 computational molecule is employed to enforce
the immersed structures boundary conditions into the flow field. The Moving
Least Squares (MLS) reconstruction [38] is used to accurately interpolate the
velocity and forcing fields between the Eulerian and Lagrangian points. Both
rigid bodies and soft membranes are considered. The former ones are
transported by solving the Euler-Newton equations for the center of mass and
the revolution angles; the latter ones are transported as a collection of
Lagrangian points, connected by means of a triangular surface mesh, responding
to a bending resistance, an elastic potential, and the enclosed volume
conservation constraint.

The proposed model is validated against well known benchmark data for rigid
and deformable capsules. First, the second order accuracy of the present
formulation is demonstrated through a grid-refinement study with five
different Eulerian meshes on a lid-driven square cavity. Then, the position
and velocity of a sphere settling under gravity are computed and compared
against the experimental findings by Ten Cate et al. [39]. The falling
dynamics is analyzed for different values of Reynolds number, Froude number
and solid/fluid density ratio. The tumbling motion of inertial rigid particles
with different shape is also considered. On one side, the Jefferey orbit of a
prolate spheroid is recovered for small values of the Reynolds number (Re);
while, one the other side, the revolution period for an oblate spheroid and
for a disk-like particle is obtained as a function of Re and the existence of
a critical value is demonstrated above which rotation is inhibited.
Furthermore, the validation of the deformable particle constitutive model is
carried out by computing the steady deformation of an elastic capsule under
shear [40]. The tumbling period of a weakly-deformable prolate spheroid under
shear is systematically analyzed as a function of: strain stiffness, regulated
by the non-dimensional shear rate (G = 0.01, 0.025, 0.05, 0.075, 0.1); bending
resistance (Eb = 0.001, 0.005, 0.01); and membrane mass ($\rho s/\rho f=$ 1,
1.2, 1.5, 2, 5). Interestingly, the revolution period seems to be only
influenced by the bending resistance. On the contrary, membrane stiffness and
mass are responsible for the larger or smaller surface relative strain without
affecting the revolution period.

##  1 Method

###  1.1 The 3D-BGK lattice-Boltzmann method

The fluid evolution is modeled on a 3D computational molecule in the velocity
space with 19 directions, the D3Q19 lattice. These reticular directions are
identified by 19 vectors in $\mathbb{R}^{3}$ corresponding to discrete
velocities ${\bf e_{i}}({\bf x})$ ($i=0,\ldots,18$), where ${\bf x}$ are the
spatial coordinates, whose Cartesian coordinates are the columns of the
following matrix:

${\bf
e}=\setcounter{MaxMatrixCols}{19}\begin{pmatrix}0&1&-1&0&0&0&0&1&1&1&1&-1&-1&-1&-1&0&0&0&0\\\
0&0&0&1&-1&0&0&1&-1&0&0&1&-1&0&0&1&1&-1&-1\\\
0&0&0&0&0&1&-1&0&0&1&-1&0&0&1&-1&1&-1&1&-1\\\ \end{pmatrix}$ (1)

On this cubic lattice, 19 discrete distribution functions $\\{f_{i}\\}$
($i=0,\ldots,18$) are transported responding to the forced Boltzmann equation,

${f_{i}({\bf x}+{\bf e}_{i}\Delta t,t+\Delta t)-f_{i}({\bf x},t)=-\frac{\Delta
t}{\tau}[f_{i}({\bf x},t)-f_{i}^{eq}({\bf x},t)]+\Delta t{\cal F}_{i}}\,,$ (2)

with $\Delta t$ the time step, $\tau$ the unique relaxation time in the BGK-
approximation [41], $\\{f^{eq}_{i}\\}$ ($i=0,\ldots,18$) the local equilibrium
distribution functions, and ${\cal F}_{i}$ ($i=0,\ldots,18$) the components of
the forcing term accounting for the body boundary condition. The kinematic
viscosity of the flow is related to $\tau$ as $\nu=c_{s}^{2}(\tau-\frac{\Delta
t}{2})$ being $c_{s}=1/\sqrt{3}$ the reticular sound speed. The moments of the
distribution functions define the fluid density $\rho=\sum_{i}f_{i}$, the
momentum $\rho{\bf u}=\sum_{i}f_{i}{\bf e}_{i}+\frac{\Delta t}{2}{\bf f}_{lb}$
and the pressure $p=c_{s}^{2}\sum_{i}f_{i}=c_{s}^{2}\rho$ [42].
The local equilibrium distribution functions $\\{f_{i}^{eq}\\}$
($i=0,\ldots,18$) are chosen to be representative of the Maxwell-Boltzmann
distribution projected onto the lattice as follows:

$f_{i}^{eq}({\bf x},t)=\omega_{i}\rho\left[1+\frac{1}{c_{s}^{2}}({\bf
e}_{i}\cdot{\bf u})+\frac{1}{2c_{s}^{4}}({\bf e}_{i}\cdot{\bf
u})^{2}-\frac{1}{2c_{s}^{2}}{\bf u}^{2}\right]\,,$ (3)

with $\omega_{0}=1/3$, $\omega_{1-6}=1/18$, and $\omega_{7-18}=1/36$.
Following the argument by Guo et al. [43], ${\cal F}_{i}$ reads:

${{\cal F}_{i}=\Bigl{(}1-\frac{1}{2\,\tau}\Bigr{)}\omega_{i}\Bigl{[}\frac{{\bf
e_{i}}-{\bf u}}{c_{s}^{2}}+\frac{{\bf e_{i}}\cdot{\bf u}}{c_{s}^{4}}{\bf
e_{i}}\Bigr{]}\cdot{\bf f}_{lb}}\,$ (4)

where ${\bf f}_{lb}({\bf x},t)$ is a the body force term accounting for the
presence of immersed objects. Dirichlet boundary conditions are imposed on
external boundaries treated with the known–velocity bounce back procedure by
Zou and He [36]. It is proved that in such a framework one can recover the
forced Navier–Stokes equation with second order accuracy [44, 45].

Although, for each $\tau$ in the stability range $\tau>\frac{1}{2}$ the scheme
converge to the same solution, higher value of $\tau$ lead to higher value of
the reference velocity $u_{ref}$ and, consequently, higher non-dimensional
time interval $L_{ref}/u_{ref}$ corresponding to a single time step ($L_{ref}$
being the reference length) [42]. As the matter of fact, $\tau$ regulates the
viscosity and $\nu$ regulates $u_{ref}$ through the Reynolds number Re =
$\frac{u_{ref}L_{ref}}{\nu}$. Due to these reasons, all of the simulations in
the present work are computed with $\tau=0.75$.

###  1.2 Immersed Boundary treatment

In the Immersed-Boundary technique [25], an obstacle in the flow can be
considered as a collection of Lagrangian markers superimposed to the Eulerian
withstanding fluid lattice (Figure.1a). Here, the surface of the immersed
geometry is discretized by means of a triangulated mesh, the Lagrangian
markers coinciding with the centroids of each triangular element. The
Lagrangian grid and Eulerian lattice interact through a transfer function
built on a moving–least-squares (MLS) reconstruction [38], exchanging the
distribution functions, while the body force term in Eq.(4), ${\bf f}_{ib}$,
is evaluated through the formulation by Favier et al. [46].

Figure 1: Immersed-boundary technique schematic. a. Geometry description by
means of triangular elements, whose centroids coincide with the Lagrangian
markers. b. Definition of the support domain (blue dashed square) centered on
the selected Lagrangian marker (black circle) and the nine Eulerian points
(per plane) involved in the forcing (blue squares). c. Sketch of the two nodes
(vertices of an edge) involved in the elastic spring model. Right: four nodes
(belonging to the two triangles sharing an edge) involved in the bending
spring model.

Given a Lagrangian marker (with index $l$), 27 Eulerian points are identified,
falling into the support domain, defined as a cube with side equal to
$r_{w}=1.3\,\Delta x$ and centered on the Lagrangian marker (see Figure.1b).
Given the solution at time level $n$, the velocity component of Lagrangian
marker l is evaluated as,

${U_{i}({\bf x})=\sum_{k=1}^{27}\,\phi^{l}_{k}({\bf x})u^{k}_{i}}\,,$ (5)

where $u^{k}_{i}$ indicates the $i$-th component of the velocity at the $k$-th
Eulerian point associated with the marker $l$ and $\phi$ is the transfer
operator obtained minimizing with respect to ${\bf a(x)}$ the following
weighted L2-norm:

${J=\sum_{k=1}^{27}W({\bf x}-{\bf x}^{k}})[{\bf p}^{T}({\bf x}^{k}){\bf
a(x)}-u^{k}_{i}]^{2}\,.$ (6)

In the above equation, ${\bf p}^{T}$ is a basis function vector and ${\bf
a(x)}$ a vector of coefficients such that $\sum_{k=1}^{27}\,\phi^{l}_{k}({\bf
x})u^{k}_{i}\,={\bf p}^{T}({\bf x}^{k}){\bf a(x)}$, ${\bf x}$ the Lagrangian
marker position, and $W({\bf x}-{\bf x}^{k})$ a weight function. In this work,
a linear basis function, ${\bf p}^{T}=(1,x,y,z)$, is considered along with an
exponential weight function:

$W({\bf x}-{\bf x}^{k})=\begin{cases}e^{-(r_{k}/\alpha)^{2}}\,,&\mbox{if
}r_{k}\leq 1\\\ 0\,,&\mbox{if }r_{k}>1\ \,,\\\ \end{cases}$ (7)

where, $\alpha=0.3$ and $r_{k}$ is the distance between the Lagrangian point
and the associated $k$-th Eulerian point normalized over the size of the
support domain, $r_{w}$, $r_{k}=\frac{|{\bf x}-{\bf x}^{k}|}{r_{w}}$. For each
Lagrangian marker the volume force ${\bf F}_{l}$ required to impose the
boundary condition can be evaluated and then transferred to the Eulerian
points into the support domain using the same functions used before. Being
${\bf U}_{b}^{l}({\bf x})$ the desired velocity at the boundary related to the
$l$-th triangle, the Lagrangian volume force is:

${{\bf F}_{l}{\bf(x)}=\frac{{\bf U}_{b}^{l}{\bf(x)}-{\bf U(x)}}{\Delta t}}\,,$
(8)

and then the body force for the $k$-th Eulerian point is finally,

${{\bf f}_{lb}}^{k}=\sum_{l}c_{l}\phi^{l}_{k}{\bf F}_{l}\,,$ (9)

where the sum is extended to all Lagrangian markers having the $k$-th point in
their support domain. The scale coefficient $c_{l}$ is conveniently obtained
by imposing the conservation between the Lagrangian and Eulerian grids of the
total force acting on the fluid [47]. The proposed algorithm concludes by
inserting ${\bf f}_{ib}$ in equation (4) and starting the next step by
transporting the distribution functions ${\bf f}_{i}$ through equation (2). In
this way, the information given by the forcing term evaluated at the time
level $n$ is transferred to the $n+1$-th time level.

###  1.3 Fluid-Structure Interaction

Immersed bodies motion is determined by dynamics IB technique described in
[26], solving the Newton equation for each mass point, accounting for both
internal and external stresses, namely, elastic boundary response and
hydrodynamics forces, respectively. Here, a weak coupling approach is adopted.
[33]

Network model for deforming bodies. Immersed bodies internal stresses are
computed as elastic strain response, bending resistance, and enclosed volume
conservation (if applicable). The stretching elastic potential acting on the
two vertices sharing the $l$-th edge (Figure.1c) is given as

${V_{l}^{s}=\frac{1}{2}k_{s,l}(l_{l}-l_{l,0})^{2}}\,,$ (10)

being $k_{s,l}$ [N/m] the elastic constant associated to the edge $l$, $l_{l}$
its current length, and $l_{l,0}$ its length in the stress-free configuration.
The nodal forces corresponding to the elastic energy for vertices 1 and 2
connected by l reads:

${\begin{cases}{\bf F}_{1}^{int,s}=-k_{s,l}(l-l_{0})\frac{{\bf
r}_{1,2}}{l}\,,\\\ {\bf F}_{2}^{int,s}=-k_{s,l}(l-l_{0})\frac{{\bf
r}_{2,1}}{l}\,,\\\ \end{cases}}$ (11)

where ${\bf r}_{i,j}={\bf r}_{i}-{\bf r}_{j}$ and ${\bf r}_{i}$ is the
position vector of i with respect to j. Here, the elastic constant for the
$l$-th edge is computed trough the capillary number Ca$=\frac{\rho\nu
u_{ref}}{k_{s,0}}$ and then distributed with a surface-weighted-averaging
process (Van Gelder model [48]),
$k_{s,l}=\frac{k_{s,0}T\sum_{i}A^{l}_{i}}{l^{2}}$. In the above relationships
$\rho$ and $\nu$ correspond to surrounding fluid density and kinematic
viscosity while $u_{ref}$ is a reference velocity, $k_{s,0}$ the body global
elastic constant measured in [N/m], $T$ the immersed body thickness and
$A^{l}_{i}$ the area of the two triangular elements sharing the $l$-th edge.

The out-of-plane deformation of two elements sharing an edge is computed as a
bending resistance involving the four vertices of this couple of triangles.
Given the initial (stress-free) and current local curvature, $k_{v,0}$ and
$k_{v}$, respectively, the bending potential reads:

${V_{v}^{b}=\frac{1}{2}k_{b}(k_{v}-k_{v,0})^{2}}\,,$ (12)

with $k_{b}$ the bending constant. The local curvature is computed by
measuring the variation of the angle between the unit normal vectors of two
adjacent elements, ${\bf n}_{1}$ and ${\bf n}_{2}$, ($\theta-\theta_{0}$),
with $\theta_{0}$ the angle in the stress free configuration (Figure.1c). The
bending constant associated to the two selected triangles is related to
$k_{s,l}$ through the flexural rigidity $Eb=\frac{k_{b}}{k_{s,l}r^{2}}$. [49]
The resulting forces acting on the four nodes read:

${\begin{cases}{\bf F}^{int,b}_{1}=\beta_{b}\bigl{[}b_{11}({\bf
n}_{1}\times{\bf r}_{32})+b_{12}({\bf n}_{2}\times{\bf r}_{32})\bigr{]}\,\\\
{\bf F}^{int,b}_{2}=\beta_{b}\bigl{[}b_{11}({\bf n}_{1}\times{\bf
r}_{13})+b_{12}({\bf n}_{1}\times{\bf r}_{34}+{\bf n}_{2}\times{\bf
r}_{13})+b_{22}({\bf n}_{2}\times{\bf r}_{34})\bigr{]}\,\\\ {\bf
F}^{int,b}_{3}=\beta_{b}\bigl{[}b_{11}({\bf n}_{1}\times{\bf
r}_{21})+b_{12}({\bf n}_{1}\times{\bf r}_{42}+{\bf n}_{2}\times{\bf
r}_{21})+b_{22}({\bf n}_{2}\times{\bf r}_{42})\bigr{]}\,\\\ {\bf
F}^{int,b}_{4}=\beta_{b}\bigl{[}b_{11}({\bf n}_{1}\times{\bf
r}_{23})+b_{22}({\bf n}_{2}\times{\bf r}_{23})\bigr{]}\,\\\ \end{cases}}\,,$
(13)

with $b_{ij}=\frac{{\bf n}_{i}\cdot{\bf n}_{j}}{|{\bf n}_{i}||{\bf n}_{j}|}$
and
$\beta_{b}=k_{b}\frac{sin(\theta)cos(\theta_{0})-cos(\theta)sin(\theta_{0})}{\sqrt{1-cos^{2}(\theta)}}$.

The fluid volume enclosed by of the incompressible membrane is constrained
with a penalty force acting on the ith-element centroid:

${{\bf F}_{i}^{int,v}=-k_{v}\left(1-\frac{V}{V_{0}}\right){\bf
n}_{i}A_{i}}\,.$ (14)

Here, $k_{v}$ is the incompressibility constant, $V$ and $V_{0}$ the current
and initial enclosed volume and $A_{i}$ the current area of the $i$-th
element.

Hydrodynamics Stresses. Pressure and viscous stresses exerted by the fluid on
the $i$-th triangular element are:

${{\bf F}_{i}^{ext,p}(t)=-(p^{+}_{i}-p^{-}_{i})\,{\bf n}^{+}_{i}\,A_{i}}\,,\\\
$ (15) ${{\bf
F}_{i}^{ext,\tau}(t)=(\bar{\bf\tau}^{+}_{i}-\bar{\bf\tau}^{-}_{i})\cdot{\bf
n}^{+}_{i}\,A_{i}}\,,$ (16)

where $\bar{\bf\tau}^{+}_{i}$, $\bar{\bf\tau}^{-}_{i}$ and $p^{+}_{i}$,
$p^{-}_{i}$ are the viscous stress tensor and the pressure on the centroid of
the $i$-th element taken from the external (+) and internal fluid (-),
respectively. ${\bf n^{+}}_{i}$ is the normal outward unit verctor. The
pressure and velocity derivatives in Eq.s (15) and (16) are evaluated using a
probe in the normal positive and negative direction of each element, being the
probe length $1.2\,\Delta x$, and using the presented moving least squares
recontruction [47]. In this framework, the velocity derivatives evaluated at
the probe are considered equal to the ones on the triangular element centroid
and all force contributions are computed with respect to the centroid of each
elements and then transferred to the vertices. [32]

Elastic Membranes. The total force ${\bf F}_{v}^{tot}(t)$ acting on the $v$-th
element of the immersed body is evaluated in time and the position of the
vertices is updated at each time step considering the membrane mass uniformly
distributed over the vertices,

${m_{v}\dot{{\bf u}}_{v}={\bf F}_{v}^{tot}(t)={\bf F}_{v}^{int,s}(t)+{\bf
F}_{v}^{int,b}(t)+{\bf F}_{v}^{int,v}(t)+{\bf F}_{v}^{ext,p}(t)+{\bf
F}_{v}^{ext,\tau}(t)}\,,$ (17)

where $m_{v}$ and $\dot{{\bf u}}_{v}$ indicate the mass and acceleration of
vertex $v$, respectively.

The Newton equation of motion is integrated with the Verlet algorithm using,
as first tentative velocity, the value obtained by interpolating the fluid
velocity from the surrounding lattice nodes, ${\bf u}_{v,0}(t+\Delta t)$,

${{\bf x}_{v}(t+\Delta t)={\bf x}_{v}(t)+{\bf u}_{v,0}(t+\Delta t)\Delta
t+\frac{1}{2}\frac{{\bf F}_{v}^{tot}(t)}{m_{v}}\Delta t^{2}+O(\Delta
t^{3})}\,;$ (18)

then, the velocity at the time level $t+\Delta t$ is computed as:

${{\bf u}_{v}(t+\Delta t)=\frac{\frac{3}{2}{\bf x}_{v}(t+\Delta t)-2{\bf
x}_{v}(t)+\frac{1}{2}{\bf x}_{v}(t-\Delta t)}{\Delta t}+O(\Delta t^{2})}\,.$
(19)

Finally, the position at the time level $t+\Delta t$ is computed employing
again eq. (18) with the velocity obtained by eq. (19).

Rigid objects. Rigid motion is readily obtained integrating all stress
contributions, eqs. (15) and (16), over the boundary of the immersed body and
then updating both linear and angular velocity: ${\bf F}^{tot}(t)=m\,\dot{{\bf
u}}(t)$ and ${\bf M}^{tot}(t)={\bf I}\,\dot{{\bf\omega}}(t)$. Here $m$ is the
particle mass, ${\bf F}^{tot}(t)$ is the total force exerted by the fluid on
the particle; ${\bf M}^{tot}(t)$ is the total moment acting on the rigid body,
and ${\bf I}$ is the tensor of inertia. Finally, ${\bf u}(t)$ and
${\bf\omega}(t)$ are computed as:

${{\bf u}(t)=\frac{2}{3}(2{\bf u}(t-\Delta t)-\frac{1}{2}{\bf u}(t-2\Delta
t)+\dot{{\bf u}}(t)\Delta t)+O(\Delta t^{2})}\,,$ (20)
${{\bf\omega}(t)=\frac{2}{3}(2{\bf\omega}(t-\Delta
t)-\frac{1}{2}{\bf\omega}(t-2\Delta t)+\dot{{\bf\omega}}(t)\Delta t)+O(\Delta
t^{2})}\,.$ (21)

###  1.4 Parallel Computing Framework

Recently, a number of parallelization techniques tailored on the BGK–LB scheme
have been proposed [50, 51, 52]. The D3Q19 computational molecule is composed
of a set of nine points per plane arranged along three computational rows thus
resulting in a square with side length equal to 2 $\Delta x$. Due to this
compact nature of the D3Q19 computational molecule, one of the most suitable
choice is to split the computational lattice into Cartesian blocks and entrust
the solution of Eq.(2) to different threads; each thread is appointed for the
computational nodes within its boundaries and for its communicating boundary
planes (back and forward) to adjacent threads. This logic is readily
implemented by using the Message Passing Interface (MPI) Cartesian CPU-grid
generator function with dimensions $N^{cpu}_{x},\,N^{cpu}_{y},\,N^{cpu}_{z}$
in the three spatial directions, respectively. So that, the total number of
threads evaluating the computational lattice is given by $N^{cpu}_{x}\times
N^{cpu}_{y}\times N^{cpu}_{z}$. At each time step the distribution functions
corresponding to boundaries (either planes, edges or angular points) between
blocks are transferred from one to the adjacent block forward and backward in
the three coordinate directions. Two ghost-node lines are defined for each
spatial direction of each Cartesian block to have that the support domain of a
Lagrangian marker falling into a block would be completely defined within the
same block. The immersed structure dynamics is handled using the MPI Cartesian
block subdivision, each CPU computing the dynamics of the Lagrangian markers
falling into its boundaries. Specifically, each thread reads the immersed
structure local positions and velocities while only the thread owing the
selected Lagrangian marker computes its internal and external stresses thus
resulting in non–null values of the forces defined in Eq.s (11), (13), (14),
(15), and (16). In this way the forces exerted by the structure derive from
the sum of those calculated by all of the CPUs. Indeed, the sum is done
through an MPI collective communicator.

##  2 Results and discussion

###  2.1 Scheme validation and accuracy: Lid-driven cavity test case

Figure 2: Lid-driven cavity. a. Sketch of the physical problem. (b,c)
Comparison between the distributions of $u_{x}$, $u_{y}$, and $u_{z}$ taken at
$x=z=0.5L$ (b) and $x=y=0.5L$ (c) and the benchmark data by Napolitano and
Pascazio [53] for the lid-driven cavity at Re=100 computed with L=50, 100,
150, 200 and 250 $\Delta x$. Figure 3: Lid-driven cavity with a spherical
obstacle immersed in. a. Contour plot of $u_{y}$ (left) and $u_{z}$ (right)
taken at $x=0.5\,L$ for L = 250 $\Delta x$. (b,c) Distribution of $u_{x}$,
$u_{y}$, and $u_{z}$ taken at $x=z=0.5L$ (b) and $x=y=0.5L$ (c) for the lid-
driven cavity at Re=100 computed with L=50, 100, 150, 200 and 250 $\Delta x$
with a spherical obstacle in the center by diameter $d=0.2\,L$.

To validate the computational method for the fluid evolution and determine its
accuracy a lid-driven cavity test case is computed at first. Specifically, a
box with side L is considered, discretized with five meshes, containing 50,
100, 150, 200 and 250 $\Delta x$ cells. The wall at $y=y_{max}$ (see
Figure.2a) moves with velocity ${\bf u}_{w}=(0,0,u_{max})$. The upper wall
velocity is regulated by the Reynolds number Re (=$\frac{u_{max}L}{\nu}$),
$\nu$ being the fluid kinematic viscosity, related to the relaxation parameter
$\tau$ equal to one for all of these computations. The resulting velocity
distributions taken at x=z=0.5 L (Figure.2b) and x=y=0.5 L (Figure.2c) are
compared with the benchmark solutions by Napolitano and Pascazio [53] thus
validating the quality of the present predictions. In particular, for the
distribution at x=z=0.5 L, the largest difference for $u_{z}$ with respect to
the reference solution is found for L=50 $\Delta x$ at 0.48 L, resulting in a
relative error of 0.015. For the distribution taken at x=y=0.5 L, the largest
difference for $u_{y}$ with respect to the reference solution is found for
L=50 $\Delta x$ at 0.75 L and the relative error is 0.018. The $L^{2}$-norm of
the $u_{z}$ relative error with respect to the reference solution (obtained
with L=250 $\Delta x$), $\epsilon=||\frac{u-u_{L=250\Delta x}}{u_{L=250\Delta
x}}||^{2}$, is depicted in Figure.4. The distribution of $\epsilon$ as a
function the mesh size (parameterized with 1/N) shows the second-order
accuracy of the method, with the error already largely smaller than 1$\%$ for
L=50 $\Delta x$.

Figure 4: Grid refinement study. Mesh-refinement study on $u_{z}$
distributions for the lid-driven cavity with and without the spherical
obstacle immersed in.

Moreover, in order to analyze the overall scheme accuracy within the coupling
technique between fluid and the immersed boundary, a solid spherical obstacle
with diameter 0.2 L is placed in the center of the computational domain (see
Figure.3a). The solid body is discretized with a number of triangular elements
corresponding to about five Lagrangian points for each Eulerian cell. The
computed flow field is reported in Figure.3a depicting the velocity components
along z and y axis for the reference Eulerian grid L = 250 $\Delta x$.
Moreover, the distributions of $u_{x}$, $u_{y}$, and $u_{z}$ taken at x=z=0.5
L and x=y=0.5 L are documented for all of the five computed Eulerian grids in
Figure.3b and Figure.3c, respectively. The $L^{2}$-norm of the $u_{z}$
relative errors depicts that overall, while preserving the second-order space
accuracy, scheme errors are slightly higher when considering the coupled
fluid-structure system although well below O(10-2) when N$\geq 100$ (see
Figure.4).

###  2.2 Rigid spheres settling under gravity

Figure 5: Rigid spheres settling under gravity. a. Sketch of the physical
problem with characteristic dimensions and non-dimensional groups used for the
computations. b. Comparison between the computed settling velocity
distributions and benchmark experimental data by Ten Cate et al [39]. The
dashed line corresponds to the velocity distribution obtained with d = 50
$\Delta x$. c. Distribution of the particles’ centroid positions obtained
during the settling process. d. Measure of the speedup performances computed
over two different computational meshes, namely with $L_{z}$ = 160 and 400
$\Delta x$ over different number of threads running in parallel. Solid lines
refer to present model computations while dashed lines correspond to
theoretical scaling.

The simulation of a rigid sphere settling under gravity in a square cylinder
filled with fluid at rest is presented here against the experiments performed
by A. ten Cate et al. [39] used as benchmark. These experiments analyzed both
centroid vertical velocity and position using Particle Image Velocimetry (PIV)
technique, thus providing accurate measures of the falling spheres from the
initial condition at rest until the stop at the bottom of the channel. The
computational domain is a box with size
$(L_{x},L_{y},L_{z})=(6.67\,d,6.67\,d,10.67\,d)$, $d$ corresponding to the
diameter of the settling sphere discretized here with 30 $\Delta x$. The
particle is a rigid plain sphere placed in the center of the channel at
$L_{p,z}=8\,d$; its surface is described by a triangular mesh of about 10’000
elements. As previously analyzed by the authors, the ratio between the
Lagrangian mesh and the Eulerian grid characteristic lengths cannot be
arbitrary into this specific immersed boundary framework [26, 33]. A
convenient choice would be to keep this ratio between 0.3 and 0.5
corresponding to three-to-five Lagrangian points per each Eulerian cell. The
sphere density is fixed by the ratio between the solid and the fluid density
$\gamma=\rho_{s}/\rho_{f}$ while the Reynolds number
(Re=$\frac{du_{inf}}{\nu}$) regulate the fluid viscosity $\nu$ and is based on
the settling velocity $u_{inf}$ obtained by the relation for the drag
coefficient $C_{d}$ by Abraham [54]:

${C_{d}=C_{0}\left(1+\frac{\delta_{0}}{\sqrt{Re}}\right)^{2}}\,,$ (22)

where $C_{0}\delta_{0}^{2}=24$ and $\delta_{0}=9.06$ and reads,

${u_{inf}=\sqrt{\frac{4gd}{3C_{d}}(\gamma-1)}}\,.$ (23)

Finally, the Froude number (Fr =$\frac{u_{inf}}{\sqrt{gd}}$) returns the
gravity acceleration $g$. The values for Re, Fr, and $\gamma$ used in the
computations are reported in Figure.5a. The agreement between present model
predictions and the experimental data is depicted in term of the distributions
of the particles’ centroid velocities during the settling process (Figure.5b).
The relative error between the two set of data, based on the displacement
achieved in the peak z-velocity,
$\epsilon=\frac{max(|u_{LBM}|)-max(|u_{exp}|)}{max(|u_{exp}|)}$ ranges from
0.0023 to 0.0039 for Re = 1.5 and 31.9, respectively. For Re = 31.9 and Fr =
0.334 the ascendant part of the velocity distribution shows a slight
displacement from the experimental data while the descendant part and the peak
value are well reproduced. In fact, due to the increasing of both the Froude
and the Reynolds numbers, the particle passes through multiple computational
cells in a single Newtonian time step. The proposed model for enforcing the
boundary conditions on the immersed body needs to gradually reconstruct the
pressure field based on the corrections made on the velocity field avoiding
any refilling procedure for solid-to-fluid cells transition [33]. This gives a
slight inaccuracy into the ascendant part of the velocity distribution that is
essentially determined by the restraining of the particle because of the
pressure field near the southern wall [39]. By increasing the resolution of
the Eulerian grid (d = 50 $\Delta x$) the pressure field is accurately
reconstructed and the displacement between the two velocity distributions
disappears (see the dashed line in Figure.5b).

The trajectory of the settling spheres are shown in Figure.5c. To assess the
performance of the Cartesian parallel framework used and presented in the
Method section, the test case with the higher Re is re-computed on two
different meshes with $L_{z}$ = 160 and 400 $\Delta x$ over a growing number
of threads running in parallel. As depicted in Figure.5d, for the larger mesh
the weak-speedup performance is quite similar to the theoretical expectations
while for the smaller mesh the communications between the CPUs dominate the
wall time increasing the overall computational burden. Then, a strong-speedup
test is performed by computing three different meshes and varying with the
same growth factor the number of CPUs running in parallel for each Cartesian
directions thus keeping constant the number of computational cells entrusted
to each tread. This test returns similar wall times for the three meshes
documented in Table.1.

CPU | Points in $L_{z}$ | Total Wall Time
---|---|---
$1^{3}$ | $100$ | $1680\,s$
$2^{3}$ | $200$ | $1721\,s$
$4^{3}$ | $400$ | $1778\,s$
Table 1: Parallel performance for different lattice dimensions and number of
CPUs used.

###  2.3 Tumbling motion of rigid spheroids under shear

Figure 6: Rigid spheroids under shear. a. Sketch of the physical problem. b.
Comparison between the angular velocity obtained with the analytical solution
by Jefferey [55] and present model predictions for a prolate plain spheroid
freely moving under shear. (c,d) Angular velocity distribution for an oblate
spheroid (c) and a discoidal particle (d) under shear obtained at different
Reynolds number. e. Revolution period as a function of the Reynolds number for
an oblate spheroid and a discoidal particle. The solid lines correspond to the
analytical predictions.

To further validate the proposed model, the revolution of a plain
neutrally–buoiant rigid spheroid in a shear flow established by two
countermoving walls is studied (Figure.6a). First, the revolution period of an
initially resting prolate spheroid placed in the center of a cubic domain is
computed for a shear flow at Re = 10-2. The particle Reynolds number follows
the definition by Huang et al [56] and is based on the moving walls velocity,
$u_{max}$, the particle major axis length, d, and the kinematic viscosity of
the surrounding fluid, Re=$\frac{4\gamma d^{2}}{\nu}$ being
$\gamma=\frac{u_{max}}{1/2L_{y}}$ the shear rate. Such particle is a spheroid
with aspect ratio 2 (ratio of the major diameter over the two minor ones), the
length of its major axis is d ($=60\Delta x$) and its surface is discretized
with 20’500 triangular elements. The particle is placed with the major axis
parallel to x in the center of the cubic box with side
$L_{x}=L_{y}=L_{z}=4\,d$. The distribution of the rotational velocity around x
is compared with the analytical solution obtained by Jefferey [55] based on
the hypothesis of minimum dissipation energy thus finding an excellent
agreement considering both the revolution period and the angular velocity
distribution (Figure.6b). Specifically, the present model prediction for the
revolution period is 7.853 against the analytical value of 7.854 while the
largest displacement for $\omega d/u_{max}$ is found at 15.67 $tu_{max}/d$
returning a relative error of 0.000832.

Then, the revolution period of an oblate spheroid rotating under shear at
increasing values of the Reynolds number is measured. This spheroid is an
oblate particle with aspect ratio 2 (ratio of the two major diameters over the
minor one), d corresponding to the length of the two major diameters
discretized with 60 $\Delta x$, again the particle surface is discretized with
20’500 triangles. Initially at the rest, the particle is placed in the center
of a box with dimension $L_{x}=2.5\,d$, $L_{y}=5.0\,d$, and $L_{z}=6.8\,$.
Ding and Aidun [57] and Zettner and Yoda [58] found the existence of a
critical Reynolds number (Rec) for the angular tumbling velocity; for
$Re>Re_{c}$ the particle stops its revolution motion assuming a fixed position
with the plane containing the major axis forming an angle $\theta_{f}$ with
respect to z direction. Precisely, the tumbling period Tt as a function of the
Reynolds number is predicted by the scaling law, $T_{t}=c(Re_{c}-Re)^{-1/2}$.
Huang et al [56] observed that for an oblate particle with aspect ratio 2 the
critical Reynolds number is 80 and c = 180 while Ding and Aidun [57] found
Re${}_{c}=81$ and $c=200$. In the present formulation the tumbling angular
velocity recovers the described phenomenology, as demonstrated in Figure.6c;
the scaling law is recovered by using $c=179.5$ and Re${}_{c}=81.3$ (see
Figure.6e) and $\theta_{f}=0.29\,\pi$ for $Re=90$. Lastly, the tumbling motion
of a circular cylinder with aspect ratio 2 (ratio between the circle diameter
d and its height) is analyzed. Specifically, the discoidal particle presents a
diameter d corresponding to 60 $\Delta x$ and its surface covered by about
30’000 triangles. Interestingly, the discoidal particle tumbling presents a
critical Reynolds number of 192.5; above this value the revolution of the disc
does not occur, as depicted in Figure.6d, and $\theta_{f}=0.16\,\pi$ for
$Re=200$. The distribution of tumbling period as a function of the Reynolds
number is documented in Figure.6e. Such distribution does not follow a power
law in the form $T_{t}=c(Re_{c}-Re)^{-1/2}$ as seen in the oblate spheroid
tumbling. However, a scaling law may be established in the form
$T_{t}=c[(A-Re)(Re_{c}-Re)^{-\frac{1}{50}}]$ with $c=\frac{1}{10}$, A = 182,
and Rec = 195.

###  2.4 Deforming capsules under linear laminar flow

Figure 7: Capsules deforming under shear. a. Schematic representation of a
spherical capsule deforming under linear laminar flow. b Comparison between
the distributions of the Taylor parameter over time obtained for G = 0.0375,
0.075, and 0.15 with the proposed model and the benchmark data by Lac et al.
[40]. c Relative error obtained for the three investigated values of G with
respect to the benchmark data by Lac et al. [40]. d Variation of the capsule
relative volume over time as a function of G. e Configurations of the steady-
state deformed capsules for G = 0.0375, 0.075, and 0.15 predicted by the
proposed Dynamic-IB scheme as compared to the Kinematic-IB scheme and to the
computations by Lac et al. [40] as well as with the experimental data by
Walter et al [59].

The dynamics of the shear-induced deformation of an initially spherical
capsule is studied in this section. Referring to Figure.7a, a spherical
capsule with diameter d is placed at the center of a cubic box with side
length H = 10 d. The capsule diameter is discretized by 10 $\Delta x$ while
the capsule surface is discretized with $n_{L}=$ 10000 triangular elements.
The imposed shear rate is $\gamma=\frac{u_{max}}{H/2}$, with $u_{max}$ the
velocity magnitude of the top and bottom walls; the resulting Reynolds number
equals 0.01 ($=\frac{\gamma(d/2)^{2}}{\nu}$). The mechanical behavior of the
capsule is regulated by the elastic constant $k_{s}$, computed using the
dimensionless shear rate G $=\frac{\nu\rho\gamma d/2}{k_{s}}$. The mass of the
capsule is defined considering a unitary solid-to-fluid density ratio,
$\frac{\rho s}{\rho f}=1$, being $\rho s$ and $\rho f$ the solid and fluid
densities, respectively; the volume constraint constant is chosen as
$k_{v}=10$.

The quality of the present model predictions is measured against well known
published data used as benchmark [40]. Specifically, the evolution of the
Taylor deformation parameter in the y-z plane, Dyz, at x = 0.5 H is analyzed,
Dyz being defined as the weighted ratio between the major (M) and the minor
(m) diameter in the selected plane, namely D${}_{yz}=\frac{M-m}{M+m}$. The
distribution of Dyz for G = 0.0375, 0.075, and 0.15 is documented in Figure.7b
against the data by Lac et al. [40]. An excellent agreement is found for all
of the three values of the shear rate. The relative error with respect to
benchmark data ($\epsilon=\frac{D_{yz}-D_{yz,Lac}}{D_{yz,Lac}}$) is smaller
than 1% in all of the computations (see Figure.7c). Moreover, the distribution
of the relative variation of the enclosed volume ($\delta
v=\frac{v-v_{0}}{v_{0}}$) with respect to the undeformed configuration is
depicted in Figure.7d. Overall, $\delta v$ is well confined within the order
of magnitude of 10-6.

As previously observed by Walter et al [59] for an organosiloxane membrane, a
capsule with a nearly null Poisson ratio immersed in a simple shear flow will
present corrugations near its equator region due to the compression generated
as response to shear-induced stretching. Specifically, the membrane section
laying on the vorticity plane will stretch finding the steady equilibrium
between hydrodynamics and internal stresses; while the off-vorticity section
responds to this stretching with a compression. This mechanism brings to a
phenomenon largely studied in literature known as membrane buckling [60, 61,
1]. This shape-instability has a numerical origin, consistent with the null
bending resistance employed model, while the wrinkles location are physical
insofar as they occur in the compressed portions of such membranes. [2, 62]
Interestingly, the present dynamic-IB model is able to reproduce the membrane
buckling in agreement with the computations performed by Lac et al. [40] and
the experimental findings by Walter et al [59] (see Figure.7e). On the
contrary, the kinematic-IB scheme, in which the immersed body is deformed
advecting the Lagrangian points velocity with the fluid velocity, reads a
smooth velocity distribution in all planes and directions and responds with a
perfectly smooth membrane configuration (Figure.7e). Note that the kinematic-
IB and the dynamic-IB only differ by the technique used to displace Lagrangian
points, as discussed by Coclite et al. [26].

###  2.5 Membrane tumbling motion as a function of the bending resistance

Figure 8: Tumbling motion of a deformable prolate membrane. a. Schematic of
the physical problem. (b,c,d) Distribution of the relative variation of the
membrane surface as a function of G for Eb = 0.001 (b), 0.005 (c), and 0.01
(d). e. Relative variation of the membrane volume for different values of G
and Eb. f. Membrane angular velocity over time for Eb = 0.001, 0.005, 0.01.

The tumbling motion of an inertial membrane in a linear laminar flow is
considered. As depicted in Figure.8a a prolate capsule with major axis d is
placed in the center of a cubic box with side length H = 10 d. The minor axis
dimension of the capsule is equal to half the major axis length. d corresponds
to 10 $\Delta x$ while the capsule surface is discretized with $n_{L}=$ 10000
triangular elements. The fluid boundary conditions are kept identical to the
previous test case, so that: the imposed shear rate is
$\gamma=\frac{u_{max}}{H/2}$, with $u_{max}$ the velocity magnitude of moving
walls; the Reynolds number Re$=\frac{\gamma(d/2)^{2}}{\nu}$ = 0.01. The
membrane elastic constant $k_{s}$ is computed by the dimensionless shear rate
G $=\frac{\nu\rho\gamma d/2}{k_{s}}$. The bending resistance is related to the
strain stiffness and is regulated by Eb$=\frac{k_{b}}{k_{s,l}(d/2)^{2}}$. The
volume constraint constant is chosen as $k_{v}=10$. The tumbling motion is
systematically assessed in term of the membrane stiffness (G = 0.01, 0.025,
0.05, 0.075, and 0.1); the bending resistance, Eb = 0.001, 0.005, 0.01; and
the solid/fluid density ratio, $\rho s/\rho f=$ 1, 1.2, 1.5, 2, and 5.

Figure 9: Particle configurations. Membrane configurations for G = 0.01 and Eb
= 0.001 taken at $\theta/\pi=$ 0, 0.25, 0.5, 0.75, and 1 with superimposed the
contour of the module of the internal forces vectors indicating the exerted
hydrodynamics stresses.

First, a freely moving neutrally buoyant ($\rho s/\rho f=$ 1) membrane
immersed in a linear shear flow is considered. As depicted in Figures.8b-d
such elastic capsule stretches in a periodic fashion with a period that
depends on Eb. As clearly shown, the amplitude of the relative surface
stretching $\delta s=\frac{s-s_{0}}{s_{0}}$ depends on the stiffness.
Precisely, the peak surface strain grows with G while is retrained by higher
values of Eb. Such membranes strain the most for $\gamma t=$ 1.40, 1.35, and
1.30 with an half period of 3.705, 3.472, and 3.419 for Eb = 0.001, 0.005, and
0.01, respectively; indeed, they are about in the unperturbed configuration
(null relative surface deviation) at $\gamma t=$ 3.705, 3.472, and 3.419.
Interestingly, the distribution of the relative volume deviation $\delta v$
shows a different behavior. All of the distributions characterized by the same
value of the bending resistance overlap, as demonstrated in Figures.8e.
Moreover, we computed the tumbling velocity in the yz plane around their body
center. The distributions of $\omega=\omega_{LB}d/u_{max}$ in time are
reported in Figure.8f and demonstrate that there is a specific revolution
period for each value of the bending stiffness while within these values the
variation of G is uninfluential. Morevoer, such membranes experience the
largest strain for multiple of $\frac{\pi}{4}$ corresponding to the points
with the larger angular acceleration. Five configurations of the capsule
revolution are sketched in Figures.9 for G = 0.01 and Eb = 0.001. Note that,
the higher internal constitutive forces are registered far from the particle
equator due to larger hydrodynamics stresses (vectors in Figures.9). On the
contrary, in the equator region the membrane is responding to the
hydrodynamics stresses with almost null internal forces.

Figure 10: Effect of the membrane mass on the tumbling period. (a, b, c)
Distribution of the relative variation of the membrane surface for G = 0.01
and Eb = 0.001 (a), 0.005 (b), and 0.01 (c) as a function of the membrane
mass. b. Membrane angular velocity for G = 0.01 and Eb = 0.001, 0.005, 0.01
obtained varying the membrane mass.

Lastly, the effect of the membrane mass on the particle tumbling is considered
by varying the ratio between the solid and the fluid densities for G = 0.01
and Eb = 0.001, 0.005, and 0.01. Specifically, the membrane mass is
distributed considering a unitary thickness uniformly on the membrane surface.
Five different values of the density ratio between the solid shell and the
withstanding fluid are considered, namely $\rho_{s}/\rho_{f}=$ 1, 1.2, 1.5, 2,
and 5. Indeed, higher values of the solid shell mass imply higher membranes
inertia and consequently larger relative surface stretches (see Figure.10a-c).
Interestingly, for the investigated range of parameters the tumbling period is
not affected by particles inertia and it is only dependent on the bending
resistance coefficient. As the matter of fact, Figure.10d demonstrates that
particles with different inertia present the same angular velocity
distribution in time.

## Conclusions and future work

A combined dynamic-Immersed Boundary (IB)–Lattice Boltzmann (LB) method is
designed and employed originally to predict the interaction between low
Reynolds number flows and deformable or rigid objects with arbitrary shape and
buoyancy. A Moving Least Squares reconstruction was implemented to accurately
interpolate the velocity and force fields between Eulerian grid describing the
fluid evolution and the Lagrangian mesh defining the immersed object dynamics.

The quality of the present model predictions is measured by several validation
tests for both rigid and deforming particles. Firstly, the accuracy of the
present LB-IB scheme is computed by a grid-refinement study for the flow in a
lid-driven square cavity with and without a spherical obstacle immersed in.
This analysis demonstrates that the second order accuracy is conserved by the
IB treatment. Then, the settling under gravity of a sphere is assessed for
different values of the Reynolds and Froude numbers and sphere mass. Moreover,
the tumbling motion of inertial rigid particles with different shape is
considered. On one side, the Jefferey orbit of a prolate spheroid is recovered
for a prolate ellipsoid; on the other hand, the existence of a critical value
of the Reynolds number regulating the revolution period for an oblate spheroid
and a disk rotating under shear is demonstrated. The dynamics of a spherical
elastic capsule deforming under shear is computed and the obtained
distribution of the Taylor parameter for different mechanical stiffness are
validated against well known computations used as benchmark. Interestingly,
the present model is able to capture the characteristic membrane buckling
emerging in consolidated experiments and simulations for thin elastic
capsules. Weakly deformable prolate spheroids tumbling under shear are then
considered and their revolution period measured as a function of the
mechanical stiffness, the bending resistance and the membrane mass.
Interestingly, this systematic study revealed that, for the analyzed range of
parameters, the revolution period is influenced by only the bending
resistance. The membrane stretching stiffness and its mass are responsible for
the larger or smaller surface relative strain but do not seem to play any role
into the spheroids tumbling period itself.

Collectively, these data demonstrate the present formulation as a versatile
and reliable framework for the microfluidics fluid/structure interaction
problems. This approach can be readily applied to model the dynamics of
biological micro-vesicles and particles for biomedical applications as well as
to predict the transport of colloids into industrial microchips. The dynamic
behavior of capsule containing a different fluid with respect to the
suspension fluid is a challenging problem and will be the core of future
investigations.

## Acknowledgments

This investigation has been partially supported by the Italian Ministry of
Education, University and Research under the Programme “Department of
Excellence” Legge 232/2016 (Grant No. CUP - D94I18000260001)

## Competing Interests

The authors declare no competing interests.

## References

  * Barthès-Biesel [2009] Dominique Barthès-Biesel.  Capsule motion in flow: Deformation and membrane buckling.  _Comptes Rendus Physique_ , 10(8):764–774, 2009\.
  * Barthès-Biesel [2016] Dominique Barthès-Biesel.  Motion and deformation of elastic capsules and vesicles in flow.  _Ann. Rev. Fliud Mech._ , 48:25–52, 2016.
  * Freund [2014] Jonathan B. Freund.  Numerical simulation of flowing blood cells.  _Ann. Rev. Fliud Mech._ , 46:67–95, 2014.
  * Sebastian and Dittrich [2018] Bernhard Sebastian and Petra S. Dittrich.  Microfluidics to mimic blood flow in health and disease.  _Ann. Rev. Fliud Mech._ , 50:483–504, 2018.
  * Secomb [2017] Timothy W. Secomb.  Blood flow in the microcirculation.  _Ann. Rev. Fliud Mech._ , 49:443–461, 2017.
  * Guglietta et al. [2020] Fabio Guglietta, Marek Behr, Luca Biferale, Giacomo Falcucci, and Mauro Sbragaglia.  On the effects of membrane viscosity on transient red blood cell dynamics.  _Soft Matter_ , 2020.
  * Müller et al. [2014] Kathrin Müller, Dmitry A. Fedosov, and Gerhard Gompper.  Margination of micro- and nano-particles in blood flow and its effect on drug delivery.  _Scientific Reports_ , 4:4871 EP –, 05 2014.  URL https://doi.org/10.1038/srep04871.
  * Mollica et al. [2018] Hilaria Mollica, Alessandro Coclite, Marco E. Miali, Rui C. Pereira, Laura Paleari, Chiara Manneschi, Andrea DeCensi, and Paolo Decuzzi.  Deciphering the relative contribution of vascular inflammation and blood rheology in metastatic spreading.  _Biomicrofluidics_ , 12(4):042205, 2018.  doi: 10.1063/1.5022879.  URL https://doi.org/10.1063/1.5022879.
  * Coclite et al. [2018] A. Coclite, G. Pascazio, M.D. de Tullio, and P. Decuzzi.  Predicting the vascular adhesion of deformable drug carriers in narrow capillaries traversed by blood cells.  _Journal of Fluids and Structures_ , 82:638 – 650, 2018\.  ISSN 0889-9746.  doi: https://doi.org/10.1016/j.jfluidstructs.2018.08.001.  URL http://www.sciencedirect.com/science/article/pii/S0889974618301944.
  * F. et al. [2008] Gentile F., Ferrari M., and Decuzzi P.  The transport of nanoparticles in blood vessels: The effect of vessel permeability and blood rheology.  _Annals of Biomedical Engineering_ , 36(2):254–261, 2008.  doi: 10.1007/s10439-007-9423-6.
  * Giverso et al. [2010] C. Giverso, M. Scianna, L. Preziosi, N. Lo Buono, and A. Funaro.  Individual cell-based model for in-vitro mesothelial invasion of ovarian cancer.  _Mathematical Modelling of Natural Phenomena_ , 5(1):203–223, 2010.  doi: 10.1051/mmnp/20105109.
  * Krüger [2012] Timm Krüger.  _Computer simulation study of collective phenomena in dense suspensions of red blood cells under shear_.  Springer Science & Business Media, 2012.
  * Bialké et al. [2015] Julian Bialké, Thomas Speck, and Hartmut Löwen.  Active colloidal suspensions: Clustering and phase behavior.  _Journal of Non-Crystalline Solids_ , 407:367–375, 2015\.
  * Nazockdast and Morris [2016] Ehssan Nazockdast and Jeffrey F Morris.  Active microrheology of colloidal suspensions: Simulation and microstructural theory.  _Journal of Rheology_ , 60(4):733–753, 2016.
  * Wittkowski et al. [2017] Raphael Wittkowski, Joakim Stenhammar, and Michael E Cates.  Nonequilibrium dynamics of mixtures of active and passive colloidal particles.  _New Journal of Physics_ , 19(10):105003, 2017\.
  * Lauga [2016] Eric Lauga.  Bacterial hydrodynamics.  _Ann. Rev. Fliud Mech._ , 48:105–130, 2016.
  * Buzhardt and Tallapragada [2019] Jake Buzhardt and Phanindra Tallapragada.  Dynamics of groups of magnetically driven artificial microswimmers.  _Physical Review E_ , 100(3):033106, 2019.
  * Stricker [2017] Laura Stricker.  Numerical simulation of artificial microswimmers driven by Marangoni flow.  _Journal of Computational Physics_ , 347:467–489, 2017\.
  * Theers et al. [2018] Mario Theers, Elmar Westphal, Kai Qi, Roland G Winkler, and Gerhard Gompper.  Clustering of microswimmers: interplay of shape and hydrodynamics.  _Soft matter_ , 14(42):8590–8603, 2018.
  * Nitti et al. [2020] Alessandro Nitti, Josef Kiendl, Alessandro Reali, and Marco D de Tullio.  An immersed-boundary/isogeometric method for fluid–structure interaction involving thin shells.  _Computer Methods in Applied Mechanics and Engineering_ , 364:112977, 2020.
  * Patel et al. [2017] HV Patel, S Das, JAM Kuipers, JT Padding, and EAJF Peters.  A coupled volume of fluid and immersed boundary method for simulating 3d multiphase flows with contact line dynamics in complex geometries.  _Chemical Engineering Science_ , 166:28–41, 2017.
  * Lai et al. [2008] Ming-Chih Lai, Yu-Hau Tseng, and Huaxiong Huang.  An immersed boundary method for interfacial flows with insoluble surfactant.  _Journal of Computational Physics_ , 227(15):7279–7293, 2008.
  * Maiti and Nigam [2007] RN Maiti and KDP Nigam.  Gas- liquid distributors for trickle-bed reactors: a review.  _Industrial & Engineering Chemistry Research_, 46(19):6164–6182, 2007.
  * Sun and Sakai [2016] Xiaosong Sun and Mikio Sakai.  Numerical simulation of two-phase flows in complex geometries by using the volume-of-fluid/immersed-boundary method.  _Chemical Engineering Science_ , 139:221–240, 2016.
  * Peskin [2002] Charles S. Peskin.  The immersed boundary method.  _Acta Numerica_ , 11:479–517, 1 2002.  doi: 10.1017/S0962492902000077.
  * Coclite et al. [2019] A. Coclite, S. Ranaldo, M.D. de Tullio, P. Decuzzi, and G. Pascazio.  Kinematic and dynamic forcing strategies for predicting the transport of inertial capsules via a combined lattice Boltzmann immersed boundary method.  _Computers & Fluids_, 180:41–53, 2019.  ISSN 0045-7930.  doi: https://doi.org/10.1016/j.compfluid.2018.12.014.  URL http://www.sciencedirect.com/science/article/pii/S0045793018304304.
  * De Rosis et al. [2014] Alessandro De Rosis, Stefano Ubertini, and Francesco Ubertini.  A comparison between the interpolated bounce-back scheme and the immersed boundary method to treat solid boundary conditions for laminar flows in the lattice Boltzmann framework.  _Journal of Scientific Computing_ , 61(3):477–489, 2014.  ISSN 1573-7691.  doi: 10.1007/s10915-014-9834-0.  URL http://dx.doi.org/10.1007/s10915-014-9834-0.
  * Ye et al. [2014] Swe Soe Ye, Yan Cheng Ng, Justin Tan, Hwa Liang Leo, and Sangho Kim.  Two-dimensional strain-hardening membrane model for large deformation behavior of multiple red blood cells in high shear conditions.  _Theoretical Biology and Medical Modelling_ , 11(1):19, 2014.
  * Krueger et al. [2011] T. Krueger, 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.  _Computers and Mathematics with Applications_ , 61(12):3485 – 3505, 2011.  ISSN 0898-1221.  doi: http://dx.doi.org/10.1016/j.camwa.2010.03.057.
  * Melchionna [2011] Simone Melchionna.  A model for red blood cells in simulations of large-scale blood flows.  _Macromolecular Theory and Simulations_ , 20(7):548–561, 2011.
  * Zhang and Hisada [2004] Q. Zhang and T. Hisada.  Studies of the strong coupling and weak coupling methods in FSI analysis.  _International Journal for Numerical Methods in Engineering_ , 60(12):2013–2029, 2004.  doi: 10.1002/nme.1034.
  * de Tullio and Pascazio [2016] M. D. de Tullio and G. Pascazio.  A moving-least-squares immersed boundary method for simulating the fluid-structure interaction of elastic bodies with arbitrary thickness.  _Journal of Computational Physics_ , submitted, 2016.
  * Coclite et al. [2016] A. Coclite, M. D. de Tullio, G. Pascazio, and P. Decuzzi.  A combined lattice Boltzmann and immersed boundary approach for predicting the vascular transport of differently shaped particles.  _Computers & Fluids_, 136:260 – 271, 2016.  ISSN 0045-7930.  doi: http://dx.doi.org/10.1016/j.compfluid.2016.06.014.
  * Coclite et al. [2017] A. Coclite, H. Mollica, S. Ranaldo, G. Pascazio, M. D. de Tullio, and P. Decuzzi.  Predicting different adhesive regimens of circulating particles at blood capillary walls.  _Microfluidics and Nanofluidics_ , 21(11):168, 2017.  ISSN 1613-4990.  doi: 10.1007/s10404-017-2003-7.  URL https://doi.org/10.1007/s10404-017-2003-7.
  * Coclite and Gambaruto [2019] Alessandro Coclite and Alberto M. Gambaruto.  Injection of deformable capsules in a reservoir: A systematic analysis.  _Fluids_ , 4(3), 2019.  ISSN 2311-5521.  doi: 10.3390/fluids4030122.  URL https://www.mdpi.com/2311-5521/4/3/122.
  * Zou and He [1997] Qisu Zou and Xiaoyi He.  On pressure and velocity boundary conditions for the lattice Boltzmann BGK model.  _Physics of Fluids_ , 9(6):1591–1598, 1997.  doi: http://dx.doi.org/10.1063/1.869307.
  * Guo et al. [2011] Zhaoli Guo, Chuguang Zheng, and Baochang Shi.  Force imbalance in lattice Boltzmann equation for two-phase flows.  _Phys. Rev. E_ , 83:036707, Mar 2011.  doi: 10.1103/PhysRevE.83.036707.
  * Liu and Gu [2005] Gui-Rong Liu and Yuan-Tong Gu.  _An introduction to meshfree methods and their programming_.  Springer Science & Business Media, 2005.
  * Ten Cate et al. [2002] A Ten Cate, CH Nieuwstad, JJ Derksen, and HEA Van den Akker.  Particle imaging velocimetry experiments and lattice-Boltzmann simulations on a single sphere settling under gravity.  _Physics of Fluids_ , 14(11):4012–4025, 2002\.
  * Lac et al. [2004] E Lac, D Barthes-Biesel, NA Pelekasis, and J Tsamopoulos.  Spherical capsules in three-dimensional unbounded Stokes flows: effect of the membrane constitutive law and onset of buckling.  _Journal of Fluid Mechanics_ , 516:303–334, 2004.
  * Bhatnagar et al. [1954] P. L. Bhatnagar, E. P. Gross, and M. Krook.  A model for collision processes in gases. i. small amplitude processes in charged and neutral one-component systems.  _Phys. Rev._ , 94:511–525, May 1954.  doi: 10.1103/PhysRev.94.511.
  * Krüger et al. [2016] Timm Krüger, Halim Kusumaatmaja, Alexander Kuzmin, Orest Shardt, Goncalo Silva, and Erlend Magnus Viggen.  _The Lattice Boltzmann Method - Principles and Practice_.  Springer, 10 2016.  ISBN 978-3-319-44647-9.  doi: 10.1007/978-3-319-44649-3.
  * Guo et al. [2002] Zhaoli Guo, Chuguang Zheng, and Baochang Shi.  Discrete lattice effects on the forcing term in the lattice Boltzmann method.  _Phys. Rev. E_ , 65:046308, Apr 2002.  doi: 10.1103/PhysRevE.65.046308.
  * Junk and Yang [2009] Michael Junk and Zhaoxia Yang.  Pressure boundary condition for the lattice boltzmann method.  _Computers & Mathematics with Applications_, 58(5):922–929, 2009.
  * Junk and Yang [2011] Michael Junk and Zhaoxia Yang.  Asymptotic analysis of lattice boltzmann outflow treatments.  _Communications in Computational Physics_ , 9(5):1117–1127, 2011.
  * Favier et al. [2014] J. Favier, A. Revell, and A. Pinelli.  A lattice Boltzmann-immersed boundary method to simulate the fluid interaction with moving and slender flexible objects.  _Journal of Computational Physics_ , 261:145–161, 2014\.  doi: 10.1016/j.jcp.2013.12.052.
  * Vanella and Balaras [2009] Marcos Vanella and Elias Balaras.  A moving-least-squares reconstruction for embedded-boundary formulations.  _Journal of Computational Physics_ , 228(18):6617 – 6628, 2009.  ISSN 0021-9991.  doi: http://dx.doi.org/10.1016/j.jcp.2009.06.003.
  * Gelder [1998] Allen Van Gelder.  Approximate simulation of elastic membranes by triangulated spring meshes.  _Journal of graphics tools_ , 3(2):21–41, 1998\.
  * Tan et al. [2012] J. Tan, A. Thomas, and Y. Liu.  Influence of red blood cells on nanoparticle targeted delivery in microcirculation.  _Soft Matter_ , 8(6):1934–1946, 2012.  doi: 10.1039/c2sm06391c.
  * Körner et al. [2006] Carolin Körner, Thomas Pohl, Ulrich Rüde, Nils Thürey, and Thomas Zeiser.  Parallel lattice Boltzmann methods for CFD applications.  In _Numerical Solution of Partial Differential Equations on Parallel Computers_ , pages 439–466. Springer, 2006.
  * Schmieschek et al. [2017] Sebastian Schmieschek, Lev Shamardin, Stefan Frijters, Timm Krüger, Ulf D Schiller, Jens Harting, and Peter V Coveney.  Lb3d: A parallel implementation of the lattice-Boltzmann method for simulation of interacting amphiphilic fluids.  _Computer Physics Communications_ , 217:149–161, 2017.
  * Woodgate et al. [2018] Mark A Woodgate, George N Barakos, Rene Steijl, and Gavin J Pringle.  Parallel performance for a real time lattice Boltzmann code.  _Computers & Fluids_, 173:237–258, 2018.
  * Napolitano and Pascazio [1991] M Napolitano and G Pascazio.  A numerical method for the vorticity-velocity Navier-Stokes equations in two and three dimensions.  _Computers & Fluids_, 19(3-4):489–495, 1991\.
  * Abraham [1970] Farid F Abraham.  Functional dependence of drag coefficient of a sphere on Reynolds number.  _The Physics of Fluids_ , 13(8):2194–2195, 1970\.
  * Jeffery [1922] George Barker Jeffery.  The motion of ellipsoidal particles immersed in a viscous fluid.  _Proceedings of the Royal Society of London. Series A, Containing papers of a mathematical and physical character_ , 102(715):161–179, 1922.
  * Huang et al. [2012] Haibo Huang, Xin Yang, Manfred Krafczyk, and Xi-Yun Lu.  Rotation of spheroidal particles in Couette flows.  _Journal of Fluid Mechanics_ , 692:369–394, 2012.
  * Aidun et al. [1998] Cyrus K Aidun, Yannan Lu, and E-Jiang Ding.  Direct analysis of particulate suspensions with inertia using the discrete Boltzmann equation.  _Journal of Fluid Mechanics_ , 373:287–311, 1998.
  * Zettner and Yoda [2001] CM Zettner and M Yoda.  Moderate-aspect-ratio elliptical cylinders in simple shear with inertia.  _Journal of Fluid Mechanics_ , 442:241–266, 2001.
  * Walter et al. [2001] Anja Walter, Heinz Rehage, and Herbert Leonhard.  Shear induced deformation of microcapsules: shape oscillations and membrane folding.  _Colloids and Surfaces A: Physicochemical and Engineering Aspects_ , 183:123–132, 2001.
  * Li and Sarkar [2008] Xiaoyi Li and Kausik Sarkar.  Front tracking simulation of deformation and buckling instability of a liquid capsule enclosed by an elastic membrane.  _Journal of Computational Physics_ , 227(10):4998–5018, 2008.
  * Knoche and Kierfeld [2011] Sebastian Knoche and Jan Kierfeld.  Buckling of spherical capsules.  _Physical Review E_ , 84(4):046608, 2011.
  * Dupont et al. [2015] Claire Dupont, A-V Salsac, Dominique Barthes-Biesel, Marina Vidrascu, and Patrick Le Tallec.  Influence of bending resistance on the dynamics of a spherical capsule in shear flow.  _Physics of Fluids_ , 27(5):051902, 2015.