1. Introduction
By adding a small amount of solid particles to turbulent flows, we can drastically modulate turbulence. This interesting phenomenon has been well known since the last century when seminal experiments of particle-laden pipe flow (Tsuji & Morikawa Reference Tsuji and Morikawa1982) and free jet (Levy & Lockwood Reference Levy and Lockwood1981) were conducted. Gore & Crowe (Reference Gore and Crowe1989) summarised experimental data at that time to conclude that the particle size
$D$
compared with the integral length
$L$
(i.e. the size of the largest eddies) of turbulence is the important parameter and turbulence can be augmented (or attenuated) when
$D/L$
is larger (or smaller) than
$0.1$
. The target of the present study is the attenuation of turbulence by particles smaller than
$L$
.
Gore & Crowe (Reference Gore and Crowe1989) already mentioned the mechanism of the turbulence attenuation, where small particles acquire the kinetic energy of turbulence from the most energetic (i.e. the largest) eddies. However, the physical mechanism of turbulence modulation is still under debate because it is difficult to draw a concrete picture of the complex phenomenon only by experiments, and because it is also difficult to conduct parametric study in laboratory. In contrast, since the seminal study by Elghobashi (Reference Elghobashi1994), direct numerical simulations (DNS) of multi-phase turbulence have been developing to be powerful enough for effective parameter surveys. Therefore, the combination of state-of-the-art measurements and DNS is rapidly advancing the understanding of turbulence modulation by solid particles (Balachandar & Eaton Reference Balachandar and Eaton2010; Brandt & Coletti Reference Brandt and Coletti2022). Even if we only look at the literature in the simplest case, i.e. the modulation of turbulence in a periodic cube, which is also the target of the present study, a number of DNS were conducted (ten Cate et al. Reference ten Cate, Derksen, Portela and van den Akker2004; Homann & Bec Reference Homann and Bec2010; Yeo et al. Reference Yeo, Dong, Climent and Maxey2010; Lucci, Ferrante & Elghobashi Reference Lucci, Ferrante and Elghobashi2010, Reference Lucci, Ferrante and Elghobashi2011; Gao, Li & Wang Reference Gao, Li and Wang2013; Wang et al. Reference Wang, Ayala, Gao, Andersen and Mathews2014; Schneiders, Meinke & Schröder Reference Schneiders, Meinke and Schröder2017; Oka & Goto Reference Oka and Goto2022; Peng, Sun & Wang Reference Peng, Sun and Wang2023; Su et al. Reference Su, Zhang, Fu, Xiang and Wang2023; Cannon, Olivieri & Rosti Reference Cannon, Olivieri and Rosti2024; Chiarini, Cannon & Rosti Reference Chiarini, Cannon and Rosti2024) to reveal the statistics and dynamics of particle-turbulence interactions. Among these studies, it is particularly important to note the observation made by ten Cate et al. (Reference ten Cate, Derksen, Portela and van den Akker2004), Yeo et al. (Reference Yeo, Dong, Climent and Maxey2010), Gao et al. (Reference Gao, Li and Wang2013) and Wang et al. (Reference Wang, Ayala, Gao, Andersen and Mathews2014) that, when turbulence is attenuated, the energy spectrum is attenuated (and augmented) in the wavenumber range corresponding to length scales larger (and smaller) than particles, because this is consistent with the physical mechanism of turbulence attenuation examined in the present article; namely, particles acquire the energy from the largest eddies and dissipate it in the shedding particle-size vortices around them.
Thanks to parametric studies by DNS and accumulations of experimental data, the condition of turbulence attenuation (or augmentation) was proposed in terms of not only the normalised particle size
$D/L$
(Gore & Crowe Reference Gore and Crowe1989), but also other parameters such as the volume fraction
$\Lambda$
, the mass loading and the Stokes number,

where
$\tau _p$
and
$T$
denote the particle velocity relaxation time due to the Stokes drag and the turnover time of the largest vortices, respectively. For example, Tanaka & Eaton (Reference Tanaka and Eaton2008) proposed a dimensionless number and even recently, Peng et al. (Reference Peng, Sun and Wang2023) also suggested another parametrisation for the turbulence modulation. Note that in the present study, we define the Stokes number (1.1) by
$T$
because it quantifies the particles’ ability to follow the most energetic motion in turbulence, but other definitions of the Stokes number are sometimes used (see e.g. Bordoloi & Variano Reference Bordoloi and Variano2017). Note also that the target of the present study is turbulence modulation by particles much heavier than carrier fluid, where added mass effects are negligible.
We also investigated this issue by DNS of the attenuation of turbulence in a periodic cube due to spherical particles (Oka & Goto Reference Oka and Goto2022) to show that the turbulent kinetic energy
$K'$
can be reduced by approximately half, even if the volume fraction
$\Lambda$
is as small as
$8\times 10^{-3}$
, if
$St$
is sufficiently larger than
$1$
and
$D$
is smaller than
$L$
. This implies that both
$D/L$
and
$St$
are important parameters. We also demonstrated that when
$St$
is large enough, particle-size vortices are shed from particles because of large-enough relative velocity between fluid and particles, and then the energy which the particles acquire from the largest eddies is dissipated in these shedding vortices. Then, we showed that we can estimate the attenuation rate,

where
$K^{\prime}_0$
is the average turbulent kinetic energy of the single-phase flow, by considering the average energy dissipation rate
$\epsilon _p$
in the wake. More precisely, the average energy input rate
$\epsilon _0$
, which is assumed to be the same as in the single-phase turbulence, is the sum of the average flux
$\epsilon _c$
of the energy cascade in the particle-laden turbulence and the average dissipation rate
$\epsilon _p$
around particles:

Since
$\epsilon _0$
and
$\epsilon _c$
are the energy flux in the single- and multi-phase flows, respectively, we can estimate them as

by the dissipation law of Taylor (Reference Taylor1935). Here,
$K_0$
and
$L_0$
denote the kinetic energy of mean flow and integral length in the single-phase flow, respectively, which are assumed to be the same as those in the particle-laden flow (see the verification of these assumptions in Appendix A). The coefficient
$C_\epsilon$
is the dissipation coefficient, which depends on flow (Goto & Vassilicos Reference Goto and Vassilicos2009). Substituting (1.4) into (1.3), we can relate
$\mathcal {A}$
to
$\epsilon _p$
as

with
$\alpha$
being defined as the ratio
${K_{0}}/{K^{\prime}_{0}}$
of the kinetic energy of mean flow to the turbulent kinetic energy in the single-phase flow. Oka & Goto (Reference Oka and Goto2022) numerically verified (1.5) by estimating the energy dissipation rate around the particles as

where
$\Delta u$
is the relative velocity between fluid and particles,
$C_p$
is a constant of
$O(1)$
, and
$\langle \:\:\cdot \:\:\rangle _p$
denotes the average over particles and time. The superscript dagger of
$\epsilon _p^\dagger$
indicates a theoretical expression (1.6) for
$\epsilon _p$
. Although the physical picture shown here is based on the classical view described in the second paragraph of this section, the crucial point to derive (1.5) is that we consider the change in the cascading energy flux. Recently, Balachandar, Peng & Wang (Reference Balachandar, Peng and Wang2024) further developed this picture towards the subgrid modelling of particle-laden turbulence.
Although Oka & Goto (Reference Oka and Goto2022) numerically demonstrated that (1.5) can accurately estimate the attenuation rate
$\mathcal {A}$
, (1.2), this demonstration was made only for spherical particles. However, it is known that anisotropic particles sometimes more effectively modify turbulence than spherical ones (Voth & Soldati Reference Voth and Soldati2017, § 6.2). For example, Paschkewitz et al. (Reference Paschkewitz, Dubief, Dimitropoulos, Shaqfeh and Moin2004) showed by DNS with a constitutive equation that rigid fibres can significantly reduce the drag of turbulent channel flow. Many other researchers also investigated the issue of turbulence attenuation and turbulent drag reduction due to anisotropic particles by DNS (Eshghinejadfard, Hosseini & Thévenin Reference Eshghinejadfard, Hosseini and Thévenin2017; Ardekani & Brandt Reference Ardekani and Brandt2019; Schneiders et al. Reference Schneiders, Fröhlich, Meinke and Schröder2019; Wang, Xu & Zhao Reference Wang, Xu and Zhao2021; Olivieri, Cannon & Rosti Reference Olivieri, Cannon and Rosti2022; Cannon et al. Reference Cannon, Olivieri and Rosti2024) and experiments (Ljus, Johansson & Almstedt Reference Ljus, Johansson and Almstedt2002; Bellani et al. Reference Bellani, Byron, Collignon, Meyer and Variano2012; Capone, Romano & Soldati Reference Capone, Romano and Soldati2015). As a study related to the present one, Zhao, George & van Wachem (Reference Zhao, George and van Wachem2015) numerically showed that turbulence is more significantly attenuated with anisotropic particles with larger
$St$
and larger aspect ratios. Mandø (Reference Mandø2009) also examined, by experiments, the modulation of turbulent jet by particles with different shapes and they observed that for a common equivalent diameter,

where
$V$
is the volume of a particle, anisotropic particles (disks and spheroids) attenuate turbulence intensity more than spherical ones.
These previous studies motivate us to investigate the effects of particle shape in turbulence attenuation. However, to systematically investigate the effects, we must examine cases with different shapes of particles under each condition of particles (i.e.
$\Lambda$
,
$D_*/L$
and
$St$
) and flow (i.e. boundary conditions and the Reynolds number). Thus, in the present study, we restrict ourselves in the case where the particle shape is prolate spheroidal and turbulence is in a triply periodic domain without walls. More concretely, we conduct DNS with different values of
$D_*/L$
,
$St$
and the aspect ratio
$\chi$
of prolate spheroidal particles for a fixed value of the volume fraction
$\Lambda$
and the Reynolds number of turbulence. Then, we demonstrate that the above-mentioned physical mechanism of turbulence attenuation is valid even for the anisotropic particles and the estimation (1.5) of
$\mathcal {A}$
in terms of
$\epsilon _p$
holds. We also generalise (1.6) for
$\epsilon _p$
around anisotropic particles.
2. Direct numerical simulation
2.1. Numerical method
We conduct DNS of turbulence with solid prolate spheroids in a periodic cube using the same numerical code used in our previous study (Fujiki et al. Reference Fujiki, Awai, Motoori and Goto2024), in the appendix of which we showed the validation of the code. Here, we briefly describe the numerical method.
First, for fluid motion, we suppose that it is governed by the Navier–Stokes equation,

and the continuity equation,

for an incompressible fluid. Here,
$\boldsymbol {u}(\boldsymbol {x}, t)$
is the velocity field at position
$\boldsymbol {x}$
and time
$t$
,
$p(\boldsymbol {x}, t)$
is the pressure field, and
$\rho _{f}$
and
$\nu$
are the fluid mass density and kinematic viscosity, respectively. In (2.1),
$\boldsymbol {f}(\boldsymbol {x})$
is an external steady force, which is expressed as

(Goto, Saito & Kawahara Reference Goto, Saito and Kawahara2017; Oka & Goto Reference Oka and Goto2022), and
$\boldsymbol {f}^{\leftarrow p}(\boldsymbol {x}, t)$
is the force due to particles. Using the simplified mark and cell (SMAC) method, we numerically solve (2.1) and (2.2) under periodic boundary conditions in a cube with side
$L_{box}=2\pi$
). The spatial derivatives are evaluated using the second-order central difference method. The convection and viscous terms are temporally integrated using the second-order Adams–Bashforth method and Crank–Nicholson method, respectively.
Table 1 lists the parameters such as the number
$N^3$
of grid points for DNS and statistics of the simulated turbulence. Here, we have evaluated the integral length as
$ L = 3\pi \int _{0}^{\infty }{k^{-1}E(k)\mathrm {d}k}/4 \int _{0}^{\infty }{E(k)\mathrm {d}k}$
, where
$E(k)$
is the temporally averaged energy spectrum of turbulence, and the Taylor length as
$ \lambda = \sqrt {10\nu K^{\prime } / \epsilon }$
, where

is the average turbulent kinetic energy per unit mass and
$\epsilon$
is its dissipation rate. Here,
$ \langle \ \cdot \ \rangle$
and
$ \overline {\ \cdot \ }$
denote the spatial and temporal averages, respectively. Then, we have evaluated the Taylor-length-based Reynolds number as
$ R_{\lambda } = u^{\prime }\lambda /\nu$
with
$u^{\prime }=\sqrt {2K^{\prime }/3}$
and the Kolmogorov length as
$\eta =\epsilon ^{-1/4}\nu ^{3/4}$
. In table 1, we also list the Courant–Friedrich–Lewy (CFL) number, which we have set sufficiently smaller than unity for numerical stability.
Table 1. Parameters and statistics of the single-phase turbulence. We have taken temporal averages over approximately
$250T$
with
$T={L}/{u^{\prime }}$
being the integral time. Here,
$R_\lambda$
,
$L$
and
$\eta$
are the Taylor-length-based Reynolds number, integral length and Kolmogorov length, respectively; and
$c={u^{\prime }}\Delta t/\Delta x$
is the CFL number defined by the fluctuating velocity
$u^{\prime }$
, the grid width
$\Delta x = L_{box}/N$
and the time increment
$\Delta t$
for the temporal integration.

Next, we describe the numerical method to simulate the motion of spheroidal particles with mass density
$\rho _{p}$
, major radius
$a$
and minor radius
$b$
. Particle’s velocity
$\boldsymbol {v}_{p}(t)$
and angular velocity
$\boldsymbol {\omega }_{p}(t)$
obey Newton’s equations of motion,

and

respectively. Here,
$m$
(
$=4\pi \rho _{p} a b^2 /3$
) is the particle mass and
$\boldsymbol {I}$
is the moment of inertia. On the right-hand sides of these equations,
$\boldsymbol {f}^{\leftarrow f}$
and
$\boldsymbol {T}^{\leftarrow f}$
denote the force and torque acting on a particle by the surrounding fluid, respectively, and
$\boldsymbol {f}^{\leftrightarrow p}$
and
$\boldsymbol {T}^{\leftrightarrow p}$
denote the interaction force and torque between particles, respectively. In the present DNS, since we neglect gravity, the energy is injected only by the external force (2.3) on fluid.
We treat particle-fluid interaction by an immersed boundary method proposed by Uhlmann (Reference Uhlmann2005) to evaluate forces
$\boldsymbol {f}^{\leftarrow p}$
in (2.1) and
$\boldsymbol {f}^{\leftarrow f}$
in (2.5) and torque
$\boldsymbol {T}^{\leftarrow f}$
in (2.6). On the other hand, we treat particle-particle interaction by the discrete element method (DEM) to evaluate
$\boldsymbol {f}^{\leftrightarrow p}$
in (2.5) and
$\boldsymbol {T}^{\leftrightarrow p}$
in (2.6). In the present DEM, similarly to Moriche et al. (Reference Moriche, Hettmann, García-Villalba and Uhlmann2023), we consider only the normal component of the contact force due to elastic collisions, which is proportional to the squared overlap distance.
To simulate the rotational motion governed by (2.6), we employ the method used by Moriche, Uhlmann & Dušek (Reference Moriche, Uhlmann and Dušek2021), where we introduce the frame attached to each particle with its origin being the particle centre and with axes parallel to the particle’s major and minor axes. Using the angular velocity
$\boldsymbol {\omega }_{p}'$
in this frame, we can rewrite (2.6) as

Here,
$\boldsymbol {I}'$
is the principal moment of inertia,

and
$\boldsymbol {R}$
is the rotation matrix between the two frames expressed as

with quaternion
$\boldsymbol {q}= (q_{0}, q_{1}, q_{2}, q_{3} )$
. In terms of
$\boldsymbol {\omega }'_{p}=({\omega }'_{px},{\omega }'_{py},{\omega }'_{pz})$
, we can express the temporal evolution of the quaternion as

Thus, we solve (2.7) and (2.10) instead of (2.6) for the rotational motion of particles. We numerically integrate (2.5), (2.7) and (2.10) with the Euler method for the particle-particle interaction terms and the second-order Adams–Bashforth method for the other terms. We solve the governing equations for fluid and particles alternately so that we can treat coupled motion of fluid and particles.
Note that we have introduced in (2.8) an artificial coefficient
$\beta$
to express the inhomogeneous distribution of mass density in a particle;
$\beta =1$
corresponds to the homogeneous density distribution. By the introduction of
$\beta$
, we can discuss whether the rotation of particles is important or not in turbulence modulation (see § 3.3).
Table 2. Particle parameters for (
$a$
)
$D_*/{L}=0.16$
(RUN384 in table 1) and (
$b$
)
$D_*/{L}=0.24$
(RUN256).

2.2. Particle parameters
We examine interaction between prolate spheroidal particles and the statistically stationary turbulence with properties given in table 1. Then, the particles are characterised by four parameters: the equivalent diameter
$D_*=2(ab^2)^{1/3}$
, which is defined as (1.7), the aspect ratio
$\chi =a/b$
(
$\chi \geqslant 1$
), the mass density ratio
$\gamma = \rho _p/ \rho _f$
between the particle and fluid, and coefficient
$\beta$
in (2.8) determining the inertial moment. We fix the particle volume fraction
$\Lambda =(4\pi ab^2 N_{p}/3L^3_{0})$
, where
$N_{p}$
is the number of particles, at
$6.0\times 10^{-3}$
, which is slightly more dilute than the system in our previous study (Oka & Goto Reference Oka and Goto2022). Note that since
$\Lambda$
is fixed, the number
$N_p$
of particles increases as
$D_*$
decreases.
First, we describe particle parameters with homogeneous mass density (
$\beta =1$
). We consider the two cases with relatively small diameters
$D_*/{L}=0.16$
and
$0.24$
(i.e.
$D_*/{\eta }=7.5$
and
$11$
). Here, we use the values of
$L$
and
$\eta$
in the single-phase flow. Then, we change the aspect ratio (
$\chi =1$
,
$3$
,
$5$
and
$7$
) and mass density ratio (
$\gamma =8$
,
$32$
,
$128$
and
$512$
) to examine
$32$
cases in total. We show results on the basis of the Stokes number (1.1) instead of
$\gamma$
. The Stokes number
$St$
denotes the ratio between the velocity relaxation time,

where we define
$\xi =\sqrt {\chi ^2-1}$
and
$f(\chi )=\chi \ln (\chi +\xi )/\xi$
, of the spheroidal particles (Shapiro & Goldenberg Reference Shapiro and Goldenberg1993) and the integral time
$T$
of single-phase turbulence. Note that (2.11) is a reference particle relaxation time estimated by assuming the Stokes drag. We list in table 2 the parameters of particles. The number of grid points in DNS is
$N^3=384^3$
for
$D_*/{L}=0.16$
and
$N^3=256^3$
for
$D_*/{L}=0.24$
to ensure
$D_*/\Delta x=12$
, where
$\Delta x$
(
$=L_{box}/N$
) is the grid width.
Next, we explain the additional particle parameter
$\beta$
, which expresses the inhomogeneous mass density distribution inside each particle. Note that, by changing the value of
$\beta$
in (2.8), we can change the particles’ ability of rotation. When
$\beta =1$
, the mass density is homogeneously distributed; whereas when
$\beta \lt 1$
, it is larger around the particle centre and such a particle can rotate more easily than the one with
$\beta =1$
. Since
$\beta$
cannot be arbitrarily large, we examine particles with
$\beta \leqslant 1$
. To quantify the rotational ability of particles, we define the rotational Stokes number as

which denotes the ratio between the relaxation time (Jayaram et al. Reference Jayaram, Jie, Zhao and Andersson2023),

of the particle tumbling motion and the turnover time,
$\tau _{D_*}$
(
$=\epsilon _0^{-1/3}D_*^{2/3}$
), of vortices at scale
$D_*$
. In (2.13),
$c_0$
and
$c_1$
are

respectively. The rotational Stokes number
$St_{r}$
represents the ability for the particles to follow the swirling of vortices with size
$D_*$
. In the present study, in addition to the
$32$
kinds of particles with homogeneous density (
$\beta =1$
), we simulate particles with different values of
$St_{r}$
corresponding to
$\beta =1/4$
,
$1/16$
and
$1/64$
in the eight cases of the most elongated (
$\chi =7$
) (see table 2).
Before showing the results, we recall that turbulence driven by (2.3) has mean flow. In Appendix A, we show that particles with the dilute volume fraction do not significantly change the average kinetic energy
$K$
of the mean flow or the integral length
$L$
. These features are important because they are assumed in (1.4), which is required to derive the formula (1.5) of attenuation rate. In the next section, we demonstrate that (1.5) holds even for anisotropic particles.
3. Results
3.1. Attenuation of turbulent kinetic energy
In the following analysis, we evaluate statistics by taking temporal average over approximately
$250T$
. The error bars in the following figures indicate the standard deviation. We first show results for particles with homogeneous mass density distribution (
$\beta =1$
).
Figure 1 shows the
$St$
-dependence of the temporal average
$K'$
of turbulent kinetic energy normalised by the value
$K^{\prime}_{0}$
for the single-phase flow. The blue smaller symbols are results for
$D_*/{L}=0.16$
and black larger ones are for
$D_*/{L}=0.24$
. The shape of symbols with different ellipticity represents the particle shape with different aspect ratios
$\chi = 1$
,
$3$
,
$5$
and
$7$
. Note that the shape of the symbols is schematic. We see in this figure that, for given
$D_*$
and
$\chi$
, turbulent kinetic energy is attenuated more significantly for larger
$St$
. This tendency is common in both cases with spherical particles (circles in the figure) and with spheroidal ones (ellipses). We therefore conclude that the largeness of
$St$
is a necessary condition for turbulence attenuation irrespective of particle shape.

Figure 1. Average turbulent kinetic energy
$K'$
normalised by the value
$K^{\prime}_{0}$
for the single-phase flow as a function of
$St$
. Blue smaller symbols are results for
$D_*/{L}=0.16$
and black larger ones are for
$D_*/{L}=0.24$
. The shape of symbols with different ellipticity represents the aspect ratios (
$\chi =1$
,
$3$
,
$5$
and
$7$
) of particles; note that the shape of symbols is schematic and different from the examined spheroids. Error bars indicates the standard deviation.
Our previous study (Oka & Goto Reference Oka and Goto2022) for spherical particles demonstrated that shedding vortices behind particles play important roles in turbulence attenuation because they are relevant for the additional dissipation rate
$\epsilon _p$
. Similar vortices can be observed around spheroidal particles with large
$St$
. Figures 2(a) and 2(b) show vortices (blue objects) identified by isosurfaces of the enstrophy and particles (yellow ones). The particles shown in these two panels have the same size (
$D_*/{L}=0.16$
) and aspect ratio (
$\chi =7$
) but different values of the Stokes number: (
$a$
)
$St=1.3$
and (
$b$
)
$86$
. Panels (
$c$
) and (
$d$
) show a cross-section of panels (
$a$
) and (
$b$
), respectively, which is perpendicular to the four columnar vortices driven by external force (2.3). Since the most energetic large-scale flow is parallel to this cross-section, we can easily judge whether or not shedding vortices exist by the enstrophy indicated by the background colour in panels (
$c$
) and (
$d$
). We see in the visualisations of panel (
$a,c$
) for
$St=1.3$
(and supplementary movie 1 available at https://doi.org/10.1017/jfm.2025.137) that the enstrophy is distributed independently of particles. This implies that vortices are generated through energy cascade even in the presence of particles, and there is no additional energy dissipation around them in this case with small
$St$
. In contrast, vortices are created behind particles with larger Stokes number (
$St=86$
) as shown in figure 2(b,d). Figure 2(e) shows the magnification of a subdomain of figure 2(b). It is clear in figure 2(b,d,e) that elongated vortices are shed along particles. Since these shedding vortices are as small as
$10{\eta }$
, they are almost immediately dying due to viscous effects. This additional energy dissipation in the shedding vortices leads to the reduction in part of the cascading energy flux, and consequently turbulent kinetic energy is attenuated. This is consistent with the physical picture described in the introduction.
Figure 1 also shows that the attenuation rate of turbulent kinetic energy depends on not only
$St$
, but also equivalent diameter
$D_{*}$
and aspect ratio
$\chi$
. More specifically, for given
$St$
, turbulent kinetic energy is attenuated more significantly for smaller
$D_{*}$
or larger
$\chi$
.
Concerning particle size
$D_*$
, we can explain its effect on the turbulence attenuation by the same argument as for spherical particles (Oka & Goto Reference Oka and Goto2022); that is, smaller particles can produce larger additional energy dissipation rate
$\epsilon _p$
per unit mass because it is inversely proportional to
$D_*$
[see (1.6)]. We emphasise that this statement is valid under the condition that the volume fraction
$\Lambda$
is constant.

Figure 2. Visualisation of vortices (blue objects) identified by isosurfaces of enstrophy and spheroidal particles (yellow ones) with (
$a$
)
$St=1.3$
and (
$b$
)
$St=86$
. Particle size is
$D_*/{L}=0.16$
and aspect ratio is
$\chi =7$
. Panels
$(c)$
and
$(d)$
show
$z=0$
planes of panels (
$a$
) and (
$b$
), respectively. The background colour indicates the magnitude of the enstrophy (lighter colours show the larger value). Panel (
$e$
) shows a subdomain of panel (
$b$
). See also supplementary movie 1.
However, concerning aspect ratio
$\chi$
, it is unclear whether we can describe the
$\chi$
-dependence solely in terms of the additional energy dissipation rate
$\epsilon _{p}$
around particles. This is a crucial issue because for larger Stokes numbers (say,
$St \gtrsim 20$
), more elongated particles attenuate turbulence more significantly and the
$\chi$
-dependence of the attenuation rate is non-negligible (figure 1).
To consider the effect of particle shape, we therefore examine the dependence of
$\epsilon _{p}$
on
$St$
,
$D_{*}$
and
$\chi$
. Here, we numerically compute the local average of the energy dissipation rate around each particle as

where
$\Omega _{p}$
denotes the spheroidal shell between the particle surface and another spheroid with the major and minor radii being
$a+3D_*/4$
and
$b+3D_*/4$
, respectively, and
$V_{\Omega _p}$
denotes the volume of
$\Omega _{p}$
. The first term on the right-hand side of (3.1) denotes the local energy dissipation rate around particles. However, since this quantity includes the dissipation rates due to both the shedding vortices and vortices generated by energy cascade, we subtract the latter contribution expressed by the second term so that we can estimate the additional dissipation rate. Here, note that
$\epsilon$
in the second term of (3.1) is the average quantity in time and space (i.e.
$\epsilon = \overline {\langle \epsilon ({\boldsymbol x},t) \rangle }$
).
Figure 3(a) shows the thus-evaluated
$\epsilon _{p}$
, normalised by the energy dissipation rate
$\epsilon _0$
in the single-phase flow, as a function of
$St$
. We see that
$\epsilon _{p}$
gets larger for (i) larger Stokes number
$St$
, (ii) larger aspect ratio
$\chi$
and (iii) smaller equivalent diameter
$D_{*}$
. These behaviours are similar to those of turbulence attenuation rate shown in figure 1, implying that the turbulence attenuation rate relates to
$\epsilon _{p}$
. In fact, plotting in figure 3(b) all the data according to (1.5), we can confirm that (1.5) holds irrespective of
$St$
,
$D_*$
and
$\chi$
.
In summary, (1.5) describing turbulence attenuation is valid even for the anisotropic particles. For the estimation of
$\epsilon _{p}$
required for (1.5), we have numerically evaluated it by (3.1). This direct evaluation however lacks a concrete relevance of
$\epsilon _{p}$
to particle properties such as
$St$
,
$D_{*}$
and
$\chi$
. On the other hand, for spheres,
$\epsilon _{p}$
is described in terms of the translational relative velocity
$\Delta u$
and particle diameter
$D$
according to (1.6). In the next subsection, we generalise (1.6) for anisotropic particles.
3.2. Estimation of energy dissipation rate due to particles
To estimate the energy dissipation rate
$\epsilon _p$
around particles, we first assume that
$\epsilon _p$
is determined by the translational motion of particles. Then, we can estimate it by using the force
$F_p$
acting on a particle:

where
$C_{D}$
is the drag coefficient,
$A$
(
$= \pi D_{*}^2/4$
) is the cross-sectional area of a sphere with diameter of
$D_*$
and is the translational relative velocity between the particle
$\boldsymbol {v}_{p}$
and its surrounding fluid. Here, denotes the average on the outer surface of the spheroidal shell
$\Omega _{p}$
around each particle. Since a particle subjected to
$F_p$
inputs energy to flow at the rate of
$F_{p}\Delta u$
, the average energy dissipation rate per unit mass in wakes around all particles in the system is expressed as

This estimation of
$\epsilon _p^\dagger$
differs from (1.6) used for spheres by Oka & Goto (Reference Oka and Goto2022) because they assumed that
$C_{p}$
is a constant of
$O(1)$
. For anisotropic particles, it is important to consider the dependence of
$C_D$
on the shape of particles as well as flow state surrounding them. In fact, the instantaneous drag coefficient
$C_{D}$
depends on the particle Reynolds number
$Re_{p}=\Delta u D_*/\nu$
, aspect ratio
$\chi$
and angle
$\phi$
between the major axis of the particle and relative velocity. In the present study, to evaluate
$C_D$
from
$\Delta u$
and
$\phi$
, we use the expression,

proposed by Sanjeevi, Dietiker & Padding (Reference Sanjeevi, Dietiker and Padding2022) according to their numerical results. Here,
$ C_{D, \phi } = (a_1/Re_{p}+a_2/Re_{p}^{a_3} ) \exp (-a_4 Re_{p}) + a_5 \lbrace 1-\exp (-a_4 Re_{p}) \rbrace$
, with
$a_i$
(
$i=1,\ldots ,5$
) being the coefficients depending on
$\chi$
and
$\phi$
. The values of these coefficients are listed in table 2 of Sanjeevi et al. (Reference Sanjeevi, Dietiker and Padding2022). Figure 4(a) shows the
$\chi$
-dependence of
$\langle C_{D}\rangle _{p}$
for the case of the smallest particles with the largest
$St$
(i.e.
$D_*/{L}=0.16$
and
$\gamma =512$
). For larger
$\chi$
,
$\langle C_{D}\rangle _{p}$
gets larger. This implies that
$\epsilon _p^\dagger$
is larger for more elongated particles. We have also confirmed similar monotonically increasing trend with respect to
$\chi$
for the other particles. It is not trivial whether (3.4), which is derived for a particle in uniform flow (Sanjeevi et al. Reference Sanjeevi, Dietiker and Padding2022), gives a good approximation in turbulence. However, the observation in figure 4(a) is consistent with the direct evaluation of
$\epsilon _p$
shown in figure 3(a). In addition, the increasing trend of
$\langle C_D\rangle _p$
with
$\chi$
is also consistent with the result for decaying turbulence by Schneiders et al. (Reference Schneiders, Fröhlich, Meinke and Schröder2019), who explained the trend in terms of particle orientations.

Figure 4. (
$a$
) Averaged drag coefficient
$C_D$
for particles with
$\gamma =512$
and
$D_{*}/{L}=0.16$
as a function of
$\chi$
. (
$b$
) Relative velocity normalised by
$u'$
as a function of
$St$
. Symbols are the same as in figure 1.
Furthermore, figure 4(b) shows the average
$\langle \Delta u \rangle _p$
of the relative velocity for all the simulated particles. For given
$D_*$
and
$\chi$
, the relative velocity is larger for larger
$St$
, and it tends to approximately
$u^\prime$
for
$St \gg 1$
. This is reasonable because
$St$
represents the ability for particles to follow the swirling motion around the largest vortices. Incidentally, although the feature that
${\langle \Delta u \rangle _p}\rightarrow 1.5u'$
–
$2u'$
for
$St\gg 1$
is independent of the kind of external force, the prefactor depends on forcing(see Oka & Goto Reference Oka and Goto2022, figure 6).

Figure 5. (
$a$
) Comparison of the energy dissipation rates between the numerical evaluation
$\epsilon _{p}$
using (3.1) and estimation
$\epsilon _{p}^\dagger$
using (3.3). Both values are normalised by the energy dissipation rate
$\epsilon _0$
in the single-phase flow. (
$b$
) Numerical verification of (1.5), where we use the value
$\epsilon _p^\dagger$
, instead of
$\epsilon _p$
, for the energy dissipation rate. The dashed lines in the both panels show linear lines with a slope of
$1$
. Symbols are the same as in figure 1.

Figure 6. Temporal average
$K'$
of turbulent kinetic energy normalised by the value
$K^{\prime}_{0}$
for the single-phase flow as functions of (
$a$
)
$St$
and (
$b$
)
$St_r$
. The open symbols indicate results for particles with homogeneous mass density with the same equivalent diameters (
$D_*/{L}=0.16$
and
$0.24$
) and aspect ratio (
$\chi =7$
), but different Stokes numbers. The partially filled ones indicate results for particles with different values of the inertial moment; the width of the filled part of symbols is narrower for smaller
$\beta$
.
Using the numerically evaluated averages of
$C_D\Delta u^3$
, we estimate the energy dissipation rate
$\epsilon _p^\dagger$
around particles by using (3.3). Figure 5(a) shows the relation between
$\epsilon _p^\dagger$
and
$\epsilon _p$
. We can observe that
$\epsilon _p^\dagger$
well coincides with
$\epsilon _p$
. This means that the estimation (3.3) accurately quantifies the energy dissipation rate around particles. In other words, as explicitly demonstrated in figure 5(b), we can use the estimation
$\epsilon _p^\dagger$
in place of
$\epsilon _p$
in (1.5) to describe the turbulence attenuation rate in the case for spheroidal particles.
Although overall data shown in figure 5(b), as well as figure 3(b), support (1.5) and (3.3), we observe systematic deviation from the diagonal line in these figures. More concretely, for both sets of different particle diameters (blue and black symbols represent
$D_*=0.16L$
and
$0.24L$
, respectively), the deviation is larger for smaller
$St$
(i.e. smaller
$\epsilon _p$
) cases. This is because, when
$St\lesssim 10$
, the velocity difference between particles and fluid is not determined only by the translational motion as will be discussed in the next subsection.
We emphasise the importance of the estimation (3.3) of
$\epsilon _p$
because it tells us that the energy dissipation rate around particles gets larger for larger
$\Delta u$
, larger
$C_D$
or smaller
$D_*$
in a quantitative manner. In particular, since
$\Delta u=O(u')$
for
$St\gg 1$
(figure 4
b), we may predict
$\epsilon _p^\dagger$
, and therefore
$\mathcal {A}$
, as a function of particle shape parameters
$D_*$
and
$\chi$
. We also emphasise that the estimation (3.3) is non-trivial because we estimate
$\epsilon _p$
by considering only translational motion of particles. In other words, we have neglected contribution from particle rotation in (3.3). As quantitatively demonstrated by Schneiders et al. (Reference Schneiders, Fröhlich, Meinke and Schröder2019) for decaying turbulence, in fact, the contribution to
$\epsilon _p$
from translational motion of particles dominates those from their rotations and fluid acceleration. We further discuss effects of particle rotation in the next subsection.
3.3. Effects of the rotation of particles
As demonstrated in figure 5(b), the degree of turbulence attenuation is described by additional energy dissipation rate
$\epsilon _p^\dagger$
around particles. We have also shown in figure 5(a) that
$\epsilon _p^\dagger$
can be quantified by the translational relative velocity between particles and fluid. In this subsection, we verify that the rotation of particles is less important on the turbulence attenuation when
$St\gg 1$
. To this end, we simulate different kinds of particles, which have inhomogeneous distribution of mass density inside, with different values of
$\beta$
(
$=1/4$
,
$1/16$
and
$1/64$
) for the inertial moment (see the second row from the right in table 2). The rotational Stokes number
$St_{r}$
of these particles is much smaller than particles with homogeneous density, which means that the particles are more likely to rotate together with circulation in vortices of size
$D_*$
. In such a case, the relative velocity, and therefore the energy dissipation rate
$\epsilon _p$
around particles, might be affected by difference in the swirling of vortices and tumbling of particles.
Figure 6 shows the ratio
${K'}/{K^{\prime}_{0}}$
of turbulent kinetic energy in single-phase flow and turbulence with particles with smaller values of
$\beta$
as functions of
$St$
in panel (
$a$
) and
$St_r$
in panel (b). The partially filled symbols indicate results for particles with inhomogeneous mass density (
$\beta \neq 1 $
), while open symbols show those for particles with
$ \beta =1 $
. This figure shows that when
$ St $
is sufficiently large (say,
$ St \gtrsim 10 $
), even if
$St_r $
is significantly different, the degree of the attenuation is almost identical when
$ D_*/L $
,
$St $
and
$ \chi $
are the same. This implies that the effect of the rotation for
$ \epsilon _p $
is negligible and we can estimate the energy dissipation rate by (3.3) for large-
$ St $
particles. This result is consistent with those by Schneiders et al. (Reference Schneiders, Fröhlich, Meinke and Schröder2019) and explained as follows. Since the rotational relative velocity
$ \Delta u_r $
is determined by swirls of vortices with size
$D_* $
, it is estimated by
$ \epsilon^{1/3} D_*^{1/3} $
in the inertial range of turbulence. If
$ St $
is sufficiently large, this cannot be larger than the translational relative velocity
$ \Delta u (\approx u^{\prime}) $
because
$ u' \approx \epsilon^{1/3} L^{1/3} $
and, therefore,
$ \Delta u_r / \Delta u \approx (D_*/L)^{1/3} \lt 1 $
.
It is also interesting to observe in figure 6 that the attenuation rate depends on
$St_r$
for the case with the smallest value of
$St$
(
$=5.4$
) indicated by the small blue symbols. This implies that when
$1\lesssim St\lesssim 10$
, the attenuation occurs due to energy dissipation in shedding vortices created by not only translational but also rotational motion. Therefore, the assumption for the estimation (3.3) of
$\epsilon _p$
is violated, which may lead to deviations in figures 2(
$b$
) and 4(
$b$
) of the data for
$St\lesssim 10$
.
3.4. Modulation of energy spectrum
Before concluding this article, we investigate the scale-dependent modulation of turbulent kinetic energy. In physical space, as shown in figure 2, a key ingredient of the turbulence attenuation is shedding vortices in the wake behind particles because they dissipate the energy at particle scale
$D_*$
. More precisely, because of the additional energy dissipation, such vortices reduce the cascading energy flux in scales between integral length
$L$
and particle size
$D_*$
, whereas they augment turbulence at scales smaller than
$D_*$
. This is quantitatively confirmed in figure 7, which shows the energy spectrum
${E}(k)$
, normalised by the value
${E_0}(k)$
in the single-phase flow, for spherical particles (
$\chi =1$
) in panel (
$a$
) and most elongated spheroidal ones (
$\chi =7$
) in panel (b). The darker lines are results with larger
$St$
. The wavenumber
$k_{D_*}$
(
$=2\pi /D_*$
) corresponding to participle size
$D_*$
is indicated by the red line. In both cases of spheres and spheroids, we observe that the energy spectrum is indeed attenuated in the range
$k\lesssim k_{D_*}$
and augmented for
$k \gtrsim k_{D_*}$
. This observation is similar to previous results (ten Cate et al. Reference ten Cate, Derksen, Portela and van den Akker2004; Yeo et al. Reference Yeo, Dong, Climent and Maxey2010; Gao et al. Reference Gao, Li and Wang2013; Wang et al. Reference Wang, Ayala, Gao, Andersen and Mathews2014) for spherical particles and consistent with the above-mentioned turbulence attenuation mechanism. Incidentally, the cross-over (
${E}/{E_0} \gtrless 1$
) of the modulation for spheroids (figure 7
$b$
) is obscurer than that for spheres owing to the effect of the length scale separation between major and minor radii. We can confirm this in the inset of figure 7(
$b$
), which shows the dependence of
$E/E_0$
on the aspect ratio
$\chi$
(
$=1$
,
$3$
,
$5$
and
$7$
). The cross-over becomes indeed obscurer for larger
$\chi$
.

Figure 7. Energy spectra
$E$
normalised by the value
$E_0$
in the single-phase flow as a function of the wavenumber
$k$
for (
$a$
) spheres (
$\chi =1$
) and (
$b$
) spheroids for
$\chi = 7$
with
$D_*/L=0.24$
. Darker lines show larger
$St$
. The red vertical line represents the wavenumber
$k_{D_*}=2\pi /D_*$
corresponding to the particle diameter. The inset in panel (
$b$
) shows results (close-up around
${k_D}_*$
) for particles with
$\gamma =512$
and
$D_*/L=0.24$
but different aspect ratios
$\chi$
. Darker blue lines show larger
$\chi$
.
4. Conclusions
We can describe the attenuation of turbulence by solid particles in a periodic cube in terms of the additional energy dissipation rate
$\epsilon _p$
in the wake behind particles. More precisely, by considering that the cascading energy flux is reduced by
$\epsilon _p$
, we can derive (1.5) for the attenuation rate
$\mathcal {A}$
defined as (1.2). Although (1.5) was verified for turbulence attenuation due to spherical particles (Oka & Goto Reference Oka and Goto2022), we have shown in the present study the numerical evidence that the formula holds even for spheroidal particles and the method how to estimate
$\epsilon _p$
for such anisotropic particles. In the following, we summarise concrete conclusions.
As demonstrated in figure 1, turbulent kinetic energy is more reduced by particles with larger Stokes number
$St$
, (1.1), smaller size
$D_*$
, (1.7), or larger aspect ratio
$\chi$
because
$\epsilon _p$
gets larger for particles with these properties (figure 3
a) when the volume fraction
$\Lambda$
of particles is kept constant. Incidentally, although we speculate that
$\epsilon _p$
is proportional to
$\Lambda$
for fixed
$D_*$
because of (1.6), its numerical verification is left for a future study. Figure 3(b) shows that once we estimate
$\epsilon _p$
, we can estimate
$\mathcal {A}$
through (1.5). This is the first main conclusion and it implies that, similarly to the case with spherical particles, the key quantity is
$\epsilon _p$
.
Next, we have shown that
$\epsilon _p$
can be estimated by (3.3), which depends on
$St$
,
$D_*$
and
$\chi$
. For given
$D_*$
and
$\chi$
,
$\epsilon _p$
is larger for larger
$St$
because the average relative velocity
$\Delta u$
is larger (figure 4
b). Therefore, according to (3.3),
$\epsilon _p$
takes larger values in the shedding vortices visualised in figure 2(b,d,e) in the case with large
$St$
, which leads to larger
$\mathcal {A}$
because of (1.5). Equation (3.3) also explains that
$\epsilon _p$
, and therefore
$\mathcal {A}$
, become larger for smaller
$D_*$
. Here, it is important that the estimation
$\epsilon _p^\dagger$
of the dissipation rate, (3.3), depends on the drag coefficient
$C_D$
, which is a function of
$\chi$
,
$\Delta u$
and the angle
$\phi$
between the major axis of particle and relative velocity. We have numerically estimated the average of
$C_D$
(figure 4
a) to evaluate
$\epsilon _p$
by (3.3) to show that the obtained value of
$\epsilon _p^\dagger$
almost perfectly coincides with the direct estimation of
$\epsilon _p$
(figure 5
a). Therefore, as demonstrated in figure 5(b), we can use
$\epsilon _p^\dagger$
in place of
$\epsilon _p$
in (1.5).
These results lead to the second conclusion that we can predict the attenuation rate by (1.5) through (3.3). In fact, when
$St$
is sufficiently large, particles acquire their kinetic energy from the largest vortices, and the cascading energy flux is reduced. This explains the reduction of the averaged energy spectrum
${E}(k)$
in the wavenumber range
$k\lesssim k_{D_*}$
(figure 7). On the other hand, the shedding vortices augment
${E}(k)$
for
$k\gtrsim k_{D_*}$
corresponding to length scales smaller than
$D_*$
. Since these behaviours of
${E}(k)$
for the anisotropic particles are similar to those for spherical ones (ten Cate et al. Reference ten Cate, Derksen, Portela and van den Akker2004; Yeo et al. Reference Yeo, Dong, Climent and Maxey2010; Gao et al. Reference Gao, Li and Wang2013; Wang et al. Reference Wang, Ayala, Gao, Andersen and Mathews2014), we conclude that the physical mechanism of the attenuation is also common.
We emphasise that we neglect effects of rotational motion of particles for
$\epsilon _p$
in (3.3). This is justified because the relative velocity between fluid and particles due to their rotation cannot be larger than the swirling velocity
$\epsilon ^{1/3} D_*^{1/3}$
of vortices with particle size
$D_*$
(
$\lt L$
), which is smaller than the translational relative velocity of
$O(u')$
(
$\approx \epsilon ^{1/3} L^{1/3}$
) for large-
$St$
particles. Although for small-
$St$
particles, the rotation of particles can affect the relative velocity, this cannot be pronounced in the examined turbulence at the moderate Reynolds number (
$R_\lambda \approx 50$
). This is explicitly verified in figure 6, which shows that the attenuation rate is independent of the inertial moment of particles for large-
$St$
particles. We therefore conclude for the examined large-
$St$
cases that the reason why more elongated particles attenuate turbulence more significantly is not their rotational motion, but it is due to the fact that
$\epsilon _p$
gets larger because smaller vortices are more likely to be shed from more elongated particles.
Funding.
This study was partly supported by the JSPS Grants-in-Aid for Scientific Research 20H02068 and 23K13253. The DNS were conducted on the supercomputer Fugaku through the HPCI System Research Projects (hp220232 and hp230288) and the Plasma Simulator under the auspices of the NIFS Collaboration Research Program (NIFS22KISS010).
Supplementary movies.
Supplementary movies are available at https://doi.org/10.1017/jfm.2025.137.
Declaration of interests.
The authors report no conflict of interest.
Appendix A. Verification of the assumptions for (1.4)
For the estimation of the energy flux
$\epsilon _c$
in (1.4), which is required for (1.5), we have assumed
$K\approx K_0$
and
$L\approx L_0$
. To numerically verify these assumptions, we show in figure 8(
$a$
) the average kinetic energy
$K$
of the mean flow and in figure 8(
$b$
) the integral length normalised by those (
$K_0$
and
$L_0$
) in the single-phase flow, respectively. We can confirm that
$K/K_0$
and
$L/L_0$
are approximately unity regardless of the particle properties. Although Chiarini et al. (Reference Chiarini, Cannon and Rosti2024) showed that the mean velocity field can be significantly modulated in the case with a large volume fraction of particles, this is not the case in the dilute regime investigated in the present study. It is also worth mentioning that heavy particles with
$St \gtrsim 1$
can be swept out by the largest-scale stationary vortices (Oka & Goto Reference Oka and Goto2021). Therefore, if we consider a denser case, particle clustering can affect the statistics of turbulence modulation.

Figure 8. Average of (
$a$
) kinetic energy
$K$
of the mean flow and (
$b$
) integral length
$L$
normalised by the values
$K_0$
and
$L_0$
for the single-phase flow, respectively, as a function of
$St$
. The symbols are the same as in figure 1.