1. Introduction
Interactions between turbulent flows and dispersed particles exist in many natural phenomena and engineering applications. When the particle size or inertia is comparable to or larger than those of the smallest eddies in a turbulent flow, dispersed particles would inevitably modify the local flow structures or even the overall properties of turbulent flow. This phenomenon is known as turbulence modulation. Specifically, how the presence of particles changes the overall intensity of turbulent flows and what physical parameters govern this change have yet to be made fully clear (Balachandar & Eaton Reference Balachandar and Eaton2010; Brandt & Coletti Reference Brandt and Coletti2022).
Many previous efforts were devoted to revealing the mechanisms responsible for the modulation of turbulent intensity due to the presence of particles. Through direct numerical simulations, Squires & Eaton (Reference Squires and Eaton1990) found that particles with sizes much smaller than the Kolmogorov length (labelled as small-size particles) could significantly attenuate the turbulent intensity of homogeneous isotropic turbulence (HIT) due to the additional dissipation resulting from particle presence. Paris (Reference Paris2001) used particle image velocimetry (PIV) to investigate the impact of small-size particles on the intensity of turbulent channel flow. The measurements showed only slightly enhanced dissipation rates, but the turbulent intensity was significantly attenuated, which indicated that introducing additional dissipation was not the only mechanism of turbulence attenuation. Opposite to the previous two studies, Yang & Shy (Reference Yang and Shy2005) observed augmentation of turbulent intensity in HIT after releasing microglass, copper and lead beads into HIT. The transfer of potential energy to turbulent kinetic energy (TKE) was responsible for these augmentations. Tanaka & Eaton (Reference Tanaka and Eaton2010) measured the flow around particles in HIT with high-resolution PIV. Their results showed significantly reduced turbulent intensity around particles compared with the far field. This observation also qualitatively agreed with results from direct numerical simulations (DNS) with static particles (Burton & Eaton Reference Burton and Eaton2005; Botto & Prosperetti Reference Botto and Prosperetti2012), indicating that the enhanced flow inertia due to the addition of heavy particles could be an essential mechanism for turbulence attenuation. Kajishima et al. (Reference Kajishima, Takiguchi, Hamasaki and Miyake2001) investigated the effect of settling particles on the intensity of a vertical channel flow. The shedding vortices from settling particles were identified as an important reason for the enhancements of TKE. Uhlmann (Reference Uhlmann2008) further showed that even without vortex shedding, the flow instabilities generated in the particle wake could still trigger turbulence augmentations.
Besides the overall turbulent intensity, particles could also modify the homogeneity of turbulence and the distribution of TKE on different scales. Ten Cate et al. (Reference Ten Cate, Derksen, Portela and Van Den Akker2004) numerically investigated the energy spectra in HIT laden with finite-size particles and found that TKE on the scales above the diameter of the particles was suppressed, and TKE on the small scales was enhanced. Vreman (Reference Vreman2016) showed the existence of TKE flux in particle-laden HIT from the far field to the particle surface due to the low TKE and high dissipation rate regions formed around particles. Turbulence modulation becomes more complex for background flows that are intrinsically inhomogeneous and anisotropic. With small-size glass and copper beads, the experiments of Kulick, Fessler & Eaton (Reference Kulick, Fessler and Eaton1994) showed that turbulence modulation occurred inhomogeneously and anisotropically in a turbulent channel flow. The most intense modulation of TKE was observed for the wall-normal velocity component at the channel centre. For finite-size particles, some recent fully resolved DNS demonstrated that particles resulted in a more homogeneous and isotropic distribution of TKE in the turbulent channel and pipe flows (Shao, Wu & Yu Reference Shao, Wu and Yu2012; Picano, Breugem & Brandt Reference Picano, Breugem and Brandt2015; Eshghinejadfard et al. Reference Eshghinejadfard, Abdelsamie, Hosseini and Thévenin2017). Through a complete TKE budget analysis, Peng, Ayala & Wang (Reference Peng, Ayala and Wang2019c) attributed this modification to the exhibition of TKE production in the buffer layer and the enhancements of TKE transport and intercomponent transfer, but what parameters control these mechanisms are yet to be revealed.
Aside from revealing modulation mechanisms, it is equally essential to summarize the controlling parameters and build models to quantify turbulence modulation. By examining the experimental results of particle-laden turbulent jet and pipe flows from 1960 to 1980, Gore & Crowe (Reference Gore and Crowe1989) found that particles with sizes exceeding 1/10 of the size of energetic eddies tended to enhance turbulence, whereas smaller size particles attenuated turbulence. The particle density and flow conditions that played a role in determining turbulence modulation were not considered. Based on the effects of shedding vortices on turbulence modulation, Hetsroni (Reference Hetsroni1989) proposed to use the particle Reynolds number to predict whether turbulence would be augmented or attenuated by particles. When the particle Reynolds number exceeds 400, turbulence augmentation would be observed; otherwise, turbulence attenuation is expected. Righetti & Romano (Reference Righetti and Romano2004) treated the particle Stokes number defined with the characteristic time of large eddies as the criterion for turbulence modulation. Turbulence augmentation and attenuation are expected when this parameter is above and lower than unity, respectively. Combining the influence of particle Stokes number, flow Reynolds number and the departure between Kolmogorov scale and integral scale, Tanaka & Eaton (Reference Tanaka and Eaton2008) introduced a new non-dimensional parameter called particle momentum number to predict the direction of turbulence modulation. At large or small particle momentum numbers, particles would enhance turbulence, whereas the opposite could happen at intermediate particle momentum numbers. A similar criterion was introduced by Luo, Luo & Fan (Reference Luo, Luo and Fan2016), which included the same aspects of the particle momentum number, but differed in the detailed definition. Recently, Yu et al. (Reference Yu, Xia, Guo and Lin2021) proposed two new criteria to qualitatively predict the modulation of the overall turbulence intensity in the turbulent channel flow and the local effect near the channel centre. Unlike the previous ones, the new criteria considered the wall effects by including the channel width rather than only making judgments based on turbulent integral scales. For HIT, Oka & Goto (Reference Oka and Goto2022) refined the criterion by Gore & Crowe (Reference Gore and Crowe1989), they argued that turbulent attenuation required the particle size to be not only below the integral length scale but also significantly higher than the ratio between the Taylor microscale and the square root of particle-to-fluid density ratio.
Although solid progress was made in understanding turbulence modulation, two areas still need improvement. First, previous studies only partially investigated how the dimensional or dimensionless parameters governed turbulence modulation. For example, the effect of particle volume fraction was investigated by Ten Cate et al. (Reference Ten Cate, Derksen, Portela and Van Den Akker2004), Lucci, Ferrante & Elghobashi (Reference Lucci, Ferrante and Elghobashi2010), Cisse et al. (Reference Cisse, Saw, Gibert, Bodenschatz and Bec2015) and Yousefi, Ardekani & Brandt (Reference Yousefi, Ardekani and Brandt2020), the effect of particle Stokes number was studied by Lucci, Ferrante & Elghobashi (Reference Lucci, Ferrante and Elghobashi2011), Abdelsamie & Lee (Reference Abdelsamie and Lee2013) and Oka & Goto (Reference Oka and Goto2022), and the effect of particle size was reported in Shao et al. (Reference Shao, Wu and Yu2012), Uhlmann & Chouippe (Reference Uhlmann and Chouippe2017) and Peng et al. (Reference Peng, Ayala and Wang2019c) among many others, but systematic parameterization studies are rare. Second, most of the abovementioned criteria can only qualitatively predict the overall modulation effect. Models that quantitatively predict the modulation levels are yet to be built. In the recent work of Oka & Goto (Reference Oka and Goto2022), the relative change of TKE between the single-phase and particle-laden HIT was found to be linearly dependent on the ratio between the particle-induced dissipation rate and the dissipation rate of single-phase HIT. However, this model cannot predict the modulation level since the particle-induced dissipation rate is not an a priori property.
In order to enhance the understanding of turbulence modulation, the present study conducts fully resolved DNS for forced particle-laden HIT with particles of size around the Taylor microscale. The simulation results serve two purposes, to deliver a more complete parameterization study on the modulation of TKE and dissipation rate in HIT and to build models to predict the modulation levels of TKE and dissipation rate based on only the information of single-phase HIT and dispersed particles. The choice of forced HIT as carrier turbulence brings convenience for measuring the modulation of overall turbulent intensity and dissipation rate, as time averaging can be taken in addition to spatial averaging. From the net energy transfer perspective, it will be shown that the coupling between the forcing scales and the particle scale is weak.
The remainder of the paper is structured as follows. In § 2, dimensional analysis is conducted to identify the non-dimensional parameters that affect the turbulence modulation and determine the simulation set-ups accordingly. The numerical approach, a brief validation on the choice of grid resolution, and a discussion on the large-scale forcing to sustain turbulence are also presented. The parameterization study on the turbulence modulation levels and the models to make quantitative predictions of the modulation effects are given in § 3. Finally, the main conclusions of the present study are summarized in § 4.
2. Parameter setting and simulation details
Simulations in this work are conducted by the lattice Boltzmann method (LBM) with the interpolated bounce-back schemes to enforce the no-slip boundary condition on the particle surfaces. This numerical approach has been validated thoroughly in many cases of both laminar and turbulent flows, with and without particles; thus, no further validation is performed here. For more information on the numerical approach and its validation tests, readers may refer to our previous work in Peng (Reference Peng2018) and Brändle de Motta et al. (Reference Brändle de Motta2019).
We consider a HIT laden with non-settling, monodisperse, spherical solid particles. The physical parameters in the simulations are prescribed through a dimensional analysis. The modulation levels of the TKE and dissipation rate are defined as 
$K_r = (K_{sp} - K_{pl})/K_{sp}$ and 
$D_r = (\varepsilon _{sp} - \varepsilon _{pl})/\varepsilon _{sp}$, respectively, where 
$K_{sp}$, 
$\varepsilon _{sp}$ are TKE and dissipation rate in the single-phase HIT, respectively, and 
$K_{pl}$, 
$\varepsilon _{pl}$ are their corresponding quantities in particle-laden cases. The greater 
$K_r$ and 
$D_r$ are, the particles introduce a more significant modulation. For HIT laden with finite-size spherical particles, the modulation levels are affected by TKE of the background turbulence (BT) 
$E_{sp}$, the dissipation rate of the BT 
$\varepsilon _{sp}$, the viscosity 
$\nu$, the fluid density 
$\rho _f$, the particle density 
$\rho _p$, the particle diameter 
$d_p$ and the particle volume fraction 
$\phi _{v}$ – seven parameters in total. The first four parameters identify a background HIT, and the last three uniquely define a dispersed phase. The seven parameters contain only three fundamental units; thus, the modulation levels are controlled by four non-dimensional parameters. There are multiple ways to choose these four non-dimensional parameters. Here we choose the following ones that have clear physical meanings. They are (1) 
$Re_{\lambda }$, the Reynolds number defined by the Taylor microscale in the unladen background turbulence that measures the scale separation of the carrier flow, (2) 
$\phi _{v}$, the particle volume fraction, (3) 
$d_p/\eta$, the ratio between the particle diameter and the Kolmogorov length that measures the relative size of particles respect to turbulent eddies and (4) 
$\rho _p/\rho _f$, the density ratio between the particle and fluid phases. The other non-dimensional parameters, such as particle Stokes number 
$St = \tau _p/\tau _f$, and particle mass fraction 
$\phi _m$, can be derived from the four selected non-dimensional parameters.
The third and last non-dimensional parameters in the following equations may be combined to yield the Stokes number 
$St = \tau _p/\tau _{K} = (\rho_p/\rho_f)(d_p/\eta)^2/18$, which is the non-dimensional parameter to measure the inertia of particles compared with the inertia of the smallest turbulent eddies:
\begin{equation} K_r = f\left( Re_{\lambda}, \phi_{v},\frac{d_p}{\eta},\frac{\rho_p}{\rho_f}\right),\quad \varepsilon_r = g\left( Re_{\lambda}, \phi_{v},\frac{d_p}{\eta},\frac{\rho_p}{\rho_f}\right). \end{equation}
Early studies of HIT laden with inertial particles had revealed that 
$St$ was a reasonable indicator to measure the turbulence modulation by particles, together with the particle volume fraction and mass loading. For finite-size particles, whether 
$St$ is still a good indicator for the level of turbulence modulation was investigated by Lucci et al. (Reference Lucci, Ferrante and Elghobashi2011), and numerical results obtained from particle-laden decaying HIT led to an unfavourable conclusion. However, under the numerical setting of Lucci et al. (Reference Lucci, Ferrante and Elghobashi2011), the majority effect of turbulence modulation was due to the sudden dissipation rate jump at the moment of releasing particles into the decaying HIT, which did not provide sufficient evidence to answer the question for forced HIT.
A complete parameterization for the modulation effect requires varying all four non-dimensional parameters. The most rigorous way to achieve this is to conduct an orthogonal experimental design. However, it would result in too many cases and lead to unaffordable computational costs. In order to maintain a reasonable computational cost, we set up 13 particle-laden cases in total, whose parameters are tabulated in table 1. Case 1-A (where ‘1’ is the index for the BT, and ‘A’ labels a group of particle parameters) to 1-C share the same 
$Re_{\lambda }$, 
$St$ and 
$\phi _{v}$, and these cases serve the purpose to answer whether 
$St$ is the indicator of turbulence modulation in forced particle-laden HIT. Case 1-B, 1-D and 1-E are different only in the volume fraction 
$\phi _{v}$, with the other three non-dimensional parameters all identical. Case 1-E, 1-F and 1-G are different only in the relative particle diameter 
$d_p/\eta$. Case 1-F, 1-H and 1-I differ only in the particle-to-fluid density ratio 
$\rho _p/\rho _f$. Cases 1-E, 2-E and 3-E examine the impact of BT on turbulence modulation. Another set, Cases 1-G, 2-G and 3-G, also serves this purpose.
Table 1. Parameter settings of the examined cases. Quantities from the second to the last column are particle-to-fluid density ratio, relative particle diameter, particle mesh size, particle Stokes number, particle volume fraction, particle mass fraction, the number of particles.
The mesh sizes to represent a particle, or particle mesh sizes, are 
$d_p/\Delta x > 16$ in the present study. The resolution requirement to fully resolve the particle boundary layer mainly depends on the particle Reynolds number 
$Re_p$. It is often required that 
$d_p/\Delta x \gg \sqrt {Re_p}$. In our study, the averaged particle Reynolds numbers are below 30 in most of the simulated cases, where the slip velocity of a finite-size particle can be quantified by the method proposed in Kidanemariam et al. (Reference Kidanemariam, Chan-Braun, Doychev and Uhlmann2013). This requirement is roughly satisfied in our simulations where 
$16 \le d_p/\Delta x \le 32$. In our previous study, we also compared the results of particle-laden decaying HIT simulations with two particle mesh sizes, 
$d_p/\Delta x = 12$ and 
$d_p/\Delta x = 24$, at a higher initial 
$Re_{\lambda } = 87.6$ than those adopted in the present study. The comparisons of flow and particle statistics showed that 
$d_p/\Delta x = 12$ already captured flow TKE, dissipation rate, and the fluctuation levels of particle velocity and angular velocity quite well (Peng et al. Reference Peng, Ayala, Brándle de Motta and Wang2019a). Therefore, the particle mesh size 
$d_p/\Delta x = 16$ is sufficient for the present study. This mesh size was adopted only for two cases out of the 13 cases in total. In the other 11 cases, finer grid resolutions are used.
To investigate the effect of 
$Re_{\lambda }$, three unladen BT, labelled as 1-BT, 2-BT and 3-BT from strong to weak, were created using the stochastic forcing scheme of Eswaran & Pope (Reference Eswaran and Pope1988). The estimated rate of the energy input 
$E_{in}$ is given as
\begin{equation} E_{in} = \frac{4N_f\sigma_f^2T_f}{1+T_f(\sigma_f^2T_f N_f k_0^2)^{1/3}\beta}, \end{equation}
where 
$N_f$ is the number of total force modes, 
$\sigma _f^2$ is the forcing magnitude, 
$T_f$ is the forcing time scale, 
$k_0 = 1$ is the lowest wavenumber in spectral units and 
$\beta = 0.8$ is a fitting parameter in the scheme. For all the simulated cases in this study, the identical 80 large-scale modes whose wavenumbers satisfying 
$0 < |{\boldsymbol k}| < \sqrt {8}$ are forced. The forcing scheme of Eswaran & Pope (Reference Eswaran and Pope1988) was chosen because of two reasons. First, (2.2) gives excellent predictions of the energy input when the time scale in the scheme is small compared with the Kolmogorov time (Rosa et al. Reference Rosa, Parishani, Ayala and Wang2015). Second, it is convenient to use, given that the scheme does not require predefined TKE spectra as input. After a simulation reaches its statistically stationary state, the ensemble averages of the turbulent statistics are computed and tabulated in table 2. The uncertainties after ‘
$\pm$’ for the averaged properties are computed as
\begin{equation} \sigma_{\bar{A}} = \sigma_{A}\sqrt{\frac{2T_c}{T_{ave}}}, \end{equation}
where 
$\sigma _A$ is the standard deviation of time-dependent quantity 
$A(t)$, 
$T_{ave}$ is the duration of the statistical averaging and 
$T_c$ is the correlation time of 
$A(t)$, which is obtained from the correlation coefficient 
$R(T_c) \equiv \langle A(t_1)A(t_1 + T_c)\rangle /\sigma _A^2 = 0.5$ (Wang et al. Reference Wang, Ayala, Gao, Andersen and Mathews2014). The present study uses a consistent grid mesh of 
$512^3$ for all simulated cases. This mesh generates a spatial resolution of 
$k_{max}\eta > 7$, which fully resolves the smallest turbulent eddies. The benchmark results from our in-house pseudospectral method (PSM) code for Case 1-BT, i.e. the strongest turbulence, with a grid mesh of 
$256^3$ are in the table for comparison. Excellent matches are observed between LBM and PSM results. The normalized energy and dissipation rate spectra of Case 1-BT are compared in figure 1. Once again, almost perfect agreements are obtained, which validates the choice of the grid resolution in the present study. The lack of clear inertial subranges in the energy and dissipation rate spectra indicates that the Reynolds numbers investigated in the present study are not high enough to separate the large and small scales fully. Therefore, caution should be taken when generalizing the conclusions of the present study to flows at higher Reynolds numbers.
Table 2. Statistics of the background HIT. Quantities from the second to the last columns are turbulent kinetic energy 
$E$, fluctuation velocity 
$u'$, dissipation rate 
$\varepsilon$, Reynolds number based on the Taylor microscale 
$Re_\lambda$, Kolmogorov length 
$\eta$, longitudinal integral length scale 
$L$, skewness 
$S$ and flatness 
$F$ of the longitudinal velocity gradient, and the eddy turnover time 
$T_e$. Unless otherwise specified, results throughout the article are presented in the spectral units.
Figure 1. Normalized spectra of TKE and dissipation rate, i.e. 
$E(k)/(\sum _k E(k))$ and 
$\varepsilon (k)/(\sum _k \varepsilon (k))$ for the background turbulence in Case 1-BT with the highest Reynolds number 
$Re_\lambda \approx 63$ examined in the present study.
After the simulations of the BT reach statistically stationary states, particles are released into the flow at randomly picked locations. Then, the flows laden with particles evolve until the new statistical stationary states are reached. Each particle-laden case's statistics are collected over 3000 time-frames covering 30 eddy turnover times, the same as in the single-phase cases.
Conceptually, when finite-size particles are present, the continuity of the flow field is disrupted, and there was concern that the large-scale forcing should no longer be used to create sustained turbulent fields (see the argument in Lucci et al. (Reference Lucci, Ferrante and Elghobashi2010)). However, the impact of the finite-size particles on the forced HIT simulations with large-scale forcing still needs to be quantified. In order to assess this influence, we examine the spectra of the energy input, i.e.
\begin{equation} W(k) = \frac{1}{T_{ave}}\int_{t_0}^{t_0 + T_{ave}}\sum_{ k - 0.5 <\vert \boldsymbol{\kappa}\vert\le k + 0.5}(\,\widetilde{{\boldsymbol f}}({\boldsymbol{\kappa},t})\boldsymbol{\cdot}\widetilde{{\boldsymbol u}}({\boldsymbol{\kappa},t}))\,{\rm d}t, \end{equation}
for the three cases with the highest particle volume fractions, i.e. Cases 1-A, 1-B and 1-C. As seen from figure 2, in the particle-laden cases, most of the energy is still introduced to flow fields in the form of large scales. The percentage of the energy input under the desired large-scale modes, i.e. 
$0 < |{\boldsymbol k}| < \sqrt {8}$ is 
$92.6\,\%$ 
$91.9\,\%$ and 
$91.7\,\%$ in Case 1-A, 1-B and 1-C, respectively. Physically speaking, when finite-size particles are present in a sustained turbulent field, a certain amount of energy would be introduced in the form of the work done by the particles, which occurs around the particle scales (see the TKE budget analyses in Peng et al. (Reference Peng, Ayala and Wang2019c)). From the energy input perspective, using large-scale stochastic forcing to sustain HIT would not likely sabotage the physical foundation of the present study, at least not qualitatively. Since the regions occupied by particles cannot be forced, the total amounts of forces added in the particle-laden cases would be smaller than their counterparts in the corresponding unladen cases. In order to fairly assess the levels of turbulence modulation introduced by particles, we slightly increase the forcing magnitude 
$\sigma _f^2$ in the particle-laden cases, so that the 
$E_{in}^{pl}$ estimated by (2.2) equals to 
$E_{in}^{sp}/(1-\phi _{v})$, where 
$E_{in}^{pl}$ and 
$E_{in}^{sp}$ are the estimated rates of energy input in the particle-laden and single-phase cases, respectively.
Figure 2. Comparison of the normalized spectra of energy input 
$W(k)/(\sum _{k}W(k))$ in the three particle-laden cases with the highest particle volume fraction 
$\phi _v = 0.12$ examined in the present study and their corresponding single-phase case. Energy in the single-phase case was purely introduced under large scales, i.e. 
$W(k\ge 3) =0$.
3. Results and discussions
The statistics of the 13 particle-laden cases and the three unladen cases are compared systematically in this section. Specific attention is given to the modulations of TKE and dissipation rate, the two essential quantities that define a HIT. In order to enhance the availability of the data, the statistics are presented in both tables and figures. Table 3 shows the statistics among Cases 1-BT, 1-A, 1-B and 1-C. Aside from those statistics listed in table 2, the statistic averages of particle translational and angular velocity fluctuations, i.e. 
$u_p'$ and 
$\omega _p'$, are also given. The three particle-laden cases have the same BT, particle volume fraction and particle Stokes number but their density ratios and particle relative sizes are varied. With the same particle Stokes numbers, the modulation levels of TKE and dissipation rate are more significant in the cases with smaller but heavier particles. This observation confirms that the particle Stokes number is not a good indicator of the turbulence modulation levels with finite-size particles, consistent with the observations of Lucci et al. (Reference Lucci, Ferrante and Elghobashi2010) in decaying HIT, and Shen et al. (Reference Shen, Peng, Wu, Chong, Lu and Wang2022) in forced HIT. It should be noted that for particles with moderate density ratios, the particle Stokes number is more frequently defined as 
$(1+ 2\rho _p/\rho _f)(d_p/\eta )^2/36$. The modulation of TKE and dissipation rate in the cases with the same particle volume fraction and the background HIT are plotted against both definitions of the particle Stokes number in figure 3. The conclusion that the particle Stokes number is not an appropriate quantifier for turbulence modulation is unaffected.
Table 3. Flow and particle statistics of cases with the same 
$Re_{\lambda }$, 
$\phi _v$ and 
$St$, but different in 
$d_p/\eta$ and 
$\rho _p/\rho _f$.
Figure 3. Comparison of turbulence modulation with the particle Stokes number using two different definitions: (a) 
$St_1 = (\rho _p/\rho _f)(d_p/ \eta )^2/18$; (b) 
$St_2 = (1+ 2\rho _p/\rho _f)(d_p/ \eta )^2/36$.
The higher modulation levels with smaller but heavier particles are associated with at least two mechanisms, the increased system inertia due to the addition of heavy particles that reduce the energy input and the formation of low TKE but high dissipation rate regions around particle surfaces (Balachandar & Eaton Reference Balachandar and Eaton2010). The presence of heavy particles increases the inertia of the fluid–particle system and lowers the fluctuation velocities created by the large-scale forcing. As a result, smaller amounts of energy are input into the system under the same magnitude of forcing. This mechanism could imply a strong dependency of turbulence modulation on the mass fraction of particles, as shown in figure 4. As shall be seen later in this section, such an observation would also benefit our modelling of turbulence modulation against the non-dimensional parameters summarized in (2.1a,b). Moreover, high dissipation rate regions would form around particles due to larger relative motion (Burton & Eaton Reference Burton and Eaton2005; Botto & Prosperetti Reference Botto and Prosperetti2012). Since cases with smaller particles have larger total surface areas when the particle volume fractions are identical, their attenuation effects on TKE and dissipation rates are more profound.
Figure 4. Dependencies of TKE and dissipation rate modulation on the particle mass fraction in Case 1-A, 1-B and 1-C with the same BT, particle volume fraction and particle Stokes number.
Since the particle Stokes number is not a good indicator for turbulence modulation, the particle size and particle-to-fluid density effects must be examined separately. The statistics of the three cases with different particle sizes are compared in table 4. The modulations of TKE and dissipation rate are also presented in figure 5, and both inversely depend on the particle size. As discussed earlier, the more significant modulations with smaller particles come from the larger total particle surface areas that create more low energy and high dissipation rate regions. On the other hand, a larger individual particle not only leads to smaller fluctuation velocities but also follows the ambient fluid motion less tightly, resulting in a more significant dissipation rate enhancement around its surface (Oka & Goto Reference Oka and Goto2022). For finite-size particles with moderate 
$St$, the latter effect could not compete with the effect of surface area increase.
Figure 5. Comparison of TKE and dissipation rate modulation between Case 1-E, 1-F and 1-G with the same background turbulence, density ratio and particle volume fraction but different relative particle sizes.
Table 4. Flow and particle statistics of cases with the same 
$Re_{\lambda }$, 
$\rho _p/\rho _f$ and 
$\phi _v$, but different in 
$d_p/\eta$.
The statistics of the three cases with different particle-to-fluid density ratios are tabulated in table 5. As also shown in figure 6, while the modulation of TKE amplifies with the density ratio, the dissipation rates are almost independent. This result implies that the amplification of the dissipation rate is dominated by the increased particle surfaces and only marginally influenced by particle inertia. For forced HIT, the dissipation rates are statistically identical to energy input rates. Although the three cases have almost the same energy input, their turbulence intensities differ. This observation indicates that heavy particles also change the energy distribution among different scales besides modifying the overall energy input. As shown by the comparison of energy spectra in figure 7, particles decrease TKE in large scales and enhance TKE in small scales. The ‘pivot’ wavenumber of this modulation corresponds to the particle diameter 
$k_p = 16$ under the current setting. This observation has been reported in many previous studies of turbulence modulation by finite-size particles in HIT (Ten Cate et al. Reference Ten Cate, Derksen, Portela and Van Den Akker2004; Lucci et al. Reference Lucci, Ferrante and Elghobashi2010). Two mechanisms could be responsible for this modulation of energy spectra. On the one hand, both experiments (Tanaka & Eaton Reference Tanaka and Eaton2010) and numerical simulations (Burton & Eaton Reference Burton and Eaton2005; Botto & Prosperetti Reference Botto and Prosperetti2012) had observed the formation of high dissipation rate regions around particles. Part of the TKE driving the relative motion between the fluid and a particle would dissipate directly in the viscous boundary layer near the particle surface, which forms a unique cross-scale energy transfer from large scales to particle–surface boundary scales in addition to the usual energy cascade process. On the other hand, the cascade of TKE from large to small scales in a single-phase HIT is through large vortices stretching and breaking into smaller ones. The disturbances introduced by particles could stimulate this process, which is shown by the comparison of vortex structures in figure 8, where many more small-scale vortices were created in the particle-laden case compared with the single-phase case. Compared with less dense particles, heavier particles yield more significant slip motions relative to the ambient fluid motions, which break down large scales more effectively and form more dissipative regions. As a result, lower TKE levels were reached in the cases with heavier particles under similar energy input rates.
Table 5. Flow and particle statistics of cases with the same 
$Re_{\lambda }$, 
$\phi _v$ and 
$d_p/\eta$, but different in 
$\rho _p/\rho _f$.
Figure 6. The comparison of TKE and dissipation rate modulation between Case 1-F, 1-H and 1-I with the same BT, particle volume fraction, relative particle sizes, but different particle-to-fluid density ratios.
Figure 7. Comparison of the energy spectra between Case 1-F, 1-H and 1-I with different particle-to-fluid density ratios: (a) TKE spectra; (b) zoom-in plot of TKE spectra around the pivot wavenumber 
$k_{p} = 16$.
Figure 8. Vortex structures of (a) the single-phase (Case 1-BT) and (b) particle-laden (Case 1-F) HIT. The vortices are visualized by the Q-criterion with 
$Q = 5000$. The colour on the isosurfaces represent the magnitude of the flow speed 
$\vert {\boldsymbol u}\vert$.
Table 6 shows the turbulent statistics of the three cases with different particle volume fractions and the corresponding single-phase case. There is no doubt that the modulation levels of TKE and dissipation rate increase with the particle volume fraction. As shown in figure 9, the modulation levels of TKE and dissipation rate scale linearly with 
$\phi _v^{(2/3)}$, which are consistent with the experimental measurements of Cisse et al. (Reference Cisse, Saw, Gibert, Bodenschatz and Bec2015) for a von Kármán flow laden with neutrally buoyant particles. This again confirms that the low TKE and high dissipation rate regions around the particles play essential roles in turbulence modulation since 
$\phi _v^{(2/3)}$ represents the total surface area of dispersed particles. In the present study, the modulation levels of TKE increase twice as fast with the total particle surface area than the modulation levels of dissipation rate, whereas an identical increasing rate was reported in the work of Cisse et al. (Reference Cisse, Saw, Gibert, Bodenschatz and Bec2015). This difference is probably because heavy particles are used in the present study, while neutrally buoyant particles were used in the experiments of Cisse et al. (Reference Cisse, Saw, Gibert, Bodenschatz and Bec2015). As shown earlier in figure 6, the increase of particle-to-fluid density significantly escalates the modulation of TKE but introduces little impact on the modulation of the dissipation rate.
Figure 9. Comparison of TKE and dissipation rate modulation between Case 1-B, 1-D and 1-E with the same background turbulence, relative particle sizes, particle-to-fluid density ratio, but different particle volume fractions.
Table 6. Flow and particle statistics of cases with the same 
$Re_{\lambda }$, 
$d_p/\eta$ and 
$\rho _p/\rho _f$, but different in 
$\phi _v$.
Finally, the flow and particle statistics for cases with different BT are shown in table 7. The corresponding modulation levels of TKE and dissipation rate are plotted in figure 10. Increasing 
$Re_\lambda$ has opposite effects on the modulation of TKE in the two groups of cases with different particle sizes, i.e. increasing 
$Re_\lambda$ leads to greater modulation in Case E with 
$d_p/\eta = 10.134$, but opposite changes in Case G with smaller particles 
$d_p/\eta = 7.166$. Accounting for the uncertainties of the results, the impacts brought by changing 
$Re_\lambda$ on the modulation of TKE are minor. This phenomenon is probably because values of 
$Re_\lambda$ are close in the three BT. The influence of 
$Re_\lambda$ on the dissipation rate modulation is more evident than its counterpart on TKE. Overall, particles with the same values of 
$\rho _p/\rho _f$, 
$d_p/\eta$ and 
$\phi _v$ lead to greater relative increases of dissipation rate for flows with higher 
$Re_\lambda$. This trend is not hard to comprehend. As discussed earlier, the modulation levels of the dissipation rate are dominated by the total area of particle surfaces in the domain, which can be quantified as 
$(N_p{\rm \pi} d_p^2)/V$, where 
$V = (N_p{\rm \pi} d_p^3)/(6\phi _v)$ is the size of the computational domain. This dimensional quantity scales with 
$1/d_p$ or 
$1/\eta$, since 
$d_p/\eta$ is maintained as constant in each group of cases. As shown in table 2, the Kolmogorov length 
$\eta$ reduces as 
$Re_\lambda$ increases, so the particles introduce larger total surface areas in cases with higher 
$Re_\lambda$ and cause more significant modulations to the dissipation rate accordingly.
Table 7. Flow and particle statistics among cases with the same 
$\phi _v$, 
$d_p/\eta$ and 
$\rho _p/\rho _f$, but different in 
$Re_{\lambda }$.
Figure 10. Comparison of TKE and dissipation rate modulation between Case 1-E(G), 2-E(G) and 3-E(G) with the same particle-to-fluid density ratio, particle volume fraction, relative particle sizes, but different background turbulence.
The above assessment gives us a qualitative impression of how an individual flow or particle parameter affects turbulence modulation. In order to build a reliable model to quantify the modulation levels on TKE and dissipation rate, the effects of these dimensionless parameters need to be considered together. Starting from the results of the dimensional analysis in (2.1a,b), and further assuming the dependencies of 
$K_r$ and 
$D_r$ on the four non-dimensional parameters, i.e. 
$\phi _v$, 
$Re_\lambda$, 
$d_p/\eta$ and 
$\rho _p/\rho _f$, follow a power law, we have
\begin{equation} K_r = C_K \phi_v^{a}\,Re_\lambda^{b}\left(\frac{d_p}{\eta}\right)^{c}\left(\frac{\rho_p}{\rho_f}\right)^{d},\quad D_r = C_\varepsilon\phi_v^{e}\,Re_\lambda^{f}\left(\frac{d_p}{\eta}\right)^{g}\left(\frac{\rho_p}{\rho_f}\right)^{h}. \end{equation}The power dependencies of the modulation on the non-dimensional parameters are inspired by some classical turbulence theories, and they are rather empirical and lack a rigorous physical foundation. Different sets of coefficients are tested in the above models, but they could not fit the data points well. As analysed earlier, one of the primary mechanisms responsible for turbulence modulation is the enhanced system inertia that reduces the rate of energy input under large-scale forces. This implies that the mass fraction of the particles, i.e.
\begin{equation} \phi_m = \frac{\rho_p \phi_v}{\rho_p\phi_v + (1-\phi_v)\rho_f}, \end{equation}
could be a better parameter to model the modulation effects than 
$\phi _v$. After replacing 
$\phi _v$ with 
$\phi _m$ in (3.1a,b), the numerical data can be matched well, as shown in figure 11, and the following models are obtained to evaluate the modulation levels of TKE and dissipation rates in forced HIT. They are
\begin{gather} K_r=5.2\phi_m^{1.1}\left(\frac{d_p}{\eta}\right)^{{-}0.4}\left( \frac{\rho_p}{\rho_f}\right)^{{-}0.3}, \end{gather}
\begin{gather}D_r=0.105\,Re_\lambda^{1.0}\phi_m^{1.0}\left(\frac{d_p}{\eta}\right)^{{-}0.8}\left( \frac{\rho_p}{\rho_f}\right)^{{-}0.35}. \end{gather}
The coefficients in these models are qualitatively consistent with the observations made earlier in this study. The absence of 
$Re_\lambda$ in 
$K_r$ indicates that the modulation of TKE is not sensitive to the flow Reynolds number 
$Re_\lambda$. In contrast, the positive dependency of 
$Re_\lambda$ in 
$D_r$ indicates that the modulation of dissipation rate is more significant in more energetic turbulence. The negative powers for 
$d_p/\eta$ in the models of 
$K_r$ and 
$D_r$ show that small-size particles introduce more significant modulation to both TKE and dissipation rate. Since the dependencies of 
$K_r$ and 
$D_r$ on the particle size are 
$-0.4$ and 
$-0.8$, the particle size has a more significant impact on the modulation of dissipation rate than the counterpart of TKE, which can also be seen in figure 5. The negative powers on 
$\rho _p/\rho _f$ do not mean that less dense particles bring more significant modulations since the mass fraction 
$\phi _m$ positively scales with 
$\rho _p/\rho _f$. Overall, the modulation levels of TKE and dissipation rate still positively scale with the particle-to-fluid density ratio, especially when the density ratio and the particle volume fraction are small, as shown in figure 12. The positive dependency of the modulation levels on the particle volume fraction 
$\phi _v$ has been absorbed into the corresponding dependency on 
$\phi _m$. Although 
$K_r\propto \phi _m^{1.1}$ and 
$D_r\propto \phi _m^{1.0}$ do not precisely reflect 
$K_r\propto \phi _v^{2/3}$ and 
$D_r\propto \phi _v^{2/3}$ presented earlier in figure 9, the overall nonlinear dependencies of the turbulence modulation on the particle volume fraction are still preserved.
Figure 11. Empirical models for quantitative predictions of the turbulence modulation of (a) TKE and (b) dissipation rate in forced HIT due to finite-size particles. The coefficients in these models are obtained by fitting in the simulation results obtained in the present study.
Figure 12. Predicted turbulence modulation as functions of the particle-to-fluid density ratio: (a) modulation of TKE 
$K_r$; (b) modulation of the dissipation rate 
$D_r$.
Finally, relevant results in the literature are collected and compared with the proposed models in figure 13 to examine their versatility. These results come from the following studies on forced HIT laden with finite-size particles.
(i) Ten Cate et al. (Reference Ten Cate, Derksen, Portela and Van Den Akker2004). This work used a LBM coupled with the immersed boundary method to simulate the particle-flow system. A similar stochastic force was used to sustain the turbulence. Results from five particle-laden cases in this work are extracted for comparison. The Reynolds number of the BT was
$Re_\lambda = 60.98$. The particle-to-fluid $\rho _p/\rho _f$ ranged from $1.146$ to $1.728$, the relative particle size $d_p/\eta$ was fixed at $16$, and the particle volume fraction $\phi _v$ ranged from $2\,\%$ to $10\,\%$. The flow was resolved with a $256^3$ mesh, and the particle mesh size was $d_p/\Delta x = 8$. The results of $K_r$ and $D_r$ used for comparison were not directly given but calculated from table 4 of this work.(ii) Yeo et al. (Reference Yeo, Dong, Climent and Maxey2010). This work used pseudospectral codes coupled with the force coupling method to study the particle and bubble-laden HIT. The turbulence was driven by the stochastic forcing scheme of Eswaran & Pope (Reference Eswaran and Pope1988), the same as in the present study. Four cases where particles being the dispersed phase are used for comparison (the ‘S2-M’ case where only force monopole terms were used has been discarded), two of which were with neutrally buoyant particles
$\rho _p/\rho _f = 1$, and the other two cases were with solid particles $\rho _p/\rho _f = 1.4$. The $Re_\lambda$ of the BT was fixed at $58.7$. Two particle volume fractions $5.7\,\%$ and $6.0\,\%$, and two relative particle sizes, $d_p/\eta = 11.06$ and $7.72$, were used. The flow was resolved with $128^3$ and $192^3$ meshes, and the particle mesh sizes were $d_p/\Delta x = 5.32-5.56$.(iii) Bellani et al. (Reference Bellani, Byron, Collignon, Meyer and Variano2012). This work conducted PIV measurements to compare turbulence modulation with spherical and ellipsoidal particles. The case with spherical particles is used for comparison. Here
$Re_\lambda$ of the BT is $269$. The particles are nearly neutrally buoyant, with $\rho _p/\rho _f\approx 1.02$ (the fluid is water at $5\,^\circ {\rm C}$, and the density of particles is 1020 kg m$^{-3}$). The particle volume fraction is $0.14\,\%$ and the relative particle size $d_p/\eta$ is $21$.(iv) Uhlmann & Chouippe (Reference Uhlmann and Chouippe2017). This work used a second-order finite-volume method coupled with the immersed boundary method to simulate the particle-flow system. The forcing scheme of Eswaran & Pope (Reference Eswaran and Pope1988) was also used to generate the turbulence. Two particle-laden cases are used for comparison. These two cases had
$Re_\lambda$ of $119$ and $142.8$, $d_p/\eta = 5.5$ and $11$, and identical $\rho _p/\rho _f = 1.5$ and $\phi _v = 0.5\,\%$. The HIT was resolved with large meshes of $1024^3$ and $2048^3$, and the particle mesh size was $d_p/\Delta x \approx 16$. It should be noted that this work's main purpose was to investigate particles’ acceleration and clustering phenomenon in HIT rather than turbulence modulation. Thus, the modulation effects of the two particle-laden cases were mentioned briefly in the text rather than presented formally in tables or figures. A data point of this kind was also found in the work of Chouippe & Uhlmann (Reference Chouippe and Uhlmann2015) (Case A-G0), whose main purpose was to investigate the effect of gravity on the turbulence–particle interactions.(v) Shen et al. (Reference Shen, Peng, Wu, Chong, Lu and Wang2022). This work used our earlier version of lattice Boltzmann codes to investigate the effect of particle Stokes number on the turbulence modulation in HIT. The numerical approach is very similar to the one used in the present study, with only two differences. First, the flow solver in the work of Shen et al. (Reference Shen, Peng, Wu, Chong, Lu and Wang2022) was based on the D3Q19 lattice, whereas the flow solver in the present work used the D3Q27 lattice. The latter lattice preserves better isotropy in the simulations of circular pipe flows and is numerically more stable than the former one (Peng et al. Reference Peng, Geneva, Guo and Wang2018). Second, when particles were present, Shen et al. (Reference Shen, Peng, Wu, Chong, Lu and Wang2022) did not change the magnitude of the stochastic forcing, while in the present study, the magnitude of the forcing term was amplified to account for the volume occupied by the solid particles. Nine cases in Shen et al. (Reference Shen, Peng, Wu, Chong, Lu and Wang2022) are used for comparison. Those cases have the same
$Re_\lambda = 73.8$ and the same $\phi _v - 9\,\%$, but different $\rho _p/\rho _f$ ranging from 5 to 20 and $d_p/\eta$ ranging from 8.89 to 17.77. The flow was simulated with a grid resolution of $256^3$. Particle mesh sizes between $d_p/\Delta x = 8$ to $16$ were used to represent particles with different diameters.(vi) Oka & Goto (Reference Oka and Goto2022). This work used a second-order finite-difference scheme coupled with the immersed boundary method to investigate turbulence modulation in a forced HIT. Two forcing schemes were used to generate turbulence. The first scheme used a time-independent forcing field defined as sine/cosine waves, and the second scheme was a large-scale deterministic force defined in the Fourier space. Twelve particle-laden cases with the latter forcing scheme are adopted for comparison. These 12 cases had a fixed Reynolds number
$Re_\lambda = 94$ in the BT, a fixed particle volume fraction $\phi _v = 0.82\,\%$, four density ratios $\rho _p/\rho _f = 2$, $8$, $32$ and $128$, and three particle sizes $d_p/\eta = 8$, $16$ and $32$, where each combination of those parameters formed a separate case. The flows in those cases were resolved with a $256^3$ mesh, with a particle mesh size ranging from $d_p/\Delta x = 8$ to $32$, corresponding to the three particle sizes. It should be noted that the four cases with $d_p/\eta = 32$ had only eight particles in the domain, so the uncertainty of the results could be significant. This work also reported four cases with a single large particle of $d_p/\eta = 64$ in HIT. Considering the uncertainty of the results and the fact that the diameter of the particle is $1/4$ of the size of the periodic domain, these four cases are not included here.
Figure 13. Comparisons between the established empirical models in (3.3) and the relevant results reported in the literature: (a) TKE; (b) dissipation rate.
Given that most of the above investigations used particles with relatively small density ratios, different from the heavy particles chosen in the present study, three additional simulations with less dense particles are conducted. Their parameters are tabulated in table 8.
Table 8. Parameters of the three additional particle-laden simulations.
As shown in figure 13, except for the results in Shen et al. (Reference Shen, Peng, Wu, Chong, Lu and Wang2022), most data points from the literature concentrate around the bottom left-hand corner of the plots. This is mainly because the particle-to-fluid density ratios in the corresponding studies are close to unity, which leads to small 
$\phi _m$ values. Overall, the predictions of TKE modulation by the model in (3.3) are reasonable, but the model for the dissipation rate modulation does not fit all the literature data points. However, caution must be given when assessing the credibility of the data points. The turbulent statistics are obtained from time averaging over different time frames, but the uncertainties for the averaging were omitted for most data points. It is also questionable whether some studies have provided sufficient spatial resolution in the simulations, especially when representing particles. Per our experience and some benchmark studies in the literature (e.g. Breugem Reference Breugem2012; Peng, Ayala & Wang Reference Peng, Ayala and Wang2019b) for particle-laden flows, a particle mesh size of 
$d_p/\Delta x > 12$ is often required to resolve the fluid–particle interactions to the satisfactory accuracy.
The implementation details of the large-scale forcing may also play a role in the deviation of the model predictions from the data points. The modulation levels reported by Shen et al. (Reference Shen, Peng, Wu, Chong, Lu and Wang2022) are slightly higher than the counterparts in the present study. This difference is likely due to the amplification of the forcing magnitude that enlarges the energy input and compensates for turbulence attenuation in the present study. Some data points in figure 13 indicate enhancements of dissipation rate, even though the corresponding TKE are attenuated. Since for a forced HIT the dissipation rate should statistically balance with the rate of energy input, i.e. 
$\langle \varepsilon \rangle = \langle\, {\boldsymbol f} \boldsymbol {\cdot } {\boldsymbol u}\rangle$, where 
${\boldsymbol f}$ is the driving force of HIT and 
${\boldsymbol u}$ is the fluid velocity. When particles attenuate 
${\boldsymbol u}$, 
$\langle \varepsilon \rangle$ is unlikely to be augmented if 
${\boldsymbol f}$ stays the same. Unfortunately, those implementation details were often not given, which brings difficulties in assessing the credibility of the data points. Before more evidence confirms the versatility of the established models, it is safer to restrict the use of these models for moderate particle-to-fluid density ratios and Reynolds numbers. Beyond the forced HIT investigated in this study, the importance of interpreting turbulence modulation observations based on specific numerical or experimental settings should also be emphasized. For example, in particle-laden turbulent channel/pipe flows, maintaining a fixed flow flux (e.g. Picano et al. Reference Picano, Breugem and Brandt2015) or a fixed pressure drop (e.g. Peng et al. Reference Peng, Ayala and Wang2019c) among the unladen and laden cases would lead to opposite turbulence modulation observations. The angle between gravity and flow directions, i.e. horizontal or vertical channel, could also generate a significant impact on how particles modify the TKE distribution in different velocity components (Kulick et al. Reference Kulick, Fessler and Eaton1994; Kussin & Sommerfeld Reference Kussin and Sommerfeld2002).
4. Discussion and conclusions
The present work conducted a parameterization study to investigate the dependencies of turbulence modulation effects on the flow and particle parameters in particle-laden forced HIT and built models to quantitatively predict the modulation of TKE and dissipation rate due to finite-size heavy particles. The main findings are summarized as follows.
Previously, there have been concerns that the presence of finite-size particles would sabotage the reliability of using large-scale forcing to study the two-way interactions between particles and the HIT. We found that over 
$91\,\%$ of the energy input from the applied forcing scheme was still associated with the large scales, even for particles with the diameter around the Taylor microscale and a volume fraction of 
$0.12$.
Turbulence modulation in HIT by spherical particles depends on four non-dimensional parameters: 
$Re_\lambda$ of the BT; the particle volume fraction 
$\phi _v$; the particle-to-fluid density ratio 
$\rho _p/\rho _f$; and the relative particle size 
$d_p/\eta$. The particle Stokes number is not a good indicator to measure the modulation of TKE and dissipation rate by finite-size particles.
Overall, the presence of heavy particles attenuates turbulence in HIT. This attenuation is mainly associated with two mechanisms: the increase of system inertia that reduces the rate of energy input and the formation of low TKE; and high dissipation rate regions around the particle surface. Some other mechanisms leading to turbulence augmentation could also exist, but their contributions are less dominating in the investigated problem. For example, the flow instability generated by the particle wake is an essential mechanism of turbulence modulation leading to the generation of TKE (Balachandar & Eaton Reference Balachandar and Eaton2010). However, this mechanism only becomes dominating when the particle Reynolds numbers are large enough (Kajishima et al. Reference Kajishima, Takiguchi, Hamasaki and Miyake2001; Uhlmann Reference Uhlmann2008). A recent study by Yu et al. (Reference Yu, Xia, Guo and Lin2021) concluded that turbulence augmentation in a vertical channel flow could be observed when the non-dimensional parameter
\begin{equation} \chi = \frac{Re_p}{Re_b^{0.33}(d_p/H)^{0.61}\rho_r^{0.05}} > 20, \end{equation}
where 
$Re_p$ is the particle Reynolds number, 
$Re_b$ is the bulk Reynolds number of the channel flow, 
$d_p$ is the particle diameter, 
$H$ is the half-channel width and 
$\rho _r$ is the particle-to-fluid density ratio. In the present study of particle-laden forced HIT, the particle Reynolds numbers due to the slip motion between a particle and its ambient fluids are small in most cases. For example, the averaged 
$Re_p$ in Case 1-A, 1-B and 1-C are 29.6, 19.2 and 11.3, respectively, where the slip velocity of a finite-size particle was calculated using the method proposed in Kidanemariam et al. (Reference Kidanemariam, Chan-Braun, Doychev and Uhlmann2013). Due to dominating attenuation mechanisms, the particle wakes under these small particle Reynolds numbers cannot generate sufficient TKE to compensate for the TKE reduction.
The attenuation of TKE and dissipation rate is more profound with smaller particles than with larger particles. This trend is mainly because the former case has larger total particle surface areas and creates more regions with low TKE and high dissipation rates. The modulation of TKE and dissipation rate is found to scale linearly with the total particle surface area, i.e. 
$\phi _v^{(2/3)}$. The modulation of TKE increases with the particle-to-fluid density ratio, but the modulation of the dissipation rate is not sensitive to the density ratio. On the contrary, the increase of Reynolds number in the BT brings significant modulation to the dissipation rate, but it has only minor effects on the modulation of TKE. Compared with the dissipation rate, the presence of heavy particles leads more significant attenuation effect on TKE. Aside from reducing the energy input by enhancing the system inertia, heavy particles also modify the energy cascade from large to small scales. By breaking up flow scales larger than the particle diameter and creating smaller ones, particles make the energy cascade from large to small scales more efficiently, escalating the attenuation of TKE.
Empirical models have been developed to quantitatively predict the modulation of TKE and dissipation rate in HIT based on non-dimensional parameters. These models are compared with the results of turbulence modulation in the literature, showing a certain level of agreement. Certain deviations are also found between the model predictions and the literature data, perhaps partially due to the limited reliability of the literature data.
Acknowledgements
C.P. acknowledges the financial support from Shandong University through the Qilu Young Scholar Program. Computing resources are provided by the National Supercomputing Center at Jinan and the Center for Computational Science and Engineering of the Southern University of Science.
Funding
This work has been supported by the National Natural Science Foundation of China (C.P., grant numbers 12102232, U2006221, 52176040), (L.-P. W., grant numbers 91852205, 11961131006, U2241269); the Natural Science Foundation of Shandong Province, China (C.P., grant numbers 2022HWYQ-023, ZR2021QA010); the Taizhou-Shenzhen Innovation Center, Guangdong Provincial Key Laboratory of Turbulence Research and Applications (L.-P. W., 2019B21203001); Guangdong-Hong Kong-Macao Joint Laboratory for Data-Driven Fluid Mechanics and Engineering Applications (L.-P. W., grant number 2020B1212030001); and Shenzhen Science & Technology Program (L.-P. W., grant number KQTD20180411143441009).
Data availability
Data will be made available on request.
Declaration of interests
The authors report no conflict of interest.
