3.1. Decay rates of TKE in single-phase DHIT

Before analysing the decay rate in particle-laden HIT, we quantify the corresponding decay rate in the unladen single-phase cases for later comparisons. As mentioned in the last section, with the first set-up, we start the decaying flows from eight different time frames in a sustainable HIT simulation. The interval between arbitrary two consecutive frames exceeds $20$ eddy turnover times to ensure sample independence of each realization. For the second set-up, there is no corresponding unladen DHIT for each specific particle-laden case since the flow has already been entirely disturbed. However, we may still create unladen DHIT cases with similar large-scale forcing for the purpose of comparison. Like the first set-up, these unladen DHIT cases are initialized with flow fields taken at different time frames of a sustainable HIT simulation. Sufficient time intervals are also kept between those time frames to ensure mutual independence among various realizations.

Figure 1 shows the temporal evolution of the normalized TKE, $E/E_0$, in some selective single-phase DHIT, where $E_0 = 3u_0'^2/2$ denotes the value of the initial TKE and $u_0'$ is the initial root-mean-square velocity. After the driving force is turned off, the flow still takes some time to get rid of the memory effects of the large-scale forcing and reaches a fully developed stage where the power-law decay is established. To quantify the exponent $n$ more accurately, it is important to determine the beginning point of the power-law decay. As shown in figure 1, the power-law decay starts around $(t-t_0)/\tau _0 = 3$ for different initial Reynolds numbers, where $\tau _0 = {u_0^{\prime }}^2/\epsilon _0$ is the initial large eddy turnover time, $\epsilon _0$ represents the initial dissipation rate. Other than visually identifying the starting point of power-law decay, the Taylor microscale $\lambda \propto t^{1/2}$ is another commonly used criterion (George Reference George2013). This criterion validates since a power-law decaying TKE $E(t) = E_0 t^{-n}$ would result in

(3.1)\begin{equation} \frac{{\rm d} E}{{\rm d} t} ={-}n E_0 t^{{-}n-1} ={-}n\frac{E}{t} ={-}\epsilon ={-}\frac{10\nu E}{\lambda^{2}}, \end{equation}

where $\nu$ is the kinematic viscosity, which leads to $\lambda \propto t^{1/2}$ for arbitrary exponent $n$. We also choose the $\lambda \propto t^{1/2}$ criterion to identify the starting point of the power-law decay.

Figure 1. Decay of the fluid TKE in single-phase DHIT with various initial $Re_\lambda$. The time is normalized by the initial large eddy turnover time $\tau _0 ={u_0^{\prime }}^2/\epsilon _0$, where $u_0'$ and $\epsilon _0$ are the root-mean-square velocity and dissipation rate at the initial time $t_0$, respectively. The black vertical dashed line at $t - t_0 = 3\tau _0$ shows roughly the starting point of the power-law decay.

The temporal evolution of $\lambda$ normalized by its initial value, i.e. $\lambda /\lambda _{0}$, in the single-phase cases is shown in figure 2 for some representative cases. The green dash-dotted line is a reference line for $\lambda \propto t^{1/2}$. For cases with both set-ups, the power-law decay starts roughly from $t - t_0 = 4\tau _0$. We used the least squares to fit the data points after $t - t_0 = 4\tau _0$ onto a straight line of $\lambda \propto t^{1/2}$ and obtained the $R^2$ value no less than $0.98$ in all cases. Other than the criterion mentioned above, some references, e.g. Huang & Leonard (Reference Huang and Leonard1994) and Djenidi et al. (Reference Djenidi, Kamruzzaman and Antonia2015), identify the power-law decay stage when the changing rate of the Taylor Reynolds number, $Re_\lambda$, becomes small. We also show $Re_\lambda /Re_{\lambda 0}$ of the representative cases in figure 3 and compute the relative change of $Re_\lambda$ after $t - t_0 = 4\tau _0$. It is confirmed that $\Delta Re_\lambda /Re_{\lambda 0}$ is less than 10 % for all cases, where $\Delta Re_\lambda$ is the change of $Re_\lambda$ from $t - t_0 = 4\tau _0$ to the end of the simulation, which is approximately $t_0+8\tau _0$. Based on this evidence, we believe using $t - t_0 = 4\tau _0$ as the starting time of the power-law decay is appropriate for the present study. The same criteria are used to determine the starting points of power-law decay in the particle-laden cases.

Figure 2. Temporal development of Taylor microscale $\lambda$ normalized by its initial value, $\lambda /\lambda _{0}$, of single-phase flow for several selected cases with different $Re_\lambda$. Here, $\lambda _{0}$ denotes the value of Taylor microscale at $t = 0$. The black vertical dotted line indicates the time, $t - t_0 = 4\tau _0$, after which the power-law decay is established.

Figure 3. Taylor Reynolds number $Re_\lambda$ of single-phase flow for several sets of cases with different $Re_\lambda$ in both set-ups.

After determining the starting points, the power-law decay exponent $n_{sp}$ for the single-phase DHIT can be quantified from the numerical results using the least square fitting. To enhance the accessibility of the data, we compile those data in table 2 for the simulated single-phase DHIT. Most of the exponents fall between $1.5$ and $1.8$, which is higher than the counterparts observed in the experiments of grid turbulence, as summarized by Yoffe & McComb (Reference Yoffe and McComb2018). However, such deviation could be mainly because the Reynolds numbers in our simulations are relatively small. It is well known that the decay exponent strongly depends on $Re_{\lambda }$ and increases as the turbulence weakens, e.g. from $10/7$ to $5/2$ as analysed by Batchelor (Reference Batchelor1948). In a range of Reynolds number $Re_{\lambda,0}$ between 25.8 and 51.6, Yoffe & McComb (Reference Yoffe and McComb2018) observed $1.55\le n \le 1.61$ in their direct numerical simulations, which is similar to the values we obtained. It is also notable from table 2 that for similar Reynolds numbers, $n$ can vary significantly with the specific realizations used for the flow initialization.

Table 2. Reynolds number and power-law decay exponent of single-phase DHIT.

3.2. Decay rates of TKE in particle-laden DHIT

We now move to the particle-laden cases. Figure 4 shows the temporal variation TKE in our particle-laden DHIT simulations, together with the data acquired from some previous studies (Lucci et al. Reference Lucci, Ferrante and Elghobashi2010; Luo et al. Reference Luo, Wang, Li, Tan and Fan2017; Schneiders et al. Reference Schneiders, Meinke and Schröder2017). All those previous studies generated particle-laden DHIT with the first set-up. Although these previous studies only had data points in limited time intervals as they mainly focused on studying the turbulence modulation at the initial decay stage, we can still observe that power-law decay starts to be established. One should note that the results of Luo et al. (Reference Luo, Wang, Li, Tan and Fan2017) have been rescaled in figure 4 since this work may have used a different definition of $\tau _0$, i.e. $0.5u_0'^2/\epsilon _0$ instead of $u_0'^2/\epsilon _0$ as in other studies. Together with the results from the present cases, it is safe to conclude that the power-law decay exists even when finite-size particles are present.

Figure 4. Decay of the fluid TKE in particle-laden DHIT for some representative cases: (a) cases with the first set-up, the symbols are the datasets abstracted from previous studies (Lucci et al. Reference Lucci, Ferrante and Elghobashi2010; Luo et al. Reference Luo, Wang, Li, Tan and Fan2017; Schneiders et al. Reference Schneiders, Meinke and Schröder2017); (b) cases with the second set-up.

The starting points of the power law are also determined using the $\lambda \propto t^{1/2}$ criterion. As shown in figure 5, although different ways of adding particles do affect the temporal evolution of $\lambda$ to some extent, all particle-laden cases establish the power-law decay after $t/\tau _0 \approx 4$. The $R^2$ values of the data points to fit the ideal case $\lambda \propto t^{1/2}$ are no less than 0.85 for all cases.

Figure 5. Temporal development of Taylor microscale $\lambda$ normalized by its initial value in particle-laden DHIT: (a) cases with the first set-up; (b) cases with the second set-up.

To further understand the power-law decay in particle-laden DHIT, we investigate the time evolution of total kinetic energy $K$ for the particle phase, which is the summation of translational and rotational kinetic energy, $K = 0.5m_p{\langle v^{\prime }\rangle }^2+0.5I{{\langle \omega ^{\prime }\rangle }}^2$, where $m_p$, $I$, $\langle v^{\prime }\rangle$ and ${\langle \omega ^{\prime }\rangle }$ stands for the particle mass, the moment of inertia, the averaged particle translation and angular velocity, respectively. In figure 6, the normalized particle kinetic energy change, $K/K_0$, for a few selected cases with both set-ups are presented. Like the fluid TKE, the power-law decay of particle kinetic energy also establishes after the system evolves for a certain amount of time. Compared with cases with the first set-up, cases with the second set-up generally require longer periods to reach the power-law decay stage, probably because the flow has been pre-mixed, so the memory effects of the large-scale forcing require more time to disappear. At the power-law decay stage, the exponents for the fluid and particle kinetic energy, i.e. $n_{pl}$ and $n_{pl,p}$, can be computed for each case. To enhance the data availability, the original data for these cases are provided in the Supplementary Materials.

Figure 6. Decaying particle kinetic energy in particle-laden DHIT. Here, $K_0$ denotes the values of particle kinetic energy at $t = t_0$.

Since the present study aims to investigate how particles affect the decay of TKE, we would mainly focus on the cases initialized with the first set-up, as their corresponding single-phase DHIT can be clearly defined. Figure 7 shows the correlation between $n_{pl}$ and $n_{pl,p}$ for cases belonging to the first set-up. An almost linear correlation between the two exponents is observed, which implies that the particle–fluid system can be approximately regarded as a single-phase flow at the power-law decay stage. The linearly correlated decay exponents between the fluid and particle phases could also bring convenience to modelling $n_{pl}$, as one shall see shortly in § 3.3.

Figure 7. Correlation between $n_{pl}$ and $n_{pl,p}$ for cases with the first set-up.

To find out how the power-law exponents $n_{pl}$ in the particle-laden DHIT change with the particle parameters, we classify the simulated cases according to their parameters and plot the probability distribution functions (p.d.f.s) of $n_{pl}$ in figure 8. Cases with lower particle-to-fluid density ratios $\rho$ have overall higher $n_{pl}$ than the counterparts with higher $\rho$. For $\rho = 1$, the peak of p.d.f. occurs between 1.6 and 1.7, and approximately 54.2 % of the simulated cases have $n_{pl}>1.6$. However, less than 20 % of cases with $\rho = 10$ have $n_{pl}>1.6$. The particle volume fraction $\phi _p$ also significantly impacts the p.d.f. of $n_{pl}$, where smaller values of $n_{pl}$ are more likely to be obtained with larger $\phi _p$. As can be seen in figure 8(b), for the cases $\phi _p=0.05$, $0.10$ and $0.20$, the peaks of the corresponding p.d.f.s fall at $n_{pl} = 1.5\unicode{x2013}1.6$, 1.4–1.5 and 1.3–1.4, and the percentage of cases showing $n_{pl}>1.6$ are 44.8 %, 35.4 % and 14.9 %, respectively. The particle size, however, does not show any clear preference for the distribution of $n_{pl}$, as shown in figure 8(c).

Figure 8. Histogram of the probability density of $n_{pl}$. Cases with different (a) particle-to-fluid density ratio $\rho$, (b) particle volume fraction $\phi _p$, (c) relative particle diameter $d_p/\langle \eta _{0}\rangle$ and (d) particle mass fraction $\phi _m$ distinguished by different colours.

As discussed in our recent work (Peng et al. Reference Peng, Sun and Wang2023), the contribution from $\rho$ and $\phi _p$ may be combined as the influence of particle mass fraction $\phi _m$, indicating the impact of enhanced system inertia on turbulence modulation. In figure 8(d), we classify the simulated cases into three groups with different $\phi _m$, i.e. $\phi _m \leq 0.1$, $0.1<\phi _m\leq 0.5$ and $\phi _m>0.5$, and plot their p.d.f.s of $n_{pl}$. For the cases with $\phi _m >0.5$, the peak of the corresponding p.d.f. falls between $n_{pl}=1.3\text { and }1.4$, while the other two groups peak beyond $n_{pl} = 1.5$. The percentage of cases with $n_{pl}>1.6$ are 55.6 %, 34.0 % and 4.3 % for $\phi _m \leq 0.1$, $0.1<\phi _m\leq 0.5$ and $\phi _m >0.5$, respectively. These results show a clear trend for lower $n_{pl}$ with higher $\phi _m$. To verify this trend, we select three particle mass fractions $\phi _m = 0.05$, 0.182 and 0.714, and plot $n_{pl}$ in the corresponding cases in figure 9. Although the data points are quite scattered, cases with higher $\phi _m$ are likely to have smaller exponents, indicating slower decay rates at the power-law decay stage.

Figure 9. $n_{pl}$ and $n_{pl,p}$ for three representative cases with $\phi _m = 0.050, 0.182, 0.714$ in the first set-up.

We use $n_{pl}/n_{sp}$ to quantify the modulation of the power-law decay exponent $n$ due to the presence of particles. The addition of particles increases $n$ when $n_{pl}/n_{sp} > 1$, and vice versa. In all 286 simulated cases with the first set-up, only 22.0 % have $n_{pl}/n_{sp} > 1$. In figure 10, we classify these cases into three groups according to their particle mass fractions, i.e. $\phi _m \leq 0.1$, $0.1<\phi _m\leq 0.5$ and $\phi _m >0.5$, and plot their probability distribution on $n_{pl}/n_{sp}$. The number of cases with $n_{pl}/n_{sp} > 1$ are 37.5 %, 23.6 % and 2.9 % for the three groups. The addition of particles in most of the cases makes TKE decay slower than the unladen flow during the power-law decay stage, which is opposite to the particle effect at the initial stage right after the release of particles, as reported in previous studies (Lucci et al. Reference Lucci, Ferrante and Elghobashi2010; Gao et al. Reference Gao, Li and Wang2013; Schneiders et al. Reference Schneiders, Meinke and Schröder2017).

Figure 10. Histogram of the probability density of $n_{pl}/n_{sp}$. Cases with different particle mass fraction $\phi _m$ distinguished by different colours.

In DHIT, particles not only transfer kinetic energy into the fluid phase through interfacial work that strengthens turbulence, but also enhance the dissipation rate, leading to faster decay of TKE. Both mechanisms exist throughout the whole process of decay, but their relative importance changes with time. At the initial stage, particularly right after particles are released, the mechanism of enhancing the dissipation rate dominates as the flow around the particles has not adapted to the presence of particles. As time evolves, the dissipation enhancement reduces and the interfacial transfer of TKE becomes the dominating mechanism since the kinetic energy of heavy particles decays relatively slower due to their large inertia. The same phenomenon was also reported in previous studies. In the study of Lucci et al. (Reference Lucci, Ferrante and Elghobashi2010), the time variations of the dissipation rate and particle work were both computed. Although they did not show the direct comparisons between the enhancement of dissipation rate $\Delta \epsilon = \langle \epsilon \rangle - \langle \epsilon _{sp}\rangle$ and the particle work $\psi$ after the addition of particles, one can further process their data and see the trend of $\psi$ surpasses $\Delta \epsilon$ as the flow evolves. Schneiders et al. (Reference Schneiders, Meinke and Schröder2017) directly showed the difference between $\psi$ and $\Delta \epsilon$ as a function of time. Although particles in this work are on the Kolmogorov scale compared with the Taylor scale particles in the present study, $\psi - \Delta \epsilon$ was also found to be negative at the initial stage and turned positive at a later time. The simulation results of the present study provide another piece of evidence for this observation, as shown for two selected cases in figure 11. One should also note that when a pair of single-phase and particle-laden DHIT enters the power-law decay stage, the initial amounts of TKE are different. Therefore, the faster or slower decay does not mean the amount of TKE reduction over a certain amount of time, but rather how long the power-law decay may last.

Figure 11. Temporal variation of the particle-induced dissipation rate $\Delta \epsilon$ and the particle work $\psi$ normalized by the initial dissipation rate of single-phase flow $\epsilon _{0,sp}$ in two selected cases of the present study. The particle work $\psi$ is computed from the TKE balance as ${\rm d} E/{\rm d} t + \epsilon$.

3.3. Modelling of the decay rates

In this section, we develop a model to predict the decay rate change due to the presence of particles. The model starts from the TKE budget equation, which can be derived from the method of volume averaging (Crowe et al. Reference Crowe, Schwarzkopf, Sommerfeld and Tsuji2011). Considering the isotropy and homogeneity of HIT, the TKE budget equation of the fluid phase without the mean velocity reads

(3.2)\begin{equation} \rho_f\phi_f V\frac{{\rm d}}{{\rm d} t}\left\langle\frac{1}{2}u'_\alpha u'_\alpha\right\rangle ={-}\phi_f V\rho_f \langle\epsilon\rangle + \oint_{SI}n_\beta({-}p\delta_{\alpha\beta}u_\alpha + \tau_{\alpha\beta}u_\alpha)\,{\rm d} S, \end{equation}

where $V$ is the averaging volume, which can be the volume of the whole domain $L^3$ and $\phi _f = 1 - \phi _p$ is the fluid volume fraction. Here, $\langle \cdots \rangle$ indicates that a quantity has been phase averaged, i.e. averaged over the volume occupied by a specific phase, and $\langle \epsilon \rangle$ is the phase-averaged dissipation rate. The last term of the above equation is the interfacial energy flux integrated over the entire particle surface submerged in $V$, i.e. $SI$, and $n_\beta$ is the outer normal direction pointing to the particle phase. If we ignore the particle rotation and note that on the particle surfaces, the fluid velocity must satisfy a no-slip condition, we have

(3.3)\begin{equation} \oint_{SI} n_\beta({-}p\delta_{\alpha\beta}u_\alpha + \tau_{\alpha\beta}u_\alpha)\,{\rm d} S = \sum_{i} v_{\alpha,i}\oint_{\delta S} n_\beta({-}p\delta_{\alpha\beta} + \tau_{\alpha\beta})\,{\rm d} S, \end{equation}

where $v_{\alpha,i}$ is the translational velocity of the $i$th particle immersed in $V$, and $\delta S$ is the particle surface area. For the particle phase, if we ignore the particle–particle interactions and external body forces, the equation of particle motion reads

(3.4)\begin{equation} m_p\frac{{\rm d} v_{\alpha,i}}{{\rm d} t} ={-}\oint_{\delta S} n_{\beta} ({-}p\delta_{\alpha\beta} + \tau_{\alpha\beta})\,{\rm d} S. \end{equation}

Multiplying the above equation by the particle velocity and summing over all $N_p$ particles, we have

(3.5)\begin{equation} m_p\frac{{\rm d}}{{\rm d} t}\sum_{i}\left(\frac{1}{2} v_{\alpha,i}v_{\alpha,i}\right) = N_p m_p \frac{{\rm d}}{{\rm d} t}\left\langle\frac{1}{2}v_{\alpha,i}v_{\alpha,i}\right\rangle ={-}\sum_i v_{\alpha,i}\oint_{\delta S} n_{\beta}({-}p\delta_{\alpha\beta} + \tau_{\alpha\beta})\,{\rm d} S. \end{equation}

Applying the short-hand notation, $k_f = \langle 0.5u'_\alpha u'_\alpha \rangle$ and $k_p = \langle 0.5 v_{\alpha,i}v_{\alpha,i}\rangle$, and combining (3.2) and (3.5) together, the interfacial energy flux terms are cancelled out and we obtain

(3.6)\begin{equation} (1 - \phi_p)\frac{{\rm d} k_f}{{\rm d} t} + \phi_{p}\rho \frac{{\rm d} k_p}{{\rm d} t} ={-}(1-\phi_p)\langle\epsilon\rangle. \end{equation}

In figure 12, we compare the amount of particle translational and rotational kinetic energy in our simulations. The results show that the assumption in the above derivation to ignore the particle rotation applies safely to the present simulations, given that the particle-to-fluid density ratio is considerably large. However, for more general cases where the rotational kinetic energy occupies a significant part of total kinetic energy, ignoring the particle rotation in the above derivations could lead to certain imbalances, so we may add a correction coefficient $C_{\omega }$ ($C_{\omega } \ge 1$) to the second term in (3.6) to ensure balance.

Figure 12. Correlations between the averaged translational and rotational kinetic energy of particles: (a) correlations for all cases with the first set-up at the beginning of the power-law decay; (b) temporal evolution of the correlation in a few selected cases.

Since the decay of both the fluid and particle kinetic energy follows a power law, and the exponents are roughly the same (see figure 7), we can replace ${\rm d} k_p/{\rm d} t$ as $(k_{p0}/k_{f0})\,{\rm d} k_f/{\rm d} t$ in the balance equation, where $k_{p0}$ and $k_{f0}$ are the initial kinetic energy of the particle and fluid phase when the power-law decay begins, respectively:

(3.7)\begin{equation} \left[(1 - \phi_p) + C_{\omega}\phi_p\rho\frac{k_{p0}}{k_{f0}}\right] \frac{{\rm d} k_f}{{\rm d} t} ={-}(1-\phi_p)\langle\epsilon\rangle. \end{equation}

The remaining task is to determine the averaged dissipation rate $\langle \epsilon \rangle$ on the right-hand side. It is well known that particles can create highly dissipative regions around their surface in particle-laden turbulent flows. The dissipation rate in the highly dissipative regions increases with the magnitude of the slip motion between the particle and its ambient fluid, while the volume of the region is related to the thickness of the boundary layer around a particle, which is essentially a function of the particle Reynolds number. Outside the high dissipation rate regions, the dissipation rate is close to the value in the unladen case (Burton & Eaton Reference Burton and Eaton2005; Wang et al. Reference Wang, Ardila, Ayala, Gao and Peng2016) if the particle volume fractions are low, and the flow outside the boundary layer can be regarded as undisturbed. Based on this understanding, we may formulate the dissipation rate in the particle-laden HIT as

(3.8)\begin{equation} (1 - \phi_p)\langle \epsilon\rangle = h\phi_p \left(C_{\epsilon}\frac{\nu\Delta u^2}{d_p^2}\right) + (1 - \phi_p)\langle \epsilon_{ul}\rangle, \end{equation}

where $h\phi _p$ is the total volume fraction of the highly dissipative region, $C_\epsilon$ is the coefficient for dissipation rate enhancement in the region, $\Delta u$ is the magnitude of the slip velocity and $\langle \epsilon _{ul}\rangle$ is the phase-averaged dissipation rate in the unladen case. Substituting (3.8) into (3.7) leads to

(3.9)\begin{equation} \frac{{\rm d} k_f}{{\rm d} t} ={-}\frac{h\phi_p\left(C_\epsilon \dfrac{\nu\Delta u^2}{\langle\epsilon_{ul}\rangle d_p^2}\right) + (1-\phi_p)}{(1 - \phi_p) + C_{\omega}\phi_p\rho \dfrac{k_{p0}}{k_{f0}}}\langle \epsilon_{ul}\rangle. \end{equation}

Obviously, in the limiting case of $\phi _p \rightarrow 0$, the above equation would regress to the TKE budget equation of single-phase DHIT.

The above equation indicates the two competing mechanisms on the decay rate of TKE induced by particles. The first mechanism is the increased system inertia appearing in the denominator, which reduces the decay rate. This mechanism was also identified in the point-particle simulations, where heavy particles with large inertia were found to transfer energy back to the fluid phases at the late stage of decay (Ferrante & Elghobashi Reference Ferrante and Elghobashi2003). The second mechanism is the enhanced dissipation in the vicinity of particle surfaces, which increases the decay rate. The competition between the two mechanisms could lead to either faster or slower decay of TKE.

Since the decay rates of both single-phase and particle-laden DHIT follow a certain power-law, i.e.

(3.10)\begin{equation} k_{f}(t) = k_{f0}(t - t_{0})^{n}, \end{equation}

where $k_{f0}$ is the initial TKE when the power-law decay starts, the above equation can be approximately written as

(3.11)\begin{equation} \frac{{\rm d} k_f}{{\rm d} t} ={-}n\frac{k_{f}}{t - t_0} \approx{-}n\frac{k_{f0}}{t - t_0} \end{equation}

in a short period of time after $t_0$. Given (3.9), we have

(3.12)\begin{equation} \frac{n_{pl}}{n_{sp}} = \frac{h\phi_p\left(C_\epsilon\dfrac{\nu\Delta u^2}{\langle \epsilon_{ul}\rangle d_p^2}\right) + (1-\phi_p)}{(1 - \phi_p) + C_{\omega}\phi_p\rho \dfrac{k_{p0}}{k_{f0}}}\frac{k_{f0,sp}}{k_{f0,pl}}. \end{equation}

Ideally, we may have $k_{f0,sp} = k_{f0,pl}$ if the corresponding particle-laden and single-phase cases have the same initial TKE. Although we will examine (3.12) with the simulation results shortly, it should be noted that it is difficult to establish the ideal match-up between the laden and unladen cases that satisfy all the assumptions made above. The main reason for the mismatch is that DHIT needs some time to reach the power-law decay stage, which brings difficulty to the cancellation of $t - t_0$ and $\langle \epsilon _{ul}\rangle (t)$ in the particle-laden and single-phase cases. The primary use of (3.12) is to provide some theoretical guidance on when to expect faster or slower decay of TKE in a particle-laden HIT compared with the unladen flow.

When the particle response time is greater than the Kolmogorov time but smaller than the characteristic time of the large eddies, we may expect $\Delta u/u_k = C_{\Delta u}St_k^{0.5}$, where $St_k = (2\rho +1)(d_p/\eta )^2/36$ is the particle Stokes number and $u_k$ is the Kolmogorov velocity (Balachandar Reference Balachandar2009; Ling, Parmar & Balachandar Reference Ling, Parmar and Balachandar2013). Substituting this empirical relationship into (3.12), we eventually obtain

(3.13)\begin{equation} \frac{n_{pl}}{n_{sp}} = \frac{h\phi_pC_\epsilon C_{\Delta u} \dfrac{(2\rho + 1)}{36} + (1-\phi_p)}{(1 - \phi_p) + C_{\omega}\phi_p\rho \dfrac{k_{p0}}{k_{f0}}}\frac{k_{f0,sp}}{k_{f0,pl}}. \end{equation}

To examine (3.13) with the simulation results, we need to determine the fitting parameters in (3.13), i.e. $h$, $C_\epsilon$, $C_{\Delta u}$, $C_\omega$, $k_{p0}/k_{f0}$ and $k_{f0,sp}/k_{f0,pl}$. Those parameters depend on specific simulation settings, so their optimal values could vary with cases. The rotational kinetic energy $0.5I\omega '_i\omega '_i$ is usually minor compared with the translational kinetic energy (see figure 12), so $C_\omega = 1$ can be a reasonable guess. For simplicity, $k_{p0}/k_{f0}$ and $k_{f0,sp}/k_{f0,pl}$ can be also set to one.

Here, $h$ is a coefficient indicating the size of the volume where the dissipation rate is enhanced, which can be estimated as

(3.14)\begin{equation} h = \frac{(a + \delta_{BL})^3 - a^3}{a^3}, \end{equation}

where $a = 0.5d_p$ is the particle radius. Additionally, $\delta _{BL} \sim a/\sqrt {Re_p}$ is the thickness of the particle boundary layer, which is a function of the particle Reynolds number. From the previous numerical results of radial distribution of the dissipation rate, e.g. Burton & Eaton (Reference Burton and Eaton2005), Lucci et al. (Reference Lucci, Ferrante and Elghobashi2010), Brändle de Motta et al. (Reference Brändle de Motta, Estivalezes, Climent and Vincent2016), Wang et al. (Reference Wang, Ardila, Ayala, Gao and Peng2016) and Peng et al. (Reference Peng, Ayala, Brändle de Motta and Wang2019), and the trials to determine the slip velocity between a particle in the turbulent field and its ambient fluids, e.g. Naso & Prosperetti (Reference Naso and Prosperetti2010) and Kidanemariam et al. (Reference Kidanemariam, Chan-Braun, Doychev and Uhlmann2013), with various flow and particle conditions, it is safe to set $\delta _{BL} = d_p$, leading to $h = 26$. Here, $C_{\epsilon }$ represents the relative enhancement of the dissipation rate around particles. Assuming a Stokes disturbance flow around the particle, the local dissipation rate reads (Wang et al. Reference Wang, Ardila, Ayala, Gao and Peng2016)

(3.15)\begin{equation} \epsilon = \frac{\nu\Delta u^2}{r^2}\left(\frac{3}{2}\right)^2 \left\{\cos^2\theta \left[3\left(\frac{a}{r}\right)^2 - 6\left(\frac{a}{r}\right)^4 + 2\left(\frac{a}{r}\right)^6\right] + \left(\frac{a}{r}\right)^6\right\}, \end{equation}

which means

(3.16)\begin{equation} C_{\epsilon} = \frac{d_p^2}{\nu\Delta u^2} \frac{\displaystyle\int_{0}^{2{\rm \pi}}{\rm d}\phi\int_{0}^{\rm \pi}\sin\theta \,{\rm d}\theta \int_{a}^{3a} r^2\epsilon \,{\rm d} r}{\displaystyle\int_{0}^{2{\rm \pi}}{\rm d}\phi\int_{0}^{\rm \pi} \sin\theta \,{\rm d}\theta\int_{a}^{3a} r^2 \,{\rm d} r} = 0.3704. \end{equation}

The dissipation rate around the particle is more significantly enhanced for non-Stokes disturbance flow. Since there is no analytic solution to estimate $C_\epsilon$, we may amplify $C_\epsilon$ by a factor of 2.5, which comes from the conditionally averaged dissipation rate profiles published by Burton & Eaton (Reference Burton and Eaton2005) and Wang et al. (Reference Wang, Ardila, Ayala, Gao and Peng2016).

Finally, for particles whose sizes are much smaller than the Kolmogorov length $\eta$, $C_{\Delta u}$ can be chosen as $(1-3/(2\rho + 1))^2\in (0,4)$ (Balachandar Reference Balachandar2009). However, this formula implies that the enhanced dissipation rate term vanishes for neutrally buoyant particles with $\rho = 1$, which conflicts with our observation from PR-DNS. Here, we keep $C_{\Delta u}$ as a free parameter whose optimal value is obtained by fitting the numerical results. As one shall see in the next section, the optimal value for $C_{\Delta u}$ that fits our simulation data is 0.52.

3.4. Model validations and its implication of turbulence modulation

Following the parametrization in the last section, (3.13) can be essentially simplified as

(3.17)\begin{equation} \frac{n_{pl}}{n_{sp}} = \frac{0.3478\phi_p(2\rho + 1) + (1-\phi_p)}{(1 - \phi_p) + \phi_p\rho}. \end{equation}

One can directly observe that the relative particle size $d_p/\eta$ does not appear in the equation, which indicates its minor effects on the power-law exponent. This observation is consistent with the results shown in figure 8(c), where the data points associated with each particle size almost cover the whole range of $n_{pl}$, and no explicit dependency on $d_p/\eta$ could be seen. The previous study of Lucci et al. (Reference Lucci, Ferrante and Elghobashi2010) reported a similar observation for the turbulence modulation at the early stage, where particles with fixed volume and mass fraction showed only minor size effects on the change of dissipation rate.

The comparison between the numerical simulation results and the model predictions is shown in figure 13. Overall, the match is reasonable, but the numerical data contain significant uncertainties, which creates difficulties for more precise validations. Unlike in sustained turbulence, where time averaging can significantly reduce the uncertainty of quantifying turbulence modulation, in time-evolving DHIT, the relative change of the decay exponent, i.e. $n_{pl}/n_{sp}$, is strongly affected by the specific flow realization used for the initialization. When different initial flow fields are created at different time frames of the same sustainable HIT simulation, the decay rate of the unladen DHIT could change significantly, i.e. $n_{sp,max}/n_{sp,min} > 1.33$. For a given set of particle parameters, the Reynolds number when the power-law decay starts, $Re_{\lambda,d}$, and the decay exponent in particle-laden DHIT, $n_{pl}$, also vary with the initial conditions, as demonstrated in figure 14. Although we try multiple realizations for each set of particle parameters, the uncertainties of $n_{pl}/n_{sp}$ are still significant. To enhance data availability, the value of $n_{pl}/n_{sp}$ for each specific simulation is provided in the Supplementary Materials.

Figure 13. Comparison between the model predicted $n_{pl}/n_{sp}$ and the simulation results.

Figure 14. Variation of $n_{pl}/n_{sp}$ for cases with the same set of particle parameters.

With certain manipulation and approximation, (3.17) can be recast as

(3.18)\begin{equation} \frac{n_{pl}}{n_{sp}} = \frac{0.3478\phi_p(2\rho + 1) + (1-\phi_p)}{(1 - \phi_p) + \phi_p\rho} \approx 1 - \frac{(0.3 - 0.35\rho)\phi_p}{(1 - \phi_p) + \phi_p\rho}. \end{equation}

The last term in this equation is close to the definition of the mass increase due to the addition of particles. If more mass is added to enhance the system's inertia, a smaller decay rate is achieved compared with the corresponding unladen case. From (3.18), we have $0.35 - 0.3\rho > 0$, or $\rho < 1.17$, to result in a faster decay rate in $n_{np}/n_{sp} > 1$. The value of $1.17$ may not be accurate because of the various assumptions and estimations we made leading to (3.17), but the overall trend that heavier particles result in slower decay of TKE at the power-law decay stage is supported by the numerical results. In figure 15(a), we classify different cases according to the density ratio and plot their probability distribution on $n_{np}/n_{sp}$. For $\rho =2, 5, 10$, 76.4 %, 88.9 % and 88.6 % of the simulated cases have $n_{np}/n_{sp} < 1$, respectively. However, only $58.3\,\%$ of the neutrally buoyant particle cases have $n_{np}/n_{sp} < 1$.

Figure 15. (a) Histogram of the probability density of $n_{pl}/n_{sp}$. Cases with different particle-to-fluid density ratio $\rho$ distinguished by different colours. (b) Average values of ${n_{pl}}/{n_{sp}}$ for cases in set-up 1. The dotted line represents ${n_{pl}}/{n_{sp}} = 1$.

The importance of density ratio on turbulence modulation was widely recognized. In the initial stages of turbulent decay, most previous studies observed that the addition of particles accelerated TKE decay, and such acceleration amplifies with the density ratio $\rho$ and particle volume fraction $\phi _p$ (Lucci et al. Reference Lucci, Ferrante and Elghobashi2010; Gao et al. Reference Gao, Li and Wang2013; Schneiders et al. Reference Schneiders, Meinke and Schröder2017). However, the current study shows that when the turbulence decay enters the later power-law decay stage, the presence of particles tends to slow down the turbulence decay. Schneiders et al. (Reference Schneiders, Meinke and Schröder2017) found that even at the early stage, particles could release kinetic energy to the fluid by doing work, compensating for the dissipation enhancement around the particle surface. We believe that this mechanism is more pronounced in the power-law decay stage, resulting in a slower decay rate.

Despite most cases with heavy particles, e.g. $\rho = 5,10$, showing lower decay rates than the single-phase case, some exceptions still show $n_{np}/n_{sp} > 1$. These exceptions are exclusively with the small volume fraction $\phi _p = 0.05$. Based on this observation, it is reasonable to deduce that a sufficiently large particle volume fraction is a necessary condition to trigger slower decay rates in the power-law decay stage.

Finally, previous studies based on the Lagrangian point-particle simulations reported that particles resulted in faster TKE decay when the particle Stokes number $\tau _p/\tau _k = (1+2\rho )(d_p/\eta )^2/36 > 1$. This criterion does not apply to DHIT laden with finite-size particles presented in the present study. Even at the power-law decay stage where turbulence is relatively weak and $\eta$ is large, most of our simulations still have $\tau _p/\tau _k > 1$. However, slower decay of TKE was more frequently observed. This conflict is partially because the particle Stokes number is no longer a suitable indicator for turbulence modulations resulting from finite-size particles (Lucci, Ferrante & Elghobashi Reference Lucci, Ferrante and Elghobashi2011; Shen et al. Reference Shen, Peng, Wu, Chong, Lu and Wang2022; Peng et al. Reference Peng, Sun and Wang2023), although some researchers believed the particle Stokes number was still applicable but the fluid characteristic time should change from Kolmogorov time to turnover time of large scales (Oka & Goto Reference Oka and Goto2022). More importantly, it indicates a very different picture of turbulence modulation at the early and late stages of decay. At the early decay stage, the enhanced dissipation rate around particle surfaces overwhelms the interface transfer of kinetic energy, while the interface transfer of kinetic energy becomes dominant as the flow further evolves.