# Particle chirality does not matter in strong turbulence

G<EMAIL_ADDRESS>M.P.A. Assen1 D. Lohse1,2 &
R<EMAIL_ADDRESS>1Physics of Fluids Group
and Max Planck Center for Complex Fluid Dynamics, Department of Science and
Technology, J.M. Burgers Center for Fluid Dynamics, and MESA+ Institute,
University of Twente, P.O. Box 217, 7500AE Enschede, The Netherlands;
2Max Planck Institute for Dynamics and Self-Organization, Am Fassberg 17,
37077 Göttingen, Germany;
3Dipartimento di Ingegneria Industriale, University of Rome “Tor Vergata”, Via
del Politecnico 1, 00133 Roma, Italy;
4Gran Sasso Science Institute, Viale F. Crispi 7, 67100 L’Aquila, Italy.

###### Abstract

We use three–dimensional direct numerical simulations of homogeneous isotropic
turbulence in a cubic domain to investigate the dynamics of heavy, chiral,
finite–size inertial particles and their effects on the flow. Using an
immersed–boundary method and a complex collision model, four–way coupled
simulations have been performed and the effects of particle–to–fluid density
ratio, turbulence strength, and particle volume fraction have been analysed.
We find that freely falling particles on the one hand add energy to the
turbulent flow but, on the other hand, they also enhance the flow dissipation:
depending on the combination of flow parameters, the former or the latter
mechanism prevails, thus yielding enhanced or weakened turbulence.
Furthermore, particle chirality entails a preferential angular velocity which
induces a net vorticity in the fluid phase. As turbulence strengthens, the
energy introduced by the falling particles becomes less relevant and stronger
velocity fluctuations alter the solid phase dynamics, making the effect of
chirality irrelevant. Moreover, comparing the time–history of collision events
for chiral particles and spheres (at the same volume fraction) suggests that
the former tend to entangle, in contrast to the latter which rebound
impulsively.

##  1 Introduction

Turbulent, inertial particle–laden flows are ubiquitous in natural and
engineering environments (Balmforth & Provenzale, 2001). Aerosols in clouds
(Shaw, 2003), particle–driven gravity currents (Necker et al., 2005),
sediments in oceans or rivers (Hargrave, 1985), sand storms in the atmosphere
(Voilland, 2018) as well as fluidisation of pipe or open–channel flows in
industrial processes (Lashgari et al., 2014) are just a few examples among
many. Owing to the high practical interest, investigation of these problems
has increased in the last decade, especially thanks to the massive
improvements achieved in terms of measurement and simulation performance.
There are many aspects that make this topic an extraordinary challenge: the
different properties of both particles and the carrier fluid, the multi-way
coupling between continuous and dispersed phases and the relevant difference
in spatial and time scales (Brandt & Coletti, 2022; Toschi & Bodenschatz,
2009). The interaction between the two phases can affect the flow both at
large and small scales depending on factors such as particle shape, size, and
volume/mass fraction. Depending on their relative size to the viscous scale,
heavy particles can attenuate or enhance turbulence (Yu et al., 2021), in
contrast to small, settling particles which always lead to turbulence damping
(Zhao et al., 2010). Moreover, particles of greater mass density than the
surrounding fluid can detach from the flow, leading to preferential sampling,
small-scale fractal clustering, and significant relative velocities (Bec et
al., 2024).

In recent years, studies have predominantly focused on particles with
vanishing inertia (behaving as passive tracers) or highly symmetric shapes
such as spheres or ellipsoids (Shaw, 2003; Voilland, 2018; Sozza et al.,
2022). However, investigations into anisotropic particles have gained traction
only recently. Elastic fibers, for instance, induce similar turbulence damping
as spheres but can also modify the energy scale-by-scale distribution
(Olivieri et al., 2022; Dotto & Marchioli, 2019). Inertial rods also possess
the potential to influence turbulence-regeneration mechanisms, with rigid,
inertia-less fibers generating additional stresses that dampen counter-
rotating vortices in the wall region, leading to drag reduction. Inertial
fibers may preferentially segregate into low-speed streaks, thereby delaying
fluid flow (Voth & Soldati, 2017).

Despite the abundance of literature in this field, most studies have focused
on isotropic or highly symmetric particle shapes (Brandt & Coletti, 2022),
overlooking a significant portion of the parameter space concerning strongly
asymmetric large particles in turbulence (Voth & Soldati, 2017). The present
paper is a first attempt aimed at filling this gap. Among several
possibilities we have considered chiral particles which are distinct from
their mirror image (Roos & Roos, 2015). We have chosen chiral particles
because they break spatial reflection symmetry and couple translational and
rotational degrees of freedom. A particle forced to translate in a fluid will
also rotate and vice versa; this feature is also present in nature and can be
appreciated by watching the falling of maple seeds. Questions we ask are: how
do chiral particles behave in a turbulent flow? Are they able to modify some
features of turbulence? What is the influence of parameters like
solid–to–fluid density ratio, turbulence strength and solid phase volume
fraction?

The paper is organised as follows: in $\S$ 2 we outline the numerical model
and the set-up of the simulations. The findings are discussed in $\S$ 3 and
the complex interplay between particle and flow dynamics is illustrated by
taking into account the effect of multiple parameters separately. Closing
remarks and an outlook are given in $\S$ 4. The paper is completed with
additional material in Appendix A, B, and C.

##  2 Problem and methods

###  2.1 Setup

We consider a three–periodic cubic domain filled with a Newtonian fluid of
density $\rho_{f}$ laden by $N_{p}$ identical chiral particles of density
$\rho_{p}$. Each particle is made of three perpendicular cylinders, arranged
as in Figure 1a, with the ends rounded by spherical caps. Let $V_{p}$ be the
volume of the particle and $D_{eq}$ the diameter of the equivalent sphere
(volumwise) which is assumed as unit length throughout the paper. The particle
can be then enclosed in a cubic bounding-box of edges
$\ell_{x}$=$\ell_{y}$=$\ell_{z}$=$1.5D_{eq}$ while the diameter of the
cylindrical legs is $d\simeq 0.42D_{eq}$. The particle is convex, as its
centroid falls outside its volume, with coordinates $G$=$(0.75,0.42,1.08)$ in
$D_{eq}$ units measured from the origin $O$. The asymmetric shape implies
principal axes not aligned with any particle legs yielding for their unit
vectors $\widehat{\bf x}_{p}$=$(-0.774,0.447,0.447)$, $\widehat{\bf
y}_{p}$=$(0.000,-0.707,0.707)$ and $\widehat{\bf z}_{p}$=$(0.633,0.548,0.548)$
and for the principal moments of inertia $I_{xp}=0.199$, $I_{yp}=0.049$ and
$I_{zp}=0.027$ in $\rho_{p}D_{eq}^{5}$ units. It is worthwhile pointing out
that the present particles are equivalent to a sphere of diameter $D_{eq}$
only volumewise; in fact, the wet surface of the former is about $50\%$ bigger
than the equivalent sphere while the encumbrance of the boundingbox $\approx
3.4$ times that of a cube enclosing the sphere. These differences are relevant
for the particle/particle and particle/fluid interaction as will be discussed
in the section of Results.

The flow domain is cubic with edges of length
$L_{x}$=$L_{y}$=$L_{z}$=$10D_{eq}$ in each direction and a constant
gravitational acceleration $\mathbb{g}$ anti–parallel to the $z$ direction.
Particles are initially distributed randomly, both in position and
orientation, within the computational domain, the only constraint being that
they cannot intersect spatially: a schematic of the system is shown in Figure
1b.

After an initial transient, the system attains a statistical steady state over
which data are collected for the analysis; depending on the problem
parameters, statistic convergence is achieved for different time windows. We
define the large–eddy–turnover time $T_{e}=u_{rms}^{2}/\varepsilon$ with the
velocity fluctuations $u_{rms}$ and the kinetic energy dissipation rate
$\varepsilon$ of the turbulent flow and the ‘observation time’ $T_{o}$ the
integration time window (in $D_{eq}/U$ units, $U$ being the flow velocity
scale): simulations of HIT without particles have been run for $33\leq
T_{o}/T_{e}\leq 10^{3}$, while for the simulations with a single particle we
had $510\leq T_{o}/T_{e}\leq 860$. Finally the cases with $N_{p}=20$ particles
have been integrated for $75\leq T_{o}/T_{e}\leq 840$. For each simulation the
duration of the integration has been decided by monitoring the quantities of
interest and verifying that their values did not change appreciably during the
final part of the simulation.

Before concluding this section we wish to point out that in some studies
dealing with settling particles in turbulence (see for example Chouippe &
Uhlmann (2015)) an elongated domain in the direction of gravity has been
employed in order to avoid possible correlation effects produced by the
preferential velocity direction. In this study we have preliminarily verified
that, at least for the flow quantities considered in the present analysis, the
differences in the results obtained in a cubic and in an elongated domain are
negligible; selected results from these tests are given in Appendix C.

Figure 1: (a) Representation of a chiral particle with global (red), body
centroid (black) and principal inertia (green) reference frames. The green
bullet (G) outside the particle is its center of mass. (b) Perspective view of
the system with randomly distributed particles in homogeneous isotropic
turbulence for $Re_{\lambda}=75$, volume fraction $\phi=1\%$ and density ratio
$\rho_{p}/\rho_{f}=2$. Colours of the background are nondimensional vorticity
magnitude contours ranging from $0$ (blue) to $20$ (yellow).

###  2.2 Carrier phase

The governing relations for the fluid phase are the incompressible Navier-
Stokes equations and the continuity equation which in nondimensional form read

$\frac{\partial\mathbb{u}}{\partial
t}+(\mathbb{u}\cdot\nabla)\mathbb{u}=-\nabla
p+\frac{1}{Re}\nabla^{2}\mathbb{u}+\boldsymbol{{\rm f}},\ \ \ \ \
\nabla\cdot\mathbb{u}=0.$ (1)

Here $\mathbb{u}$ and $p$ are fluid velocity vector and kinematic pressure,
respectively, while $Re=D_{eq}U/\nu$ is the Reynolds number based on the
length $D_{eq}$ and the flow velocity $U$ scales. The last term in the first
of equations (1) is a vector containing two volume forces $\boldsymbol{{\rm
f}}=\boldsymbol{f}_{{}_{T}}+\boldsymbol{f}_{{}_{P}}$:
$\boldsymbol{f}_{{}_{T}}$ sustains the HIT flow while
$\boldsymbol{f}_{{}_{P}}$ is the IBM forcing which accounts for the presence
of solid particles and allows for the two–way coupling between solid and fluid
phases.

We numerically solve the governing equations (1) with our in–house advanced
finite difference code (“AFiD”) which is extensively described and validated
in van der Poel et al. (2015); Spandan et al. (2020). The main features are:
spatial derivatives are approximated by conservative, second–order accurate
finite–differences discretised on a staggered mesh which is uniform and
homogeneous ($\Delta x\\!=\\!\Delta y\\!=\\!\Delta z\\!=\\!\Delta$). A
combination of Crank–Nicolson and low–storage third–order Runge–Kutta schemes
is used to integrate, respectively, the viscous terms implicitly and all other
terms explicitly. Finally pressure and momentum are strongly coupled through a
fractional–step method as described in Rai & Moin (1991); Verzicco & Orlandi
(1996).

###  2.3 Homogeneous Isotropic Turbulence forcing method

The HIT forcing $\boldsymbol{f}_{{}_{T}}$ is prescribed, similarly to Chouippe
& Uhlmann (2015), adopting the method introduced by Eswaran & Pope (1988). It
allows to attain a statistically stationary velocity field by forcing large
scale (small wavenumber) components. This is achieved through a random
Uhlenbeck–Ornstein process via a vector
$\boldsymbol{\hat{b}}(\boldsymbol{k},t)$ whose components are prescribed in
time according to

$\hat{b}_{i}(\boldsymbol{k},t+\Delta
t)=\hat{b}_{i}(\boldsymbol{k},t)\Bigg{(}1-\frac{\Delta
t}{T_{L}}\Bigg{)}+e_{i}(\boldsymbol{k},t)\Bigg{(}2\sigma^{2}\frac{\Delta
t}{T_{L}}\Bigg{)}^{\frac{1}{2}},\hskip 14.22636pti=1,2,3.$ (2)

Here $\boldsymbol{k}=(l,m,n)$ is the wavenumber vector, $\Delta t$ the
numerical integration time step and $e_{i}(\boldsymbol{k},t)$ a complex random
number drawn from a standardised Gaussian distribution, $T_{L}$ and
$\sigma^{2}$ are time scale and variance of the random process.

The forcing term $\boldsymbol{\hat{f}}_{{}_{T}}(\boldsymbol{k},t)$, in Fourier
space, is non–zero only in the band $|\boldsymbol{k}|\leq k_{f}=2.3$ and the
final forcing is obtained by an orthogonal projection of
$\boldsymbol{\hat{b}}(\boldsymbol{k},t)$

$\boldsymbol{\hat{f}}_{{}_{T}}(\boldsymbol{k},t)=\boldsymbol{\hat{b}}(\boldsymbol{k},t)-\boldsymbol{k}(\boldsymbol{k}\cdot\boldsymbol{\hat{b}}(\boldsymbol{k},t+\Delta
t))/(\boldsymbol{k}\cdot\boldsymbol{k}),$ (3)

which is needed to ensure the divergence–free condition. More details about
the HIT forcing can be found in Eswaran & Pope (1988); Chouippe & Uhlmann
(2015).

Note that equation (3) defines the forcing in Fourier space while equations
(1) are solved in the physical space. Therefore the former expression must be
transformed into values to be prescribed at the Eulerian grid points. To
determine the field in physical space
$\boldsymbol{f}_{{}_{T}}(\boldsymbol{x},t)$, we compute the Fourier
coefficients from equation (3), with
$\boldsymbol{\hat{f}}_{lmn}(t)=\boldsymbol{\hat{f}}(l,m,n,t)$, from which the
value at the discrete position
$\boldsymbol{x}_{i,j,k}=(i\Delta,j\Delta,k\Delta)$ is obtained by

$\boldsymbol{{f}}_{{}_{T}}(\boldsymbol{x}_{ijk},t)=\sum_{l=-N_{f}}^{N_{f}}\sum_{m=-N_{f}}^{N_{f}}\sum_{n=-N_{f}}^{N_{f}}\boldsymbol{\hat{f}}_{lmn}(t)\exp(\j\boldsymbol{k}\cdot\boldsymbol{x}_{i,j,k}),~{}~{}\forall
i,j,k$ (4)

where $N^{3}_{f}=27$ is the number of forced wavenumbers and $\j$ the
imaginary unit.

The HIT forcing and additional simulation parameters are reported in Table 1
where it is also shown that the spatial resolution is always smaller than the
Kolmogorov scale.

A validation of the flow solver is provided by computing the one–dimensional
longitudinal energy spectrum $E_{11}(k_{1})$ and the isotropy parameter

$I(k_{1})=\frac{E_{11}(k_{1})-k_{1}\partial E_{11}(k_{1})/\partial
k_{1}}{2E_{22}(k_{1})},$ (5)

for a case at $Re_{\lambda}\approx 140$, the latter being the Reynolds number
based on the Taylor micro–scale $\lambda=(15\nu u^{2}_{rms})^{\frac{1}{2}}$
and on the rms value of the velocity fluctuations $u_{rms}$.

The comparison with the data by Jiménez et al. (1993) and Chouippe & Uhlmann
(2015) is given in Figure 2, showing a very good agreement for both
quantities.

As an aside, we note that results obtained with a resolution $512^{3}$
($\eta/\Delta=0.92$) and $1024^{3}$ ($\eta/\Delta=1.83$) perfectly agree and,
in both cases, the viscous range shows an exponential decay, evidenced by the
linear slope when plotted in semi–log axes (Jiménez et al., 1993). The reason
for the sharp decay of the benchmark data in the viscous range is that they
had been obtained using the spectral code developed by Rogallo (1981), which
is very strict on isotropy. In fact, to approximate the flow as isotropically
as possible, it applies every time step a mask in $k$–space, such that the
energy outside the sphere $|k\eta|=2$ is set to zero (J. Jiménez, Personal
Communication).

Figure 2: Single phase HIT at $Re_{\lambda}=140$: (a) isotropy parameter $I$, as defined in equation (5), compared with Chouippe & Uhlmann (2015). (b) One–dimensional energy spectrum, compared with Jiménez et al. (1993). RUN | $T_{L}$ | $\sigma$ | $Re_{\lambda}$ | $K$ | $\varepsilon$ | $\lambda/D_{eq}$ | $\eta/\Delta$ | $N^{3}$
---|---|---|---|---|---|---|---|---
Re15 | $5.0\times 10^{-1}$ | $1.7\times 10^{-2}$ | $13.0$ | $8.4\times 10^{-3}$ | $2.7\times 10^{-4}$ | $1.74$ | $6.30$ | $256^{3}$
Re30 | $8.5\times 10^{-1}$ | $6.2\times 10^{-3}$ | $33.9$ | $1.1\times 10^{-1}$ | $7.5\times 10^{-3}$ | $1.22$ | $2.75$ | $256^{3}$
Re60 | $8.5\times 10^{-1}$ | $4.4\times 10^{-1}$ | $66.7$ | $1.3\times 10^{0}~{}~{}$ | $2.4\times 10^{-1}$ | $0.72$ | $2.31$ | $512^{3}$
Re80 | $8.2\times 10^{-2}$ | $2.6\times 10^{0}~{}~{}$ | $78.5$ | $2.2\times 10^{0}~{}~{}$ | $5.3\times 10^{-1}$ | $0.64$ | $1.90$ | $512^{3}$
Re140 | $8.5\times 10^{-2}$ | $8.3\times 10^{0}~{}~{}$ | $134.0$ | $1.7\times 10^{+1}$ | $1.0\times 10^{+1}$ | $0.40$ | $1.83$ | $1024^{3}$
Table 1: Run parameters of the single–phase simulations: $\lambda=(15\nu
u^{2}_{rms}/\varepsilon)^{1/2}$ is the Taylor micro–scale while $Re_{\lambda}$
is the corresponding Reynolds number, $\eta=(\nu^{3}/\varepsilon)^{1/4}$ the
Kolmogorov length scale, $N=N_{x}=N_{y}=N_{z}$ is the number of gridpoints of
the Eulerian mesh in each direction, $K$ is the kinetic energy, and $T_{L}$
and $\sigma$ are the input parameters for the HIT forcing.

###  2.4 Immersed boundary method

As mentioned above, the term $\boldsymbol{{\rm f}}$ on the right hand side of
the first equation (1) contains also the IBM part $\boldsymbol{f}_{{}_{P}}$
which allows for two–way coupling between particles and fluid. Here we have
employed the method proposed by de Tullio & Pascazio (2016) in combination
with the regularised Dirac delta function of Roma et al. (1999). The surface
of each solid object is discretised by $N_{l}$ equilateral triangular elements
(see Figure 1a) and Lagrangian markers are placed at their centroids. The IBM
forcing term is computed on these Lagrangian markers, to satisfy velocity
boundary conditions, and then transferred back to the Eulerian mesh to impose
the effect of the solid phase on the flow. In order to avoid unforced Eulerian
nodes, which would result in surface ‘holes’, the size of triangles edges must
not exceed $0.7\Delta$ (de Tullio & Pascazio, 2016) thus implying that, for
increasing $Re_{\lambda}$, not only the Eulerian mesh but also the Lagrangian
triangulation must be refined to maintain the correct phase coupling.
Accordingly, the values of $N_{l}$ were $4912$, $10048$ or $21262$, depending
on the flow resolution.

It is worth mentioning that $N_{p}$ particles settling under their own weight
impose on the fluid a force constantly aligned with gravity which, in a
tri–periodic domain, would produce a continuous vertical acceleration of the
flow and this would prevent the attainment of a statistical steady state.
Following Chouippe & Uhlmann (2015), this is avoided by subtracting, at every
time step, the mean of the IBM terms applied by the particles on the flow,
thus imposing a force with zero average in all directions: this implies that
the net momentum of the system must be on average zero. Thus the downward
motion of the falling particles is compensated by an upward motion of the
fluid such to give zero net vertical momentum (this is the analogous of the
time dependent pressure gradient in turbulent channel flows used to ensure a
constant flow rate (Kim et al., 1987) in time). An analogous observation
applies to angular momentum since free falling chiral particles rotate in a
preferred direction and this would generate a net circulation. In a
tri–periodic cube, however, the Stokes theorem imposes zero circulation in all
directions. Therefore the angular velocity of particles is compensated by
fluid vorticity in the opposite direction so to maintain the circulation zero
at every instant.

These observations must be kept in mind when analysing the results since while
the resulting system can be efficiently simulated within homogeneous and
isotropic conditions, an analogous laboratory experiment might not be
immediately possible to set up.

###  2.5 Particles

Figure 3: Representative quantities for the Stokes dynamics of a free–falling
sphere (a) and chiral particles (b) in a stagnant fluid. (a) Vertical velocity
(red) and drag force (blue) are compared with the reference Stokes values
(black). (b) Angular velocities in the body reference frame (green axes in
Figure 1) for a left–handed (blue) and right–handed (red) chiral particles.
Note that here, since the external turbulence forcing is absent, the velocity
scale $U$ is undefined; all results are therefore scaled using
$U_{g}=(|\rho_{p}/\rho_{f}-1||\boldsymbol{g}|D_{eq})^{1/2}$.

The dynamics of the dispersed solid phase is obtained by solving, for each
particle, the six degrees of freedom of the Newton–Euler equations

$\hskip-2.84544pt\Bigg{(}\rho_{p}V_{p}\frac{d\boldsymbol{u}_{c}}{dt}\\!=\\!\oint_{A_{p}}\\!\\!\\!\\!\boldsymbol{\tau}\cdot\boldsymbol{n}dA+(\rho_{p}-\rho_{f})V_{p}\boldsymbol{g}+\boldsymbol{F}_{c}\Bigg{)}^{*}\\!\\!,\
\
\Bigg{(}\boldsymbol{I}_{p}\frac{d\boldsymbol{\omega}_{c}}{dt}+\boldsymbol{\omega}_{c}\times(\boldsymbol{I}_{p}\boldsymbol{\omega}_{c})\\!=\\!\oint_{A_{p}}\\!\\!\\!\\!\boldsymbol{r}\times(\boldsymbol{\tau}\cdot\boldsymbol{n})dA+\boldsymbol{T}_{c}\Bigg{)}^{*}\\!\\!,$
(6)

where the notation $(\cdot)^{*}$ implies a content in dimensional units. Here
$A_{p}$ is the particle wet surface,
$\boldsymbol{\tau}=-p\boldsymbol{I}+\mu_{f}(\nabla\boldsymbol{u}+\nabla\boldsymbol{u}^{T})$
the Cauchy stress tensor for a Newtonian fluid with $\boldsymbol{I}$ the unit
tensor, $\boldsymbol{n}$ is the outward normal of the particle surface, and
$\boldsymbol{I}_{p}$ the inertial tensor of the particle which, in the
principal inertial reference frame (Figure 1a), has only the diagonal
components. $\boldsymbol{F}_{c}$ and $\boldsymbol{T}_{c}$ represent,
respectively, force and torque acting on the particle as result of
collisions/physical contact with other particles.

Following Breugem (2012) and relying on the third Newton’s law of motion, we
replace the integral terms on the right–hand–side of equations (6) with the
volume forcing of the fluid phase, obtaining in nondimensional form:

$\frac{d\boldsymbol{u}_{c}}{dt}\approx\frac{\rho_{f}}{\rho_{p}}\frac{6}{\pi}\Bigg{(}-\sum_{i=1}^{N_{l}}\boldsymbol{f}_{i}\Delta
V_{i}+\frac{d}{dt}\int_{V_{p}}\boldsymbol{u}dV\Bigg{)}-\frac{\widehat{{\bf
k}}}{Fr}+\frac{\rho_{f}}{\rho_{p}}\frac{6}{\pi}\boldsymbol{F}_{c},$
$\boldsymbol{I}_{p}\frac{d\boldsymbol{\omega}_{c}}{dt}+\boldsymbol{\omega}_{c}\times(\boldsymbol{I}_{p}\boldsymbol{\omega}_{c})\approx\frac{\rho_{f}}{\rho_{p}}\frac{6}{\pi}\Bigg{(}-\sum_{i=1}^{N_{l}}\boldsymbol{r}_{i}^{n}\times\boldsymbol{f}_{i}\Delta
V_{i}+\frac{d}{dt}\int_{V_{p}}\boldsymbol{r}\times\boldsymbol{u}dV\Bigg{)}+\frac{\rho_{f}}{\rho_{p}}\frac{6}{\pi}\boldsymbol{T}_{c},$
(7)

in which $\widehat{{\bf k}}$ is the vertical unit vector and

$Fr=\frac{U^{2}\rho_{p}}{D|\boldsymbol{g}|(\rho_{p}-\rho_{f})}$ (8)

the Froude number. In the above expression we note that surfaces are
triangulated and each triangle is tagged by a Lagrangian marker at its
centroid; thus the index $i$ indicates the Lagrangian point over the surface
and $\Delta V_{i}$ is the volume of the Eulerian cell intersected.

It is important to note that equation (7) clearly shows the particle dynamics
to depend on two nondimensional parameters, namely $\rho_{p}/\rho_{f}$ and
$Fr$, which in principle could be varied independently. Distinguishing between
$\rho_{p}/\rho_{f}$ and $Fr$, however, might look pleonastic since, in
standard laboratory experiments within the same gravity field,
$\rho_{p}/\rho_{f}$ and $Fr$ are not independent and equation (7) would in
fact depend only on a single parameter. Accordingly, for the vast majority of
our simulations we have fixed the magnitude of gravity ($|\boldsymbol{g}|=2$)
and varied only $\rho_{p}/\rho_{f}$ as independent parameter. This value of
$|\boldsymbol{g}|$ was initially selected since, for $\rho_{p}/\rho_{f}=2$, it
yielded $Fr=1$.

However, centrifuge setups (Jiang et al., 2020) or microgravity environments
(Futterer et al., 2013) could allow the independent variation of
$\boldsymbol{g}$, thus making relevant studying the separate effect of $Fr$
and $\rho_{p}/\rho_{f}$. In a few simulations, either of individual particles
or of particle crowds, $\rho_{p}/\rho_{f}$ and $Fr$ have been varied
independently, and a negligible effect of the former parameter compared to the
latter has been observed; Appendix A reports some cases to confirm that indeed
the results, for a fixed Froude number, hardly depend on density ratio
variations.

The above model has been validated by considering a free falling sphere in an
initially quiescent viscous fluid. With a density ratio $\rho_{p}/\rho_{f}=2$
($Fr=1$) and a fluid viscosity such to yield $Re\approx 1$ in the stationary
state, the Stokes solution is expected to apply and this has been used for
comparison. Domains up to $15\times 15\times 15$ in $D_{eq}$ units with a
resolution of $288^{3}$ nodes have been employed and the results are given in
Figure 3a for the vertical velocity $V$ and the drag force $F_{D}$; the
agreement between the present data and the analytical Stokes solution for
analogous quantities looks excellent, once the stationary state is attained.
Similar calculations for the torque always yielded values within the round–off
error as expected from the symmetry of the system geometry.

Figure 4: Initial transient dynamics of a chiral particle falling in a
stagnant fluid, for different initial orientations, at $Re_{p}\approx 10$. The
snapshots are taken at time lugs of $\Delta tD_{eq}/U_{g}=10$. For the results
of panels (b) and (c) the same scaling velocity as in Figure 3 has been used.

Unfortunately, exact results in the low–$Re$ regime, as in Figure 1a, are not
available for chiral particles. Therefore we could only run consistency
checks, showing opposite steady rotation for left– and right–handed particles.
More in detail, starting from any initial orientation, a free falling particle
undergoes a transient dynamics during which rotations occur about all three
axes until the mass–centroid attains the lowest possible position. Once this
stable attitude is gained, only a constant rotation about the $y_{p}$
(free–falling) axis remains which has opposite sign for different chirality of
the particle (Figure 3b). Additional simulations have revealed that the
long–term dynamics is independent of the initial orientation which, instead,
affects only the duration of the transient dynamics (Figure 4); this has
important consequences on the particle interaction with turbulence, which will
be discussed in the next Section.

Figure 5: Possible configurations for close particles; the volumwise
equivalent spheres are represented by transparent solids: (a) Particles with a
contact point but with non–intersecting spheres. (b) Particles not in contact
with intersecting spheres. The white segments in panel (a) are the axes of the
particles legs used to compute proximity and contact.

An important aspect of the present problem is the evaluation of
particle/particle collisions. In fact, we will see that, even if the volume
fraction of the solid phase never exceeds $\phi=2\%$, the complex shape of
particles involves collisions with multiple contacts and even entanglement.
Therefore the correct interaction dynamics is relevant for the statistics of
the whole system. Given the complex particle shape, determining interactions
is not trivial as two bodies can have their centroids farther than $D_{eq}$
even if there is contact (Figure 5a) or vice versa (Figure 5b).

In fact, the collision model consists of two main ingredients: the proximity
detection and the contact loads computation. Concerning the former, each
particle position and orientation is hard–coded through the coordinates of its
nodal points $P$, $Q$, $R$ and $S$ (Figure 5a), the diameter $d$ of each leg
(Figure 1a), and the radius $d/2$ of each distal spherical cap. For any couple
among the $N_{p}$ particles the centroids distance is computed and for those
whose value is compatible with a contact, the minimum distance among the
segments $h$ is evaluated. Knowing the particle shape it is easy to compute
the thickness of the fluid gap between the two surfaces $d_{c}$ and, if it is
smaller than two Eulerian grid cells, the contact loads are computed. The
reason for using $d_{c}\leq 2\Delta$ as threshold value is that, in IBM
methods, the flow at the first external point is not obtained from the
Navier–Stokes equations but from a model which enforces the boundary
conditions; if the gap between two immersed boundaries is thinner than
$2\Delta$ only a single Eulerian point is available to satisfy two boundary
conditions which is not a well–posed problem. Following Breugem (2012), the
short–range repulsive force, which the $j$–th particle exerts on the $i$–th,
is prescribed according to
$\boldsymbol{F}_{c,ij}=-\alpha_{c}[(|\boldsymbol{d}_{ij}|-d_{c})/d_{c}]^{2}\boldsymbol{d}_{ij}/|\boldsymbol{d}_{ij}|$
with $\boldsymbol{d}_{ij}$ the oriented minimum surface distance during the
collision. $\alpha_{c}=10^{-4}$ as in Breugem (2012) and
$\boldsymbol{F}_{c,ji}=-\boldsymbol{F}_{c,ij}$ as prescribed by the third
Newton’s law of motion. Finally, the contact force produces also a contact
moment which is immediately computed from
$\boldsymbol{T}_{c,ij}=\boldsymbol{r}_{c}\times\boldsymbol{F}_{c,ij}$, with
$\boldsymbol{r}_{c}$ being the vector connecting the particle centroid with
the contact point (Figure 5a).

We refer to Assen et al. (2022) for additional details and a thorough
validation of the collision model.

##  3 Results

###  3.1 Turbulence intensity effect

Figure 6: Single chiral particle with $\rho_{p}/\rho_{f}=2$ ($Fr=1$) in a
turbulent flow for the case Re15 of Table 1. (a) Time evolution of the centre
of mass velocity components. (b) Time evolution of angular velocity components
expressed in the particle body frame (green axes in figure 1).

In order to understand how chiral bodies interact with the surrounding flow,
we preliminarily consider the case of a single particle (volume fraction
$\phi=0.05\%$) in a turbulent environment of increasing intensity. We begin
with a particle at $\rho_{p}/\rho_{f}=2$ ($Fr=1$) in weak turbulence, as for
the case Re15 of Table 1. Figure 6 shows also in this case a net rotation
about the $y_{p}$ body axis sustained by the vertical falling through the
coupling of rotational and translational degrees of freedom.

Figure 7: Single chiral particle with $\rho_{p}/\rho_{f}=2$ ($Fr=1$) in a
turbulent flow for different turbulence strengths (see Table 1). (a) Time
evolution of the centre of mass vertical falling velocity. (b) Time evolution
of the angular velocity about the $y_{p}$ particle body axis (green axes in
Figure 1).

Here the dynamics does not differ much from that of Figure 3b, except for the
unsteadiness whose origin is twofold: particle settling velocity and
interaction with turbulence. In fact, $|\boldsymbol{u}_{c}|$ is proportional
to $U_{g}=(|\rho_{p}/\rho_{f}-1||\boldsymbol{g}|D_{eq})^{1/2}$, in turn
setting the particle Reynolds number $Re_{p}=U_{g}D_{eq}/\nu$ and, for high
enough values, this yields an unsteady falling even in a quiescent fluid. For
the present flow parameters it results $Re_{p}\approx 140$ and auxiliary
simulations have shown that indeed the settling dynamics is unsteady even in
absence of external perturbations. The interaction with turbulence, however,
introduces further unsteadiness as the particle dynamics is coupled with
velocity fluctuations whose strength is quantified by $Re_{\lambda}$. In this
case it results $Re_{\lambda}\approx 14.7$ with a velocity ratio
$u_{rms}/<|v_{z}|>\approx 0.08$, suggesting that turbulence acts only as a
perturbation on the particle dynamics
($v_{z}=\boldsymbol{u}_{c}\cdot\widehat{{\bf k}}$ is the falling particle
velocity).

It is interesting to note that for the parameters of case Re15, single–phase
turbulence yielded $Re_{\lambda}\approx 13$ while the introduction of a
particle raises the value by more than $10\%$; this is due to the constant
energy injection caused by the particle free fall (at the expenses of its
potential energy) which adds to the external forcing and increases the
turbulence strength. Though the increase is small, the phenomenon underlines
the two–way interaction between fluid and solid phases whose dynamics is
strongly connected.

As $Re_{\lambda}$ increases, the energy introduced by the particle becomes
less relevant to turbulence while the stronger velocity fluctuations alter the
solid body dynamics. In fact, for the parameters of case Re80 the effective
$Re_{\lambda}$ and $u_{rms}$ are indistinguishable from the single–phase
values while particle mean vertical and angular velocities become negligible
when compared to their fluctuations.

The results are summarised in Figure 7, which shows the mean vertical and
angular velocities, both normalised by their standard deviations. It is worth
pointing out that the ratio $v_{z}/\sigma_{v_{z}}$ of Figure 7a decreases with
$Re_{\lambda}$ both because the standard deviation $\sigma_{v_{z}}$ increases
and because $|v_{z}|$ itself decreases. On the other hand,
$|\dot{\theta}_{y}|$ remains approximately constant within this range of
$Re_{\lambda}$ although $\sigma_{\dot{\theta}_{y}}$ increases and therefore
their ratio decreases.

An important effect of particle chirality is its mean angular velocity which,
on account of the zero total circulation imposed by the domain boundary
conditions, induces an opposite mean vorticity in the flow field. This is
indeed observed from the probability density functions of the vertical fluid
vorticity component (Figure 8a) which are biased towards negative values and
become more symmetric as $Re_{\lambda}$ increases and particle rotation
becomes less relevant.

These results are consistent with those of Figure 7 suggesting that as
turbulence strengthens the particle tends to be passively advected by
turbulent fluctuations rather than altering turbulence itself through its
helical free fall.

A partial confirmation of this conjecture comes from Figure 8b, which reports
the ratio between $u_{rms}$ and $<|v_{z}|>$: when this ratio is small particle
dynamics dominates over turbulence which, therefore, is affected by particle
chirality. In contrast, large $Re_{\lambda}$ values results in
$u_{rms}/|v_{z}|>1$ and turbulent fluctuations dominate over particle forcing,
yielding a more isotropic dynamics.

Figure 8: Single chiral particle with $\rho_{p}/\rho_{f}=2$ ($Fr=1$) in a
turbulent flow at different $Re_{\lambda}$. (a) Probability density function
of the fluid vertical vorticity component. A preferred vertical vorticity, due
to the immersed chiral particles, is only obtained for weak turbulence. (b)
Ratio of turbulent fluctuation intensity and mean vertical particle velocity
as function of the effective flow $Re_{\lambda}$.

###  3.2 Density ratio effect

Given the role played by the relative magnitude of turbulence fluctuations
($u_{rms}$) and the particle falling velocity
$<|v_{z}|>\sim(|\rho_{p}/\rho_{f}-1||\boldsymbol{g}|D_{eq})^{1/2}$, it should
be clear that the system is sensitive not only to the strength of turbulence
forcing but also to the density ratio $\rho_{p}/\rho_{f}$ (or $Fr$). In fact,
the immediate consequence of increasing $\rho_{p}/\rho_{f}$ is to enhance the
falling velocity and therefore the mean angular velocity. However, translation
and rotation in a chiral particle are not rigidly coupled (as a nut and bolt
assembly), as both depend on hydrodynamic loads. Therefore they vary by
different amounts. Moreover, the potential energy transfer of the falling
particle grows with $\rho_{p}/\rho_{f}$, thus further increasing the
turbulence level which reacts back on the particle dynamics. Finally,
$\rho_{p}/\rho_{f}$ changes particle inertia to translation and rotation, in
turn altering the trajectory sensitivity to turbulent fluctuations.

Figure 9: Single chiral particle with $\rho_{p}/\rho_{f}$ from $2$ to $10$
($Fr$ from $1$ to $5/9$) in a turbulent flow with the forcing of Re15 case of
Table 1: (a) Time evolution of the centre of mass velocity components. (b)
Time evolution of angular velocity components expressed in the particle body
frame. In the inset the same curves are plotted with a running average of
$20D_{eq}/U$ time units to evidence their mean values.

Examples of particle dynamics are given in Figure 9 where results at
$\rho_{p}/\rho_{f}$ from $2$ to $10$ ($Fr\in[1,5/9]$) are compared for a
turbulent forcing as the Re15 case of Table 1: since for these flows the
interplay between mean values and standard deviations is more complex, we
complement Figure 9 with Table 2 where additional data are provided. A clear
trend can be observed for the mean vales of $<v_{z}>$ and $<\dot{\theta}_{y}>$
which both increase in magnitude with the ratio $\rho_{p}/\rho_{f}$. Also the
effect of energy transfer, from particle to turbulence, is clearly evidenced
from the growth of $Re_{\lambda}$ and $u_{rms}$. A similar monotonic increase
is found for the mean vertical component of the fluid vorticity which is a
direct consequence of the chirality induced angular velocity of the particle.

A closer look at the data, however, reveals that the fivefold increase of
$\rho_{p}/\rho_{f}$ produces growths of $<v_{z}>$ of $\approx 3.6$ and
$<\dot{\theta}_{y}>$ of $\approx 2.7$, thus confirming the loose coupling of
particle rotation and translation. Concerning the effect on turbulence,
$<\omega_{z}>/\sigma_{\omega_{z}}$ grows by $\approx 1.6$. This is confirmed
by the histograms of Figure 10a which drift in the direction of negative
vorticity consistently with the monotonic increase of $<\dot{\theta}_{y}>$. On
the other hand, despite the growth of $Re_{\lambda}$, the ratio
$u_{rms}/<|v_{z}|>$ decreases with $\rho_{p}/\rho_{f}$, consistently with the
observation that particles of larger density are less sensitive to turbulent
fluctuations and their chirality makes them behave as impellers which impart
rotation to the surrounding fluid.

$\rho_{p}/\rho_{f}$ | $<v_{z}>$ | $<\dot{\theta}_{y}>$ | $Re_{\lambda}$ | $u_{rms}$ | $<\omega_{z}>/\sigma_{\omega_{z}}$
---|---|---|---|---|---
$2$ | $-1.41$ | $-0.17$ | $14.7$ | $0.12$ | $-0.74$
$4$ | $-2.80$ | $-0.29$ | $20.4$ | $0.21$ | $-0.67$
$6$ | $-3.73$ | $-0.36$ | $23.4$ | $0.28$ | $-0.93$
$10$ | $-5.16$ | $-0.46$ | $27.4$ | $0.38$ | $-1.17$
Table 2: Flow parameters of single chiral particle simulations with variable
$\rho_{p}/\rho_{f}$ and turbulent forcing as case Re15 of Table 1.  Figure 10:
Single chiral particle with $\rho_{p}/\rho_{f}$ from $2$ to $10$ ($Fr$ from
$1$ to $5/9$) in a turbulent flow with the forcing of Re15 case of Table 1:
(a) Probability density function of the fluid vertical vorticity component.
(b) Ratio of turbulent fluctuation intensity and mean vertical particle
velocity as function of $\rho_{p}/\rho_{f}$.

We have repeated the same simulations as before, increasing the turbulence
forcing as for the case Re30 of Table 1. For the sake of conciseness, we
report only the results of Figure 11, which shows an effect of stronger
turbulence only on the vorticity histograms of the lighter particles.
Accordingly the ratio $u_{rms}/<|v_{z}|>$ shows a substantial increase for
$\rho_{p}/\rho_{f}=2$ and $4$ ($Fr=1$ and $2/3$) while it is essentially
unchanged for the higher density ratios.

Figure 11:  Single chiral particle with $\rho_{p}/\rho_{f}$ from $2$ to $10$
($Fr$ from $1$ to $5/9$) in a turbulent flow with the forcing of Re30 case of
Table 1: (a) Probability density function of the fluid vertical vorticity
component. (b) Ratio of turbulent fluctuation intensity and mean vertical
particle velocity as function of $\rho_{p}/\rho_{f}$. Note the very different
vertical scale in Figures 11b and 10b.

The picture emerging from the simulation campaigns with a single chiral
particle is that its fall injects additional energy in the fluid whose
turbulence is altered; in turn, particle falling and rotational velocities are
affected by the strength of velocity fluctuations. Therefore the system
dynamics emerges out of the balance of a complex two–way coupling between
solid and fluid phases. Nevertheless, a clear trend is observed such that
stronger velocity fluctuations reduce the effects of particle and its
chirality on turbulence, while higher particle/fluid density ratios increase
them.

###  3.3 Volume fraction effect

In real systems, such as helically swimming active matter (Liebchen & Levis,
2022) or plants aiming at maximising seeds transport and dispersal (Mazzolai
et al., 2021), multiple chiral particles interact simultaneously among them
and react back on the surrounding fluid. Accordingly, we have investigated how
the system dynamics change when several particles are placed in a turbulent
environment and what is the effect of their volume fraction $\phi$.

The immediate effect of increasing the number of particles is to enhance the
forcing on the turbulence, owing to growing energy injection in the fluid
phase. However, the presence of multiple particles increases also energy
dissipation, by friction with solid surfaces and the small scale of their
wakes. Therefore turbulence can be strengthened or suppressed, depending on
which of the two mechanisms prevails.

Furthermore, the presence of obstacles in the bulk flow reduces the size
attainable by the largest flow scales resulting in further modifications of
the energy spectrum. Finally, the collision among complex shape particles
favours prolongated interactions rather than impulsive rebounds, as in the
case of spheres or other compact objects, and this changes the solid body
dynamics and the interaction with the flow.

$N_{p}$ | $\phi\ (\%)$ | $<v_{z}>$ | $<\dot{\theta}_{y}>$ | $Re_{\lambda}$ | $u_{rms}$ | $\varepsilon$ | $<\omega_{z}>/\sigma_{\omega_{z}}$
---|---|---|---|---|---|---|---
$0$ | $0$ | – | – | $33.9$ | $0.27$ | $7.5\times 10^{-3}$ | $0.00$
$1$ | $0.05$ | $-1.18$ | $-0.17$ | $31.9$ | $0.27$ | $8.0\times 10^{-3}$ | $-0.43$
$5$ | $0.25$ | $-1.17$ | $-0.17$ | $29.5$ | $0.29$ | $1.2\times 10^{-2}$ | $-0.86$
$10$ | $0.50$ | $-1.14$ | $-0.17$ | $28.1$ | $0.30$ | $1.6\times 10^{-2}$ | $-1.21$
$20$ | $1.00$ | $-1.10$ | $-0.16$ | $24.4$ | $0.31$ | $2.5\times 10^{-2}$ | $-1.59$
$40$ | $2.00$ | $-1.03$ | $-0.16$ | $22.9$ | $0.34$ | $4.0\times 10^{-2}$ | $-1.67$
Table 3: Flow parameters for multiple particle simulations with variable
volume fraction $\phi$, $\rho_{p}/\rho_{f}=2$ ($Fr=1$) and turbulent forcing
as case Re30 of Table 1.  Figure 12:  Multiple chiral particles with $\phi$
from $0.05\%$ to $2\%$ in a turbulent flow with the forcing of Re30 case of
Table 1 and $\rho_{p}/\rho_{f}=2$: (a) Kinetic energy dissipation rate as
function of the volume fraction $\phi$. (b) Comparison of one–dimensional
energy spectra for single–phase, one chiral particle, $20$ chiral particles
and $20$ spheres.

Table 3 reports some representative data for flows with a turbulent forcing as
in case Re30 of Table 1 and an increasing number of particles with
$\rho_{p}/\rho_{f}=2$ ($Fr=1$). It is indeed confirmed that a larger number of
particles evolving within the same domain enhances velocity fluctuations,
consistently with the augmented energy injection; however, also viscous
dissipation increases and this happens at a faster rate than the increase of
$u_{rms}$. The dependence of $\varepsilon$ on $\phi$ is displayed in Figure
12a, which shows a linear behaviour which matches the findings of Fornari et
al. (2019) for finite size spheres in turbulence. The different growths of
$\varepsilon$ and $u_{rms}$ with $\phi$ result in a decreasing $Re_{\lambda}$
indicating that, for this setup, more particles weaken turbulence by
transferring more energy to small scales. This mechanism is quantitatively
confirmed by the spectra of Figure 12b, which shows an increased energy
content at the smallest scales and a reduced content at larger scales, which
is similar to what was observed by Gao et al. (2013) in HIT flows with spheres
of diameter $\approx\lambda$. The spectra of Figure 12b also show that, for
given $\phi$, the present chiral particles are much more effective than
spheres in modifying turbulence. This might be due to the angular momentum of
the dispersed phase which in turn induces a net rotation of the fluid
turbulence. Indeed, Biferale et al. (2013) have already shown that, in HIT,
altering the nonlinear terms so to have a sign–definite helicity content
heavily affects the energy cascade and turbulence spectra.

An additional reason for the marked difference between spheres and chiral
particles could is the collision dynamics which is illustrated in Figure 13. A
set of $N_{p}$ particles yields $N_{p}(N_{p}-1)/2$ distinct couples, if pairs
$\\#i$–$\\#j$ and $\\#j$–$\\#i$ are equivalent. Whenever particles $\\#i$ and
$\\#j$ are close enough to activate the collision model, we assemble a
four–digit integer $\\#ij$ from particles id (for example, particles $\\#i=12$
and $\\#j=19$ give $\\#ij=1219$. $\\#ij$ is always formed with $\\#i<\\#j$.)
and record the time of interaction. Figure 13 shows, in a representative time
window, the interactions for $N_{p}=20$ spheres and chiral particles with the
same volume fraction $\phi=1\%$ for two flows with different turbulent forcing
and density ratios $\rho_{p}/\rho_{f}$. It is evident that spheres tend to
have short duration impulsive rebounds while chiral particles, with their
complex shape, favour entanglement which entail multiple contacts over several
time units. Analysis of the flow fields has evidenced an energy dissipation
enhancement caused by particle collision which is consistent with the energy
spectra of Figure 12b.

Clearly, frequency and duration of the collisions depend strongly on
turbulence strength and particle density (Figure 13). The comprehension of
their statistics deserves a dedicated study which will be presented in a
forthcoming paper.

Figure 13:  Collision events versus time for $N_{p}=20$ particles ($\phi=1\%$)
in a turbulent flow: a) forcing as Re30 case of Table 1 and
$\rho_{p}/\rho_{f}=2$; b) forcing as Re60 case of Table 1 and
$\rho_{p}/\rho_{f}=7$; the four–digit number on the y–axis identifies the two
interacting particles: for example, $\\#i,j=0412$ indicates that particles
$\\#i=4$ and $\\#j=12$ have collided. Red bullets for spherical particles,
blue open circles for chiral particles.

Concerning the particle angular velocity, Table 3 shows that it is essentially
constant with $\phi$ while the amplitude of the mean fluid vorticity increases
(Figure 14a); this can be explained by considering that $<\dot{\theta}_{y}>$
is the mean angular velocity of each particle. When their number increases the
fluid must rotate faster to compensate for the total angular momentum which
yields a global zero circulation.

The values of $<v_{z}>$ slightly decrease with $\phi$. This result agrees with
the findings of Fornari et al. (2016) for spheres at similar volume fractions.
The ratio $u_{rms}/<|v_{z}|>$ is however well below unity for all volume
fractions (Figure 14b) and the vorticity histograms of Figure 14a behave
consistently, showing always asymmetric vorticity distributions.

Figure 14:  Multiple chiral particles with $\phi$ from $0.05\%$ to $2\%$ in a
turbulent flow with the forcing of Re30 case of Table 1 and
$\rho_{p}/\rho_{f}=2$: (a) Probability density function of the fluid vertical
vorticity component. (b) Ratio of turbulent fluctuation intensity and mean
vertical particle velocity as function of the volume fraction $\phi$.

After having investigated the effects of the single governing parameters
separately, we can now vary them simultaneously and verify whether the system
behaves according to the gained understanding. The picture emerging from the
previous simulations is that particles fall because of gravity and increase
turbulent fluctuations, by agitating the surrounding fluid at the expenses of
their potential energy. At the same time, also turbulent energy dissipation is
enhanced. Thus the overall turbulence strength ($Re_{\lambda}$) results from
the balance of these two competing effects. On the other hand, particles
settling velocity decreases as the turbulence strengthens, thus giving a
two–way coupled dynamics. The chiral shape of the present particles couples
rotational and translational degrees of freedom which entails a net angular
momentum of the solid phase and, correspondingly, a mean vorticity of the
fluid turbulence. Finally, increasing the volume fractions $\phi$ of the
chiral particles enhances the previous effects but introduces also a complex
collision dynamics which further complicates the interactions.

On account of the above arguments, the relative strength of particle forcing
and turbulence intensity can be quantified by the ratio $u_{rms}/<|v_{z}|>$,
while the chirality–induced flow rotation through the mean vorticity is
normalised by its standard deviation $<\omega_{z}>/\sigma_{\omega_{z}}$.
Several simulation campaigns for different values of $Re_{\lambda}$ and
$\rho_{p}/\rho_{f}$ have been performed. The results shown in Figure 15
confirm that indeed the behaviour is well described by these quantities for a
wide range of governing parameters. Some additional ad hoc cases have been run
using different values of $\boldsymbol{g}$ in order to vary independently $Fr$
and $\rho_{p}/\rho_{f}$: also these data are reported in Figure 15. They
collapse with the rest of the simulations in which only $Fr$ has been varied
as main parameter.

It is important to note that most of the flows of Figure 15 have been computed
with $N_{p}=20$ particles ($\phi=1\%$) although we have verified that the same
plot for all cases obtained using $N_{p}=1$ ($\phi=0.05\%$) behaves similarly
even if data collapse onto a different curve. On the other hand, Table 3 and
the related discussion have evidenced that $\phi$ affects the flow features in
several ways and parametrising its mechanisms is not trivial. However, we can
attempt to estimate the effect of $\phi$ from Table 3, which has been obtained
for a single value of $Re_{\lambda}$ and $Fr$, and verify whether the same
correction works for all other flows. More in detail, from Table 3 the
quantities $u_{rms}/<|v_{z}|>$ and $<\omega_{z}>/\sigma_{\omega_{z}}$ at
$N_{p}=1$ and $N_{p}=20$ have been employed to compute their ratios to be used
as correction coefficients for all cases with $N_{p}=1$ computed for different
values of $Re_{\lambda}$ and $Fr$. Indeed, the open symbols in Figure 15
confirm that the behaviour is the same and all data now collapse on the same
curve.

Figure 15:  Mean flow vorticity normalised by its standard deviation versus
ratio of turbulent velocity fluctuations with mean settling particle velocity
in linear (a) and logarithmic scale (b). The black solid line is the curve
$0.45(u_{rms}/v_{z})^{-1}$ which best fits all the simulations.

A final comment is in order concerning the huge error bar of the rightmost
points in Figure 15a: in this case the strength of turbulence (produced by the
forcing $\boldsymbol{f}_{T}$ alone) is such to make negligible the effect of
the particle on the flow whose dynamics tend to that of a passive tracer.
Accordingly, the settling velocity and the induced flow rotation should vanish
while turbulent fluctuations strengthen. We have found that achieving
statistical convergence for these cases requires very long simulation times
(see also the large variations of the ratio $T_{o}/T_{e}$ reported in section
2.1) and for $u_{rms}/<|v_{z}|>>2.5$ obtaining convergence becomes
unpractical. We have however verified that, regardless of turbulence
intensity, simulations with $\rho_{p}/\rho_{f}=1$ always yield
$<|v_{z}|>\rightarrow 0$ and $<\omega_{z}>\rightarrow 0$, thus confirming that
the data of Figure 15 approach the correct limit.

##  4 Conclusions and outlook

Using direct numerical simulations of homogeneous isotropic turbulence with
finite–size, heavy chiral particles, we have investigated the dynamics of the
system under the influence of the main governing parameters which are:
turbulence strength, solid–to–fluid density ratio (or Froude number) and
volume fraction of the solid phase.

A particle, falling under the effect of an external gravity field, converts
its potential energy into a kinetic counterpart which is transferred to the
fluid through viscous and pressure forces. On the one hand this mechanism
enhances the intensity of velocity fluctuations but, on the other hand, it
also increases the energy dissipation through viscous friction and the
generation of small scales in the wake. Depending on the specific flow
parameters, the former or the latter mechanisms can prevail, thus boosting or
suppressing turbulence intensity, respectively.

Turbulence of increasing strength, in turn, affects particle dynamics by
reducing its mean falling velocity and decreasing the relevance of the
additional energy input with respect to that introduced by the HIT forcing
alone.

Particle chirality, by coupling translational and rotational degrees of
freedom, produces a mean flow vorticity which is evidenced through its
probability density function with a non–zero peak position.

Our results suggest that the ratio of turbulent velocity fluctuations to the
mean particle settling velocity $u_{rms}/v_{z}$ can be used to characterise
the intensity of the chirality induced mean rotation of the flow which
decreases as the previous ratio increases: the implication of this finding is
that particle chirality does not matter much in strong turbulence.

The influence of the particle volume fraction $\phi$ is more complex than the
other parameters since, in addition to the above mentioned mechanisms, there
are the interactions among particles and their collisions. The latter have
shown to be characterised by long entanglement times in contrast to spheres
which instead rebound impulsively. These differences are evident also from the
turbulent energy spectra which, for a given volume fraction ($\phi=1\%$), show
a more pronounced energy reduction caused by chiral particles than by spheres.
Despite the more complex dynamics, yet the quantity $u_{rms}/v_{z}$ describes
correctly the induced flow rotation as confirmed by Figure 15.

This paper is only a first attempt to investigate the interaction of
relatively complex shape chiral particles with homogeneous isotropic
turbulence and many things could be done to extend the study in different
ways: more detailed flow statistics should be computed in order to understand
how turbulence is modified not only with respect to rotation. Heavier and
lighter spheres in HIT (such as droplets in gases and droplets in liquids,
respectively), owing to centrifugal effects, are known to collect in different
flow regions. Thus chiral particles lighter than the surrounding fluid might
behave differently from the present ones. In this paper we have verified that
left– and right–handed particles behave identically except for the sign of the
induced flow rotation. A $50\%$ mixture of left– and right–handed particles is
expected to produce zero net rotation of the flow although it is likely to
alter other turbulence properties very differently from a set of spheres with
the same $\phi$. Finally, the statistics of particle collisions have shown to
strongly depend on the various flow parameters (Figure 13). A deeper
understanding of their features would be very useful to characterise the
dynamics of real suspensions.

All these points are currently under investigation or will be dealt with in
the near future. They will be the subject of forthcoming papers.

[Acknowledgements]We wish to thank Mr. Xander M. de Wit and Prof. Federico
Toschi for many fruitful clarifying discussions about the dynamics of chiral
bodies and their interaction with turbulence.

[Funding]This research has received funding from the Dutch Research Council
under the project “Shaping turbulence with smart particles” (file number
OCENW.GROOT.2019.03). We acknowledge the access to several computational
resources, all of which were used for this work: European High Performance
Computing Joint Undertaking for awarding us access to Discoverer under the
project EHPC-REG-2022R03-208 and the national e-infrastructure of SURFsara, a
subsidiary of SURF cooperation.

[Declaration of interests]The authors report no conflict of interest.

[Data availability statement]The data that support the findings of this study
are available upon reasonable request.

##  Appendix A Density ratio independence if $Fr=1$

In this appendix, we provide evidence of the reduced sensitivity of equations
(7) to density ratio compared to the Froude number. Figure 16 reports some
cases run for $2\leq\rho_{p}/\rho_{f}\leq 20$ and $Fr=1$, showing only limited
differences. Indeed, some effects are evident for the settling velocity,
especially for the lowest density ratio. These however have to be compared
with the results of Figure 9 where a smaller variation of $Fr$ produces much
larger effects.

As an aside we note that if
$U_{g}=\sqrt{D(1-\rho_{f}/\rho_{p})|\boldsymbol{g}|}$ is used as scaling
velocity for all equations (rather than the flow velocity $U$) the Froude
number turns out to be always $Fr=1$ and the resulting system show only
limited sensitivity to the density ratio $\rho_{p}/\rho_{f}$.

Figure 16:  Single chiral particle falling in a stagnant fluid at
$Re_{p}\approx 100$ for different $\rho_{p}/\rho_{f}$ and fixed $Fr=1$. (a)
Time evolution of the centre of mass vertical velocity component. (b) Time
evolution of angular velocity $y-$component expressed in the particle body
frame (green axes in figure 1).

##  Appendix B Kinetic energy balance

As a further consistency check of the numerics, we consider the energy budget
obtained from equations (1) after a scalar product with $\mathbb{u}$ and
integrating over the entire domain $\Omega$. We obtain, for the kinetic energy
$E_{k}=\mathbb{u\cdot u}/2$, the following relation:

$\frac{d}{dt}\langle
E_{k}\rangle_{{}_{\Omega}}=-\varepsilon_{\Omega}+\Pi_{{}_{T}}+\Pi_{{}_{P}}$
(9)

with the dissipation rate

$\varepsilon_{\Omega}=\frac{1}{Re}\Biggl{\langle}\Bigg{|}\frac{\partial\mathbb{u}}{\partial\mathbb{x}}\Bigg{|}^{2}\Biggr{\rangle}_{\Omega},$
(10)

the work done by the turbulence forcing $\Pi_{{}_{T}}$

$\Pi_{{}_{T}}=\langle\mathbb{u\cdot}\boldsymbol{f}_{{}_{T}}\rangle_{{}_{\Omega}},$
(11)

and the fluid–particle coupling term $\Pi_{{}_{P}}$,

$\Pi_{{}_{P}}=\langle\mathbb{u\cdot}\boldsymbol{f}_{{}_{P}}\rangle_{{}_{\Omega}}.$
(12)

Figure17 shows the time evolution of the left– and right–hand side of equation
(9), confirming the balance.

Figure 17: Time evolution of left– and right–hand side of equation (9) for
$N_{p}=20$ chiral particles of $\rho_{p}/\rho_{f}=2$ in the forcing case Re30
of Table 1.The right–hand side terms are shown both individually (a) and in
sum (b).

##  Appendix C Computational domain size in gravity direction

Here we provide the results of additional simulations run in a computational
domain whose size has been doubled in the gravity direction
($L_{z}=20D_{eq}$). These tests were run to ascertain a possible influence of
the vertical domain size on account of particles settling velocity. Figure 18
shows the time evolution settling and angular velocity for a single particle
at $\rho_{p}/\rho_{f}=2$ in a turbulent flow with the forcing as case Re30 of
Table 1: their mean values are $v_{z}/\sigma_{v_{z}}=-6.7\pm 1.0$, $-2.3\pm
0.3$ for $L_{z}=20D_{eq}$ and $v_{z}/\sigma_{v_{z}}=-5.8\pm 1.1$, $-2.1\pm
0.3$ for $L_{z}=10D_{eq}$, thus confirming that the cubic domain is indeed
adequate for our analyses.

In Figure 19 we report the results of a similar test performed for a crowd of
particles at $\rho_{p}/\rho_{f}=2$ in a turbulent flow with the forcing as
case Re30 of Table 1 and $\phi=1\%$ ($N_{p}=40$); once again the results show
very good agreement with the analogous ones obtained in a cubic computational
domain.

Figure 18:  Single chiral particle of $\rho_{p}/\rho_{f}=2$ in a turbulent
flow with the forcing of Re30 case of Table1: (a) Time evolution of the centre
of mass vertical velocity component. (b) Time evolution of angular velocity
$y-$component expressed in the particle body frame (green axes in Figure 1).
Figure 19:  Crowd of chiral particles of $\rho_{p}/\rho_{f}=2$ in a turbulent
flow with the forcing of Re30 case of Table1: (a) Time evolution of the centre
of mass vertical velocity component. (b) Time evolution of angular velocity
$y-$component expressed in the particle body frame (green lines in Figure 1),
(c) Probability density function of the fluid vertical vorticity component,
(d) streamlines of the nondimensional velocity field coloured according to the
vertical component ranging from $-0.3$ (blue) to $0.3$ (yellow).

## References

  * Assen et al. (2022) Assen, M.P.A., Ng, C.S., Will, J.B., Stevens, R.J.A.M., Lohse, D. & Verzicco, R. 2022 Strong alignment of prolate ellipsoids in Taylor-Couette flow. Journal of Fluid Mechanics 935, A7.
  * Balmforth & Provenzale (2001) Balmforth, N. & Provenzale, A. 2001 Geomorphological fluid mechanics, 1st edn. Springer Berlin, Heidelberg.
  * Bec et al. (2024) Bec, J., Gustavsson, K. & Mehlig, B. 2024 Statistical models for the dynamics of heavy particles in turbulence. Annual Review of Fluid Mechanics 56 (1), 189–213.
  * Biferale et al. (2013) Biferale, L., Musacchio, S. & Toschi, F. 2013 Split energy–helicity cascades in three–dimensional homogeneous and isotropic turbulence. Journal of Fluid Mechanics 730, 309–327.
  * Brandt & Coletti (2022) Brandt, Luca & Coletti, Filippo 2022 Particle-laden turbulence: Progress and perspectives. Annual Review of Fluid Mechanics 54 (1), 159–189, arXiv: https://doi.org/10.1146/annurev-fluid-030121-021103.
  * Breugem (2012) Breugem, W.-P. 2012 A second-order accurate immersed boundary method for fully resolved simulations of particle-laden flows. Journal of Computational Physics 231 (13), 4469–4498.
  * Chouippe & Uhlmann (2015) Chouippe, A. & Uhlmann, M. 2015 Forcing homogeneous turbulence in direct numerical simulation of particulate flow with interface resolution and gravity. Physics of Fluids 27 (12), 123301, arXiv: https://doi.org/10.1063/1.4936274.
  * de Tullio & Pascazio (2016) de Tullio, M.D. & Pascazio, G. 2016 A moving-least-squares immersed boundary method for simulating the fluid–structure interaction of elastic bodies with arbitrary thickness. Journal of Computational Physics 325, 201–225.
  * Dotto & Marchioli (2019) Dotto, L & Marchioli, C 2019 Orientation, distribution, and deformation of inertial flexible fibers in turbulent channel flow. Acta Mech 230, 597–621.
  * Eswaran & Pope (1988) Eswaran, V. & Pope, S.B. 1988 An examination of forcing in direct numerical simulations of turbulence. Computers & Fluids 16 (3), 257–278.
  * Fornari et al. (2016) Fornari, W., Picano, F. & Brandt, L. 2016 Sedimentation of finite–size spheres in quiescent and turbulent environments. Journal of Fluid Mechanics 788, 640–669.
  * Fornari et al. (2019) Fornari, W., Zade, S., Brandt, L. & Picano, F. 2019 Settling of finite-size particles in turbulence at different volume fractions. Acta Mechanica 230, 413–430.
  * Futterer et al. (2013) Futterer, B., Krebs, A., Plesa, A.-C., Zaussinger, F., Hollerbach, R., Breuer, D. & Egbers, C. 2013 Sheet-like and plume-like thermal flow in a spherical convection experiment performed under microgravity. Journal of Fluid Mechanics 735, 647–683.
  * Gao et al. (2013) Gao, H., Li, H. & Wang, L.-P. 2013 Lattice boltzmann simulation of turbulent flow laden with finite-size particles. Computer and Mathematics with Applications 65, 194–210.
  * Hargrave (1985) Hargrave, B.T. 1985 Particle sedimentation in the ocean. Ecological Modelling 30 (3), 229–246.
  * Jiang et al. (2020) Jiang, H., Zu, X., Wang, D., Huisman, S.G. & Sun, C. 2020 Supergravitational turbulent thermal convection. Science Advances 6, eabb8676.
  * Jiménez et al. (1993) Jiménez, J., Wray, A., Saffman, P.G. & Rogallo, R.S. 1993 The structure of intense vorticity in isotropic turbulence. Journal of Fluid Mechanics 255, 65–90.
  * Kim et al. (1987) Kim, J., Moin, P. & Moser, R. 1987 Turbulence statistics in fully developed channel flow at low reynolds number. Journal of Fluid Mechanics 177, 133–166.
  * Lashgari et al. (2014) Lashgari, I., Picano, F., Breugem, W.-P. & Brandt, L. 2014 Laminar, turbulent, and inertial shear-thickening regimes in channel flow of neutrally buoyant particle suspensions. Physical Review Letters 113, 254502.
  * Liebchen & Levis (2022) Liebchen, B. & Levis, D. 2022 Chiral active matter. Europhysics Letters 65, 67001.
  * Mazzolai et al. (2021) Mazzolai, B., Mariani, S., Ronzan, M., Cecchini, L., Fiorello, I., Cikalleshi, K. & Margheri, L. 2021 Morphological computation in plant seeds for a new generation of self-burial and flying soft robots. Frontiers 8.
  * Necker et al. (2005) Necker, F., Härtel, C., Kleiser, L. & Meiburg, E. 2005 Mixing and dissipation in particle-driven gravity currents. Journal of Fluid Mechanics 545, 339–372.
  * Olivieri et al. (2022) Olivieri, S., Cannon, I. & Rosti, M. E. 2022 The effect of particle anisotropy on the modulation of turbulent flows. Journal of Fluid Mechanics 950, R2.
  * van der Poel et al. (2015) van der Poel, E. P., Ostilla-Mónico, R., Donners, J. & Verzicco, R. 2015 A pencil distributed finite difference code for strongly turbulent wall-bounded flows. Computers & Fluids 116, 10–16.
  * Rai & Moin (1991) Rai, M. M. & Moin, P. 1991 Direct simulations of turbulent flow using finite-difference schemes. Journal of Computational Physics 96 (1), 15–53.
  * Rogallo (1981) Rogallo, R.S. 1981 Numerical experiments in homogeneous turbulence. NASA TM-81315 https://ntrs.nasa.gov/citations/19810022965.
  * Roma et al. (1999) Roma, A.M., Peskin, C.S. & Berger, M.J. 1999 An adaptive version of the immersed boundary method. Journal of Computational Physics 153 (2).
  * Roos & Roos (2015) Roos, G. & Roos, C. 2015 Chapter 3 - isomers and stereochemistry. In Organic Chemistry Concepts (ed. G. Roos & C. Roos), pp. 43–54. Boston: Academic Press.
  * Shaw (2003) Shaw, R.A. 2003 Particle-turbulence interactions in atmospheric clouds. Annual Review of Fluid Mechanics 35 (1), 183–227.
  * Sozza et al. (2022) Sozza, A., Cencini, M., Musacchio, S. & Boffetta, G. 2022 Instability of a dusty kolmogorov flow. Journal of Fluid Mechanics 931, A26.
  * Spandan et al. (2020) Spandan, V., Putt, D., Ostilla-Mónico, R. & Lee, A. A. 2020 Fluctuation-induced force in homogeneous isotropic turbulence. Science Advances 6, eaba0461.
  * Toschi & Bodenschatz (2009) Toschi, Federico & Bodenschatz, Eberhard 2009 Lagrangian properties of particles in turbulence. Annual Review of Fluid Mechanics 41 (1), 375–404, arXiv: https://doi.org/10.1146/annurev.fluid.010908.165210.
  * Verzicco & Orlandi (1996) Verzicco, R. & Orlandi, P. 1996 A finite-difference scheme for three-dimensional incompressible flows in cylindrical coordinates. Journal of Computational Physics 123 (2), 402–414.
  * Voilland (2018) Voilland, A. 2018 Here comes the saharan dust. Https://earthobservatory.nasa.gov/images/92358/here-comes-the-saharan-dust.
  * Voth & Soldati (2017) Voth, Greg A. & Soldati, Alfredo 2017 Anisotropic particles in turbulence. Annual Review of Fluid Mechanics 49 (1), 249–276, arXiv: https://doi.org/10.1146/annurev-fluid-010816-060135.
  * Yu et al. (2021) Yu, Z., Xia, Y., Guo, Y. & Lin, J. 2021 Modulation of turbulence intensity by heavy finite-size particles in upward channel flow. Journal of Fluid Mechanics 913, A3.
  * Zhao et al. (2010) Zhao, L.H., Andersson, H.I. & Gillissen, J.J.J. 2010 Turbulence modulation and drag reduction by spherical particles. Physics of Fluids 22 (8), 081702, arXiv: https://doi.org/10.1063/1.3478308.