# Dynamics of finite-size spheroids in turbulent flow: the roles of flow
structures and particle boundary layers

Linfeng Jiang1  Cheng Wang1  Shuang Liu1
Chao Sun1,2<EMAIL_ADDRESS>Enrico Calzavarini3
<EMAIL_ADDRESS>1Center for Combustion Energy, Key
laboratory for Thermal Science and Power Engineering of Ministry of Education,
Department of Energy and Power Engineering, Tsinghua University, Beijing
100084, China 2Department of Engineering Mechanics, School of Aerospace
Engineering, Tsinghua University, Beijing 100084, China 3Univ. Lille, Unité de
Mécanique de Lille - J. Boussinesq - ULR 7512, F-59000 Lille, France

###### Abstract

We study the translational and rotational dynamics of neutrally-buoyant
finite-size spheroids in hydrodynamic turbulence by means of fully resolved
numerical simulations. We examine axisymmetric shapes, from oblate to prolate,
and the particle volume dependences. We show that the accelerations and
rotations experienced by non-spherical inertial-scale particles result from
volume filtered fluid forces and torques, similar to spherical particles.
However, the particle orientations carry signatures of preferential alignments
with the surrounding flow structures, which is reflected in distinct axial and
lateral fluctuations for accelerations and rotation rates. The randomization
of orientations does not occur even for particles with volume equivalent
diameter size in the inertial range, here up to $60\eta$ at
$Re_{\lambda}=120$. Additionally, we demonstrate that the role of fluid
boundary layers around the particles cannot be neglected to reach a
quantitative understanding of particle statistical dynamics, as they affect
the intensities of angular velocities, and the relative importance of tumbling
with respect to spinning rotations. This study brings to the fore the
importance of inertial-scale flow structures in homogeneous and isotropic
turbulence and their impacts on the transport of neutrally-buoyant bodies with
size in the inertial range.

###### keywords:

Inertial particles, Finite-size, Accelerations, Rotations, Turbulent flows

##  1 Introduction

The motion of neutrally-buoyant finite-size material particles in a turbulent
flow is typically regarded as the result of turbulent flow fluctuations
occurring at the scale of the particle. In other words, the particles are seen
as probes of a coarse-grained turbulent field. This implies that both the
scaling trends in the fluctuations of accelerations and rotational velocities
as a function of the particle size are interpreted in terms of the known
scaling of structure functions of turbulence in the framework of Kolmogorov
1941 (K41) phenomenological theory (Porta et al., 2001; Voth et al., 2002;
Qureshi et al., 2007; Brown et al., 2009; Calzavarini et al., 2009; Homann &
Bec, 2010). Such a picture has even been extended to the case of non-spherical
axisymmetric particles with various aspect ratios to elucidate the observed
different scaling trends in tumbling and spinning rotation rates (Parsa &
Voth, 2014; Bordoloi & Variano, 2017; Bounoua et al., 2018; Kuperman et al.,
2019; Pujara et al., 2019; Oehmke et al., 2021). This interpretation assumes
that neutrally-buoyant material particles are only weakly coupled to the flow,
so that the feedback effects (also referred as two-way coupling) stemming from
the formation of boundary layers and wakes around the particle do not have a
significant role, at least in statistical terms.
This line of reasoning is markedly distinct from the one that is used to grasp
the dynamics of small, i.e. sub- or Kolmogorov-scale ($\eta$) sized particles,
where instead it is crucial to take into account the local properties of the
flow, notably the fluid acceleration and the velocity gradient and its
topological properties along fluid Lagrangian trajectories (Chong et al.,
1990; Benzi et al., 2009) . It is key here to recognize that distinctive flow
structures exist at small scales, such as the filament-like vortices where
fluid tracers get trapped (Biferale et al., 2005; Bentkamp et al., 2019) or
from which tiny inertial particles are ejected (Bec et al., 2007) and to which
spheroidal particles align (Pumir & Wilkinson, 2011; Gustavsson et al., 2014;
Ni et al., 2014). These two pictures, the dissipative and the inertial-scale
ones, are somehow conceptually disconnected, with the only bridging idea being
the call into play of so-called Faxén corrections which account for effects
induced by the local curvature of the flow although at scales where the flow
field is smooth, typically $\leq 10\eta$, hence not yet in the inertial range
(Calzavarini et al., 2009) (We note that while simple for spherical particles
the Faxén corrections are more complex for other particle shapes, see Dolata &
Zia (2021)).
In this paper, we aim at ameliorating the above described conceptual models by
demonstrating that the dynamics of neutrally-buoyant inertial-sized particles:
i) is affected by inertial-range coherent flow structures with similar effects
as the ones produced by dissipative-scale turbulent structures on small
particles, and ii) is significantly impacted by the particle feedback on the
flow, especially if the considered observable is the angular rotation. This is
here achieved by means of a novel series of numerical experiments. First, we
perform extensive numerical simulations of fully-resolved homogeneous
isotropic turbulence (HIT) seeded with neutrally-buoyant spheroidal particles
with sizes ranging from the upper bound of the dissipative range ($\sim
10\eta$) to the far inertial range ($60\eta$). Second, we conduct a sequence
of numerical experiments in the same fluid flow setting but with virtual
particles, i.e. particle toy models, where we vary the type or/and level of
coupling between the flow and the particles. As we will show such virtual
particles reveals to be a crucial tool for an original physical interpretation
of the statistical properties of translation and rotation of the real material
particles in interaction with the turbulent flow environment.
The article is organized as follows. In the next section, we will describe in
detail the dynamical governing equations for the model system of particles in
turbulence and its numerical implementation. We will dedicate a section to the
description of the numerical experiments performed both with realistic and
virtual particles. We later present the results on the particle acceleration
and rotation rate statistics. In each of these two parts, we will guide the
reader through a physical interpretation of the results that is mainly done by
contrasting the simulations of real particles with virtual ones. Finally, the
conclusions and perspectives opened-up by the study are given.

$N^{3}$ |  $\eta/\Delta x$ |  $\tau_{\eta}/\Delta t$ |  $L/\eta$ |  $T_{L}/\tau_{\eta}$ |  $\lambda/\eta$ |  $Re_{\lambda}$
---|---|---|---|---|---|---
$512^{3}$ |  $1.5$ |  230 |  175 |  47 |  22 |  120
Table 1: Parameter of the NSE numerical simulations and relevant turbulence
scales. $N^{3}$: number of spatial grid points,
$\eta=(\nu^{3}/\epsilon)^{1/4}$: Kolmogorov dissipation length scale in grid
space units $\Delta x$, $\tau_{\eta}$: Kolmogorov time scale in time-step
units $\Delta t$, $L=u^{\prime 3}/\epsilon$: integral scale,
$T_{L}=L/u^{\prime}$: large-eddy turnover time, $\lambda=(15\nu
u^{\prime}/\epsilon)^{1/2}$: Taylor micro-scale, $Re_{\lambda}$: Taylor-
Reynolds number.

$\alpha$ $\frac{l}{\eta}$ $\frac{d}{\eta}$ $\frac{D_{v}}{\eta}$ 1/2 9.95 19.9
15.8 1/2 15.8 31.6 25.1 1/2 25 50.1 39.8 1/2 39.7 79.4 63 1 10.2 10.2 10.2 1
19.7 19.7 19.7 1 39.8 39.8 39.8 1 63 63 63 2 25.1 12.5 15.8 2 39.8 19.9 25.1 2
63.2 31.6 39.8 2 100 50 63 4 39.8 9.95 15.8 4 63.25 15.8 25.1 4 100.3 25.1
39.8 4 158.75 39.7 63

Figure 1: Parameter space of particle shapes and sizes in turbulent scale
units at $Re_{\lambda}=120$. The symbols indicate the numerically simulated
particle types. All particles are axisymmetric of diameter $d$ and length $l$,
aspect ratio $\alpha=l/d$, with $\alpha=1/2,1,2,4$, the volume-equivalent-
diameter is $D_{v}=(ld^{2})^{1/3}$. The lengths are here all expressed in
terms of the dissipative scale $\eta$. The Taylor micro-scale $\lambda$
(dashed lines) and the integral scale $L$ (dashed-dotted lines) are also
reported.

##  2 Methods

We first describe the basic physical model system and the equation of motion
ruling the dynamics of particles in a fluid flow. The numerical methods
adopted and details on their implementations are very briefly outlined.

###  2.1 The particles in turbulence model system and its numerical
implementation

The fluid flow by which particles are transported is modeled by the
incompressible Navier-Stokes equations (NSE) driven by an external random
large-scale statistically homogeneous and isotropic force with constant in-
time global energy input. This reads:

$\displaystyle\partial_{t}\textbf{u}+\textbf{u}\cdot{\nabla}\textbf{u}$
$\displaystyle=$ $\displaystyle-\rho^{-1}{\nabla}p+\nu\
\nabla^{2}\textbf{u}+\textbf{f},$ (1) $\displaystyle{\nabla}\cdot\textbf{u}$
$\displaystyle=$ $\displaystyle 0,$ (2)

where $\textbf{u}(\textbf{x},t)$ denotes the fluid velocity vector field,
$p(\textbf{x},t)$ is the hydrodynamic pressure, and the parameters are the
kinematic viscosity $\nu$ and the reference liquid density $\rho$. The vector
field f refers to the external force sustaining the HIT flow. The particle-
free turbulent flow intensity is here identified by a single dimensionless
parameter, the Reynolds number based on the Taylor microscale,
$Re_{\lambda}=\lambda u^{\prime}/\nu$, where $u^{\prime}=\sqrt{\langle
u_{i}u_{i}\rangle^{v,t}/3}$ is the root-mean-square (r.m.s.) velocity of the
turbulent flow ($\langle\ldots\rangle^{v,t}$ denotes here volume and time
average), $\lambda=\sqrt{15\nu u^{\prime 2}/\epsilon}$ is the Taylor length
scale and
$\epsilon=(\nu/2)\Sigma_{i,j}\langle(\nabla_{i}u_{j}+\nabla_{j}u_{i})^{2}\rangle^{v,t}$
is the mean global energy dissipation rate.

The translation and rotation of rigid-body material particle is governed by
the Newton-Euler equations (NEE):

$\displaystyle m_{p}\frac{\rm{d}\textbf{v}}{{\rm d}t}$ $\displaystyle=$
$\displaystyle\textbf{F}+\textbf{F}_{c},$ (3) $\displaystyle\frac{{\rm
d}\bm{\mathcal{I}}\bm{\Omega}}{{\rm d}t}$ $\displaystyle=$
$\displaystyle\textbf{T}+\textbf{T}_{c},$ (4)

where $\textbf{v}(t)=d\textbf{r}/dt$ and $\bm{\Omega}(t)$ are the particle
velocity and angular velocity vectors of a particle at position
$\textbf{r}(t)$ with mass $m_{p}=\rho_{p}V_{p}$ ($\rho_{p}$ the particle
density and $V_{p}$ the volume) and $\bm{\mathcal{I}}$ the moment of inertia
tensor. Note that since we consider neutrally buoyant and homogeneous
particles, $\rho_{p}=\rho$, this allows neglecting the buoyancy force in the
particle equation of motion as well any gravity-induced torque. Hence F and T
in (3)-(4) denote here the hydrodynamic force and torque acting on the
particle, which are formally written as:

F $\displaystyle=$ $\displaystyle\oint_{S_{p}}\bm{\sigma}\cdot\textbf{n}\ dS$
(5) T $\displaystyle=$
$\displaystyle\oint_{S_{p}}(\textbf{x}-\textbf{r})\times(\bm{\sigma}\cdot\textbf{n})\
dS,$ (6)

where $\bm{\sigma}=-p\bm{I}+\rho\nu(\nabla\textbf{u}+\nabla\textbf{u}^{T})$ is
the fluid stress tensor, $\textbf{x}-\textbf{r}$ the position vector relative
to the particle center and n the outward-pointing normal to the particle
surface $S_{p}$. While $\textbf{F}_{c}$ and $\textbf{T}_{c}$ are the
additional impulsive forces and torques related to the particle-particle
collisions (the so-called four-way interaction).The two-way coupling, meaning
the coupling between the fluid and the particles is provided by the
requirement of no-slip boundary condition on the particle surfaces (see below
for its numerical implementation).

For an axisymmetric particle with symmetry axis identified by a unit vector p,
as the ones we consider in this study, the angular velocity $\bm{\Omega}$ can
be decomposed into the tumbling rotation rate,
${\bm{\dot{}}{\textbf{p}}}={\bm{\Omega}}\times\textbf{p}$, the spinning
rotation rate, $\bm{\Omega}^{s}=({\bm{\Omega}}\cdot\textbf{p})\textbf{p}$
(Voth & Soldati, 2017). These will be key quantities in our analysis of
particle angular dynamics. Similarly, it will reveal useful in our analysis to
distinguish between axial and lateral accelerations,
$\bm{a}=\rm{d}\textbf{v}/\rm{d}t$, denoted
$\bm{a}^{\parallel}=(\bm{a}\cdot\textbf{p})\textbf{p}$ and
$\bm{a}^{\perp}=\bm{a}\times\textbf{p}$ respectively. The shape of an
axisymmetric spheroid is characterized by the aspect ratio, $\alpha=l/d$,
where $l$ and $d$ are the sizes of the symmetry axis and of the one
perpendicular to it. We will use the volume equivalent diameter
$D_{v}=(d^{2}l)^{1/3}=d\ \alpha^{1/3}$ as a parameter to compare the size of
particles with different shapes.

On the computational side: the NSE turbulent dynamics is here numerically
simulated by means of a Lattice Boltzmann Method (LBM) code, the ch4-project
(Calzavarini, 2019) which has been already extensively employed in studies of
Lagrangian tracers and point-like particle dynamics in turbulence (Mathai et
al., 2016; Calzavarini et al., 2020; Jiang et al., 2020, 2021). The code uses
a tri-linear scheme for the Lagrangian-Eulerian frame interpolations. The
computational domain is cubic with equispaced grid sizes, $128^{3}$ and
$512^{3}$, corresponding to $Re_{\lambda}=32,120$ (see Table 1 for details on
numerical accuracy properties and physical characteristics of the simulated
turbulent flow at with largest $Re_{\lambda}$). The NEE are numerically
integrated with a second-order Adams-Bashforth time-stepping scheme. The
particle-fluid two-way coupling is implemented by means of the immersed
boundary method (IBM), which has been widely used in particle-laden flows
(Peskin, 2002; Mittal & Iaccarino, 2005; Uhlmann, 2005; Ardekani et al., 2016;
Shen et al., 2021) also in combination to LBM e.g. in (Suzuki & Inamuro, 2011;
Do-Quang et al., 2014). The IBM enforce the no-penetration and no-slip
boundary conditions at the fluid-particle interface by means of a localized
feedback force, $\textbf{f}_{p}$, added to the NSE (1). Such $\textbf{f}_{p}$
term is also denoted as two-way coupling. In order to ensure high-accuracy for
the implementation of the no-slip fluid boundary condition at the particle
surface, we adopt the so-called IBM multi-forcing scheme with 5 step
iterations, see (Luo et al., 2007; Wang et al., 2021). The non-spherical
particle-particle interactions are implemented by means of soft-sphere
collision forces (Costa et al., 2015; Ardekani et al., 2016) and lubrication
corrections (Brenner, 1961; Cooley & O’Neill, 1969; Costa et al., 2015;
Ardekani et al., 2016) In order to limit the effect of the pair interactions
on the statistics of particle dynamics, the concentration of particles is kept
very low (the maximal volume fraction is $1\%$) and only up to 4 particles are
seeded in the simulation for the smallest particle size. All inter-particle
collision events are removed from the data sets employed for the statistical
analysis performed in the present study. To achieve this, the particle
acceleration/rotation data around the collisions are removed by a time window
whose width is of the order of the correlation time of the particle
acceleration (i.e. $\gtrsim\tau_{\eta}$, see also Fig. 10(b)). Thus, the
effect of pair interaction is considered negligible.

Prolate ($\alpha=2,4$), spherical ($\alpha=1$), and oblate ($\alpha=0.5$)
particles are investigated in the simulations, see Figure 1 for a
representation of the explored parameter space. For all considered cases, the
minimal linear size of the particle, either $d$ and $l$, is $>15\Delta x$
(grid units). Such choice guarantees a reasonable resolution of the particle
boundary layers whose inner (or viscous) thickness is expected to be of order
$\eta$ (although it might become thinner for particles with $D_{v}\gg\eta$ as
pointed out in (Cisse et al., 2013)). The particle Stokes number, defined as
the ratio of the particle viscous response time ($\tau_{p}=D_{v}^{2}/(12\nu)$)
with respect to the characteristic flow scale at the particle scale (here the
turbulent eddy-turnover time $\tau_{\ell}=\ell^{2/3}/\epsilon^{1/3}$), can be
expressed as $St_{\ell}=\frac{1}{12}(\frac{D_{v}}{\eta})^{4/3}$ (Xu &
Bodenschatz, 2008; Fiabane et al., 2012). It varies in the range 1.8 to 20.9.
Finally, for each studied particle type-case the duration simulation is more
than 100 large-eddy turnover times.

We have validated the numerical results by comparing the translational dynamic
of spheres ($\alpha=1$) at $Re_{\lambda}=32$ with a previous reference study
(Homann & Bec, 2010). This check validates our acceleration statistics (see
Appendix A for details). The rotational dynamics of the spheroids is instead
checked by solving the rotation rate of a prolate spheroid in a pure shear
flow, a condition for which the Jeffery solution (Jeffery, 1922) is available,
and by reproducing the numerical experiments in (Suzuki & Inamuro, 2011),
which include settling of ellipsoidal particles in a quiescent flow at finite
Reynolds numbers.

Figure 2: Sketches of the six numerical experiments performed: (a) real
particles RP, (b) real particles with fixed locations RPFL, (c) virtual
particles VP, (d) virtual particles with fixed locations VPFL, (e) Jeffery
fluid tracers JFT, and (f) volume average of Jeffery fluid tracers VA . The
light blue region denotes the fluid. The black shaded areas in panels (a) and
(b) represent the two-way coupling region. The points in panels (e) and (f)
represent the Jeffery fluid tracers. The tracers marked in black in panel (f)
are used to calculate the volume averages.

###  2.2 Numerical experiments with real and virtual particles

Taking advantage of the numerical approach, which allows to easily modify the
equations of motions of the advected particles and so to create toy-models, we
carry out several additional numerical experiments. This series of numerical
experiments is described in the following.

The first numerical experiment is the simulation of real particles (RP) which
obeys the model system introduced in the previous section. Their dynamics is
governed by the NSE and NEE with fluid-solid coupling computed by the IBM
algorithm, see Tab. 2, as also sketched in the cartoon of Figure 2a.

The second numerical experiment consists of spheroidal volume averages (VA) of
fluid flow, as shown in Figure 2f. We average the fluid acceleration and
vorticity over many ($N$) spheroids with the same diameter $D_{v}$ and aspect
ratio $\alpha$ as the ones used for RP. This reads:

$\displaystyle\langle a_{i}^{2}\rangle^{VA}$ $\displaystyle=$
$\displaystyle\frac{1}{N}\sum_{n=0}^{N}\left(\frac{1}{V_{p}}\int_{V_{p}}\frac{Du_{i}}{Dt}\
d^{3}x\right)^{2}$ (7) $\displaystyle\langle\Omega_{i}^{2}\rangle^{VA}$
$\displaystyle=$
$\displaystyle\frac{1}{N}\sum_{n=0}^{N}\left(\frac{1}{V_{p}}\int_{V_{p}}\omega_{i}\
d^{3}x\right)^{2}$ (8)

where $Du_{i}/Dt=\partial_{t}u_{i}+u_{j}\partial_{j}u_{i}$ is the i-th
acceleration component the fluid inside a spheroid of volume $V_{p}$. We note
that the orientation and location of VA are fixed. However, in a statistical
homogeneous and isotropic flow, as it is here, this is identical to performing
the averages at locations random-uniformly distributed in space locations and
over randomly oriented spheroids.

In the third numerical experiment, we seed virtual particles (VP) in the HIT
and evolve the trajectory and rotation by solving the NEE, as shown in Figure
2c. For these particles, the IBM algorithm in not activated and no feedback is
implemented on the flow by the virtual particles ($\textbf{f}_{p}=0$). Thus,
the translation and rotation of the particles are driven by the time
derivative of the momentum and angular momentum of the fluid inside the
virtual particles, which means the dynamics of VP is controlled by the fluid
at the scale of VP. More explicitly the force, $\textbf{F}_{VP}$, and torque,
$\textbf{T}_{VP}$, on the virtual particles are:

$\displaystyle\textbf{F}_{VP}$ $\displaystyle=$ $\displaystyle\rho
V_{p}\left<\frac{D\textbf{u}}{Dt}\right>^{v}\equiv\int_{V_{p}}\rho\frac{D\textbf{u}}{Dt}~{}d^{3}x$
(9) $\displaystyle\textbf{T}_{VP}$ $\displaystyle=$ $\displaystyle
V_{p}\left<\frac{d\textbf{L}}{dt}\right>^{v}\equiv\int_{V_{p}}(\bm{x}-\bm{r})\times\rho\frac{D\textbf{u}}{Dt}~{}d^{3}x$
(10)

where $D\textbf{u}/Dt$ is the fluid acceleration, L is the angular momentum,
and $\langle\rangle^{v}$ denotes the average over the fluid volume of the
virtual particle.

In the fourth numerical experiment, we adopt real particles with fixed
locations (RPFL). We fix the spatial position of the real particles and only
evolve the rotation by solving the Eq. (4). In this case, the IBM is still
implemented and the particles are two-way coupled with the flow, as shown in
Figure 2b.

The fifth numerical experiment contains the virtual particles with fixed
locations (VPFL), as shown in Figure 2d. We fixed the location of VP in the
turbulence and only evolve the rotation by solving the Eq. (4). No feedback on
the flow is imposed and the rotation of VPFL is driven by the time derivative
of the angular momentum of the fluid inside the virtual particles.

Finally, to better contrast the dynamics of finite-size particles with respect
to sub-Kolmogorov scale particles, we carry out a sixth numerical experiment
with point-like spheroidal tracers, as shown in Figure 2e. The governing
equations for the point-like spheroidal tracers with position,
$\textbf{r}(t)$, and orientation, $\textbf{p}(t)$, are given by

$\displaystyle\dot{\textbf{r}}$ $\displaystyle=$
$\displaystyle\textbf{u}(\textbf{r}(t),t),$ (11)
$\displaystyle\dot{\textbf{p}}$ $\displaystyle=$
$\displaystyle\bm{\Omega}\times\textbf{p}$ (12) $\displaystyle{\bm{\Omega}}$
$\displaystyle=$
$\displaystyle\frac{1}{2}{\bm{\omega}}(\textbf{r}(t),t)+\Lambda\
\textbf{p}\times\mathcal{S}(\textbf{r}(t),t)\textbf{p}$ (13)

Here, $\textbf{u}(\textbf{r}(t),t)$ is the fluid velocity at the particle
position at time $t$. The vector
${\bm{\omega}}(\textbf{r}(t),t)=\nabla\times\textbf{u}$ is the fluid vorticity
and $\mathcal{S}=({\nabla}\textbf{u}+{\nabla}\textbf{u}^{T})/2$ is the strain-
rate matrix of the fluid velocity gradient tensor, ${\nabla}\textbf{u}$, at
the particle position, and $\Lambda=\frac{\alpha^{2}-1}{\alpha^{2}+1}$ is the
shape parameter of particles. The equations of rotation Eq. (12)(13), called
the Jeffery angular equations (Jeffery, 1922; Byron et al., 2015), have been
extensively used for the study of the dynamics of tiny spheroidal tracers
(Parsa et al., 2012; Gustavsson et al., 2014; Calzavarini et al., 2020; Jiang
et al., 2021, 2020). The aspect ratios of the point-like spheroidal tracers
(hereafter called Jeffery tracers) are chosen equal to the ones of RP. For
each aspect ratio, we seed $10^{5}$ Jeffery tracers in the flow to have
converged statistics. The numerical condition of the flow is identical to the
simulation of RP. We use $\langle\rangle^{tracer}$ to denote the temporal and
ensemble average of Jeffery tracers. Note that the translational dynamics of
the particles is the one of perfect Lagrangian tracers, Eq. (11),
independently of their geometric aspect ratio (i.e. the drag force is here
neglected). For this reason, we denote the particles as Jeffery fluid tracers
(JFT). For better clarity all the model systems described above are summarized
in Table 2.

n. | Numerical Experiments | Model Equations
---|---|---
1 | Real Particles (RP) | NSE + NEE + IBM
2 | Volume Averages (VA) | NSE
3 | Virtual Particles (VP) | NSE + NEE
4 | Real Particles with Fixed Locations (RPFL) | NSE + EE + IBM
5 | Virtual Particles with Fixed Locations (VPFL) | NSE + EE
6 | Jeffery Fluid Tracers (JFT) | NSE + tracer eq.+ Jeffery eq.
Table 2: Summary of the six performed numerical experiments with the
corresponding involved dynamical equations: Navier-Stokes equation (NSE)
(1)-(2), Newton-Euler equation (NEE) (3)-(4); Euler equation (EE) (4) only
rotation for particles with fixed location; No-slip boundary conditions (i.e.
two-way coupling) at particle-fluid interface implemented via the immersed
boundary method (IBM); tracer equation (11), Jeffery equation (12)-(13).
Figure 3:  (a) Example of trajectory for a spheroid with an equivalent
diameter $D_{v}\approx 25~{}\eta$ and aspect ratio $\alpha=2$ at
$Re_{\lambda}=120$. The slices show the local dissipation rate
$\epsilon_{\textbf{x},t}$, normalized by $\epsilon$, at the planes through the
spheroid center at time $t=0,8.7,17.4,26.1,39.1$ in $\tau_{\eta}$. (b)
Azimuthally and temporally averaged dissipation rate $\overline{\epsilon}$,
normalized by $\epsilon$, around the spheroid with a diameter $D_{v}\approx
25~{}\eta$ and aspect ratio $\alpha=2$ at $Re_{\lambda}=120$. The black line
denotes the contour line of $\rm log_{10}(\overline{\epsilon}/\epsilon)=0$.
(c) Azimuthally and temporally averaged square of $\bm{u}-\bm{v}$, normalized
by the mean square of $\bm{u}$, around the spheroid with diameter
$D_{v}=25~{}\eta$ and aspect ratio $\alpha=2$ at $Re_{\lambda}=120$. The black
line denotes the contour line of
$\overline{(\bm{u}-\bm{v})^{2}}/\langle\bm{u}^{2}\rangle^{v,t}=1$. $x_{p}$ and
$y_{p}$ represent the coordinates in the particle frame.

##  3 Results

###  3.1 Particle trajectories and fields visualizations

In Figure 3a a visualization example of the trajectory of a spheroid with
$\alpha=2$ and an equivalent diameter $D_{v}=25\eta$ over a relatively short
time span, $\sim 40$ dissipative times ($\tau_{\eta}$), obtained from the
simulation. The energy dissipation rate field in the fluid, $\epsilon$, is
represented in colour over planes intersecting the spheroid. Such a quantity
that depends on the square of velocity gradients allows to foreground the
effect the fluid-particle coupling. We observe that just a thin layer of large
dissipation appears around the particle (see the zoom-in region of the slice
at $t_{1}=8.7\tau_{\eta}$), while no wakes are observed at larger distances.
This mild effect induced by the particle presence in the flow is better
quantified by computing temporal fluid averages in the vicinity of the
particle. This is numerically performed by seeding the flow with point-like
Lagrangian probes whose positions are fixed in the particle frame of
reference. Figure 3b shows the temporally and azimuthally averaged local
dissipation rate around the spheroid ($\alpha=2$ and $D_{v}/\eta=25$). It
appears that the dissipation rate close to the particle is enhanced by one
order of magnitude with respect to the far-field value, furthermore such an
increase is uniformly distributed around the particle. Indeed, the contour
line in the figure (black colour) marking where the particle influence reaches
a saturation, indicates that the finite-size particle modifies the flow over
the spheroidal region of about 1.5 times of the linear size of the particle.
Similar results for finite-size spheres have been observed in turbulence at
$Re_{\lambda}=160$ (Cisse et al., 2013). The presence of the particle
significantly modifies also the kinetic energy around the particle, the
kinetic energy of the fluid is null at the particle interface and increases to
the unperturbed value at about twice the particle size $D_{v}$, see Figure 3c.
The same figure suggests that the slip velocity of neutrally buoyant particles
is small, in particular $|\bm{u}-\bm{v}|\sim\bm{u}_{rms}$. This implies that
the particle Reynolds number $Re_{p}=\langle|\bm{u}-\bm{v}|\rangle
D_{v}/\nu\sim(D_{v}/\lambda)Re_{\lambda}\sim O(10^{2})$. If the flow
surrounding the particle were steady and laminar this would lead to a steady
or periodically oscillating wake. In a turbulent flow, where the direction of
motion changes frequently, such an effect is weakened. For the cases of
spherical particles, this agrees with detailed study by Cisse et al. (2013).
Overall the presence of a neutrally-buoyant finite-sized non-spherical
particle has weak feedback on the turbulent flow field. This is certainly very
specific to the case of inclusions with the same density, the situation is
very different in the case of heavy particles (Brändle de Motta et al., 2016)
or large bubbles (Mathai et al., 2018). This suggests that the picture of
neutral inertial particles as essentially passive objects might be sound. We
will now address this in more detail by looking at two quantities that are
most sensitive to the fluid-particle one-way coupling, i.e., the fluid
acceleration and later the angular velocity.

###  3.2 Particle Acceleration statistics:
the VP model explains the size and shape dependencies in RP

We begin discussing the dependency of the particle acceleration intensity of
large spheroids as a function of their shape and size and the physical
mechanisms responsible for the observed dependencies. Figure 4a reports the
temporal and ensemble average (denoted in the following
$\langle\ldots\rangle$) of the single-cartesian-component acceleration
variance $\langle a_{i}^{2}\rangle$ of real particles (RP), normalized by the
acceleration variance of fluid tracers $\langle a_{i}^{2}\rangle^{tracer}$, as
a function of the equivalent diameter in turbulent dissipative scale units
$D_{v}/\eta$. Two observations are here in order. First, the acceleration
variance decreases as the size of the spheroids increases. In the range
$D_{v}/\eta\in[10,63]$ we measure power-law $\sim(D_{v}/\eta)^{-1.0\pm 0.1}$
(at $Re_{\lambda}=120$). Second, we remark a good collapse onto a single curve
of all the data points corresponding to different aspect ratios (with
$\alpha=[1/2,1,2,4]$). The observed scaling behaviour is consistent with the
one observed for spherical particles at similar $Re_{\lambda}$ values (Volk et
al., 2011). Such a scaling deviates from the $(D_{v}/\eta)^{-2/3}$ that can be
guessed on the basis of dimensional reasoning relying on the Kolmogorov (K41)
theory of turbulence (Voth et al., 2002; Qureshi et al., 2007). This reasoning
assumes that the particles are only one-way coupled to the fluid, i.e., they
do not affect the fluid flow significantly through particle-fluid
interactions. However, the existence of such a deviation has been already
reported in studies on spherical particles. In particular, it is known that
the scaling exponent has a weak Reynolds number dependence (Voth et al., 2002;
Qureshi et al., 2007). At lower $Re_{\lambda}$ than the one used in this study
the observed exponent is closer to the value $-4/3$ as observed in Ref. (Bec
et al., 2007). The role of the pressure scaling dependence on $Re_{\lambda}$
has been conjectured to be responsible for such a behaviour (Bec et al.,
2007). Our results at low Reynolds number, $Re_{\lambda}=32$ (see Figure 10 in
Appendix A) confirms such a dependence. On the other hand, the collapse of
acceleration variance for different particle shapes (Figure 4a) is a novel
feature. It is remarkable because it implies that the hydrodynamical forces
involved in the translational dynamics of spheroids are overall determined by
the particle volume.

Figure 4: (a) Normalized single component acceleration variance as a function
of the normalized particle size $D_{v}/\eta$ at $Re_{\lambda}=120$ for real
particles (RP: solid symbols) and volume averages (VA: open symbols). The
dashed lines are the fitting lines for spheres of RP (exponent $-1.0\pm 0.1$).
Inset in (a) shows the ratio of acceleration variance of VA to acceleration
variance of RP, $\langle a_{i}^{2}\rangle^{VA}/\langle a_{i}^{2}\rangle$, as a
function of $D_{v}/\eta$. (b) Normalized componentwise acceleration variance
as a function of the normalized particle size $D_{v}/\eta$ at
$Re_{\lambda}=120$ for real particles (RP: solid symbols) and virtual
particles (VP: open symbols). The dashed lines are the fitting lines for
spheres of RP. Inset in (b) shows the ratio of acceleration variance of VP to
acceleration variance of RP, $\langle a_{i}^{2}\rangle^{VP}/\langle
a_{i}^{2}\rangle$, as a function of $D_{v}/\eta$. (c) The ratio of axial
variance $\langle(\bm{a}\cdot\bm{p})^{2}\rangle$ to lateral variance
$\langle|\bm{a}\times\bm{p}|^{2}\rangle$ of particle acceleration as a
function of aspect ratio at $Re_{\lambda}=120$ for RP (triangles), VA
(square), and VP (circles), respectively. The black dash-dot line denotes the
results of Jeffery tracers. The color represents the normalized particle size
$D_{v}/\eta$ in the logarithmic scale.

In order to further substantiate these findings, we show the results of the
numerical experiment with fluid volume averages (VA) and the so-called virtual
particles (VP) in Figure 4(a and b) (empty symbols). Similar scaling relations
of acceleration variances are observed both for VA and VP in the inertial
regime, whereas the magnitudes of the acceleration variances of VP and VA are
not equal. In particular, the amplitude of accelerations fluctuations of VP
particles is closer to the ones of RP particles (see insets in Figure 4(a and
b) where their ratio is plotted). Note that the acceleration variances of VA
and VP represent respectively the volume averaging of the fluid flow
acceleration in the fixed frame (Eulerian perspective) and in an advected
frame (Lagrangian perspective). Therefore, these results stress that the
scaling behaviour of the translational dynamics is dominated by the volume
averaging instead of other surface forces such as e.g. the drag. However, an
interesting question is: Why VP are capable to better approximate the
statistics of real particles than VA? This might be due to preferential
sampling of fluid structures along their Lagrangian trajectories, a fact that
we tend to exclude because we know that these particles do not form clusters,
or this might come from a preferential orientation in space (that is dynamical
for the RP and VP cases while absent for VA). This is at odds with the
expectations that particles of sizes well into the inertial range, shall
become randomly oriented due to the averaging of turbulent fluid fluctuations
that results from the spatial volume filtering.

In order to understand if any correlation exists between the particle
acceleration and its orientation, we look at the ratio of the particle axial
acceleration variance $\langle(\bm{a}\cdot\bm{p})^{2}\rangle$ to the lateral
acceleration variance $\langle|\bm{a}\times\bm{p}|^{2}\rangle$ as a function
of their aspect ratio, Figure 4c. For spherical particles ($\alpha=1$), no
preferential alignment is expected between a and p, even for various sizes,
due to the full rotational symmetry of the body shape. For the non-spherical
case, a value of the ratio different from 0.5 marks the existence of
preferential orientation.

Figure 5: PDFs of $|\rm cos(\theta_{\textbf{a},\textbf{p}})|$ for RP spheroids
with different aspect ratios: (a) $\alpha=0.5$, (b) $\alpha=1$, (c)
$\alpha=2$, and (d) $\alpha=4$ and for Jeffery fluid tracers.

We find that the oblate particles preferentially align with their acceleration
($\langle(\bm{a}\cdot\bm{p})^{2}\rangle/\langle|\bm{a}\times\bm{p}|^{2}\rangle>0.5$)
and the prolate particles are predominantly perpendicular to the acceleration
($\langle(\bm{a}\cdot\bm{p})^{2}\rangle/\langle|\bm{a}\times\bm{p}|^{2}\rangle<0.5$).
It also turns out that the alignment of large particles is identical to that
of Jeffery fluid tracers (JFT), which indicates that the alignment is
insensitive to the particle’s size. This is further confirmed by the
measurements of the PDFs of $|cos(\theta_{\bm{a},\bm{p}})|$ for various sizes
and shapes (see Figure 5), where $\theta_{ap}$ defines the angle between the
acceleration vector $\bm{a}$ and orientation vector $\bm{p}$. The particular
trend of
$\langle(\bm{a}\cdot\bm{p})^{2}\rangle/\langle|\bm{a}\times\bm{p}|^{2}\rangle$
as a function of $\alpha$ can then be understood in terms of the knowledge of
Jeffery fluid tracers. It has been shown experimentally that prolate fluid
tracers have accelerations preferentially perpendicular to the local vorticity
(Liberzon et al., 2012). It has also been theoretically demonstrated that
prolate particles preferentially align with vorticity while oblate particles
are preferentially perpendicular to it (Pumir & Wilkinson, 2011). The
combination of these two pieces of evidence explain qualitatively the observed
alignment, depending on the aspect ratio, between the particle orientation and
the acceleration. The existence of these alignments indicates that the
translation of finite-size neutrally-buoyant particles in a flow is, in this
respect, surprisingly similar to that of point-like inertialess axisymmetric
fluid tracers, or conversely that geometrical structure of flow at the
inertial scale shows a similar structure with the one of the velocity
gradient, which occurs at the dissipation scale.

We now turn the attention to the results for VA and VP. For VA, the values of
the ratio
$\langle(\bm{a}\cdot\bm{p})^{2}\rangle/\langle|\bm{a}\times\bm{p}|^{2}\rangle$
are close to 0.5 for all considered sizes and aspect ratios, which is expected
because VA is not aligned with the corresponding acceleration vector.
Remarkably, it is found that the ratio
$\langle(\bm{a}\cdot\bm{p})^{2}\rangle/\langle|\bm{a}\times\bm{p}|^{2}\rangle$
of VP shows excellent agreement with that of RP (Figure 4c). This suggests
that the effect induced by the two-way coupling on the particle preferential
alignment, if present, is weak. The similarity between RP and VP implies that
the neutrally-buoyant finite-size particle behaves roughly like a mass of
fluid with the same particle size and shape. Furthermore, the size
independence of the alignment indicate that the Lagrangian behavior of a fluid
mass at the inertial scale is similar to the flow at the dissipative scale.
The neutrally-buoyant finite-size particles can then be regarded as probes of
the Lagrangian coarse-grained field of turbulence.

In the next section, we will check if the similarity between real and virtual
particles still holds when the rotational properties of particles are
considered.

Figure 6: (a) Mean squared angular velocity, normalized by the results of
Jeffery tracers with corresponding aspect ratio, as a function of $D_{v}/\eta$
at $Re_{\lambda}=120$. The open symbols denote the results of VA for spheres.
The dashed line denotes the fitting line for spheres with scaling exponent
-4/3. Inset shows the ratio of angular velocity variance of VA to that of RP,
$\langle\Omega_{i}\Omega_{i}\rangle^{VA}/\langle\Omega_{i}\Omega_{i}\rangle$,
as a function of $D_{v}/\eta$ for spheres. (b) Mean squared angular velocity,
normalized by the results of Jeffery tracers with corresponding aspect ratio,
as a function of $D_{v}/\eta$ at $Re_{\lambda}=120$. The open symbols denote
the results of VP. The dashed line denotes the fitting line for spheres with
scaling exponent -4/3. Inset shows the ratio of angular velocity variance of
VP to that of RP,
$\langle\Omega_{i}\Omega_{i}\rangle^{VP}/\langle\Omega_{i}\Omega_{i}\rangle$,
as a function of $D_{v}/\eta$.

###  3.3 Particle rotation statistics:
the VP model fails to describe the shape dependence in RP

As already mentioned in Sec. 1, the rotations of tiny, $D_{v}\lesssim\eta$,
neutrally-buoyant particles in turbulent flows are determined by the evolution
of the fluid velocity gradient along the Lagrangian particle trajectories. The
Jeffery equation (12)- (13) well captures the angular dynamics of such
particles. This has been recently proven in a series of combined experimental
and numerical studies (Parsa et al., 2012; Byron et al., 2015; Calzavarini et
al., 2020; Jiang et al., 2021). When it comes to larger particles, in
particular to inertial scale particles, the physics is more complicated due
to: i) the non-smooth character of the flow field at the particle scale, ii)
the spatio-temporal filtering of flow fluctuations associated with the
particle size and inertia and iii) the two-way coupling with the fluid flow.
This means that not only the shape (parametrized by $\alpha$) but also the
particle sizes become relevant parameters of the problem. This is easily
grasped in Figure 6, where we show the mean square angular velocity,
normalized by the value of the mean square angular velocity of point-like
Jeffery fluid tracers (JFT) with corresponding aspect ratio, as a function of
the normalized diameter $D_{v}/\eta$ for RP (solid symbols). As the size
increases, we observe that the mean square angular velocity decreases and
scales approximately as
$\langle{\Omega}_{i}{\Omega}_{i}\rangle\sim(D_{v}/\eta)^{-4/3}$ for $D_{v}\geq
10\eta$. Since the angular velocity of JFT is nearly constant at changing the
aspect ratio (Byron et al., 2015), the decrease of angular velocity is
exclusively an inertial term. Even if an offset between spherical and non-
spherical particles is observed, the good collapse for prolate and oblate
spheroids with different shapes substantiates the equivalent diameter $D_{v}$
as the relevant characteristic scale also for the rotational dynamics. We
recall that the scaling exponent -4/3 can be dimensionally deduced using the
relations of the K41 turbulence theory. Specifically by assuming that the
rotation rate is related to the eddy turnover time at the particle scale, as
confirmed in experiments with long slender fibers (Parsa & Voth, 2014; Bounoua
et al., 2018) and large cylinders (Bordoloi & Variano, 2017; Pujara et al.,
2018).

The VA and VP numerical experiments exhibit also a similar behaviour as RP,
i.e. a decreasing trend with $D_{v}$ consistent with the $-4/3$ power-law.
However, the magnitude of the mean squared angular velocity is larger for VA
and VP than that of RP (see again Figure 6(a,b)). This suggests that
rotational dynamics of finite-size neutrally-buoyant spheroids (RP) is also
dominated by the volume filtering effect, as proven by its scaling-behavior
with respect to $D_{v}$, while its amplitude must be affected by other
physical effects (to be discussed later). From the insets of Figure 6, we note
that the angular velocity variance discrepancy between VA and RP and VP and RP
is larger than the discrepancy of acceleration variances. In conclusion, as
far as angular velocity is concerned we can not say that VPs are a good
approximation of RP particles. This discrepancy is even more important when
one looks at the ratio between the intensity of axial as compared to lateral
rotational components, that is to say, tumbling and spinning, as we will
discuss in the following.

Figure 7:  (a) The ratio of spinning rate squared to tumbling rate squared,
$\langle\Omega_{i}^{s}\Omega_{i}^{s}\rangle/\langle\dot{p}_{i}\dot{p}_{i}\rangle$,
as a function of aspect ratio $\alpha$ for RP. The color bar represents the
normalized particle size $D_{v}/\eta$ in the logarithmic scale. Spheres denote
the results for Jeffery tracers. Leftward triangles denote the results at
$Re_{\lambda}=32$ and rightward triangles denote the results at
$Re_{\lambda}=120$. The arrows denote the direction that the size of particles
increases. The dashed line shows the results of tracers at $Re_{\lambda}=433$
in Ref. (Byron et al., 2015). (b) The ratio of mean spinning rate squared to
mean tumbling rate squared,
$\langle{\Omega}_{i}^{s}{\Omega}_{i}^{s}\rangle/\langle\dot{p}_{i}\dot{p}_{i}\rangle$,
as a function of size $D_{v}/\eta$.(c) The cartoon illustrates that the weight
of the tumbling rate in the angular velocity of prolate spheroids becomes
larger as the particle size increases. (d) The ratio of spinning rate squared
to tumbling rate squared,
$\langle\Omega_{i}^{s}\Omega_{i}^{s}\rangle/\langle\dot{p}_{i}\dot{p}_{i}\rangle$,
as a function of aspect ratio $\alpha$ for VA and VP. The color bar represents
the normalized particle size $D_{v}/\eta$ in the logarithmic scale.

In Figure 7a, we show the ratio
$\langle{\Omega}_{i}^{s}{\Omega}_{i}^{s}\rangle/\langle\dot{p}_{i}\dot{p}_{i}\rangle$
as a function of aspect ratio $\alpha$ for particles of different sizes. We
first discuss the behaviour of sub-Kolomogorov scale particles, whose dynamics
is well described by the Jeffery equation. In this case, if the particle
orientation vector were uncorrelated with the fluid velocity gradients, and p
is an isotropic vector, upon time and ensemble averaging Eq. (13), one
obtains:

$\frac{\langle{\Omega}_{i}^{s}{\Omega}_{i}^{s}\rangle}{\langle\dot{p}_{i}\dot{p}_{i}\rangle}=\frac{\frac{1}{12}\langle\omega^{2}\rangle}{\frac{1}{6}\langle\omega^{2}\rangle+\frac{1}{5}\Lambda^{2}\langle\mathcal{S}:\mathcal{S}\rangle}=\frac{5}{10+6\Lambda^{2}}.$
(14)

Note that the last equality follows from the statistical properties of the HIT
flow, where
$\frac{1}{2}\langle\omega^{2}\rangle=2\langle\mathcal{S}:\mathcal{S}\rangle=\epsilon/\nu$
(Parsa et al., 2012; Byron et al., 2015). This would produce only a mild
variation of the spinning to tumbling ratio across the whole range of aspect
ratios, as it is shown in Figure 7a (blue solid curve). However, it is known
that small particles develop strong correlations with velocity gradients that
are responsible for a sharp deviation from this prediction. The present
Jeffery’ models simulations clearly show this behaviour in Figure 7a (spheres
symbols). Incidentally, we note that our results at $Re_{\lambda}=120$ agree
with previous numerical study (dashed line) at $Re_{\lambda}=433$ (Byron et
al., 2015), underlying the weak $Re_{\lambda}$ dependence of this phenomenon.
For spheres, the angular velocity is isotropic and the ratio
$\langle{\Omega}_{i}^{s}{\Omega}_{i}^{s}\rangle/\langle\dot{p}_{i}\dot{p}_{i}\rangle$
stays around 0.5 for all cases. A marked asymmetry exists between oblate and
prolate particles, while the first tumble intensely the seconds are spin
dominated, as they tend to align with the vorticity vector (Voth & Soldati,
2017; Parsa et al., 2012). Futhermore, a saturation of these behaviours is
observed in both cases for extreme aspect ratios, i.e. $\alpha\lesssim 1/4$
and $\alpha\gtrsim 1/4$.

What happens when particles have a finite size? A naive K41-based prediction
would suggest that the tumbling scales with the length
$\langle\dot{p}_{i}\dot{p}_{i}\rangle\sim l^{-4/3}$ and that the spinning
scale with the $\langle{\Omega}_{i}^{s}{\Omega}_{i}^{s}\rangle\sim
0.5d^{-4/3}$ (where the 0.5 factor comes for the matching with the spherical
case) so that

$\frac{\langle{\Omega}_{i}^{s}{\Omega}_{i}^{s}\rangle}{\langle\dot{p}_{i}\dot{p}_{i}\rangle}\sim\frac{\alpha^{4/3}}{2}.$
(15)

However, this prediction falls largely off the result of the simulations,
where we see a marked decreasing $D_{v}$ dependence of the ratio of the mean
square spinning to tumbling rate (see Figure 7a and b). Larger particles in
proportion spin less and tumble more. This is indeed true for all aspect
ratios except for spheres (see Figure 7 d). It has been pointed out by (Oehmke
et al., 2021), that these types of dimensional reasoning ignore the
decorrelation of the velocity field along the particle surface. The rotation
rates stem from the integral of the velocity gradients along with the
particles and as such e.g. the spinning rate can be much attenuated in long
fibers ($\alpha\gg 1$). This integral effect is possibly more important when
both the diameter and the length are much larger than the dissipation scale.
As a consequence, the scaling behaviors of the spinning and tumbling rate can
deviate from K41-based scalings as another length scale is involved, and the
ratio
$\langle{\Omega}_{i}^{s}{\Omega}_{i}^{s}\rangle/\langle\dot{p}_{i}\dot{p}_{i}\rangle$
becomes to be a function of the particle size.

The experiments by Oehmke et al. (2021), the only to date to have measured
both the spinning and tumbling of finite-size particles, find that for prolate
fibers the variance of the spinning rate is always larger than the variance of
the tumbling rate. These authors comment that the observed behaviour
“surprisingly resembles that observed previously for sub-Kolmogorov fibers”.
These observations are fully confirmed by our simulations (see Appendix B for
a more detailed comparison). However, we additionally observe a new feature:

The ratio
$\langle{\Omega}_{i}^{s}{\Omega}_{i}^{s}\rangle/\langle\dot{p}_{i}\dot{p}_{i}\rangle$
for oblate and prolate spheroids shows systematic decreasing trends as the
particle equivalent size increases, and this for both $Re_{\lambda}=32$ and
$Re_{\lambda}=120$.

For oblate spheroids ($\alpha=0.5$), the ratio becomes vanishingly small at
increasing their size, which indicates that the oblate spheroids hardly spin
when their size is large compared to the dissipation scale. Similarly, prolate
spheroids of large size reduce their spinning to the point that it becomes
smaller than tumbling ( Figure 7a).

Figure 7b shows the ratio
$\langle{\Omega}_{i}^{s}{\Omega}_{i}^{s}\rangle/\langle\dot{p}_{i}\dot{p}_{i}\rangle$
as a function of the equivalent particle size $D_{v}/\eta$. It is found that
the ratio
$\langle{\Omega}_{i}^{s}{\Omega}_{i}^{s}\rangle/\langle\dot{p}_{i}\dot{p}_{i}\rangle\sim(D_{v}/\eta)^{-0.6}$
for all considered shapes. We remark that the existence of the observed
scaling behavior of the ratio
$\langle{\Omega}_{i}^{s}{\Omega}_{i}^{s}\rangle/\langle\dot{p}_{i}\dot{p}_{i}\rangle$
is not certain. However, this emphasizes that the scaling exponents of the
mean square tumbling rate and spinning rate can not be both equal to -4/3. To
our knowledge, no presently available model accounts for the observed
behaviour.

In summary, the main feature connected to the finite size is that larger
particles tend to tumble much more than spinning (see cartoon Figure 7c).
Further, we observe that the VA and VP measurements do not have this property.
While VA nearly keeps the value of the Jeffery randomly oriented prediction,
VPs do not sensibly vary the relative weight of spinning with respect to
tumbling and apparently approximately follow the simplistic K41 prediction, eq
(15), and (red line) 7d.

Figure 8:  Mean squared angular velocity, normalized by the results of tracers
with corresponding aspect ratio, as a function of $D_{v}/\eta$ for the
numerical experiments of (a) RPFL and (b) VPFL, compared with the results of
RP. Insets show the ratio of the mean squared angular velocity of the
numerical experiments to the results of RP as a function of $D_{v}/\eta$. The
ratio of the mean spinning rate squared to the mean tumbling rate squared,
$\langle\Omega_{i}^{s}\Omega_{i}^{s}\rangle/\langle\dot{p}_{i}\dot{p}_{i}\rangle$,
as a function of aspect ratio $\alpha$ for the numerical experiments of (c)
RPFL and (d) VPFL, respectively. The color bar represents the normalized
particle size $D_{v}/\eta$ in the logarithmic scale.

In order to find out the leading physical mechanism responsible for the
decrease of
$\langle{\Omega}_{i}^{s}{\Omega}_{i}^{s}\rangle/\langle\dot{p}_{i}\dot{p}_{i}\rangle$,
we carry out the numerical experiments with particles whose location is fixed,
either with the solid-fluid interface or without it, denoted respectively as
RPFL and VPFL, as shown in Figure 8. The rationale behind fixing the position
is to rule out the role of preferential sampling of the flow position (Squires
& Eaton, 1991; Maxey, 1987; Calzavarini et al., 2008). Indeed, even if this
has never been observed, for finite-sized spherical particles Fiabane et al.
(2012), we can in principle not exclude it for non-spherical large particles.
On the other hand, a direct measurement of preferential concentration in the
present simulations is not an easy task due to the limited number of particles
that we could place in the domain.

We find that the mean normalized rotation rate variance of RPFL scales
similarly RP and have a similar amplitude (see Figure 8a). Furthermore,
$\langle{\Omega}_{i}^{s}{\Omega}_{i}^{s}\rangle/\langle\dot{p}_{i}\dot{p}_{i}\rangle$
decreases as the particle size increases (see Figure 8c), which is consistent
with RP. These results indicate that the possible preferential sampling is not
responsible for the decrease of
$\langle{\Omega}_{i}^{s}{\Omega}_{i}^{s}\rangle/\langle\dot{p}_{i}\dot{p}_{i}\rangle$.
To further confirm this point, we show the rotational dynamics of VPFL in
Figure 8(b and d). Similar scaling behavior of mean rotation rate square as a
function of particle size is observed. And
$\langle{\Omega}_{i}^{s}{\Omega}_{i}^{s}\rangle/\langle\dot{p}_{i}\dot{p}_{i}\rangle$
remains constant at the value of the fixed tracers for the corresponding
aspect ratio. Since VP do not have a solid interface and do not affect the
flow, the results of these two additional numerical experiments suggest that
the two-way coupling between the particles and the flow is responsible for the
decrease of the ratio
$\langle{\Omega}_{i}^{s}{\Omega}_{i}^{s}\rangle/\langle\dot{p}_{i}\dot{p}_{i}\rangle$.

####  3.3.1 The role of particle boundary layers

At this point, given the fact that fully resolved simulations allow virtually
to measure any fluid quantity, one may wonder if it is possible to better
characterize the effect of the particle feedback on the flow. With this in
mind, we measure the fluid vorticity $\bm{\omega}$ and the rate of strain
tensor $\mathcal{S}$ in the boundary layers surrounding the particles. We then
compare the relative importance of their intensities. This is performed by
computing the variance ratio
$\tfrac{1}{4}\overline{\bm{\omega}^{2}}/\overline{\mathcal{S}:\mathcal{S}}$,
where the average is performed over time, ensemble i.e. multiple trajectories,
and over the azimuthal direction in order to improve the statistical
convergence. We note that the ratio must have the value 1 in the particle-free
HIT flow.

Figure 9(a-c) shows three typical fields of the considered ratio around the
particle up to the size of $1.5D_{v}$ It is found that
$\tfrac{1}{4}\overline{\bm{\omega}^{2}}/\overline{\mathcal{S}:\mathcal{S}}\simeq
1.1$ (vorticity dominated) in a thin shell around the particle, while at
immediately larger distances a wider region displays
$\tfrac{1}{4}\overline{\bm{\omega}^{2}}/\
\overline{\mathcal{S}:\mathcal{S}}<1$ (strain dominated), finally beyond $\sim
1.5D_{v}$ any effect has vanished.

The quantity
$\langle\tfrac{1}{4}\overline{\bm{\omega}^{2}}\rangle|_{1.5D_{v}}/\langle\overline{\mathcal{S}:\mathcal{S}}\rangle|_{1.5D_{v}}$,
where $\langle\ldots\rangle_{1.5D_{v}}$ is a volume average over shells up to
$1.5D_{v}$, for different aspect ratios and particle sizes is reported in
Figure 9(d). We can observe that the overall effect in the region affected by
the particle presence is to decrease the importance of vorticity with respect
to strain. This points to the fact that spinning is possibly reduced as a
consequence of the reduction of vorticity in the flow surrounding a particle,
in qualitative agreement with previous observations. However, we think that
more sensitive tests remain to be performed in order to soundly quantify the
effect of feedback on the flow.

Figure 9: Field of the ratio
$\tfrac{1}{4}\overline{\bm{\omega}^{2}}/\overline{\mathcal{S}:\mathcal{S}}$ at
(a) $\alpha=1$, $D_{v}=20$ (b) $\alpha=0.5$, $D_{v}=25.1$, and (c) $\alpha=2$,
$D_{v}=25.1$ at $Re_{\lambda}=120$. The colour bar indicates the value of
$\tfrac{1}{4}\overline{\bm{\omega}^{2}}/\overline{\mathcal{S}:\mathcal{S}}$.
(d) The ratio of
$\langle\tfrac{1}{4}\overline{\bm{\omega}^{2}}\rangle|_{1.5D_{v}}$ to
$\langle\overline{\mathcal{S}:\mathcal{S}}\rangle|_{1.5D_{v}}$ as a function
of $D_{v}/\eta$.

##  4 Conclusions

In summary, some statistical properties of the translational and rotational
dynamics of neutrally buoyant finite-size spheroids in HIT have been
investigated by means of interface-resolved direct numerical simulations.
Besides, the simulation of the real system, several extra numerical
experiments have been carried out in order to enhance our insight on the
essential features of the coupling of the particles with the flow. It is found
that non-spherical inertial-scale particles experience spatially volume
filtered turbulent fluctuations, similarly to spherical particles, even if the
presence of finite-size particle sensibly influences the flow around the
particle over a region that has similar aspect-ratio as the particle and
approximately doubled $D_{v}$. The variance of the acceleration of spheroids
with different shapes collapses on a single line when plotted as a function of
the equivalent diameter $D_{v}$, and at $Re_{\lambda}=120$ a power-law scaling
$~{}D_{v}^{-1}$ is observed. It turns out that the acceleration variance
versus particle size of virtual particles and fluid volume averages shows a
similar scaling behavior, which proves that the scaling behavior of the
acceleration intensity is dominated by the volume averaging of the
acceleration of the flow field. Furthermore, the orientation of particles with
size in the inertial range shows preferential alignments with flow structures.
This is a novel result in the context of large particles, because it has been
previously observed only for particles in the dissipative scale range of
turbulence. We find that the preferential alignment of the particle
acceleration shows a negligible size dependence, which suggests a similarity
of flow coherent structures from the dissipative to inertial scales.

Next, we studied the rotational dynamics of the particles as their size is
increased. The mean square angular velocity of spheroids with different shapes
show a clear scaling behavior with an exponent $~{}D_{v}^{-4/3}$, which agrees
with predictions based on the K41 theory and with the available experiments
(Parsa & Voth, 2014; Bordoloi & Variano, 2017; Oehmke et al., 2021). However,
it turns out that the ratio of the mean square spinning rate to the mean
square tumbling rate shows a systematic decrease as the particle size
increases, and this is for all considered aspect ratios. To our knowledge,
this new feature has never been observed before, and it calls for future
experimental verifications. On the basis of the analysis of virtual particles,
we interpret this phenomenon as due to the particle feedback on the flow. That
is to say to the formation of boundary layers around the particles. In such
layers the fluid strain-rate is overall increased in comparison to the
vorticity, and this favors the tumbling with respect to the spinning.

The present findings improve the understanding of translational and rotational
dynamics for neutrally-buoyant finite-size spheroids in turbulence. The
results also shed light on a connection between the flow structures at the
dissipative scales and inertial scales, which might open a new perspective for
the investigation of turbulent flow. We expect that a similar approach as the
one followed in this study will also reaveal useful understanding of the
dynamics of particles heavier/lighter than the carrying turbulent fluid flow.

## Acknowledgements

This work was supported by Natural Science Foundation of China under grant
nos. 11988102 and 91852202, and Tencent Foundation through the xplorer prize.

##  Appendix A Particle acceleration statistics at $Re_{\lambda}=32$

In order to validate the numerical results, we carried out simulations at
$Re_{\lambda}=32$ for spheres, in the same condition as in a previous study
(Homann & Bec, 2010). Figure 10a shows the single-cartesian-component
acceleration variance, $\langle a_{i}^{2}\rangle$, normalized by the
acceleration variance of tracers, $\langle a_{i}^{2}\rangle^{tracer}$, as a
function of normalized particle diameter $d/\eta$ ($d$ is the diameter of
spheres). Our simulation results show good agreement with the previous
study(Homann & Bec, 2010) and scale as $~{}(d/\eta)^{-4/3}$. At low
$Re_{\lambda}$, the scaling exponent is smaller than -2/3 which is the
exponent value expected in HIT with high $Re_{\lambda}$. The dashed-dot line
represents the prediction of the Faxén correction based on point-particle
simulation under the same $Re_{\lambda}$. When the particle is small
($d/\eta<4$), the acceleration variance is captured by the Faxén corrections.
However, the Faxén correction can not account for the scaling behavior
acceleration variance when $d/\eta>7$. The inset shows the PDFs of the
acceleration component for the variant size of spheres. The effect of particle
size on the PDF of acceleration is not significant, which is consistent with
previous studies(Xu et al., 2007; Brown et al., 2009; Homann & Bec, 2010).
Figure 10b shows the normalized integral time of the temporal correlation
function for particle acceleration, $T_{I}/\tau_{\eta}$, as a function
$d/\eta$. The integral time $T_{I}$ is defined as the integral of the
autocorrelation function from time zero to its first zero-crossing time
$T_{0}$

$T_{I}\equiv\int_{0}^{T_{0}}C(\tau)d\tau,~{}~{}C(\tau)\equiv\frac{\langle
a_{i}(t+\tau)a_{i}(t)\rangle}{\langle(a_{i}(t))^{2}\rangle}$ (16)

Our results are consistent with the correlation time in Ref. (Homann & Bec,
2010) and show the power-law behaviour $T_{I}\sim(d)^{2/3}$. The inset shows
the autocorrelation functions for various sizes of spheres, illustrating the
increasing autocorrelation time as the size of the sphere increases. In
conclusion, our translational statistics of simulation show good agreement
with the previous study(Homann & Bec, 2010), which validates our simulations.

Figure 10: a). Normalized single cartesian component of the acceleration
variance $\langle a_{i}^{2}\rangle/\langle a_{i}^{2}\rangle^{tracer}$ as a
function of the normalized particle size $d/\eta$ at $Re_{\lambda}=32$.
Squares denote the results in Ref. (Homann & Bec, 2010) at the same
$Re_{\lambda}$. The dashed-dot line denotes the Faxén correction (Calzavarini
et al., 2012). The inset shows the PDFs of acceleration at different sizes.
b). Normalized correlation time of acceleration, $T_{I}/\tau_{\eta}$, as a
function of the normalized particle size $D_{v}/\eta$. Squares denote the
correlation time in Ref. (Homann & Bec, 2010) at the same $Re_{\lambda}$. The
inset shows the autocorrelation function of a single cartesian acceleration
component.

##  Appendix B Comparison with Oehmke et al.(2021) experiments

We provide here a side-by-side comparison of the present simulations with the
experiments by Oehmke et al. (2021), which are the only to date to have
measured both spinning and tumbling rates of particles. The experimental
conditions of developed turbulence at $Re_{\lambda}\in\left[90,630\right]$ and
particle aspect ratios $\alpha=5.4,7.5,10.8$ are relatively close to the ones
explored in the present study. Figure 11(a) shows the ratio of spinning rate
variance over tumbling rate variance as a function of the aspect ratio. There
is agreement on the dominance of spinning compared to tumbling for prolate
($\alpha>1$) particles (note that the isotropic value for spheres is 0.5).
However, the decreasing trend at increasing the particle size $D_{v}$ which is
clearly observed in the simulations is only partially consistent with the
experiments. Figure 11(b) represents the same measurements in a compensated
form: the ratio of rotation rate variances is normalized by the $\alpha^{4/3}$
scaling as in the K41-based prediction Eq. (15) and plotted against the length
of the particle symmetry axis $l$ in $\lambda$ units. This choice of
representation is adapted from figure 5(b) in Oehmke et al. (2021). We observe
here an excellent agreement for the behaviour of all prolate particles.
However, the decreasing behaviour observed both for prolate and oblate
particles stresses further the failure of the (15) prediction.

Figure 11: Ratio of spinning rate variance over tumbling rate variance:
comparison of present DNS results at $Re_{\lambda}=120$ with the experimental
measurements by Oehmke et al. (2021) at $Re_{\lambda}\in\left[90,630\right]$.
(a) Particle rotation ratio as a function of aspect ratio, here the color code
maps the particle size in equivalent radius units;(b) Rotations ratio
compensated by $\alpha^{4/3}$ as a function of the particle length (i.e. the
symmetry axis) $l$ normalized by the Taylor microscale $\lambda$ (same as
figure 5(b) in Oehmke et al. PRF 2021).

## References

  * Ardekani et al. (2016) Ardekani, M. N., Costa, P., Breugem, W. P. & Brandt, L. 2016 Numerical study of the sedimentation of spheroidal particles. Int. J. Multiph. Flow 87, 16–34.
  * Bec et al. (2007) Bec, J., Biferale, L., Cencini, M., Lanotte, A., Musacchio, S. & Toschi, F. 2007 Heavy particle concentration in turbulence at dissipative and inertial scales. Phys. Rev. Lett. 98 (8), 084502.
  * Bentkamp et al. (2019) Bentkamp, L., Lalescu, C. C. & Wilczek, M. 2019 Persistent accelerations disentangle lagrangian turbulence. Nat. Commun. 10 (1), 1–8.
  * Benzi et al. (2009) Benzi, R., Biferale, L., Calzavarini, E., Lohse, D. & Toschi, F. 2009 Velocity-gradient statistics along particle trajectories in turbulent flows: The refined similarity hypothesis in the lagrangian frame. Phys. Rev. E 80, 066318.
  * Biferale et al. (2005) Biferale, L., Boffetta, G., Celani, A., Lanotte, A. & Toschi, F. 2005 Particle trapping in three-dimensional fully developed turbulence. Phys. Fluids 17 (2), 021701.
  * Bordoloi & Variano (2017) Bordoloi, A. D. & Variano, E. 2017 Rotational kinematics of large cylindrical particles in turbulence. J. Fluid Mech. 815, 199–222.
  * Bounoua et al. (2018) Bounoua, S., Bouchet, G. & Verhille, G. 2018 Tumbling of inertial fibers in turbulence. Phys. Rev. Lett. 121, 124502.
  * Brenner (1961) Brenner, H. 1961 The slow motion of a sphere through a viscous fluid towards a plane surface. Chem. Eng. Sci. 16 (3), 242–251.
  * Brown et al. (2009) Brown, R. D., Warhaft, Z. & Voth, G. A. 2009 Acceleration statistics of neutrally buoyant spherical particles in intense turbulence. Phys. Rev. Lett. 103, 194501.
  * Brändle de Motta et al. (2016) Brändle de Motta, J.C., Estivalezes, J.L., Climent, E. & Vincent, S. 2016 Local dissipation properties and collision dynamics in a sustained homogeneous turbulent suspension composed of finite size particles. International Journal of Multiphase Flow 85, 369–379.
  * Byron et al. (2015) Byron, M., Einarsson, J., Gustavsson, K., Voth, G. A., Mehlig, B. & Variano, E. 2015 Shape-dependence of particle rotation in isotropic turbulence. Phys. Fluids 27 (3), 035101.
  * Calzavarini (2019) Calzavarini, E. 2019 Eulerian-lagrangian fluid dynamics platform: The ch4-project. Software Impacts 1, 100002.
  * Calzavarini et al. (2020) Calzavarini, E., Jiang, L. & Sun, C. 2020 Anisotropic particles in two-dimensional convective turbulence. Phys. Fluids 32 (2), 023305.
  * Calzavarini et al. (2008) Calzavarini, E., Kerscher, M., Lohse, D. & Toschi, F. 2008 Dimensionality and morphology of particle and bubble clusters in turbulent flow. J. Fluid Mech. 607, 13–24.
  * Calzavarini et al. (2009) Calzavarini, E., Volk, R., Bourgoin, M., Leveque, E., Pinton, J. F. & Toschi, F. 2009 Acceleration statistics of finite-sized particles in turbulent flow: the role of faxén forces. J. Fluid Mech. 630, 179–189.
  * Calzavarini et al. (2012) Calzavarini, E., Volk, R., Leveque, E., Pinton, J. F. & Toschi, F. 2012 Impact of trailing wake drag on the statistical properties and dynamics of finite-sized particle in turbulence. Phys. D 241D (3), 237–244.
  * Chong et al. (1990) Chong, M. S., Perry, A. E. & Cantwell, B. J. 1990 A general classification of three‐dimensional flow fields. Phys. Fluids 2 (5), 765–777.
  * Cisse et al. (2013) Cisse, M., Homann, H. & Bec, J. 2013 Slipping motion of large neutrally buoyant particles in turbulence. J. Fluid Mech. 735, R1.
  * Cooley & O’Neill (1969) Cooley, M. D. A. & O’Neill, M. E. 1969 On the slow motion generated in a viscous fluid by the approach of a sphere to a plane wall or stationary sphere. Mathematika 16 (1), 37–49.
  * Costa et al. (2015) Costa, P., Boersma, B. J., Westerweel, J. & Breugem, W. P. 2015 Collision model for fully resolved simulations of flows laden with finite-size particles. Phys. Rev. E 92, 053012\.
  * Do-Quang et al. (2014) Do-Quang, M., Amberg, G., Brethouwer, G. & Johansson, A. V. 2014 Simulation of finite-size fibers in turbulent channel flows. Phys. Rev. E 89, 013006.
  * Dolata & Zia (2021) Dolata, B. E. & Zia, R. N. 2021 Faxén formulas for particles of arbitrary shape and material composition. Journal of Fluid Mechanics 910, A22.
  * Fiabane et al. (2012) Fiabane, L., Zimmermann, R., Volk, R., Pinton, J. F. & Bourgoin, M. 2012 Clustering of finite-size particles in turbulence. Phys. Rev. E 86, 035301.
  * Gustavsson et al. (2014) Gustavsson, K., Einarsson, J. & Mehlig, B. 2014 Tumbling of small axisymmetric particles in random and turbulent flows. Phys. Rev. Lett. 112, 014501.
  * Homann & Bec (2010) Homann, H. & Bec, J. 2010 Finite-size effects in the dynamics of neutrally buoyant particles in turbulent flow. J. Fluid Mech. 651, 81–91.
  * Jeffery (1922) Jeffery, G. B. 1922 The motion of ellipsoidal particles immersed in a viscous fluid. Proc. R. Soc. Lond. A 102 (715), 161–179.
  * Jiang et al. (2020) Jiang, L., Calzavarini, E. & Sun, C. 2020 Rotation of anisotropic particles in rayleigh–bénard turbulence. J. Fluid Mech. 901, A8.
  * Jiang et al. (2021) Jiang, L., Wang, C., Liu, S.and Sun, C. & Calzavarini, E. 2021 Rotational dynamics of bottom-heavy rods in turbulence from experiments and numerical simulations. Theor. App. Mech. Lett. 11, 100227.
  * Kuperman et al. (2019) Kuperman, S., Sabban, L. & van Hout, R. 2019 Inertial effects on the dynamics of rigid heavy fibers in isotropic turbulence. Phys. Rev. Fluids 4, 064301.
  * Liberzon et al. (2012) Liberzon, A., Luthi, B., Holzner, M., Ott, S., Berg, J. & Jakob, M. 2012 On the structure of acceleration in turbulence. Phys. D 241 (3), 208–215.
  * Luo et al. (2007) Luo, K., Wang, Z., Fan, J. & Cen, K. 2007 Full-scale solutions to particle-laden flows: Multidirect forcing and immersed boundary method. Phys. Rev. E 76, 066709.
  * Mathai et al. (2016) Mathai, V., Calzavarini, E., Brons, J., Sun, C. & Lohse, D. 2016 Microbubbles and microparticles are not faithful tracers of turbulent acceleration. Phys. Rev. Lett. 117 (2), 024501.
  * Mathai et al. (2018) Mathai, V., Huisman, S. G., Sun, C., Lohse, D. & Bourgoin, M. 2018 Dispersion of air bubbles in isotropic turbulence. Phys. Rev. Lett. 121 (5), 054501.
  * Maxey (1987) Maxey, M. R. 1987 The gravitational settling of aerosol particles in homogeneous turbulence and random flow fields. J. Fluid Mech. 174, 441–465.
  * Mittal & Iaccarino (2005) Mittal, R. & Iaccarino, G. 2005 Immersed boundary methods. Annu. Rev. Fluid Mech. 37, 239–261.
  * Ni et al. (2014) Ni, Rui, Ouellette, Nicholas T. & Voth, Greg A. 2014 Alignment of vorticity and rods with lagrangian fluid stretching in turbulence. J. Fluid Mech. 743, R3.
  * Oehmke et al. (2021) Oehmke, T. B., Bordoloi, A. D., Variano, E. & Verhille, G. 2021 Spinning and tumbling of long fibers in isotropic turbulence. Phys. Rev. Fluids 6, 044610.
  * Parsa et al. (2012) Parsa, S., Calzavarini, E., Toschi, F. & Voth, G. A. 2012 Rotation rate of rods in turbulent fluid flow. Phys. Rev. Lett. 109 (13), 134501.
  * Parsa & Voth (2014) Parsa, S. & Voth, G. A. 2014 Inertial range scaling in rotations of long rods in turbulence. Phys. Rev. Lett. 112, 024501\.
  * Peskin (2002) Peskin, C. S. 2002 The immersed boundary method. Acta Numer. 11, 479–517.
  * Porta et al. (2001) Porta, A. La, Voth, G. A., Crawford, A. M., Alexander, J. & Bodenschatz, E. 2001 Fluid particle accelerations in fully developed turbulence. Nature 409, 22.
  * Pujara et al. (2018) Pujara, N., Oehmke, T. B., Bordoloi, A. D. & Variano, E. A. 2018 Rotations of large inertial cubes, cuboids, cones, and cylinders in turbulence. Phys. Rev. Fluids 3, 054605.
  * Pujara et al. (2019) Pujara, Nimish, Voth, Greg A. & Variano, Evan A. 2019 Scale-dependent alignment, tumbling and stretching of slender rods in isotropic turbulence. J. Fluid Mech. 860, 465–486.
  * Pumir & Wilkinson (2011) Pumir, A. & Wilkinson, M. 2011 Orientation statistics of small particles in turbulence. New J. Phys. 13 (9), 093030.
  * Qureshi et al. (2007) Qureshi, N. M., Bourgoin, M., Baudet, C., Cartellier, A. & Gagne, Y. 2007 Turbulent transport of material particles: an experimental study of finite size effects. Phys. Rev. Lett. 99 (18), 184502.
  * Shen et al. (2021) Shen, J., Lu, Z., Wang, L. & Peng, C. 2021 Influence of particle-fluid density ratio on the dynamics of finite-size particles in homogeneous isotropic turbulent flows. Phys. Rev. E 104, 025109.
  * Squires & Eaton (1991) Squires, K. D. & Eaton, J. K. 1991 Preferential concentration of particles by turbulence. Phys. Fluids 3, 1169\.
  * Suzuki & Inamuro (2011) Suzuki, K. & Inamuro, T. 2011 Effect of internal mass in the simulation of a moving body by the immersed boundary method. Comput. Fluids 49 (1), 173–187.
  * Uhlmann (2005) Uhlmann, M. 2005 An immersed boundary method with direct forcing for the simulation of particulate flows. J. Comput. Phys. 209 (2), 448–476.
  * Volk et al. (2011) Volk, R., Calzavarini, E., Leveque, E. & Pinton, J. F. 2011 Dynamics of inertial particles in a turbulent von Kármán flow. J. Fluid Mech. 668, 223–235.
  * Voth et al. (2002) Voth, G. A., Porta, L. A., Crawford, A. M., Alexander, J. & Bodenschatz, E. 2002 Measurement of particle accelerations in fully developed turbulence. J. Fluid Mech. 469, 121–160.
  * Voth & Soldati (2017) Voth, G. A. & Soldati, A. 2017 Anisotropic particles in turbulence. Annu. Rev. Fluid Mech. 49, 249–276.
  * Wang et al. (2021) Wang, C., Jiang, L., Jiang, H., Sun, C. & Liu, S. 2021 Heat transfer and flow structure of two-dimensional thermal convection overratchet surfaces. J. Hydrodyn. p. In Press.
  * Xu & Bodenschatz (2008) Xu, H. & Bodenschatz, E. 2008 Motion of inertial particles with size larger than kolmogorov scale in turbulent flows. Phys. D 237 (14-17), 2095–2100.
  * Xu et al. (2007) Xu, H., Ouellette, N. T., Vincenzi, D. & Bodenschatz, E. 2007 Acceleration correlations and pressure structure functions in high-reynolds number turbulence. Phys. Rev. Lett. 99, 204501.