# Numerical study on turbulence modulation of finite-size particles in plane-
Couette flow

Cheng Wang1  Linfeng Jiang1,2<EMAIL_ADDRESS>Chao Sun1,3
<EMAIL_ADDRESS>1Center for Combustion Energy, Key Laboratory for
Thermal Science and Power Engineering of Ministry of Education, International
Joint Laboratory on Low Carbon Clean Energy Innovation, Department of Energy
and Power Engineering, Tsinghua University, Beijing 100084, China 2Physics of
Fluids Group, and Max Planck UT Center for Complex Fluid Dynamics, University
of Twente, 7500 AE Enschede, The Netherlands 3Department of Engineering
Mechanics, School of Aerospace Engineering, Tsinghua University, Beijing
100084, China

###### Abstract

Turbulent plane-Couette flow suspended with finite-size spheroidal particles
is studied using fully particle-resolved direct numerical simulations. The
effects of particle aspect ratio on turbulent arguments and particle
statistics are explored, leading to the same conclusions as the previous
experimental findings (Wang et al., 2022). By performing stress analysis, we
find that the presence of particles introduces extra stresses to the system
and accounts for the global drag increases. The particle-laden flow cases
exhibit spectra that are consistent with the scaling of $k^{-5/3}$ and
$k^{-3}$ in the large and small scales, respectively. While the $k^{-3}$
scaling observed in the particle-laden flow is reminiscent of bubbly flow, an
examination of the particle Reynolds number suggests that the mechanism
responsible may not be attributable to the pseudo-turbulence induced by
particles as in the case of bubbles. In the view of particle statistics, we
observe that spherical and non-spherical particles preferentially cluster in
the near-wall and the bulk region, respectively, and that the orientations of
non-spherical particles are affected by their aspect ratios, especially in the
near-wall region. The present numerical results, combined with previous
experimental findings in Wang et al. (2022), provide in-depth information on
both the fluid and the particle phase, contributing to a better understanding
of particle suspension in shear flows.

##  1 Introduction

The intriguing landscape of particles dispersed in turbulent flows has
attracted tremendous attention over the past decades due to their extensive
existence in nature and industry (Pedley & Kessler, 1992; Moffet & Prather,
2009; Lundell et al., 2011; Mittal et al., 2020). The presence of particles
could introduce an extra scale (the particle size) to the turbulent flows,
which may alter the energy cascade of the turbulence and modulate the
turbulent flow. For neutrally buoyant particles as considered in the present
work, when the suspension is dilute and/or the particles are extremely small
(Voth & Soldati, 2017; Mathai et al., 2020; Brandt & Coletti, 2022), the
particle-laden system could be simplified using the point-particle model by
taking the dispersed particles as inertialess particles which could faithfully
follow the surrounding flows. Such typical situations applying the point-
particle model include plankton in the ocean (Stocker, 2012; Qiu et al., 2022)
and pollen species in the atmosphere (Sabban & van Hout, 2011). For particles
with higher volume fractions or non-negligible inertia, the interactions
between particle and fluid can be modelled as two-way coupling using point-
particle method toward the motivation of investigating the particle dynamics
and flow modulations (see the literature in the reviews (Voth & Soldati, 2017;
Brandt & Coletti, 2022)). However, given its simplicity, the point-particle
model could not precisely account for the feedback of particles on the
surrounding flows in general situations, such as the particle boundary layer
and the wake behind them (Jiang et al., 2022). To this end, the finite-size
particles, typically with diameters larger than the dissipation length scale
$\eta$ of the surrounding turbulent flows (Voth & Soldati, 2017; Mathai et
al., 2020; Brandt & Coletti, 2022), have been a vigorous field in recent
years, both in experiments (Qureshi et al., 2007, 2008; Fiabane et al., 2012;
Will et al., 2021; Will & Krug, 2021a, b; Obligado & Bourgoin, 2022) and
numerical simulations (Calzavarini et al., 2009; Picano et al., 2015; Yousefi
et al., 2020; Peng et al., 2020; Demou et al., 2022; Li et al., 2022; Assen et
al., 2022). In particular, by numerically implementing no-slip boundary
conditions on their surface, the finite-size particles are capable of
modelling the particle dynamics (Jiang et al., 2020, 2022) and the resulting
turbulence modulation (Wang et al., 2017b; Ardekani et al., 2017; Ardekani &
Brandt, 2019).

One of the subjects of particle-laden turbulence in the non-dilute regime is
the turbulence modulation caused by particles. Due to the no-slip boundary
conditions at their surface, the finite-size particles could modulate the
surrounding flow and the entire field such as through shedding vortex in the
wake region (Risso, 2018; Mathai et al., 2020) and enhancing the dissipation
rate around them (Jiang et al., 2022). For particle-laden turbulence with
neutrally buoyant particles, the system can be characterized by four
dimensionless numbers, namely, the flow Reynolds number $Re$, the size-ratio
of particle diameter to turbulence dissipation length scale $d_{v}/\eta$, the
aspect ratio of particles $\lambda$ and the volume-fraction of particles
$\phi$. Among others, the aspect ratio of particles has been found to affect
the clustering effects and particle dynamics, and thus altering the magnitude
of turbulence modulation. For instance, Ardekani & Brandt (2019) found that
the spherical particles could result in an overall drag enhancement, while the
non-spherical particles might result in drag enhancement or reduction
depending on their aspect ratios. Similar findings on the turbulence
modulation caused by non-spherical are also reported in their previous work
(Ardekani et al., 2017). In the view of turbulence, the aspect ratios of
particles are found to affect their effects on the turbulence fluctuation
intensity, thus altering the stress contributions (Ardekani & Brandt, 2019).

In our previous experimental work in turbulent Taylor-Couette (TC) flow (Wang
et al., 2022), it has been observed that the suspended spherical particles
could result in a larger drag increase than the non-spherical particles. By
proposing a qualitative analysis of the stress balance, the drag increases
caused by the suspended particles are explained. On the other hand, it is
observed that the clustering effects of particles are affected by the aspect
ratio of particles. The spherical particles show clustering near the walls,
whereas the non-spherical particles preferentially cluster in the bulk region.
Based on the previous numerical studies of finite-size particles (Ardekani et
al., 2017; Ardekani & Brandt, 2019), it is conjectured that the preferential
clustering of particles could be responsible for their different magnitudes of
turbulence modulation. However, limited by the experimental techniques, an in-
depth and more quantitative analysis is absent then on several aspects.
Specifically, (i) how and to what extent is the basic turbulent flow modulated
by the suspended particles with different aspect ratios? (ii) how are the
particle statistics affected by their aspect ratios? (iii) what is the
difference in the turbulence modulation between the near-wall clustering of
spherical particles and the bulk clustering of non-spherical particles?

We should note that, although similar phenomena have been reported for channel
flows, the difference of shear flow from the channel flow could give rise to
new physics and observations. For example, the findings on drag modulation in
Wang et al. (2022) indicate that, regardless of the particle aspect ratio,
only drag increments is possible for particle-laden Taylor Coutte flow.
However, numerical studies in channel flows have found drag reduction in the
cases of oblate particles (e.g.Ardekani & Brandt (2019)), leaving open the
question of whether this is a physical feature of channel flow or simply a
numerical artefact. To resolve this issue, we can perform simulations in
similar shear flows to those used in experiments on Taylor-Couette flow. In
fact, unlike the pressure-driven channel flow, which has a parabolic velocity
profile and yields the maximum velocity and zero velocity gradient at the
channel centre, the Couette flow shows a non-zero gradient at the domain
centre. The difference in velocity profile could affect the velocity
fluctuation and thereby the stress balance (drag) in the fluid phase.
Additionally, the particle statistics are affected by the velocity profile of
the flow as their dynamics are related to the velocity gradient tensor of the
flow (Voth & Soldati, 2017). Therefore, investigating the shear flows laden
with particles is of great interest.

To this end, this work aims to conduct fully particle-resolved direct
numerical simulations (PR-DNS) to provide a quantitative analysis of our
previous experimental results and answer the aforementioned questions.
Performing simulation in Taylor-Couette flow is computationally expensive
because it requires not only processing cylindrical geometries but non-uniform
gird spacing to resolve the boundary layers in the radial direction. Hence, we
employ a configuration of plane-Couette flow to carry out our investigation,
which is of Cartesian coordinates and can be efficiently implemented with a
scheme based on the Lattice Boltzmann Method. Using immersed boundary method
(IBM) (Peskin, 2002), we can resolve the motion of each particle and its
surrounding flow field, allowing us to uncover the physics of both particle
dynamics and the resulting turbulence modulation. The rest of the paper is
organized as follows: Section 2 describes the flow configuration and numerical
schemes, in Section 3 we discuss the results, including the turbulent
arguments and the particle statistics, and Section 4 gives final remarks on
this work.

##  2 Numerical Methodology

###  2.1 Configurations of flow and particles

Figure 1: Sketch of the flow configuration. Two walls are moving at constant
velocities in opposite directions. For particle-laden cases, particles with
three kinds of aspect ratios are randomly suspended at the beginning and then
move freely in the domain box.

We conduct simulations with spheroidal particles suspended in a turbulent
plane-Couette flow. The flow configurations are depicted in figure 1. Two
walls are moving in the opposite direction with a constant velocity $U$. No-
slip boundary conditions are imposed on the walls, while the flow in the
streamwise and spanwise directions is periodic. The turbulent flow is governed
by the conservation equations of momentum and mass, which read as

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

where $\boldsymbol{u}$ and $p$ are the velocity vector and hydrodynamic
pressure field of the flow, $\rho$ and $\nu$ are the density and viscosity of
the fluid. Since the suspensions are composed of finite-size particles and in
the non-dilute regime, the particle can modulate the surrounding turbulence by
exerting feedback force $\boldsymbol{f}_{p}$ on the flow. The turbulence
intensity can be quantified by the bulk Reynolds number (hereafter, the
Reynolds number),

$Re=Uh/\nu,$ (3)

where $\nu$ is the kinematic viscosity of the fluid and $h$ is the half-width
of the domain in the wall-normal direction. Alternatively, the turbulence
intensity can also be related to the shear Reynolds number,

$Re_{\tau}=u_{\tau}h/\nu,$ (4)

where $u_{\tau}=\sqrt{\tau_{w}/\rho}$ is the friction velocity with $\tau_{w}$
being the shear stress at walls. The values of $Re$ and the corresponding
$Re_{\tau}$ are presented in Table 1. In this work, the values of $Re$ are set
based on that in our previous experimental work. We note that in TC flow (Wang
et al., 2022) the Reynolds number is defined as $Re_{e}=(u_{i}-u_{o})d/\nu$,
where $u_{i}$ ($u_{o}$) is the velocity of the inner (outer) wall and $d$ is
the gap width between walls. Hence, the values of $Re$ in the present work are
one-fourth of the corresponding values of $Re_{e}$ reported in Wang et al.
(2022) (see in Table 1).

$Re(=Re_{e}/4)$ | $Re_{e}$ | $Nx\times Ny\times Nz$ | $d_{v}/Ly$ | $\lambda$ | $\phi$ | $Re_{\tau}$ | $d_{v}/\eta$ | $Re^{slip}_{p}$ | $Re_{p}^{\lambda}$ | $\Delta/\delta_{\nu}$ | $\Delta/\eta$ | $u_{\tau}T/h$ | $\delta_{\tau_{w}}/\overline{\tau}_{w}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---
1600 | 6400 | $384\times 192\times 192$ | 1/10 | - | 0% | 99.8 | - | - | - | 1.04 | 0.658 | 78.0 | 0.7%
1/3 | 2% | 99.3 | 12.64 | 11.06 | 64.0 | 1.03 | 0.658 | 76.9 | 1.5%
1 | 2% | 101.3 | 12.74 | 6.38 | 64.0 | 1.06 | 0.664 | 145.1 | 0.4%
3 | 2% | 98.9 | 12.58 | 1.92 | 64.0 | 1.03 | 0.655 | 83.0 | 0.5%
1/3 | 6% | 100.4 | 13.00 | 2.15 | 64.0 | 1.05 | 0.677 | 80.4 | 3.0%
1 | 6% | 105.2 | 13.02 | 6.36 | 64.0 | 1.10 | 0.679 | 89.7 | 0.4%
3 | 6% | 100.5 | 12.78 | 4.94 | 64.0 | 1.05 | 0.666 | 81.8 | 0.4%
3200 | 12800 | $512\times 256\times 256$ | 1/10 | - | 0% | 179.2 | - | - | - | 1.40 | 0.785 | 56.9 | 0.3%
1/3 | 2% | 179.7 | 20.26 | 12.39 | 85.3 | 1.40 | 0.791 | 49.1 | 0.3%
1 | 2% | 183.5 | 20.44 | 3.76 | 85.3 | 1.43 | 0.798 | 89.6 | 0.4%
3 | 2% | 181.1 | 20.29 | 7.80 | 85.3 | 1.42 | 0.792 | 50.0 | 0.4%

Table 1: Simulation parameters. $Re$ is the bulk Reynolds number defined by
the wall velocity $U$ and half-width of the gap $h$, which is therefore one-
fourth of the corresponding values of $Re_{e}$ reported in Wang et al. (2022)
defined by the difference of wall velocity, $(u_{i}-u_{o})=2U$, and the gap
width, $d=2h$. $Nx\times Ny\times Nz$ is the number of grids in each
direction. $d_{v}/Ly$ and $d_{v}/\eta$ are the ratios of the equivalent-volume
diameter of particles to the width of the domain in the wall-normal direction
and to the dissipation length scale, respectively. $\lambda$ and $\phi$ are
the aspect ratio and the volume fraction of particles. $Re_{\tau}$ is the
friction Reynolds number defined by equation 4. $Re^{slip}_{p}$ is the
particle Reynolds number based on the global mean slip velocity and
$Re_{p}^{\lambda}=Ud_{v}^{2}/h\nu^{2}$ is the particle Reynolds number based
on local shear. $\Delta/\delta_{\nu}$ and $\Delta/\eta$ are the grid spacing
normalized by the viscous length scale $\delta_{\nu}$ and the dissipation
length scale $\eta$. $u_{\tau}T/h$ indicates the running times in terms of
turnover periods for eddies of size $h$ and velocity $u_{\tau}$.
$\delta_{\tau_{w}}/\overline{\tau}_{w}$ gives the statistic errors for each
case, where $\delta_{\tau_{w}}=\tau_{w}|_{y=2h}-\tau_{w}|_{y=o}$ is the
difference of wall stress, and
$\overline{\tau}_{w}=(\tau_{w}|_{y=2h}+\tau_{w}|_{y=o})/2$ is the mean wall
stress.

Neutrally buoyant spheroids are studied with different parameters, including
the equivalent-volume diameter $d_{v}$, the aspect ratio $\lambda$, and the
volume fraction $\phi$, chosen based on our previous experimental work (Wang
et al., 2022). Here for a spheroidal particle, we have
$d_{v}=(h_{p}d_{p}^{2})^{1/3}$ and $\lambda=h_{p}/d_{p}$, where $d_{p}$ is the
length of the symmetric axis of the particle, and $h_{p}$ is the length
perpendicular to it (see figure 1). The parameters of particles are given in
Table 1. The particle motion is governed by the Newton-Euler equations as
follows

$m_{p}\frac{\mathrm{d}\boldsymbol{v}}{\mathrm{d}t}=\boldsymbol{F}+\boldsymbol{F}_{c},$
(5)
$\frac{\mathrm{d}\boldsymbol{I}_{p}\boldsymbol{\Omega}}{\mathrm{d}t}=\boldsymbol{T}+\boldsymbol{T}_{c},$
(6)

where $\boldsymbol{v}(t)=\mathrm{d}\boldsymbol{r}/\mathrm{d}t$ and
$\boldsymbol{\Omega}(t)$ are the particle velocity and angular velocity
vectors of a particle at position $\boldsymbol{r}(t)$ with mass
$m_{p}=\rho_{p}V_{p}$ ($\rho_{p}$ is the particle density and $V_{p}$ the
volume) and $\boldsymbol{I}_{p}$ the moment of inertia tensor. In the right-
hand side (r.h.s.) of the above equations, $\boldsymbol{F}$ is the force
exerted on particles by the surrounding fluid through hydrodynamic, while
$\boldsymbol{F}_{c}$ is the collisions force accounting for the particles-
particle collisions and particles-wall collisions. $\boldsymbol{T}$ and
$\boldsymbol{T}_{c}$ are the torque defined in similar means.

At $Re=1600$, two values of $\phi$ are chosen to study its effects on
turbulence modulation. Whereas at $Re=3200$, only the case of small $\phi$
($\phi=2\%$) is conducted due to the costly computational expense at the high
volume-fraction case, which nonetheless could uncover the roles of turbulent
intensity. For each $Re$ case, we first run the cases of single-phase flow
(i.e., $\phi=0\%$) to generate a fully developed plane-Couette turbulence.
Then the particle-laden case is conducted by randomly dispersing the particle
in the entire domain with an initial velocity equal to zero. Statistical
equilibrium is ensured by reaching time-averaged statistical errors in wall
stress typically less than $1\%$ (see Table 1).

###  2.2 Numerical Schemes

In this section, we briefly describe the numerical schemes used to simulate
the particle suspension. The turbulence is numerically solved using an open-
source code with the lattice Boltzmann method (LBM), the $ch4-project$
(Calzavarini, 2019), which has been extensively validated in particle-laden
turbulence both for point-like particles (Mathai et al., 2016; Calzavarini et
al., 2020; Jiang et al., 2021) and finite-size particles (Jiang et al., 2020,
2022). The code was validated by comparing the translational dynamics of
spheres at $Re_{\lambda}=32$ with a previous reference study (Homann & Bec,
2010). The details and validations of the code can be found in our previous
work (Jiang et al., 2022).

On the fluid-phase side: The domain size is $Lx\times Ly\times Lz$ = $4h\times
2h\times 2h$, where $Lx,Ly,Lz$ are the width of the domain in the streamwise,
wall-normal and spanwise direction, respectively. Here we note that the domain
size is chosen based on producing fully-developed turbulence, which for the
current $Re$ ranges is $Lx/Ly\simeq 2$ (Owolabi & Lin, 2018). To double-check
if the domain size is long enough and the artefacts of periodic boundary
conditions affect the results, we run an addition case with $Lx\times Ly\times
Lz$ = $8h\times 2h\times 2h$. It is found that the drag coefficient ( see the
$blue$ square dot in figure 2a) and the particle statistics (not shown here)
are comparable to the cases reported in Table 1. Therefore, the domain size in
Table 1 is long enough to yield conclusive results. The computational domain
is meshed by grids distributed uniformly in each direction, and the numbers of
grids, $Nx\times Ny\times Nz$, are given in Table 1. Considering the
requirement of the numerical resolution, the maximum grid spacing is validated
to be $\Delta\leq 1.1\delta_{\nu}$ and $\Delta\leq 0.68\eta$ at $Re=1600$ (see
Table 1), where $\delta_{\nu}=\nu/u_{\tau}$ is the viscous length scale and
$\eta=(\nu^{3}/\epsilon)^{1/4}$ is the dissipation length scale with
$\epsilon$ being the global dissipation rate of the flow. In consequence,
about 5 grids are embedded in the linear viscous layer (as can be seen in
figure 5). Hence, the grid spacing adopted in all runs is sufficient to
resolve the turbulence fluctuation and flow structures, both in the boundary
layers and the bulk.

On the particle-phase side: the particle-fluid interaction is fully resolved
by employing the immersed boundary method (IBM) (Peskin, 2002). The IBM
enforces the no-penetration and no-slip boundary conditions at the
fluid–particle interface by means of a localized feedback force,
$\boldsymbol{f}_{p}$, added to equation 1. Such $\boldsymbol{f}_{p}$ term is
also denoted as two-way coupling. The particle surfaces are captured by
uniformly distributed Lagrangian nodes. Adequate grids are required to resolve
the flow around particles. The test toward this has been carried out in our
previous work (Jiang et al., 2022). It is found that at least 16 grids per
particle diameter are needed to resolve the particle boundary layers. In this
work, we surpass this criterion to attain convincing results. Specifically,
for the cases of $Re=1600$ and $Re=3200$, 19 and 26 grids per particle
diameter are respectively used to resolve the particle boundary layers and
their dynamic. 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) and Wang et
al. (2021). The processes of particle-particle collisions and particle-wall
collisions 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), where
$\boldsymbol{F}_{c}$ and $\boldsymbol{T}_{c}$ can be computed. The details of
the calculations of the hydrodynamic force and torque can be found in our
previous work (Jiang et al., 2022). In our simulations with the immersed
boundary method, the boundary of the finite-size particle is represented by
Lagrangian nodes and the particle-fluid interactions are solved by the delta
function as suggested by Peskin. Furthermore, we seed Lagrangian points inside
the finite-size particle which are fixed in the particle frame. These
Lagrangian points are used to measure the momentum time derivative of fluid
inside the finite-size particle to obtain the correct driving force. And we
use the tri-linear scheme to interpolate the velocities of the Lagrangian
points inside the finite-size particle. The Newton-Euler equations that govern
the particle dynamics are numerically integrated with a second-order
Adams–Bashforth time-stepping scheme.

Since the size-ratio $d_{p}/\eta$ is fixed, the numbers of particles $Np$ are
determined by their volume fractions, yielding $Np=76$ and $229$ for the cases
of $\phi=2\%$ and $\phi=6\%$, respectively. Therefore, the statistical
convergence of the numerical wall stress and the particle statistics are
ensured by running the simulation long enough to get sufficient data. The
running time for each case is given in Table 1. For additional information,
the time series of the fluid and particle velocity fluctuations are shown in
Appendix A.

##  3 Results and Discussion

###  3.1 Turbulence modulation: global drag, velocity field, and dissipation
rate

Figure 2: Normalized wall stress for particle-laden cases. The experimental
data (empty circles) is adapted from the TC experiments (Wang et al., 2022).
The $blue$ square in (a) denotes the result of the validation case of extended
domain size in the streamwise direction ($Lx\times Ly\times Lz$ = $8h\times
2h\times 2h$). The values of $Re$ and $\phi$ are shown in each panel.

We start by comparing the global turbulence modulation of numerical results to
that of experimental measurements (Wang et al., 2022). The global transport
quantity of plane-Couette turbulence is the shear stress at the walls

$\tau_{w,\phi}=\rho\nu\frac{\mathrm{d}U_{f}}{\mathrm{d}y}|_{y=0,2h},$ (7)

where $\frac{\mathrm{d}U_{f}}{\mathrm{d}y}|_{y=0,2h}$ is the gradient of mean
stream-wise velocity at the walls. We normalize the wall stress of particle-
laden flow by that of single-phase flow, i.e.,
$\tau_{w,\phi}/\tau_{w,\phi=0}$, and show it in figure 2, where the results
from previous experiments in Taylor-Couette flow(Wang et al., 2022) are also
depicted for comparison. The discrepancies between numerics and experiments
could be attributed to the different flow configurations. Although both flow
configurations are shear-induced flows, TC flow suffers from more complex
conditions, such as the cylindrical geometries causing asymmetric velocity
profiles in the radial direction, and the strong confinement due to the end-
plates which cause secondary flow (Grossmann et al., 2016). It is therefore
not expected to observe exact agreements between the results of numeric and
experiments. Moreover, as the particle dynamics are affected by their aspect
ratios (Voth & Soldati, 2017), the spherical and non-spherical particles could
respond in different means as the flow configuration changes, therefore
resulting in different responses in drag modulation. This could possibly
account for that in drag modulation (figure 2), compared to non-spherical
particles, the spherical particles show better agreements in TC flow and
plane-Coutte simulations. Nevertheless, concerning that we are mainly
interested in the effects of particle aspect ratio on the turbulence
modulation, the conclusion can be drawn that the spherical particle could
cause larger drag modulations than the non-spherical particles, for different
$Re$ and $\phi$ with the numerical results in figure 2. Similar findings are
also reported in our previous experimental TC flows, where the difference in
magnitude of drag modulation between spherical and non-spherical cases is not
obvious at the present $\phi$ ($\Delta_{c_{f}/c_{f,\phi=0}}\sim 1\%$, see
figure 2) but become larger at a higher value of $\phi$
($\Delta_{c_{f}/c_{f,\phi=0}}\sim 5\%$ at $\phi=10\%$).

Figure 3: Instantaneous snapshots of (a-c) stream-wise velocity indicated by
the colour bar. The overlaid arrow vectors depict flow velocity within the
plane. (d-f) dissipation rate. Data is obtained from the cases of spherical
particles ($\lambda=1$) at $Re=1600$. Columns from left to right are
$\phi=0\%$, $\phi=2\%$ and $\phi=6\%$, respectively.

For fixed $Re$ and $\lambda$ (see figures 2(a, b)), the magnitudes of drag
modulation increase with increasing $\phi$, which again agrees with the
results in our previous experimental measurements (Wang et al., 2022). To
demonstrate this $\phi$-dependence of drag modulation, in figure 3 we present
the velocity field and the corresponding fields of dissipation rate at the
centre plane in stream-wise direction (i.e. $x=Lx/2$). For clarity, only the
case of the spherical particles at $Re=1600$ is shown. Compared to the case of
single-phase (figure 3a), the presence of particles could introduce strong
disturbance both in the bulk and the boundary layers, see the emergence of the
higher value of dissipation rate coloured in red around the particles in
figures 3(b,c). On the other hand, figure 4 shows similar information as
figure 3 but for different $\lambda$ at $Re=1600$ and $\phi=6\%$. It can be
seen that the particle could modify the surrounding fields and particularly
enhance the dissipation rate near their surface, which again could be related
to the particle boundary layers (Jiang et al., 2022) and accounts for the drag
increases reported in figure 2. This enhancement of turbulence dissipation of
particles has been found in other flows, e.g. homogenous and/or isotropic
turbulence (Cisse et al., 2013; de Motta et al., 2016; Jiang et al., 2022),
and decaying isotropic turbulence (Lucci et al., 2010). Due to the no-slip
condition enforced at the particle surface, a local force arises on the fluid
surrounding the particle, which locally increases the velocity gradients close
to the particle surface and thus increases the local strain rate as well as
the dissipation rate (Lucci et al., 2010). As the distance to the particle
surface increases, this impact of the no-slip condition fades, and so does the
dissipation rate. This accounts for the observations in figures 3, 4, where
the bursts of high dissipation rate (coloured in $red$) only occur near the
particles while the rear fluid almost remains undisturbed (coloured in
$blue$). Indeed, the particles are found to enhance the dissipation rate
around them to a distance of about one radius to their surfaces (Cisse et al.,
2013; de Motta et al., 2016) for spherical particles, while for non-spherical
particles the volume of disturbed fluid is smaller in every dimension (Jiang
et al., 2022). Additionally, the enhancement of the dissipation rate around
particles distributes in a non-uniform manner, namely, the dissipation rate in
front of the particle is higher than the rear of the particles (Cisse et al.,
2013; de Motta et al., 2016). Although this can not be quantified in the
present work, one can still be hinted from figures 3(e,f) and figures 4(d-f),
where the dissipation rate around particles is asymmetric. With increasing
volume fraction of particles, the dissipation rate becomes pronounced, as
could be seen in figures 3,4. This is due to the increase of particle surface
area that affects the surrounding fluid and thus the two-way coupling force
(Lucci et al., 2010). Moreover, a closer look would find that the spherical
particles cause stronger disturbances to the boundary layers, while the non-
spherical particles mainly alter the flow in the bulk.

Figure 4: Instantaneous snapshots of (a-c) stream-wise velocity indicated by
the colour bar. The overlaid arrow vectors depict flow velocity within the
plane. and (d-f) dissipation rate. Data is obtained at $Re=1600$ and
$\phi=6\%$. Columns from left to right are for $\lambda=1/3$, $\lambda=1$ and
$\lambda=3$, respectively.

We depict the profiles in the wall-normal direction of stream-wise velocity
and dissipation rate in figure 5 and figure 6, respectively. In figure 5, the
velocity profile is shown in the wall units, where $y^{+}=y/\delta_{\nu}$ is
the distance from the lower wall ($y=0$) in units of the viscous length scale
$\delta_{\nu}$, and $u^{+}=(u-u|_{y=0})/u_{\tau}$ is the velocity difference
from the lower wall normalized by the friction velocity $u_{\tau}$. For
comparison, we also depict the viscous sublayer $u^{+}=y^{+}$ and the
logarithmic law of flat turbulent boundary layer flow, i.e.
$u^{+}=1/\kappa~{}$ln$~{}y^{+}+B$ with the typical values of $\kappa=0.40$ and
$B=5.2$ as suggested by Prandtl and von-Kármán (Pope, 2000). One can see that,
regardless of the values of $\phi$ and $\lambda$, the profile of $u^{+}=y^{+}$
is well retained with sufficient grids inside the viscous sublayer
($y^{+}<5$), which is barely affected by the suspended particles due to their
finite sizes (see the vertical dotted lines in figure 5). Outside the viscous
sublayer, the velocity profiles generally follow the logarithmic law in all
cases. The effects of particle aspect ratio become pronounced with increasing
$\phi$ (see figure 5b). Compared to single-phase flow, the cases of particle-
laden flow result in the velocity profile consisting of an extended
logarithmic layer, which is similar to the observation in channel flow
(Eshghinejadfard et al., 2018) and yields in the reduced velocity gradient in
the near-wall region. This seems to contradict the drag increases reported in
figure 2. However, as shown later, for particle-laden flow the dampening in
velocity does not imply a global drag reduction since there are extra stress
contributions from the particle phase. Moreover, in terms of the effect of
particle aspect ratio, it is found that at fixed $Re$ and $\phi$ the spherical
particles cause greater velocity reduction than the non-spherical particles in
the logarithmic layer (e.g., see figure 5b). Specifically, in figure 5b, the
value of $B$ of spherical particles ($B\simeq 2.3$) is smaller than that of
non-spherical particles ($B\simeq 3.5$) while the values of $\kappa$ are close
to $0.40$ in all particle-laden cases. This is consistent with the drag
increase found in figure 2 since the decrease in $B$ indicates a drag
enhancement and the decrease in $\kappa$ does the opposite (Ardekani et al.,
2017).

Figure 5: Velocity profiles depicted in the wall units. The vertical dotted
and dashed lines in each panel indicate, respectively, the location of the
equivalent-volume radius of particles ($r_{v}=d_{v}/2$) and the location of
the central plane in the wall-normal direction ($y=h$). Figure 6: Profiles of
dissipation rate. The insets show the zoom-in view of the dissipation rate in
the near-wall region.

As for the dissipation rate (figure 6), its magnitudes agree well with the
magnitudes of drag increases found in cases of different particle aspect
ratios. We note that, though the non-spherical particles could alter the bulk
flow stronger than the spherical ones in the bulk region (see figure 4), the
dissipation rates in the bulk are minor contributions to the total field,
resulting in negligible changes of the profiles of dissipation rate in that
region (figure 6). In contrast, spherical particles efficiently increase the
dissipation rate more than non-spherical particles near the walls (see insets
in figure 6). This is in line with the magnitudes of global drag modulation of
particles with different $\lambda$ because, as the wall stress does, the
dissipation rate in the near-wall region is positively associated with the
velocity gradient at the walls.

###  3.2 Turbulence modulation: stress analysis and turbulent energy spectrum

In the previous work (Wang et al., 2022), we decompose the total stress (i.e.,
the shear stress at walls $\tau_{w,\phi}$) into four terms, namely, the
viscous stress of the fluid phase ($\tau_{\nu}$), the turbulent stress of
fluid and particle phase ($\tau_{Tf}$ and $\tau_{Tp}$, respectively), and the
particle-induced stress ($\tau_{p}$), which can be written as

$\tau_{w,\phi}=\tau_{\nu}+\tau_{Tf}+\tau_{Tp}+\tau_{p}.$ (8)

where $\tau_{w,\phi}$ is defined by equation 7. For each term in the r.h.s,
one can explicitly express it, as suggested by (Batchelor, 1970; Zhang &
Prosperetti, 2010; Picano et al., 2015; Wang et al., 2017a, 2022), in the form
of turbulence quantities as,

$\tau_{\nu}=\rho\nu(1-\phi)\frac{\mathrm{d}U_{f}}{\mathrm{d}y}$ (9)
$\tau_{Tf}=-\rho(1-\phi)\langle u^{\prime}_{f}v^{\prime}_{f}\rangle,$ (10)
$\tau_{Tp}=-\rho\phi\langle u^{\prime}_{p}v^{\prime}_{p}\rangle,$ (11)
$\tau_{p}=\phi\langle\sigma^{p}_{xy}\rangle,$ (12)

where $u^{\prime}$ and $v^{\prime}$ are the velocity fluctuation in the
stream-wise and wall-normal direction with the subscripts $f$ and $p$ denoting
fluid and particle phase, respectively, $\sigma^{p}_{xy}$ the general stress
in the particle-phase, projected in the stream-wise plane and pointing in the
stream-wise direction. $\langle...\rangle$ means ensemble averaged in time and
over the stream-wise plane. For the case of single-phase flow, the latter two
terms are zero.

By means of numerics, we can look into the contributions of each term in r.h.s
of equation 8 to the total stress. The particle-induced stress $\tau_{p}$ is
calculated indirectly (i.e.
$\tau_{p}=\tau_{w}-\tau_{\nu}-\tau_{Tf}-\tau_{Tp}$) as the explicit expression
of $\sigma^{p}_{xy}$ remains unknow currently. In figure 7 we plot the profile
of stress balance in the wall-normal direction for each case, where we
normalize each term of stress by the total stress of that case, i.e.,
$\tau/\tau_{w,\phi}$. We first demonstrate how the particle aspect ratio,
$\lambda$, affects the contributions of each stress at given $Re$ and $\phi$
(see each row of figure 7). Taking the case of $Re=1600$ and $\phi=2\%$ as an
example (first row of figures 7), one can see that the particle aspect ratio
has slight effects on $\tau_{\nu}$ and $\tau_{Tf}$. Similarly, the turbulent
stress of the particle phase, $\tau_{Tp}$, only weakly depends on the
$\lambda$ while having a positive dependence on the $\phi$ (see the second row
of figure 7). However, the particle-induced stress, $\tau_{p}$, is affected by
the particle aspect ratios to a noticeable extent, which becomes even more
pronounced at high $\phi$ (see the second row of figure 7). Since $\tau_{p}$
is the combined results of hydrodynamic interaction and particle collisions
(Picano et al., 2015), its value could be related to both the particle
dynamics and the turbulence field. As the numerical errors from $\tau_{p}$
cannot be excluded now, we then take the case of high $\phi$ (second row of
figures 7), where the numerical errors are of less significance, as an
instance to illustrate the effects of $\lambda$ on $\tau_{p}$. First, it is
noticeable that in all cases of $\lambda$, $\tau_{p}$ reaches its peak near
the walls. Near the walls, the spherical particles cause more pronounced peaks
than the non-spherical particles in $\tau_{p}$, which validates our
conjectures in previous work (Wang et al., 2022) and is in line with the
findings in other configurations of particle-laden flow (Picano et al., 2015;
Ardekani et al., 2017). As shown by the particle statistics in the later
section, for spherical cases the occurrence of near-wall peaking in $\tau_{p}$
could originate from their near-wall preferential clustering, which could
cause the particle collisions to happen in easier and more intense ways.

Figure 7: Profiles of stress budget. Each term in r.h.s. of equation 8 is
normalized by the total shearing stress of corresponding particle-laden cases
(i.e. $\tau_{w,\phi}$ in equation 8). The two horizontal dashed lines with the
values of 1 and 0 are shown for reference.

It is noticeable that, for the case of oblate ($\lambda=1/3$, see figure 7d),
$\tau_{p}$ becomes negative in the bulk region. As the results are ensured to
be convergence by running the simulation long enough in time, the fact that it
is negative could be attributed to the hydrodynamic interaction caused by
particles. This indicates that $\tau_{p}$ could be a stress sink under
specific conditions. Similar overall drag reductions caused by particles have
also been reported in simulation results in channel flow (Ardekani & Brandt,
2019). However, we should note that the previous experimental results (Wang et
al., 2022) report that only drag enhancement is observed in particle-laden
flow under similar parameter regimes, even for the cases of oblate particles.
This is in contrast to the current simulation results, where the particle-
induced stress is found to be negative (figure 7d) and overall a drag (stress)
reduction is observed (figure 8). In this sense, the particle-induced stress
being negative in the bulk could be a numerical artefact, though there are
differences in TC and plance-Couette flows as we mentioned before. Future
study is needed to uncover the role of $\tau_{p}$. Moreover, the profiles in
figure 7 are not exactly symmetric, especially in the particle-induced stress
($\tau_{p}$). This could be due to the limitations of statistics as the
particle number determined by its volume fraction are relatively small, which
could cause inhomogeneity in the bounded wall-normal direction (i.e.
$y-$direction). Although extending the simulation time steps might improve the
data quality, as we did for the spherical cases (see the values of
$u_{\tau}T/h$ in Table 2), the asymmetries are still present in the spherical
cases. Nonetheless, the asymmetric amplitude is typically less than 5% for all
simulation cases and thereby the main conclusion on the stress analysis can be
convinced. On the other hand, as $Re$ increases (see the third row of figure
7), the profiles of stress coming from the particle phase ($\tau_{Tp}$ and
$\tau_{p}$) become almost flat. In these cases, the total stress is dominated
by the viscous stress in the boundary layers and by the fluid-phase turbulence
stress in the bulk, which is actually reminiscent of the case of single-phase
flow.

Figure 8: Contributions of each individual stress to the total shearing
stress. The total shearing stress of the corresponding cases of single-phase
(i.e., $\tau_{w,\phi=0}$ at given $Re$) is used for normalization.

The contribution ratio of each stress to the total one can be found in figure
8, where we normalize each stress by the total shear stress at the walls of
the cases of single-phase, i.e., $\tau/\tau_{w,\phi=0}$. For small $\phi$ (see
figures 8(a,c)), it can be seen that the total stress is almost dominated by
the viscous stress and the turbulent stress, while the presence of particles
negligibly affects the stress of the flow. The changes in $Re$ mainly
redistribute the contributions from the $\tau_{\nu}$ and $\tau_{Tf}$ to the
total stress. With increasing $\phi$ (see figure 8b), the role of particles
becomes pronounced as indicated by the greater contributions of $\tau_{Tp}$
and $\tau_{p}$. Indeed, one can see that the stress coming from the particle
phase could principally be responsible for the drag increases found in the
cases of particle-laden flow, even though the turbulent stress of the fluid
phase might be dampened a bit by the suspended particles at higher $\phi$.
Nevertheless, we note that, compared to the fluid phase, the contributions
from particle-phase stress ($\tau_{Tp}$ and $\tau_{p}$) are noticeable but
still minor effects due to the rather small $\phi$ here.

Figure 9: Turbulent energy spectra (a-c) and the compensated plots ((d-f) are
compensated with $k^{-5/3}$, (g-i) are compensated with $k^{-3}$). The two
dashed lines shown in the first row are for the guidance of scaling
behaviours. The vertical dashed line in each panel indicates the scale of the
equivalent-volume diameter of the particles $d_{v}$.

Besides the stress analysis, we also examine the one-dimensional spectra of
turbulent velocity fluctuations of the fluid phase for particle-laden cases as
shown in figure 9. For all cases of $Re$ and $\phi$ shown in figures 9(a-c),
one can see that the cases of particle-laden flow coincide with that of
single-phase flow at large scales while deviating above it at small scales. As
indicated by the vertical dashed line in each panel, the suspended particles
mainly enhance the energy fluctuation at scales smaller than $d_{v}$. In
contrast, the particle aspect ratios only mildly affect the spectrum. For
comparison purposes, we show the compensated plot of the power spectrum with
$k^{-5/3}$ (figures 9(d-f)) and $k^{-3}$ (figures 9(g-i)). Both the spectrum
of single-phase flow and particle-laden flow follow the well-known scaling
$k^{-5/3}$ in an inertial range as proposed by Kolmogorov (Pope, 2000). It
should be noted that the emergence of $k^{-5/3}$ scaling holds over less than
one decade, which is mainly due to the moderate turbulence intensity in the
present work. At scales smaller than $d_{v}$, however, the cases of particle-
laden flow show good agreements with the scaling $k^{-3}$. Such steeper
scaling has been robustly observed in bubbly flow, through experiments
(Mercado et al., 2010; Riboux et al., 2010; Mendez-Diaz et al., 2013; Prakash
et al., 2016; Alméras et al., 2017; Zhang et al., 2020; Dung et al., 2022) and
numerics (Pandey et al., 2020; Innocenti et al., 2021; Hidman et al., 2022;
Pandey et al., 2022), at scales approximately below the bubble diameter and
above the dissipation length scale (Hidman et al., 2022). The emergence of
$k^{-3}$ in bubbly flow could be attributed to the wake behind bubbles(Alméras
et al., 2017; Risso, 2018; Wang et al., 2019; Pandey et al., 2020) \- also
referred to as bubble-induced agitations or pseudo-turbulence - where the
velocity fluctuation produced by bubbles is directly dissipated by viscosity
(Lance & Bataille, 1991; Risso, 2018). Indeed, bubbles need to be finite-size
(i.e. large enough) to generate wakes and thereby the emergence of $k^{-3}$
scaling, typically with a bubble Reynolds number $10\leq Re_{bub}\leq 1000$
(Risso, 2018; Pandey et al., 2020). In contrast, the $k^{-3}$ scaling could
not be observed in simulations of point-like particles since the wake is
absent (Mazzitelli & Lohse, 2009).

However, this seems not to be the situation here for the particle-laden flow
even though the particles are finite-size. Unlike the cases of rising bubbles
where the relatively large slip velocity is the major contribution to the
pseudo-turbulence (Risso, 2018), the particles used in this work could only
yield in small slip velocity due to their neutral buoyancy. Here, we measure
the particle Reynolds number based on the mean slip velocity,
$Re^{slip}_{p}=(\langle u_{x}\rangle-\langle v_{x}\rangle)d_{v}/\nu$. As shown
in Table 1, the $Re_{p}^{slip}$ over the whole domain for each case is of the
order $O(1)$ to $O(10)$, which is of marginal value to observe the scaling of
$k^{-3}$ for the bubbly flow (Risso, 2018; Pandey et al., 2020). However, the
particle Reynolds number could be sensitive to the position as the flow is
inhomogeneous. To check this, we look into the profiles of the averaged
stream-wise velocity of both phases which are shown in figure 10. As it can be
seen in figure 10(a,b), the velocity of the particles differs mainly in the
near-wall region from the main flow of the fluid. Correspondingly, we measure
the profiles of the particle Reynolds number based on the horizontally
averaged slip velocity, $Re^{slip}_{p}(y)$, as shown in figure 10c. It can be
seen that when close to the walls the local particles Reynolds number is high
due to the large slip velocity in this region. Whereas in the bulk region, the
$Re^{slip}_{p}$ is rather small and could be due to their neutral buoyancy.
This being said the emergence of $k^{-3}$ scaling in the particle-laden flow
here could be originated from a different mechanism than that of bubbly flow,
which is beyond the scope of the current work. More studies are needed to
further investigate the origin of the ‘$-3$’ scaling observed in particle-
laden flow.

Figure 10: Profiles of the normalized stream-wise velocity of (a) fluid phase
$\langle u_{x}\rangle$ and (b) particle phase $\langle v_{x}\rangle$. (c)
Profile of the particle Reynolds number based on the horizontally averaged
slip velocity.

###  3.3 Particle statistics: distributions and orientations

The particle statistics, as the origins of turbulence modulation, are crucial
for understanding the underlying physics. In the section on stress analysis,
we find that the spherical particles cause stronger particle-induced stress
than the non-spherical ones in the near-wall region. To understand this, we
first examine how the particle distributions are affected by the aspect ratio.
To give a straightforward comparison to previous experimental results (Wang et
al., 2022), we note that here the particle distribution is presented by
showing its probability density function (PDF) with respect to its coordinates
in wall-normal direction (i.e., $PDF(y)\sim y$), instead of evaluating the
local volume fraction by computing the percentage of the grids inside the
particle boundary (see e.g. Wang et al. (2018)). It can be seen that, as shown
in figure 11, the aspect ratio of particles determines the particle
distribution regardless of $Re$ and $\phi$. The spherical particles (figure
11b) strongly cluster just near the walls, followed by ebbs next to the peaks.
On the contrary, the probability of non-spherical particles appearing in the
near-wall region is almost zero (figures 11(a, c)), while some above-average
peaks occur in the bulk region, suggesting mild but noticeable bulk-cluster
effects. Nevertheless, the particles are generally uniformly distributed in
most of the bulk. This could be attributed to the gentle turbulence and the
less intensive mixing caused by the flows in this region, which would result
in a weaker finite-size effect for the particles in the bulk region and
thereby weakened clustering. These clustering effects agree with our
experimental findings in Wang et al. (2022), and similar results of particle
preferential clustering have also been reported in other flows (Picano et al.,
2015; Ardekani et al., 2017; Ardekani & Brandt, 2019). On the other hand, it
is found that, in figure 11, $Re$ has minor effects on the particle
distributions. This could be due to the relatively strong finite-size effects
of particles, which allows the particles to escape from the flow structures.
It should be noted that the particle distribution could depend on $Re$ if the
turbulence intensity changes, for example, Wang et al. (2018) found that the
particles could be trapped inside the vortex and show different profiles of
distributions at smaller $Re$.

Figure 11: Particle distributions in the wall-normal direction. Panels from
left to right: $\lambda=1/3$, $\lambda=1$ and $\lambda=3$.

In addition to the particle distributions, its orientations (figure 12) are
analyzed. As expected, see figure 12b, the orientation profiles of spherical
particles show no preferential alignment to any coordinate axis since the
symmetric axis of spherical particles can be in an arbitrary direction.
However, the situation is different for non-spherical particles, particularly
in the near-wall regions. For oblate cases (figure 12a), the symmetric axes of
particles show a stronger alignment to the $y$-direction (i.e. wall-normal
direction) as the particle location approaches the walls. While for prolate
cases (figure 12c), the particles align their symmetric axes to the
$x$-direction (i.e. stream-wise direction) when close to the walls. As in the
bulk region, the particles with different aspect ratios still show
preferential alignment with the axis as in the near-wall region, i.e., the
symmetric axes of oblate (prolate) particles align with the $y$-direction
($x$-direction), while the angles between the symmetric axes of particles and
the coordinate axis are almost independent of the particle position. In
particular, we zoom in on the near-wall region and present the particle
preferential orientation in figure 13. The results show that the oblate
(prolate) particles in the near-wall region tend to orientate their symmetric
axis to the wall-normal (streamwise) direction, while in the spanwise
direction (figure 13c), particles of all explored cases show noteless
preferential orientations due to the unbounded periodic flow in $z-$direction.

Figure 12: Profiles of particle orientations in the wall-normal direction.
$\theta_{pi}$ is the angle between the symmetric axis of particles
$\boldsymbol{p}$ and the coordinate axis $\boldsymbol{i}$. The sketch in each
panel are shown to illustrate the preferential alignment of particles in the
near-wall region. Figure 13: Probability density functions (PDFs) of the
particle orientation in the near-wall region ($y^{+}\leq d_{v}^{+}$):
(a)$|cos\theta_{px}|$, (b) $|cos\theta_{py}|$, (c) $|cos\theta_{pz}|$.

##  4 Conclusion

We report a numerical study on suspensions of neutrally buoyant finite-size
particles in turbulent plane-Couette flow by employing a fully particle-
resolved method. The shape of particles varies from oblate to prolate by
changing their aspect ratio $\lambda$. The effects of particle aspect ratio on
turbulence modulation (in global arguments and local flows) and particle
statistics are studied, which agree well with previous findings in experiments
(Wang et al., 2022).

The turbulence modulation induced by suspended particles is investigated by
changing the aspect ratio $\lambda$, the volume fraction $\phi$ of the
particles and the bulk Reynolds number $Re$. For given $Re$ and $\phi$, it is
found that the spherical particles could cause the maximum drag increase
compared to non-spherical ones. We examine the field of fluid velocity and the
dissipation rate field and found that, at given $Re$ and $\phi$, the spherical
particles impact the boundary layer stronger than the non-spherical ones.
Through stress analysis, we found that the particle-induced stress $\tau_{p}$
accounts for a great contribution to the total wall stress, particularly in
the near-wall region. In addition, the magnitude of $\tau_{p}$ is larger for
the cases of spherical particles than that of non-spherical ones, which, as
shown by the particle statistics, could be attributed to the preferential
cluster of spherical particles in the near-wall regions. The contribution of
each term to the total stress in the stress balance is discussed by
integrating them within the whole domain. As a single-phase flow does, the
viscous stress and the turbulent stress of the fluid phase still account for
the major contributions to total stress in particle-laden flows. However, for
the latter cases, the global drag increases could be attributed to the
presence of particles, which give rise to the particle-phase turbulent stress
and the particle-induced stress. As $Re$ increases, the stress of the fluid
phase, particularly the turbulent stress $\tau_{Tf}$, becomes dominant, which
results in minor contributions coming from the particle phase and thereby the
smaller modulations in global drag.

The turbulent energy spectrum of particle-laden flow is discussed and the
results are found to be interesting. At the larger scale, the spectrum of the
cases of particle-laden collapse with that of single-phase flow and both
follow the scaling of $k^{-5/3}$, suggesting that the flow at large scales is
less affected by the presence of particles. While at scales smaller than the
equivalent-volume diameter of the particles ($d_{v}$), the spectrum of
particle-laden flow is found to decrease at slower rates than that of single-
phase flow. Through compensated plots, we found that the spectrum of particle-
laden flow shows good agreement with the scaling of $k^{-3}$ below $d_{v}$,
which is reminiscent of the observation in bubbly flow. However, by checking
the particle Reynolds number $Re_{p}^{slip}$, it is found that the
$Re_{p}^{slip}$ ( the order of $O(1)$ to $O(10)$) is too small for the
particles to affect the surrounding flows as the bubble does. Therefore,
future studies are required to understand the similarity in the small scale of
the turbulent energy spectrum found between the particle-laden flow and bubbly
flow, for instance, from the view of particle/bubble dynamics.

Additionally, particle statistics in different regions are found to be
affected by the aspect ratios of particles. Spherical and non-spherical
particles are found to preferentially cluster in the near-wall and bulk
regions, respectively, which are in line with findings in previous
experimental results (Wang et al., 2022). The different cluster effects, which
are associated with the intensity of particle collision rate, could account
for the magnitudes of particle-induced stress and the resulting global drag.
On the other side, spherical particles, as expected, show no preferential
alignment to any direction due to their perfect symmetry. Whereas, oblate and
prolate particles are found to preferentially align with the wall-normal and
stream-wise direction due to their anisotropy.

In collaboration with our previous experimental findings, this numerical study
provides deeper sights into the information of both fluid and particle phases,
allowing one to understand how the turbulence is modulated by the suspended
particles. Limited by computation consumption, the present work mainly focuses
on the role of particle aspect ratio, while the particle volume fractions are
kept at relatively small values. It is expected in future work to observe
greater turbulence modulation in the dense regime of particle suspensions.

Acknowledgements: We thank Enrico Calzavarini for the help with the numerical
simulations & the insightful discussions, and thank Jinghong Su, Dongpu Wang,
Lei Yi and Ning Zhu for the helpful discussions.

Funding: This work was supported by the National Natural Science Foundation of
China under grant no. 11988102, and the Tencent Foundation through the XPLORER
PRIZE.

Declaration of interests: The authors report no conflict of interest.

##  Appendix A Fluctuation of fluid and particle velocity

Since the fluid and particle velocity fluctuations contribute to the stress
decomposition in equation 8), it is obligatory to check how they converged in
time. The results are shown in figure 14 and figure 15. Only the cases at
$Re=1600$ and $\phi=2\%$ are shown here. For all cases with/without particles,
the fluid velocity fluctuates within reasonable ranges. As expected the
maximum fluctuation occurs in the streamwise direction since the flow is
driven by the horizontally moving walls. The fluctuations in wall-normal and
spanwise directions are rather small. The fluctuations of particle velocity
(figure 15) are similar to those of the fluid phase but with stronger
amplitudes. Additionally, some burst events in the particle velocity
fluctuation occur, which may be due to stronger particle collisions. To
conclude, these figures show that the velocity fluctuations of both phases
converged in time.

Figure 14: Fluid velocity fluctuation verse time at $Re=1600$: (a) Single-
phase case, (b) $\lambda=1/3$, (c) $\lambda=1$, (d) $\lambda=3$. For particle-
laden cases, $\phi=2\%$. The horizontal axis is denoted by the numerical time
steps $N_{t}$, and the vertical axis is the normalized fluid velocity
fluctuations. Figure 15: Particle velocity fluctuation verse time at $Re=1600$
and $\phi=2\%$: (a) $\lambda=1/3$, (b) $\lambda=1$, (c) $\lambda=3$. The
horizontal axis is denoted by the numerical time steps $N_{t}$, and the
vertical axis is the normalized particle velocity fluctuations.

## References

  * Alméras et al. (2017) Alméras, Elise, Mathai, Varghese, Lohse, Detlef & Sun, Chao 2017 Experimental investigation of the turbulence induced by a bubble swarm rising within incident turbulence. J. Fluid Mech. 825, 1091–1112.
  * Ardekani & Brandt (2019) Ardekani, M Niazi & Brandt, Luca 2019 Turbulence modulation in channel flow of finite-size spheroidal particles. J. Fluid Mech. 859, 887–901.
  * Ardekani et al. (2016) Ardekani, Mehdi Niazi, Costa, Pedro, Breugem, Wim Paul & Brandt, Luca 2016 Numerical study of the sedimentation of spheroidal particles. Int. J. Multiph. Flow 87, 16–34.
  * Ardekani et al. (2017) Ardekani, M Niazi, Costa, Pedro, Breugem, W-P, Picano, Francesco & Brandt, Luca 2017 Drag reduction in turbulent channel flow laden with finite-size oblate spheroids. J. Fluid Mech. 816, 43–70.
  * Assen et al. (2022) Assen, Martin PA, Ng, Chong Shen, Will, Jelle B, Stevens, Richard JAM, Lohse, Detlef & Verzicco, Roberto 2022 Strong alignment of prolate ellipsoids in taylor–couette flow. J. Fluid Mech. 935, A7.
  * Batchelor (1970) Batchelor, GK 1970 The stress system in a suspension of force-free particles. J. Fluid Mech. 41 (3), 545–570.
  * Bellani & Variano (2012) Bellani, Gabriele & Variano, Evan A 2012 Slip velocity of large neutrally buoyant particles in turbulent flows. New J. Phys 14 (12), 125009.
  * Brandt & Coletti (2022) Brandt, Luca & Coletti, Filippo 2022 Particle-laden turbulence: progress and perspectives. Annu. Rev. Fluid Mech. 54, 159–189.
  * Brenner (1961) Brenner, Howard 1961 The slow motion of a sphere through a viscous fluid towards a plane surface. Chem. Eng. Sci. 16 (3-4), 242–251.
  * Calzavarini (2019) Calzavarini, Enrico 2019 Eulerian–lagrangian fluid dynamics platform: The ch4-project. Softw. Impacts 1, 100002\.
  * Calzavarini et al. (2020) Calzavarini, Enrico, Jiang, Linfeng & Sun, Chao 2020 Anisotropic particles in two-dimensional convective turbulence. Phys. Fluids 32 (2), 023305.
  * Calzavarini et al. (2009) Calzavarini, Enrico, Volk, Romain, Bourgoin, Mickaël, Lévêque, Emmanuel, Pinton, J-F & Toschi, Federico 2009 Acceleration statistics of finite-sized particles in turbulent flow: the role of faxén forces. J. Fluid Mech. 630, 179–189.
  * Cisse et al. (2013) Cisse, Mamadou, Homann, Holger & Bec, Jérémie 2013 Slipping motion of large neutrally buoyant particles in turbulence. J. Fluid Mech. 735, R1.
  * Cooley & O’neill (1969) Cooley, MDA & O’neill, ME 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, Pedro, Boersma, Bendiks Jan, Westerweel, Jerry & Breugem, Wim-Paul 2015 Collision model for fully resolved simulations of flows laden with finite-size particles. Phys. Rev. E 92 (5), 053012.
  * de Motta et al. (2016) de Motta, JC Brändle, Estivalezes, Jean-Luc, Climent, Eric & Vincent, Stéphane 2016 Local dissipation properties and collision dynamics in a sustained homogeneous turbulent suspension composed of finite size particles. Int. J. Multiph. Flow 85, 369–379.
  * Demou et al. (2022) Demou, Andreas D, Ardekani, Mehdi Niazi, Mirbod, Parisa & Brandt, Luca 2022 Turbulent Rayleigh–Bénard convection in non-colloidal suspensions. J. Fluid Mech. 945, A6.
  * Dung et al. (2022) Dung, On-Yu, Waasdorp, Pim, Sun, Chao, Lohse, Detlef & Huisman, Sander G 2022 The emergence of bubble-induced scaling in thermal spectra in turbulence. arXiv:2207.05175.
  * Eshghinejadfard et al. (2018) Eshghinejadfard, Amir, Zhao, Lihao & Thévenin, Dominique 2018 Lattice boltzmann simulation of resolved oblate spheroids in wall turbulence. J. Fluid Mech. 849, 510–540.
  * Fiabane et al. (2012) Fiabane, Lionel, Zimmermann, Robert, Volk, Romain, Pinton, J-F & Bourgoin, Mickaël 2012 Clustering of finite-size particles in turbulence. Phys. Rev. E 86 (3), 035301\.
  * Grossmann et al. (2016) Grossmann, Siegfried, Lohse, Detlef & Sun, Chao 2016 High-Reynolds number Taylor-Couette turbulence. Annu. Rev. Fluid Mech. 48, 53–80.
  * Hidman et al. (2022) Hidman, Niklas, Ström, Henrik, Sasic, Srdjan & Sardina, Gaetano 2022 Assessing passive scalar dynamics in bubble-induced turbulence using DNS. arXiv:2211.06293.
  * Homann & Bec (2010) Homann, Holger & Bec, Jeremie 2010 Finite-size effects in the dynamics of neutrally buoyant particles in turbulent flow. J. Fluid Mech. 651, 81–91.
  * Innocenti et al. (2021) Innocenti, Alessio, Jaccod, Alice, Popinet, Stéphane & Chibbaro, Sergio 2021 Direct numerical simulation of bubble-induced turbulence. J. Fluid Mech. 918, A23.
  * Jiang et al. (2020) Jiang, Linfeng, Calzavarini, Enrico & Sun, Chao 2020 Rotation of anisotropic particles in Rayleigh–Bénard turbulence. J. Fluid Mech. 901, A8.
  * Jiang et al. (2021) Jiang, Linfeng, Wang, Cheng, Liu, Shuang, Sun, Chao & Calzavarini, Enrico 2021 Rotational dynamics of bottom-heavy rods in turbulence from experiments and numerical simulations. Theor. Appl. Mech. Lett. 11 (1), 100227.
  * Jiang et al. (2022) Jiang, Linfeng, Wang, Cheng, Liu, Shuang, Sun, Chao & Calzavarini, Enrico 2022 Dynamics of finite-size spheroids in turbulent flow: the roles of flow structures and particle boundary layers. J. Fluid Mech. 939, A22.
  * Lance & Bataille (1991) Lance, M & Bataille, J 1991 Turbulence in the liquid phase of a uniform bubbly air–water flow. J. Fluid Mech. 222, 95–118.
  * Li et al. (2022) Li, Yang, Xia, Zhenhua & Wang, Lian-Ping 2022 Inertial migration of a neutrally buoyant oblate spheroid in three-dimensional square duct poiseuille flows. Int. J. Multiph. Flow 155, 104148.
  * Lucci et al. (2010) Lucci, Francesco, Ferrante, Antonino & Elghobashi, Said 2010 Modulation of isotropic turbulence by particles of taylor length-scale size. J. Fluid Mech. 650, 5–55.
  * Lundell et al. (2011) Lundell, Fredrik, Söderberg, L Daniel & Alfredsson, P Henrik 2011 Fluid mechanics of papermaking. Annu. Rev. Fluid Mech. 43, 195–217.
  * Luo et al. (2007) Luo, Kun, Wang, Zeli, Fan, Jianren & Cen, Kefa 2007 Full-scale solutions to particle-laden flows: Multidirect forcing and immersed boundary method. Phys. Rev. E 76 (6), 066709\.
  * Mathai et al. (2016) Mathai, Varghese, Calzavarini, Enrico, Brons, Jon, Sun, Chao & Lohse, Detlef 2016 Microbubbles and microparticles are not faithful tracers of turbulent acceleration. Phys. Rev. Lett. 117 (2), 024501.
  * Mathai et al. (2020) Mathai, Varghese, Lohse, Detlef & Sun, Chao 2020 Bubbly and buoyant particle–laden turbulent flows. Annu. Rev. Cond. Matt. Phys. 11, 529–559.
  * Mazzitelli & Lohse (2009) Mazzitelli, Irene M & Lohse, Detlef 2009 Evolution of energy in flow driven by rising bubbles. Phys. Rev. E 79 (6), 066317.
  * Mendez-Diaz et al. (2013) Mendez-Diaz, S, Serrano-Garcia, JC, Zenit, R & Hernandez-Cordero, JA 2013 Power spectral distributions of pseudo-turbulent bubbly flows. Phys. Fluids 25 (4), 043303\.
  * Mercado et al. (2010) Mercado, Julian Martinez, Gomez, Daniel Chehata, Van Gils, Dennis, Sun, Chao & Lohse, Detlef 2010 On bubble clustering and energy spectra in pseudo-turbulence. J. Fluid Mech. 650, 287–306.
  * Mittal et al. (2020) Mittal, Rajat, Ni, Rui & Seo, Jung-Hee 2020 The flow physics of covid-19. J. Fluid Mech. 894, F2.
  * Moffet & Prather (2009) Moffet, Ryan C & Prather, Kimberly A 2009 In-situ measurements of the mixing state and optical properties of soot with implications for radiative forcing estimates. Proc. Natl. Acad. Sci. U.S.A. 106 (29), 11872–11877.
  * Obligado & Bourgoin (2022) Obligado, Martin & Bourgoin, Mickaël 2022 Dynamics of towed particles in a turbulent flow. J. Fluids Struct. 114, 103704.
  * Owolabi & Lin (2018) Owolabi, Bayode E & Lin, Chao-An 2018 Marginally turbulent couette flow in a spanwise confined passage of square cross section. Phys. Fluids 30 (7), 075102.
  * Pandey et al. (2022) Pandey, Vikash, Mitra, Dhrubaditya & Perlekar, Prasad 2022 Kolmogorov turbulence co-exists with pseudo-turbulence in buoyancy-driven bubbly flows. arXiv:2204.04505.
  * Pandey et al. (2020) Pandey, Vikash, Ramadugu, Rashmi & Perlekar, Prasad 2020 Liquid velocity fluctuations and energy spectra in three-dimensional buoyancy-driven bubbly flows. J. Fluid Mech. 884, R6.
  * Pedley & Kessler (1992) Pedley, TJ & Kessler, John O 1992 Hydrodynamic phenomena in suspensions of swimming microorganisms. Annu. Rev. Fluid Mech. 24 (1), 313–358.
  * Peng et al. (2020) Peng, Cheng, Ayala, Orlando M & Wang, Lian-Ping 2020 Flow modulation by a few fixed spherical particles in a turbulent channel flow. J. Fluid Mech. 884, A15.
  * Peskin (2002) Peskin, Charles S 2002 The immersed boundary method. Acta numerica 11, 479–517.
  * Picano et al. (2015) Picano, Francesco, Breugem, Wim-Paul & Brandt, Luca 2015 Turbulent channel flow of dense suspensions of neutrally buoyant spheres. J. Fluid Mech. 764, 463–487.
  * Pope (2000) Pope, Stephen B. 2000 Turbulent Flows. Cambridge University Press.
  * Prakash et al. (2016) Prakash, Vivek N., Martínez Mercado, J., van Wijngaarden, Leen, Mancilla, E., Tagawa, Y., Lohse, Detlef & Sun, Chao 2016 Energy spectra in turbulent bubbly flows. J. Fluid Mech. 791, 174–190.
  * Qiu et al. (2022) Qiu, Jingran, Mousavi, Navid, Zhao, Lihao & Gustavsson, Kristian 2022 Active gyrotactic stability of microswimmers using hydromechanical signals. Phys. Rev. Fluids 7 (1), 014311.
  * Qureshi et al. (2008) Qureshi, Nauman M, Arrieta, Unai, Baudet, Christophe, Cartellier, Alain, Gagne, Yves & Bourgoin, Mickaël 2008 Acceleration statistics of inertial particles in turbulent flow. Eur. Phys. J. B 66 (4), 531–536.
  * Qureshi et al. (2007) Qureshi, Nauman M, Bourgoin, Mickaël, Baudet, Christophe, Cartellier, Alain & Gagne, Yves 2007 Turbulent transport of material particles: an experimental study of finite size effects. Phys. Rev. Lett. 99 (18), 184502.
  * Riboux et al. (2010) Riboux, Guillaume, Risso, Frédéric & Legendre, Dominique 2010 Experimental characterization of the agitation generated by bubbles rising at high reynolds number. J. Fluid Mech. 643, 509–539.
  * Risso (2018) Risso, Frédéric 2018 Agitation, mixing, and transfers induced by bubbles. Annu. Rev. Fluid Mech. 50, 25–48.
  * Sabban & van Hout (2011) Sabban, Lilach & van Hout, René 2011 Measurements of pollen grain dispersal in still air and stationary, near homogeneous, isotropic turbulence. J. Aerosol. Sci. 42 (12), 867–882.
  * Stocker (2012) Stocker, Roman 2012 Marine microbes see a sea of gradients. Science 338 (6107), 628–633.
  * Voth & Soldati (2017) Voth, Greg A & Soldati, Alfredo 2017 Anisotropic particles in turbulence. Annu. Rev. Fluid Mech. 49 (1), 249–276.
  * Wang et al. (2021) Wang, Cheng, Jiang, Linfeng, Jiang, Hechuan, Sun, Chao & Liu, Shuang 2021 Heat transfer and flow structure of two-dimensional thermal convection over ratchet surfaces. J. Hydrodyn. 33 (5), 970–978.
  * Wang et al. (2022) Wang, Cheng, Yi, Lei, Jiang, Linfeng & Sun, Chao 2022 How do the finite-size particles modify the drag in taylor–couette turbulent flow. J. Fluid Mech. 937, A15.
  * Wang et al. (2017a) Wang, Guiquan, Abbas, Micheline & Climent, Éric 2017a Modulation of large-scale structures by neutrally buoyant and inertial finite-size particles in turbulent couette flow. Phys. Rev. Fluids 2 (8), 084302.
  * Wang et al. (2018) Wang, Guiquan, Abbas, Micheline, Yu, Zhaosheng, Pedrono, Annaig & Climent, Eric 2018 Transport of finite-size particles in a turbulent couette flow: The effect of particle shape and inertia. Int. J. Multiph. Flow 107, 168–181.
  * Wang et al. (2017b) Wang, Yayun, Sierakowski, Adam J & Prosperetti, Andrea 2017b Fully-resolved simulation of particulate flows with particles–fluid heat transfer. J. Comput. Phys. 350, 638–656.
  * Wang et al. (2019) Wang, Ziqi, Mathai, Varghese & Sun, Chao 2019 Self-sustained biphasic catalytic particle turbulence. Nat. Commun. 10 (1), 3333.
  * Will & Krug (2021a) Will, Jelle B & Krug, Dominik 2021a Dynamics of freely rising spheres: the effect of moment of inertia. J. Fluid Mech. 927, A7.
  * Will & Krug (2021b) Will, Jelle B & Krug, Dominik 2021b Rising and sinking in resonance: Mass distribution critically affects buoyancy-driven spheres via rotational dynamics. Phys. Rev. Lett. 126 (17), 174502.
  * Will et al. (2021) Will, Jelle B, Mathai, Varghese, Huisman, Sander G, Lohse, Detlef, Sun, Chao & Krug, Dominik 2021 Kinematics and dynamics of freely rising spheroids at high reynolds numbers. J. Fluid Mech. 912, A16.
  * Yousefi et al. (2020) Yousefi, Ali, Ardekani, Mehdi Niazi & Brandt, Luca 2020 Modulation of turbulence by finite-size particles in statistically steady-state homogeneous shear turbulence. J. Fluid Mech. 899, A19.
  * Zhang et al. (2020) Zhang, Lingxin, Zhou, Zecai & Shao, Xueming 2020 Numerical investigation on the drag force of a single bubble and bubble swarm. J. Fluid Mech. 32(6), 1043–1049.
  * Zhang & Prosperetti (2010) Zhang, Quan & Prosperetti, Andrea 2010 Physics-based analysis of the hydrodynamic stress in a fluid-particle system. Phys. Fluids 22 (3), 033306.